Re: [petsc-users] Issues creating DMPlex from higher order mesh generated by gmsh

[Jed Brown]

Matthew Knepley <knep...@gmail.com> writes:

> On Fri, May 5, 2023 at 10:55 AM Vilmer Dahlberg via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
>
>> Hi.
>>
>>
>> I'm trying to read a mesh of higher element order, in this example a mesh
>> consisting of 10-node tetrahedral elements, from gmsh, into PETSC. But It
>> looks like the mesh is not properly being loaded and converted into a
>> DMPlex. gmsh tells me it has generated a mesh with 7087 nodes, but when I
>> view my dm object it tells me it has 1081 0-cells. This is the printout I
>> get
>>
>
> Hi Vilmer,
>
> Plex makes a distinction between topological entities, like vertices, edges
> and cells, and the function spaces used to represent fields, like velocity
> or coordinates. When formats use "nodes", they mix the two concepts
> together.
>
> You see that if you add the number of vertices and edges, you get 7087,
> since for P2 there is a "node" on every edge. Is anything else wrong?

Note that quadratic (and higher order) tets are broken with the Gmsh reader. 
It's been on my todo list for a while.

As an example, this works when using linear elements (the projection makes them 
quadratic and visualization is correct), but is tangled when holes.msh is 
quadratic.

$ $PETSC_ARCH/tests/dm/impls/plex/tutorials/ex1 -dm_plex_filename 
~/meshes/holes.msh -dm_view cgns:s.cgns -dm_coord_petscspace_degree 2


SetFactory("OpenCASCADE");

radius = 0.3;
DefineConstant[
  nx = {1, Min 1, Max 30, Step 1, Name "Parameters/nx"}
  ny = {1, Min 1, Max 30, Step 1, Name "Parameters/ny"}
  extrude_length = {1, Min .1, Max 10, Step .1, Name "Parameters/extrusion 
length"}
  extrude_layers = {10, Min 1, Max 100, Step 1, Name "Parameters/extrusion 
layers"}
];
N = nx * ny;
Rectangle(1) = {0, 0, 0, 1, 1, 0};
Disk(10) = {0.5, 0.5, 0, radius};
BooleanDifference(100) = {Surface{1}; Delete;}{Surface{10}; Delete;};

For i In {0:nx-1}
  For j In {0:ny-1}
    If (i + j > 0)
       Translate {i, j, 0} { Duplicata { Surface{100}; } }
    EndIf
  EndFor
EndFor

Coherence;

// All the straight edges should have 8 elements
Transfinite Curve {:} = 8+1;

// Select the circles
circles = {};
For i In {0:nx-1}
  For j In {0:ny-1}
    lo = .5 - radius - 1e-4;
    hi = .5 + radius + 1e-4;
    circles() += Curve In BoundingBox {
      i + lo, j + lo, -1e-4,
      i + hi, j + hi, +1e-4
    };
  EndFor
EndFor

// Circles need 16 elements
Transfinite Curve {circles()} = 16+1;

Mesh.Algorithm = 8;

Extrude {0, 0, extrude_length} { Surface{100:100+N-1}; Layers{extrude_layers}; }

e = 1e-4;
extrude_start() = Surface In BoundingBox {-e, -e, -e, nx+e, ny+e, e};
extrude_end() = Surface In BoundingBox {
  -e, -e, extrude_length-e,
  nx+e, ny+e, extrude_length+e};

Physical Surface("start") = {extrude_start()};
Physical Surface("end") = {extrude_end()};
Physical Volume("solid") = {1:N};