<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"><!-- P {margin-top:0;margin-bottom:0;} --></style>
</head>
<body dir="ltr">
<div id="divtagdefaultwrapper" style="font-size:12pt;color:#000000;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>
</body>
</html>