On Tue, May 23, 2023 at 8:44 PM Ferrand, Jesus A. <ferra...@my.erau.edu> 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):* ferra...@my.erau.edu > > jesus.ferr...@gmail.com > -- 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/>