On Tue, May 23, 2023 at 8:44 PM Ferrand, Jesus A. <[email protected]> wrote:

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 thought I already sorted them this way. Are you adding them in a different order?

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));

I do not understand why you need the index of this value.

PetscCall(DMLabelGetStratumIS(ctype_label, StratumIdx, &pIS));

This does not look right. You lookup by value, not by value index. Thanks, Matt

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):* [email protected]

[email protected]

-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>