Gmsh unexpected behaviour - axes directions swapped

[vachanpo...@gmail.com]

Dear all,

I have encountered a strange problem while dealing with physical boundary ids. I have a simple cube mesh (annexure 1) which I am exporting to msh4 format for reading through GridIn. The mesh is constructed by extruding a square 2D mesh in x-y plane to 3D.

I was facing problems with using collect_periodic_faces() on this mesh. So I mimicked the function by writing my own test function (annexure 2). It loops over faces with ids 2*direction and 2*direction+1 (with direction=2) which lie on a boundary and prints the boundary id.

According to the construction of the mesh (see geo file) face1 must have boundary id 4 and face 2 must have boundary id 5 in the test function (coincidentally, matching with the cell-local id). However, when I execute the test, the output shows both faces have a boundary id 0.

More interestingly, if I uncomment the lines in geo file and set boundary ids for all boundaries (in the trailing lines), then the code execution shows boundary ids of y-normal and z-normal faces swapped; i.e.; face1 and face2 show boundary id 3 which was assigned to y-normal faces (see the output below).

Found face1 of local id 4 at boundary with bid 3

Found face2 of local id 5 at boundary with bid 3

If I instead define face1 and face2 with direction=1 using

const auto face_1 = cell->face(2);

const auto face_2 = cell->face(3);

in the test function and modify the local id in print statements, then the boundary ids printed are 4 and 5 which were assigned to z-normal faces.

Found face1 of local id 2 at boundary with bid 4

Found face2 of local id 3 at boundary with bid 5

The x normal faces lying on boundary have the correct boundary id. So it looks like y and z axes are getting swapped?

Am I doing something wrong? Has anyone experienced this before? I am unable to find an explanation for this behaviour.

Thanking in anticipation

Vachan

Annexure 1

------------------------------------------------------------------------

// cube.geo file

SetFactory("OpenCASCADE");

// Simple cube geometry

lc=0.1;

l=1;

// 2D geometry

p1=newp; Point(p1) = {0,0,0,lc};

p2=newp; Point(p2) = {l,0,0,lc};

p3=newp; Point(p3) = {l,l,0,lc};

p4=newp; Point(p4) = {0,l,0,lc};

l1=newl; Line(l1) = {p1,p2};

l2=newl; Line(l2) = {p2,p3};

l3=newl; Line(l3) = {p3,p4};

l4=newl; Line(l4) = {p4,p1};

Transfinite Curve{l1,l2,l3,l4} = 10;

ll1=newll; Curve Loop(ll1) = {l1,l2,l3,l4};

s1=news; Plane Surface(s1) = {ll1};

Transfinite Surface{s1} = {p1,p2,p3,p4};

Recombine Surface{s1};

// 3D extrusion

out[] = Extrude {0,0,l} {

    Surface{s1}; Layers{10}; Recombine;

};

// Physical entities

Physical Volume(1) = {out[1]};

// Physical Surface(2) = {out[3],out[5]}; // x dir faces

// Physical Surface(3) = {out[2],out[4]}; // y dir faces

Physical Surface(4) = {s1}; // z- dir face

Physical Surface(5) = {out[0]}; // z+ dir face

Annexure 2

--------------------------------------------------------------

// snippet from test function mimicking collect_periodic_faces()

for(auto cell: problem.triang.active_cell_iterators()){

    const auto face_1 = cell->face(4);

    const auto face_2 = cell->face(5);

    if(face_1->at_boundary()){

        std::cout << "Found face1 of local id 4 at boundary with bid " << face_1->boundary_id()

        << "\n";

    }

    if(face_2->at_boundary()){

        std::cout << "Found face2 of local id 5 at boundary with bid " << face_2->boundary_id()

        << "\n";

    }

}

[Daniel Arndt]

Vachan,

It is difficult to say what goes wrong here with the information you provided. How does the msh file created by gmsh look like. Is the boundary information correct there?

Which gmsh version are you using? Can you post a complete minimal example of your code?

Best,

Daniel

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/97ddd5ab-18be-4ce2-8ecf-b0d77cd70f9dn%40googlegroups.com.

[vachanpo...@gmail.com]

Daniel,

Thank you for responding.

I am not familiar with the msh file format. But I see the 'Entities' section. Definitely, some info is being assigned, but I am not sure if it is correct.

I am attaching a minimal working example containing 4 files: main.cc, cube.geo, cube.msh and CMakeLists.txt. The executable (named a.out in cmake file) takes one command line argument for the direction. I am using gmsh 4.3.0.

The code reads 'cube.msh' and loops over all faces with local ids 2*direction and 2*direction+1. If any of them is on the boundary, it prints the boundary id. The geo file is supposed to assign boundary id 2 to x direction faces, 3 to y direction faces and 4,5 to z- (normal in negative z direction) and z+ faces.

The output for x faces is as expected, but the boundary ids assigned to y normal and z normal faces are somehow swapped.

Please let me know if any additional information is required.

[vachanpo...@gmail.com]

I have just read the gmsh manual and verified that the physical tags to surfaces are indeed getting assigned properly. The relevant lines are 25-31 in the msh file provided above.

I have arranged these lines in a table and verified that the physical tags are getting assigned exactly as I prescribed in the geo file. The picture shows the first 9 space-separated entries of lines 25-31 of cube.msh. I have inserted an additional column commenting about the plane's location.

According to the msh file, both x normal surfaces have tag 2, both y normal surfaces have tag 3 and z normal surfaces have tags 4 and 5 -- exactly as in the geo file.

So, I think the issue is with the procedure I am employing to read the mesh. I am attaching the main.cc (provided in the MWE above) here again, inline. I am unable to understand why the tags for y normal and z normal surfaces get swapped when I read and print them using this code. I used a similar code for 2D meshes and it worked fine. Are there any tricky things to consider additionally in 3D?

Vachan

----------------------------------------------------------------------------

#include <fstream>

#include <string>

#include <deal.II/base/mpi.h>

#include <deal.II/base/utilities.h>

#include <deal.II/distributed/tria.h>

#include <deal.II/base/conditional_ostream.h>

#include <deal.II/grid/grid_in.h>

using namespace dealii;

int main(int argc, char **argv)

{

    Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

    constexpr int dim = 3;

    MPI_Comm mpi_comm(MPI_COMM_WORLD);

    parallel::distributed::Triangulation<dim> triang(mpi_comm);

    GridIn<dim> grid_in;

    grid_in.attach_triangulation(triang);

    std::ifstream file("cube.msh");

    grid_in.read_msh(file);

    const int periodic_dir = std::stoi(argv[1]); // first command line argument

    for(auto cell: triang.active_cell_iterators()){

        const auto face_1 = cell->face(2*periodic_dir);

        const auto face_2 = cell->face(2*periodic_dir+1);

        if(face_1->at_boundary()){

            std::cout << "Found face1 of local id " << 2*periodic_dir << " at boundary with bid "

                << face_1->boundary_id() << "\n";

        }

        if(face_2->at_boundary()){

            std::cout << "Found face2 of local id " << 2*periodic_dir+1 << " at boundary with bid "

                << face_2->boundary_id() << "\n";

        }

    }

    return 0;

}

[vachanpo...@gmail.com]

I think I have found an explanation.

Since the mesh generated is a cartesian mesh, I was assuming that the faces of physical and reference cells have normals in the same direction and the mapping between then would be a simple stretching. But when I printed out vectors connecting face centers (0 to 1, 2 to 3 and 4 to 5), I found that this was not the case. The physical and reference cells have a rotation too included in the mapping.

That was the reason why version 1 of collect_periodic_faces() was working and version 2 was not: since the direction parameter is also used to access faces in the latter version. This would produce an exception because there is a rotation of axes between physical and reference space and the boundary id wouldn't match.

This has provoked another question. Version 2 of collect_periodic_faces() says 'This compatibility version of collect_periodic_faces() only works on grids with cells in standard orientation'. But is that sufficient? In my case, the mesh, after reading through GridIn, was in standard orientation. I had verified this by printing face orientation, flip and rotation. But since the boundary ids were set according to physical axes directions and these are different from reference axes directions, the output was unexpected to me.

Any comments on whatever I concluded so far would be greatly helpful to me!

[Daniel Arndt]

Vachan,

The documentation for the second version says

"Instead of defining a 'first' and 'second' boundary with the help of two boundary_ids this function defines a 'left' boundary as all faces with local face index 2*direction and boundary indicator b_id and, similarly, a 'right' boundary consisting of all face with local face index 2*direction+1 and boundary indicator b_id. Faces with coordinates only differing in the direction component are identified."

So there is an assumption about the local face index and the direction parameter that is also used for identifying faces. This implies that a case like yours where the cell is rotated is not supported by this overload. Just as you said, it's not enough for the cell to be in standard orientation.

I'm sorry that this caused so much trouble but I'm glad that the first version (that I would probably always recommend using anyway) works for you.

Best,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/c53204dc-3519-432f-9854-f2edfb43d361n%40googlegroups.com.

[vachanpo...@gmail.com]

Daniel,

As you rightly said, the documentation very clearly states how the direction parameter would be used. But the presence of a rotation in the mapping physical and reference cell was unexpected for me considering that the mesh was Cartesian. So the tricky part here is how the mesh is being read.. May be one small comment on this too would be helpful.

Thanks for the clarification and help!