snappyHexMesh on sharp corners (cube)
 

I'm trying to mesh a *really simple* geometry with snappyHexMesh but I'm getting weird results upon viewing it in Paraview.
Before getting into details, I should say that I'm using OpenFOAM 2.0.1. I know this is a relatively old release, but right now I need to use some libraries that currently doesn't work with very recent OpenFOAM releases and I won't have time to port them in the next couple of months. Since what I'm trying to achieve is very simple, 2.0.1 might be enough for now.

This is the geometry I'm looking for: a 3000x3000x1000m domain with a cube on the floor, centered. Cube is 150m in length. The bottom part of the cube should lie on the domain "floor". The cube is just a representation of the geometry I'm going to have there (yes, the geometry is still relatively simple and will contain sharp corners).

See the attached images of what I'm obtaining. The translucent cube is where my mesh should be; it represents the 150x150x150 cube. It is like the mesh is snapping to something smaller than my geometry. And since the cube geometry I'm getting is smaller that what I want, the cube end up not touching the floor and I get cells in between the cube and floor.

These are the steps I'm taking:

1. run blockMesh on the outside domain

2. use decomposePar (scotch method)

3. use surfaceFeatureExtract by using the following command:
surfaceFeatureExtract -includedAngle 150 -writeObj constant/triSurface/cube_ig.stl cubeSFE
(in OpenFOAM 2.0.1 we don't need the dictionary)

4. run snappyHexMesh (now using ptscotch in decomposeParDict)
mpirun -np 32 snappyHexMesh -noFunctionObjects -overwrite -parallel

5. Reconstruct the mesh, remove processor* dirs, decompose again (scotch method)


Some observations:
- The STL file looks ok. I have used surfaceCheck to verify. It is so simple that I had even verified it manually.

- I'm running checkMesh as I go. Everything looks fine. The cube boundary is being added as expected.

- I'm aware that the 5th step can be replaced by using changeDictionary, but I still haven't got to that point yet (any tips are appreciated, though).

- Since my geometry can be meshed with just regular hexahedrons, I tried turning off the snap feature on snappyHexMeshDict. But that didn't really help. It doesn't matter if is on or off, still doesn't work with my current settings

- I did notice I'm getting a warning during sHM, but don't really know what it means:
I have attached the dictionaries and a simple script that goes through steps 1--5 outlined above. The output of the script is the defective mesh.

Any ideas on what I'm missing?

Best,
Regis