<div dir="ltr">Hi Anna,<div><br></div><div>I think you are looking at an output from a processor that does not have the face set loaded.</div><div><br></div><div>If I load your mesh with a single MPI, I see:</div><div><br></div><div><font face="monospace">DM Object: DM_0x84000000_0 1 MPI processes<br> type: plex<br>DM_0x84000000_0 in 3 dimensions:<br> Number of 0-cells per rank: 27499<br> Number of 1-cells per rank: 180922<br> Number of 2-cells per rank: 298443<br> Number of 3-cells per rank: 145019<br>Labels:<br> celltype: 4 strata with value/size (0 (27499), 6 (145019), 3 (298443), 1 (180922))<br> depth: 4 strata with value/size (0 (27499), 1 (180922), 2 (298443), 3 (145019))<br> Cell Sets: 1 strata with value/size (1 (145019))<br> Face Sets: 1 strata with value/size (101 (934))<br> Vertex Sets: 2 strata with value/size (101 (428), 1 (19092))<br></font></div><div><br></div><div>Which seems to match the gmsh script.</div><div><br></div><div>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.</div><div><br></div><div>Are you after having the face sets and vertex sets on all ranks?</div><div><br></div><div>--</div><div>David</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Sep 20, 2023 at 6:17 AM Anna Dalklint via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg2074183612224069869">
<div dir="ltr">
<div id="m_2074183612224069869divtagdefaultwrapper" style="font-size:12pt;color:rgb(0,0,0);font-family:Calibri,Helvetica,sans-serif" dir="ltr">
<p></p>
<div>
<p>Hello,</p>
<p><br>
</p>
<p>I am trying to load a mesh, created via the gmsh api, into a DMPlex. The code is seen below:<br>
<br>
</p>
<div> 6 gmsh::initialize(); <br>
</div>
<div> 7 elsize0 = 1.0;<br>
8 double a = 20.0; <br>
9 double b = 20.0; <br>
10 double c = 10.0; <br>
11 double d = 30.0; <br>
12 double e = 70.0; <br>
13 <br>
14 gmsh::model::geo::addPoint(0.0, 0.0, 0.0, elsize0, 1);<br>
15 gmsh::model::geo::addPoint(a, 0.0, 0.0, elsize0, 2);<br>
16 gmsh::model::geo::addPoint(0.0, d, 0.0, elsize0, 3);<br>
17 gmsh::model::geo::addPoint(a, c, 0.0, elsize0, 4);<br>
18 gmsh::model::geo::addPoint(e, d, 0.0, elsize0, 5);<br>
19 gmsh::model::geo::addPoint(e, c, 0.0, elsize0, 6);<br>
20 gmsh::model::geo::addPoint(a, 0.0, b, elsize0, 7);<br>
21 gmsh::model::geo::addPoint(0.0, 0.0, b, elsize0, 8);<br>
22 gmsh::model::geo::addPoint(0.0, d, b, elsize0, 9);<br>
23 gmsh::model::geo::addPoint(e, d, b, elsize0, 10);<br>
24 gmsh::model::geo::addPoint(e, c, b, elsize0, 11);<br>
25 gmsh::model::geo::addPoint(a, c, b, elsize0, 12);<br>
26 <br>
27 // Define lines connecting points<br>
28 gmsh::model::geo::addLine(2, 1, 1);<br>
29 gmsh::model::geo::addLine(1, 3, 2);<br>
30 gmsh::model::geo::addLine(3, 5, 3);<br>
31 gmsh::model::geo::addLine(5, 6, 4);<br>
32 gmsh::model::geo::addLine(6, 4, 5);<br>
33 gmsh::model::geo::addLine(4, 2, 6);<br>
34 gmsh::model::geo::addLine(2, 7, 7);<br>
35 gmsh::model::geo::addLine(7, 8, 8);<br>
36 gmsh::model::geo::addLine(8, 1, 9);<br>
37 gmsh::model::geo::addLine(8, 9, 10);<br>
38 gmsh::model::geo::addLine(9, 3, 11);<br>
39 gmsh::model::geo::addLine(9, 10, 12);<br>
40 gmsh::model::geo::addLine(10, 11, 13);<br>
41 gmsh::model::geo::addLine(11, 12, 14);<br>
42 gmsh::model::geo::addLine(12, 7, 15);<br>
43 gmsh::model::geo::addLine(12, 4, 16);<br>
44 gmsh::model::geo::addLine(11, 6, 17);<br>
45 gmsh::model::geo::addLine(10, 5, 18);<br>
46 <br>
47 // Define surfaces through their boundaries<br>
48 gmsh::model::geo::addCurveLoop({3, 4, 5, 6, 1, 2}, 1);<br>
49 gmsh::model::geo::addPlaneSurface({1}, 1);<br>
50 gmsh::model::geo::addCurveLoop({18, 4, -17, -13}, 2);<br>
51 gmsh::model::geo::addPlaneSurface({2}, 2);<br>
52 gmsh::model::geo::addCurveLoop({5, -16, -14, 17}, 3);<br>
53 gmsh::model::geo::addPlaneSurface({3}, 3);<br>
54 gmsh::model::geo::addCurveLoop({16, 6, 7, -15}, 4);<br>
55 gmsh::model::geo::addPlaneSurface({4}, 4);<br>
56 gmsh::model::geo::addCurveLoop({12, 13, 14, 15, 8, 10}, 5);<br>
57 gmsh::model::geo::addPlaneSurface({5}, 5);<br>
58 gmsh::model::geo::addCurveLoop({11, -2, -9, 10}, 6);<br>
59 gmsh::model::geo::addPlaneSurface({6}, 6);<br>
60 gmsh::model::geo::addCurveLoop({1, -9, -8, -7}, 7);<br>
61 gmsh::model::geo::addPlaneSurface({7}, 7);</div>
<div> 62 gmsh::model::geo::addCurveLoop({11, 3, -18, -12}, 8); <br>
63 gmsh::model::geo::addPlaneSurface({8}, 8); <br>
64 <br>
65 gmsh::model::geo::addSurfaceLoop({8, 1, 2, 3, 4, 7, 6, 5}, 1);<br>
66 gmsh::model::geo::addVolume({1}, 1);<br>
67 <br>
68 gmsh::model::geo::synchronize();<br>
70 <br>
72 gmsh::model::addPhysicalGroup(3, {1}, 1);<br>
74 gmsh::model::addPhysicalGroup(2, {7}, 101, "dir_y");</div>
<div><br>
</div>
<div><span> 98 gmsh::option::setNumber("Mesh.Algorithm", 6);</span><br>
<span> 99 gmsh::option::setNumber("Mesh.ElementOrder", 1);</span></div>
<div><span><span> 100 gmsh::model::mesh::generate(3);</span><br>
</span></div>
<div><span> 101 <span>gmsh::write("filename.msh");</span><br>
</span></div>
<p></p>
<p><br>
</p>
<p>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:</p>
<p><br>
</p>
<p></p>
<div>171 DM Object: Parallel Mesh 4 MPI processes<br>
172 type: plex<br>
173 Parallel Mesh in 3 dimensions:<br>
174 Number of 0-cells per rank: 7138 7192 7328 7171<br>
175 Number of 1-cells per rank: 45825 45992 46817 45824<br>
176 Number of 2-cells per rank: 74791 74901 76203 74568<br>
177 Number of 3-cells per rank: 36103 36100 36713 35914<br>
178 Labels:<br>
179 depth: 4 strata with value/size (0 (7138), 1 (45825), 2 (74791), 3 (36103))<br>
180 celltype: 4 strata with value/size (0 (7138), 1 (45825), 3 (74791), 6 (36103))<br>
181 Cell Sets: 1 strata with value/size (1 (36103))<br>
182 Face Sets: 0 strata with value/size ()<br>
183 Vertex Sets: 1 strata with value/size (1 (5200))</div>
<p></p>
<p><br>
</p>
<p>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. <br>
</p>
<p><br>
</p>
<p>I have tried changing the numbering of Plane Surface 7 (to e.g. 70) but the issue remains.<br>
</p>
<p><br>
</p>
<p>I run the code with the <span>-dm_plex_gmsh_mark_vertices and -dm_plex_gmsh_multiple_tags flags.</span>
<br>
</p>
<p><br>
</p>
<p>Thank you,</p>
Anna</div>
<br>
<p></p>
</div>
</div>
</div></blockquote></div>