[petsc-users] DMLabel to extract height-0 points by their DMPolytope value

[Ferrand, Jesus A.]

Dear PETSc team:

I am trying to use DMPlex and DMLabel to develop an API to write plexes to .cgns format in parallel.
To that end, I need a way to extract the height-0 points and sort them by topological type (i.e., chunk of tetrahedra, followed by chunk of pyramids, etc.).
I figured I could use the DMLabel produced by DMPlexComputeCellTypes() as follows:

** I get the "celltype" DMLabel **

PetscBool has_tetrahedra, has_hexahedra, has_pyramids, has_tri_prisms;
PetscCall(DMLabelHasStratum(ctype_label, DM_POLYTOPE_TETRAHEDRON, &has_tetrahedra));
PetscCall(DMLabelHasStratum(ctype_label, DM_POLYTOPE_HEXAHEDRON, &has_hexahedra));
PetscCall(DMLabelHasStratum(ctype_label, DM_POLYTOPE_PYRAMID, &has_pyramids));
PetscCall(DMLabelHasStratum(ctype_label, DM_POLYTOPE_TRI_PRISM, &has_tri_prisms));

PetscInt nTopology = (PetscInt)has_tetrahedra + (PetscInt)has_hexahedra + (PetscInt)has_pyramids + (PetscInt)has_tri_prisms;
PetscInt *pType, *vType, *nType;
PetscMalloc3(nTopology, &pType, nTopology, &vType, nTopology, &nType);
PetscInt counter = -1;
if(has_tetrahedra){
      counter++;
      pType[counter] = DM_POLYTOPE_TETRAHEDRON;
      vType[counter] = 4;
      PetscCall(DMLabelGetStratumSize(ctype_label, DM_POLYTOPE_TETRAHEDRON, &nType[counter]));
}

** Repeat this pattern of if-statement for the rest. **

if(has_tri_prisms) {
      counter++;
      pType[counter] = DM_POLYTOPE_TRI_PRISM;
      vType[counter] = 6;
      PetscCall(DMLabelGetStratumSize(ctype_label, DM_POLYTOPE_TRI_PRISM, &nType[counter]));
}

IS pIS;
PetscInt StratumIdx;
const PetscInt* pPoints;
for(PetscInt ii = 0; ii < nTopology; ii++){
      PetscCall(DMLabelGetValueIndex(ctype_label, pType[ii], &StratumIdx));
      PetscCall(DMLabelGetStratumIS(ctype_label, StratumIdx, &pIS));
      PetscCall(ISGetIndices(pIS, &pPoints));
      for(PetscInt jj = 0; jj < nType[ii]; jj++){
            PetscCall(DMPlexGetTransitiveClosure(dm, pPoints[ii], PETSC_TRUE, &ClosureSize, &pClosure));

      **  Assemble connectivity array for each chunk of height-0 topology **

            PetscCall(DMPlexRestoreTransitiveClosure(dm, pPoints[ii], PETSC_TRUE, &ClosureSize, &pClosure));
      }
      PetscCall(ISRestoreIndices(pS, &pPoints));
      PetscCall(ISDestroy(&pIS));
}
PetscFree3(pType, vType, nType);

I think, that in principle, this is the correct approach for my immediate objective.
However, my problem is that the DAG points returned by pIS through pPoints are outside the height-0 stratum.
I can tell because I'm printing the contents of pPoints and comparing againts my DAG's height/depth strata.

It is as if the DMLabel has a different numbering from the DAG.
Also, for DMLabel's APIs, is the "stratum value" the same thing as the "label value"?
My gutt feeling is that the former is 0 <= StratumValue < nStrata, and the latter could be potentially disjoint (e.g., LabelValue in [-1, 0, 3 , 6, 7, 8]).
Is that what DMLabelGetValueIndex() is for?

Sincerely:

J.A. Ferrand

Embry-Riddle Aeronautical University - Daytona Beach - FL
Ph.D. Candidate, Aerospace Engineering

M.Sc. Aerospace Engineering

B.Sc. Aerospace Engineering

B.Sc. Computational Mathematics


Phone: (386)-843-1829

Email(s): ferranj2 at my.erau.edu

    jesus.ferrand at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230524/c6862abd/attachment-0001.html>