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

Thu Sep 21 23:39:49 CDT 2023

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 <bourdin at mcmaster.ca>
Sent: Thursday, September 21, 2023 5:05:35 PM
To: Anna Dalklint <anna.dalklint at solid.lth.se>
Cc: petsc-users at mcs.anl.gov <petsc-users at 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 at 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

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