Hi Anna,
I think you are looking at an output from a processor that does not have
the face set loaded.
If I load your mesh with a single MPI, I see:
DM Object: DM_0x84000000_0 1 MPI processes
type: plex
DM_0x84000000_0 in 3 dimensions:
Number of 0-cells per rank: 27499
Number of 1-cells per rank: 180922
Number of 2-cells per rank: 298443
Number of 3-cells per rank: 145019
Labels:
celltype: 4 strata with value/size (0 (27499), 6 (145019), 3 (298443), 1
(180922))
depth: 4 strata with value/size (0 (27499), 1 (180922), 2 (298443), 3
(145019))
Cell Sets: 1 strata with value/size (1 (145019))
Face Sets: 1 strata with value/size (101 (934))
Vertex Sets: 2 strata with value/size (101 (428), 1 (19092))
Which seems to match the gmsh script.
If I load your mesh with more MPIs, and print the DM via `DMView(dm,
PETSC_VIEWER_STDOUT_SELF);` then I can see that some MPI ranks do not have
it. I think this is expected behavior - only ranks that own the
elements attached to the face set will have the face set (or part of it)
defined. Similarly for vertex sets. At least that's how I understand it.
Are you after having the face sets and vertex sets on all ranks?
--
David
On Wed, Sep 20, 2023 at 6:17 AM Anna Dalklint via petsc-users <
petsc-users@mcs.anl.gov> wrote:
> 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
>
>
