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

[Matthew Knepley]

On Tue, May 23, 2023 at 8:44 PM Ferrand, Jesus A. <FERRANJ2 at 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):* ferranj2 at my.erau.edu
>
>     jesus.ferrand at 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230524/9cac4d26/attachment-0001.html>