problem reading exports from GMSH

[ews-ffarella]

Edit: I will try to compile GMSH from source on the mdolab image and let you know

I am trying to read multi-block grid generatey by GMSH version 4.12.2 but somehow it fails.

The code to generate the surface is:

lc = 100;
cx = 324925.0;
cy = 5216025.0;
s = 3000.0;
r = 16000.0;

// Number of cells at inner square
Ns = 12;
// Number of cells between inner square and circle
Ni = 36;

// Horizontal gradings
blockGradi = 1.1;
blockGrads =  1.0;

// Inner square side curvature - 0.4
scr = 1.1;
sc = s * scr;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// 45 degree points angle
a0 = -45;
a1 = -135;
a2 = 135;
a3 = 45;

// Half of 45 degree points angle
ea0 = 0;
ea1 = -90;
ea2 = 180;
ea3 = 90;

ca0 = Cos((Pi/180)*a0);
ca1 = Cos((Pi/180)*a1);
ca2 = Cos((Pi/180)*a2);
ca3 = Cos((Pi/180)*a3);

sa0 = Sin((Pi/180)*a0);
sa1 = Sin((Pi/180)*a1);
sa2 = Sin((Pi/180)*a2);
sa3 = Sin((Pi/180)*a3);

cea0 = Cos((Pi/180)*ea0);
cea1 = Cos((Pi/180)*ea1);
cea2 = Cos((Pi/180)*ea2);
cea3 = Cos((Pi/180)*ea3);

sea0 = Sin((Pi/180)*ea0);
sea1 = Sin((Pi/180)*ea1);
sea2 = Sin((Pi/180)*ea2);
sea3 = Sin((Pi/180)*ea3);

// Inner square x and y position

Point(1) = {cx, cy, 0, lc};

Point(11) = {cx-s, cy-s, 0, lc};
Point(12) = {cx-s, cy+s, 0, lc};
Point(13) = {cx+s, cy+s, 0, lc};
Point(14) = {cx+s, cy-s, 0, lc};

Point(21) = {cx+r*ca1, cy+r*sa1, 0, lc};
Point(22) = {cx+r*ca2, cy+r*sa2, 0, lc};
Point(23) = {cx+r*ca3, cy+r*sa3, 0, lc};
Point(24) = {cx+r*ca0, cy+r*sa0, 0, lc};

Point(31) = {cx-sc, cy, 0, lc};
Point(32) = {cx, cy+sc, 0, lc};
Point(33) = {cx+sc, cy, 0, lc};
Point(34) = {cx, cy-sc, 0, lc};

Circle(1) = {21, 1, 22};
Circle(2) = {22, 1, 23};
Circle(3) = {23, 1, 24};
Circle(4) = {24, 1, 21};

// Gmsh provides other curve primitives than stright lines: splines,
// B-splines, circle arcs, ellipse arcs, etc. Here we define a new
// circle arc, starting at point 14 and ending at point 16, with the
// circle's center being the point 15:

Spline(21) = {11,31,12};
Spline(22) = {12,32,13};
Spline(23) = {13,33,14};
Spline(24) = {14,34,11};

Line(31) = {11, 21};
Line(32) = {12, 22};
Line(33) = {13, 23};
Line(34) = {14, 24};

Curve Loop(1) = {21, 22, 23, 24};
Surface(1) = {1};

Curve Loop(11) = {1, -32, -21, 31};
Surface(11) = {11};

Curve Loop(12) = {32, 2, -33, -22};
Surface(12) = {12};

Curve Loop(13) = {33, 3, -34, -23};
Surface(13) = {13};

Curve Loop(14) = {34, 4, -31, -24};
Surface(14) = {14};

Transfinite Curve {1, 2, 3, 4} = Ni+1 Using Progression 1.0;
Transfinite Curve {21, 22, 23, 24} = Ni+1 Using Progression 1.0;
Transfinite Curve {31, 32,33,34} = Ns+1 Using Progression blockGradi;

Transfinite Surface {1} = {11, 12, 13, 14};
Transfinite Surface {11} = {21, 22, 12, 11};
Transfinite Surface {12} = {22, 23, 13, 12};
Transfinite Surface {13} = {23, 24, 14, 13};
Transfinite Surface {14} = {24, 21, 11, 14};

Recombine Surface {1, 11, 12, 13, 14};

Physical Curve("LEFT") = {1};
Physical Curve("TOP") = {2};
Physical Curve("RIGHT") = {3};
Physical Curve("BOTTOM") = {4};
Physical Surface("TERRAIN1") = {1};
Physical Surface("TERRAIN2") = {11};
Physical Surface("TERRAIN3") = {12};
Physical Surface("TERRAIN4") = {13};
Physical Surface("TERRAIN5") = {14};

Mesh 2;

Mesh.Smoothing = 0;
Mesh.RecombineAll = 1;
Mesh.Smoothing = 0;

Mesh.CgnsExportStructured = 1;
Mesh.Format = 32;   //CGNS fails
// Mesh.Format = 28;    //P3D

So I export the data using the PLOT3D Driver, and pyhyp can load it. Somehow, the normals are inverted, but I can live with it. Now, Paraview 5.11.2 can read the CGNS output generated by GMSH (Mesh.Format = 32). So I would like to do a bit of processing, show below, modify the output of GMSH, and save it to VTM and then CGNS from paraview.

baseDir = os.path.dirname(os.path.abspath(__file__))
surfaceFile = os.path.join(baseDir, "circular_domain.p3d")
print(surfaceFile)
assert os.path.isfile(surfaceFile), surfaceFile

elevFile = os.path.join(baseDir, "grossglockner_elevation_data.vts")
print(elevFile)
assert os.path.isfile(elevFile), elevFile

import vtk
from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk

dem_reader = vtk.vtkXMLStructuredGridReader()
dem_reader.SetFileName(elevFile)
dem_reader.Update()
dem = dem_reader.GetOutput()
dem_pts = vtk_to_numpy(dem.GetPoints().GetData())
dem_pts[:, 2] = 0
pts = vtk.vtkPoints()
pts.SetData(numpy_to_vtk(dem_pts))
dem.SetPoints(pts)

if 1:   
    reader = vtk.vtkCGNSReader()
    reader.SetFileName("circular_domain.cgns")
    reader.Update()

    _org = reader.GetOutput()
    out = _org.GetBlock(0)
    for b in range(out.GetNumberOfBlocks()):
        block = out.GetBlock(b)
        pdata = vtk.vtkPolyData()
        pdata.SetPoints(block.GetPoints())
        # print(pdata.GetPoints().GetNumberOfPoints())
        probeFilter = vtk.vtkProbeFilter()
        probeFilter.SetInputData(pdata)
        probeFilter.SetSourceData(dem)
        probeFilter.Update()
        z = vtk_to_numpy(probeFilter.GetOutput().GetPointData().GetArray("Elevation"))
        print((z.min(), z.max()))
        pts3d = vtk_to_numpy(block.GetPoints().GetData())
        pts3d[:, 2] = -z
        pts = vtk.vtkPoints()
        pts.SetData(numpy_to_vtk(pts3d))
        block.SetPoints(pts)
    writer = vtk.vtkXMLMultiBlockDataWriter()
    writer.SetFileName("circular_domain_scaled.vtm")
    writer.SetInputData(_org)
    writer.Update()

# Wish i could save from VTK to CGNS
# But the we can load it in Paraview and re-export as CGNS

options = {

    # ---------------------------
    "inputFile": surfaceFile,
    "fileType": "PLOT3D",
    "unattachedEdgesAreSymmetry": False,
    "outerFaceBC": "farfield",
    "autoConnect": True,

...

but somehow, the file is still corrupted. Any idea where this comes from?

This is what my CGNS looks like in Paraview

I must point out that I am running it in a container and that the files are written on a CIFS volume. I will try to write them under '/tmp', maybe this could be help.

Thanks

[ews-ffarella]

OK, this solved it (assuming your are using the docker image)

This is taken from https://gitlab.onelab.info/gmsh/gmsh/-/blob/master/utils/docker/Dockerfile.ubuntu20.04

cd $HOME/packages
find $HOME/packages -maxdepth 1 -type d -iname "gmsh*" | xargs rm -rf

GMSH_VERSION="4.11.1"
GMSH_NAME=$(echo $GMSH_VERSION | sed 's/\./_/g')
wget "https://gitlab.onelab.info/gmsh/gmsh/-/archive/gmsh_$GMSH_NAME/gmsh-gmsh_$GMSH_NAME.tar.gz"  && tar -xf "gmsh-gmsh_$GMSH_NAME.tar.gz" && rm -f "gmsh-gmsh_$GMSH_NAME.tar.gz" && mv "gmsh-gmsh_$GMSH_NAME" "gmsh-$GMSH_VERSION"

# Or if we don't want the names being fixed
GMSH_VERSION="master"
git clone https://gitlab.onelab.info/gmsh/gmsh.git  "gmsh-$GMSH_VERSION" && cd "gmsh-$GMSH_VERSION"
git checkout 3a8640cbda19bbde95a80bdeef0525485d0f145e

cd "$HOME/packages/gmsh-$GMSH_VERSION"

mkdir build && cd build

CGNS_ROOT="$HOME/packages/CGNS-4.4.0/opt-gcc"  cmake \
    -DCMAKE_INSTALL_PREFIX="$HOME/packages/gmsh-$GMSH_VERSION/opt-gcc" \
    -DCMAKE_BUILD_TYPE=Release \
    -DENABLE_BUILD_SHARED=1 \
    -DENABLE_PRIVATE_API=1 \
    -DENABLE_CGNS=ON \
    -DENABLE_CGNS_CPEX0045=OFF \
    -DENABLE_PETSC=ON \
    -DENABLE_MPI=ON \
    -DENABLE_OPENMP=ON \
    -DENABLE_PETSC4PY=ON \
    -DOPENGL_GL_PREFERENCE=LEGACY \
    .. 

make -j8 shared && make -j8 install && cd .. && rm -rf build
cd $HOME

Now i can at least run cgnscheck on my output

Note that I reverted to version 4.11.1 because of this https://gitlab.onelab.info/gmsh/gmsh/-/issues/2652

[ews-ffarella]

Solved it! If anyone face problems with gmsh output, you can reach out to me. The trick is to use cgns_utilities, modify the BCSTANDARD dict and set the key of wall to 0 before using readGrid. Then, one can loop over the blocks and reset the cgns_type to the desired value