Hi again,
David was right, the DMView call only printed the sets on the master rank
(which I wasn't aware of). So there was no issue!
Thank you,
Anna
________________________________
From: Blaise Bourdin <bour...@mcmaster.ca>
Sent: Thursday, September 21, 2023 5:05:35 PM
To: Anna Dalklint <anna.dalkl...@solid.lth.se>
Cc: petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov>
Subject: Re: [petsc-users] Face Set not appearing when loading GMSH mesh to
DMPlex
Can you send me the geo or msh file, or instruction on how to build your
example? I am not familiar with the gmsh API
Blaise
On Sep 20, 2023, at 8:16 AM, Anna Dalklint via petsc-users
<petsc-users@mcs.anl.gov> wrote:
Caution: External email.
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
—
Canada Research Chair in Mathematical and Computational Aspects of Solid
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243
