[petsc-users] Face Set not appearing when loading GMSH mesh to DMPlex

[Anna Dalklint]

Hello,


I am trying to load a mesh, created via the gmsh api, into a DMPlex. The code is seen below:


  6     gmsh::initialize();
  7     elsize0 = 1.0;
  8     double a = 20.0;
  9     double b = 20.0;
 10     double c = 10.0;
 11     double d = 30.0;
 12     double e = 70.0;
 13
 14     gmsh::model::geo::addPoint(0.0, 0.0, 0.0, elsize0, 1);
 15     gmsh::model::geo::addPoint(a, 0.0, 0.0, elsize0, 2);
 16     gmsh::model::geo::addPoint(0.0, d, 0.0, elsize0, 3);
 17     gmsh::model::geo::addPoint(a, c, 0.0, elsize0, 4);
 18     gmsh::model::geo::addPoint(e, d, 0.0, elsize0, 5);
 19     gmsh::model::geo::addPoint(e, c, 0.0, elsize0, 6);
 20     gmsh::model::geo::addPoint(a, 0.0, b, elsize0, 7);
 21     gmsh::model::geo::addPoint(0.0, 0.0, b, elsize0, 8);
 22     gmsh::model::geo::addPoint(0.0, d, b, elsize0, 9);
 23     gmsh::model::geo::addPoint(e, d, b, elsize0, 10);
 24     gmsh::model::geo::addPoint(e, c, b, elsize0, 11);
 25     gmsh::model::geo::addPoint(a, c, b, elsize0, 12);
 26
 27     // Define lines connecting points
 28     gmsh::model::geo::addLine(2, 1, 1);
 29     gmsh::model::geo::addLine(1, 3, 2);
 30     gmsh::model::geo::addLine(3, 5, 3);
 31     gmsh::model::geo::addLine(5, 6, 4);
 32     gmsh::model::geo::addLine(6, 4, 5);
 33     gmsh::model::geo::addLine(4, 2, 6);
 34     gmsh::model::geo::addLine(2, 7, 7);
 35     gmsh::model::geo::addLine(7, 8, 8);
 36     gmsh::model::geo::addLine(8, 1, 9);
 37     gmsh::model::geo::addLine(8, 9, 10);
 38     gmsh::model::geo::addLine(9, 3, 11);
 39     gmsh::model::geo::addLine(9, 10, 12);
 40     gmsh::model::geo::addLine(10, 11, 13);
 41     gmsh::model::geo::addLine(11, 12, 14);
 42     gmsh::model::geo::addLine(12, 7, 15);
 43     gmsh::model::geo::addLine(12, 4, 16);
 44     gmsh::model::geo::addLine(11, 6, 17);
 45     gmsh::model::geo::addLine(10, 5, 18);
 46
 47     // Define surfaces through their boundaries
 48     gmsh::model::geo::addCurveLoop({3, 4, 5, 6, 1, 2}, 1);
 49     gmsh::model::geo::addPlaneSurface({1}, 1);
 50     gmsh::model::geo::addCurveLoop({18, 4, -17, -13}, 2);
 51     gmsh::model::geo::addPlaneSurface({2}, 2);
 52     gmsh::model::geo::addCurveLoop({5, -16, -14, 17}, 3);
 53     gmsh::model::geo::addPlaneSurface({3}, 3);
 54     gmsh::model::geo::addCurveLoop({16, 6, 7, -15}, 4);
 55     gmsh::model::geo::addPlaneSurface({4}, 4);
 56     gmsh::model::geo::addCurveLoop({12, 13, 14, 15, 8, 10}, 5);
 57     gmsh::model::geo::addPlaneSurface({5}, 5);
 58     gmsh::model::geo::addCurveLoop({11, -2, -9, 10}, 6);
 59     gmsh::model::geo::addPlaneSurface({6}, 6);
 60     gmsh::model::geo::addCurveLoop({1, -9, -8, -7}, 7);
 61     gmsh::model::geo::addPlaneSurface({7}, 7);
 62     gmsh::model::geo::addCurveLoop({11, 3, -18, -12}, 8);
 63     gmsh::model::geo::addPlaneSurface({8}, 8);
 64
 65     gmsh::model::geo::addSurfaceLoop({8, 1, 2, 3, 4, 7, 6, 5}, 1);
 66     gmsh::model::geo::addVolume({1}, 1);
 67
 68     gmsh::model::geo::synchronize();
 70
 72     gmsh::model::addPhysicalGroup(3, {1}, 1);
 74     gmsh::model::addPhysicalGroup(2, {7}, 101, "dir_y");

  98     gmsh::option::setNumber("Mesh.Algorithm", 6);
  99     gmsh::option::setNumber("Mesh.ElementOrder", 1);
  100   gmsh::model::mesh::generate(3);
  101   gmsh::write("filename.msh");


The mesh is generated, but the Physical Group is not correctly loaded into the DMPlex. If I run a DMView on my dm, I get:


171 DM Object: Parallel Mesh 4 MPI processes
172   type: plex
173 Parallel Mesh in 3 dimensions:
174   Number of 0-cells per rank: 7138 7192 7328 7171
175   Number of 1-cells per rank: 45825 45992 46817 45824
176   Number of 2-cells per rank: 74791 74901 76203 74568
177   Number of 3-cells per rank: 36103 36100 36713 35914
178 Labels:
179   depth: 4 strata with value/size (0 (7138), 1 (45825), 2 (74791), 3 (36103))
180   celltype: 4 strata with value/size (0 (7138), 1 (45825), 3 (74791), 6 (36103))
181   Cell Sets: 1 strata with value/size (1 (36103))
182   Face Sets: 0 strata with value/size ()
183   Vertex Sets: 1 strata with value/size (1 (5200))


Why are the elements of the physicalGroup "dir_y" (Plane Surface 7) not loaded as a face set? If I change the Plane Surface to another surface, e.g. 8, it marks the elements as the correct face set. Another weird thing is that the number of vertices in the above Vertex set is not correct.


I have tried changing the numbering of Plane Surface 7 (to e.g. 70) but the issue remains.


I run the code with the -dm_plex_gmsh_mark_vertices and -dm_plex_gmsh_multiple_tags flags.


Thank you,

Anna

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230920/73a4f65b/attachment.html>