t8code icon indicating copy to clipboard operation
t8code copied to clipboard

Feature request: support for tagged/labeled and periodic boundary information coming from gmsh files

Open jmark opened this issue 3 years ago • 5 comments

The gmsh file format allows tagged/labeled and periodic domain boundaries via so-called "physical entities". As of now, t8code discards such information and simply indicates a domain boundary via absence of any neighbor at the vertex/line/face.

jmark avatar Jan 17 '23 10:01 jmark

should be able to implement this via the cmesh attributes. Add new attribute to the list https://github.com/DLR-AMR/t8code/blob/46864dc8a1dd7ee84567fdd0412cdedfdb05d23a/src/t8_cmesh/t8_cmesh_types.h#L64

holke avatar Jan 17 '23 10:01 holke

The main problem I see sofar is, how to interpret the mesh file information. I post a minimal example of a square below with one tagged physical boundary "north" (tag no. = 9942) linked to curve 3. But I do not know how to properly link the information in "$entities" section to the "$elements" which is only read in by t8code (alongside "$nodes").

square.geo

lc = 1e-2;

Point(1) = {-1, -1, 0, lc};
Point(2) = { 1, -1, 0, lc};
Point(3) = { 1,  1, 0, lc};
Point(4) = {-1,  1, 0, lc};

Line(1) = {1, 2};
Line(2) = {3, 2};
Line(3) = {3, 4};
Line(4) = {4, 1};

Curve Loop(1) = {4, 1, -2, 3};

Plane Surface(1) = {1};

Transfinite Curve {:} = 2 Using Progression 1;

Transfinite Surface{:};
Recombine Surface{:};

Mesh 2;

Physical Curve("north", 9942) = {3};

// Mesh.MshFileVersion = 2.1;
Mesh.MshFileVersion = 4.1;
Mesh.SaveAll = 1;
Save "square.msh";

square.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
1
1 9942 "north"
$EndPhysicalNames
$Entities
4 4 1 0
1 -1 -1 0 0 
2 1 -1 0 0 
3 1 1 0 0 
4 -1 1 0 0 
1 -1 -1 0 1 -1 0 0 2 1 -2 
2 1 -1 0 1 1 0 0 2 3 -2 
3 -1 1 0 1 1 0 1 9942 2 3 -4 
4 -1 -1 0 -1 1 0 0 2 4 -1 
1 -1 -1 0 1 1 0 0 4 4 1 -2 3 
$EndEntities
$Nodes
9 4 1 4
0 1 0 1
1
-1 -1 0
0 2 0 1
2
1 -1 0
0 3 0 1
3
1 1 0
0 4 0 1
4
-1 1 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
2 1 0 0
$EndNodes
$Elements
9 9 1 9
0 1 15 1
1 1 
0 2 15 1
2 2 
0 3 15 1
3 3 
0 4 15 1
4 4 
1 1 1 1
5 1 2 
1 2 1 1
6 3 2 
1 3 1 1
7 3 4 
1 4 1 1
8 4 1 
2 1 3 1
9 4 1 2 3 
$EndElements

jmark avatar Jan 17 '23 11:01 jmark

See also this picture for better visualization. grafik

jmark avatar Jan 17 '23 11:01 jmark

ah, so the tags are linked with the geometries (curves, surfaces etc) and not the elements directly?

If so @sandro-elsweijer should be able to help. I think you can get the curve/surface etc. number from the element line in the file and could then read the tag. Or am i wrong Sandro?

holke avatar Jan 17 '23 11:01 holke

Excerpt from the msh 4 documentation:

$PhysicalNames // same as MSH version 2
  numPhysicalNames(ASCII int)
  dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
              //Here the physical name is connected to the physicalTag
  ...
$EndPhysicalNames

$Entities
  numPoints(size_t) numCurves(size_t) numSurfaces(size_t) numVolumes(size_t)
  pointTag(int) X(double) Y(double) Z(double)
    numPhysicalTags(size_t) physicalTag(int) ...
              //Here the physicalTags are connected to the pointTag
  ...
  curveTag(int) minX(double) minY(double) minZ(double)
    maxX(double) maxY(double) maxZ(double)
    numPhysicalTags(size_t) physicalTag(int) ...
              //Here the physicalTags are connected to the curveTag
    numBoundingPoints(size_t) pointTag(int; sign encodes orientation) ...
  ...
  surfaceTag(int) minX(double) minY(double) minZ(double)
    maxX(double) maxY(double) maxZ(double)
    numPhysicalTags(size_t) physicalTag(int) ...
              //Here the physicalTags are connected to the surfaceTag
    numBoundingCurves(size_t) curveTag(int; sign encodes orientation) ...
  ...
$EndEntities
...
$Nodes
  numEntityBlocks(size_t) numNodes(size_t) minNodeTag(size_t) maxNodeTag(size_t)
  entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesInBlock(size_t)
     //Here the point-/curve-/surfaceTags are connected to the nodes (entityDim : 0 --> pointTag; entityDim : 1 --> curveTag; ...)
    nodeTag(size_t)
    ...
    x(double) y(double) z(double)
       = 1) >
       = 2) >
       
    ...
  ...
$EndNodes

Though the physicalTags are connected to the entities, they are also implicitly connected to the nodes. We already read the entityDim and entityTag of a node and save it here: https://github.com/DLR-AMR/t8code/blob/441f5dceac2dca1a497c7dfc62d52fd02130166f/src/t8_cmesh/t8_cmesh_readmshfile.cxx#L1775

Each node is saved as a t8_msh_file_node_parametric_t, which holds all the information needed: https://github.com/DLR-AMR/t8code/blob/441f5dceac2dca1a497c7dfc62d52fd02130166f/src/t8_cmesh_readmshfile.h#L69:L77

But there are some problems coming with that: We only know the connected entities of the nodes. The entities of the edges and faces are not given directly. Therefore, to check if an edge has a physicalTag, we have to check, if the nodes of the edge have a physicalTag. In the case, that the edge has at least one node on the curve, we know, that the edge lies completely on that curve and that the edge has the same curveTag as the node. But if the edge covers the curve completely (the edge starts and ends together with the curve), both nodes of the edge lie on a point and hence, only the pointTags are given. The curveTag is not mentioned by any of the nodes. But this can be resolved by checking, if both points are on the curve. This information is given in the section $Entities. Therefore, we need to read the $PhysicalNames and $Entities section, to have all the information we need.

sandro-elsweijer avatar Jan 18 '23 09:01 sandro-elsweijer