1: static char help[] = "Tests DMPlexTransformCreateSplitCellLabel().\n\n";

  3: #include <petscdmplex.h>
  4: #include <petscdmplextransform.h>

  6: // Flags a single cell for refinement, which is what an error estimator would do
  7: static PetscErrorCode CreateAdaptLabel(DM dm, DMLabel *adaptLabel)
  8: {
  9:   PetscInt cStart, cEnd;

 11:   PetscFunctionBeginUser;
 12:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 13:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", adaptLabel));
 14:   PetscCall(DMLabelSetDefaultValue(*adaptLabel, DM_ADAPT_KEEP));
 15:   if (cEnd > cStart) PetscCall(DMLabelSetValue(*adaptLabel, cStart, DM_ADAPT_REFINE));
 16:   PetscFunctionReturn(PETSC_SUCCESS);
 17: }

 19: int main(int argc, char **argv)
 20: {
 21:   DM              dm, rdm;
 22:   DMPlexTransform tr;
 23:   DMLabel         adaptLabel, splitLabel;
 24:   PetscInt        cStart, cEnd, rcStart, rcEnd, numSplit;

 26:   PetscFunctionBeginUser;
 27:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 28:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 29:   PetscCall(DMSetType(dm, DMPLEX));
 30:   PetscCall(DMSetFromOptions(dm));
 31:   // Without this the refined mesh does not keep the transformation that produced it
 32:   PetscCall(DMPlexSetSaveTransform(dm, PETSC_TRUE));
 33:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));

 35:   PetscCall(CreateAdaptLabel(dm, &adaptLabel));
 36:   PetscCall(DMAdaptLabel(dm, adaptLabel, &rdm));
 37:   PetscCheck(rdm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Adaptation did not produce a mesh");
 38:   PetscCall(DMPlexGetHeightStratum(rdm, 0, &rcStart, &rcEnd));
 39:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coarse cells: %" PetscInt_FMT " Refined cells: %" PetscInt_FMT "\n", cEnd - cStart, rcEnd - rcStart));

 41:   PetscCall(DMPlexGetTransform(rdm, &tr));
 42:   PetscCheck(tr, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Refined mesh did not keep its transform");
 43:   PetscCall(DMPlexTransformCreateSplitCellLabel(tr, rdm, &splitLabel));
 44:   PetscCall(DMLabelGetStratumSize(splitLabel, 1, &numSplit));
 45:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Cells from a split parent: %" PetscInt_FMT "\n", numSplit));
 46:   // Every cell of a genuinely refined parent is marked, and cells passed through unchanged are not
 47:   PetscCheck(numSplit > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "No cell came from a split parent, but the mesh grew");
 48:   PetscCheck(numSplit <= rcEnd - rcStart, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Marked %" PetscInt_FMT " cells, but the mesh only has %" PetscInt_FMT, numSplit, rcEnd - rcStart);

 50:   // The closure of the split cells is the set of entities whose star contains one, which is what PCPATCH wants
 51:   PetscCall(DMPlexLabelComplete(rdm, splitLabel));
 52:   PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD));

 54:   PetscCall(DMLabelDestroy(&splitLabel));
 55:   PetscCall(DMLabelDestroy(&adaptLabel));
 56:   PetscCall(DMDestroy(&rdm));
 57:   PetscCall(DMDestroy(&dm));
 58:   PetscCall(PetscFinalize());
 59:   return 0;
 60: }

 62: /*TEST

 64:   # A doublet is two tetrahedra, so SBR splits the flagged cell and its neighbor to stay conforming
 65:   test:
 66:     suffix: sbr
 67:     args: -dm_adaptor cellrefiner -dm_plex_shape doublet -dm_plex_dim 3 -dm_plex_simplex 1 -dm_plex_transform_type refine_sbr

 69:   test:
 70:     suffix: regular
 71:     args: -dm_adaptor cellrefiner -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -dm_plex_transform_type refine_regular

 73: TEST*/