1: const char help[] = "Test DMPlexInsertBoundaryValues with DMPlexSetClosurePermutationTensor.\n";

  3: #include <petscdmplex.h>

  5: static PetscErrorCode bc_func(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt num_comp_u, PetscScalar *u, PetscCtx ctx)
  6: {
  7:   PetscFunctionBeginUser;
  8:   for (PetscInt i = 0; i < num_comp_u; i++) u[i] = coords[i];
  9:   PetscFunctionReturn(PETSC_SUCCESS);
 10: }

 12: int main(int argc, char **argv)
 13: {
 14:   DM        dm;
 15:   PetscFE   fe;
 16:   Vec       U_loc;
 17:   PetscInt  dim, order = 1;
 18:   PetscBool tensorCoords = PETSC_TRUE;

 20:   /* Initialize PETSc */
 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 23:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-tensor_coords", &tensorCoords, NULL));
 24:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL));
 25:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 26:   PetscCall(DMSetType(dm, DMPLEX));
 27:   PetscCall(DMSetFromOptions(dm));

 29:   PetscCall(DMGetDimension(dm, &dim));
 30:   PetscCall(PetscFECreateLagrange(PETSC_COMM_WORLD, dim, dim, PETSC_FALSE, order, order, &fe));
 31:   PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
 32:   PetscCall(PetscFEDestroy(&fe));
 33:   PetscCall(DMCreateDS(dm));

 35:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

 37:   DMLabel  label;
 38:   PetscInt marker_ids[] = {1};
 39:   PetscCall(DMGetLabel(dm, "marker", &label));
 40:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", label, 1, marker_ids, 0, 0, NULL, (PetscVoidFn *)bc_func, NULL, NULL, NULL));
 41:   PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
 42:   {
 43:     DM cdm;
 44:     PetscCall(DMGetCoordinateDM(dm, &cdm));
 45:     if (tensorCoords) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
 46:   }

 48:   PetscCall(DMCreateLocalVector(dm, &U_loc));
 49:   {
 50:     /* The tensor permutation is stored per closure size, so a wrong degree leaves the cell closure silently unpermuted */
 51:     PetscSection section;
 52:     PetscScalar *vals = NULL;
 53:     PetscInt     cStart, clSize;
 54:     IS           perm;

 56:     PetscCall(DMGetLocalSection(dm, &section));
 57:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
 58:     PetscCall(DMPlexVecGetClosure(dm, section, U_loc, cStart, &clSize, &vals));
 59:     PetscCall(DMPlexVecRestoreClosure(dm, section, U_loc, cStart, &clSize, &vals));
 60:     PetscCall(PetscSectionGetClosurePermutation(section, (PetscObject)dm, dim, clSize, &perm));
 61:     PetscCall(ISDestroy(&perm));
 62:   }
 63:   PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, U_loc, 1., NULL, NULL, NULL));
 64:   PetscCall(VecViewFromOptions(U_loc, NULL, "-u_loc_vec_view"));
 65:   PetscCall(VecDestroy(&U_loc));
 66:   PetscCall(DMDestroy(&dm));
 67:   PetscCall(PetscFinalize());
 68:   return 0;
 69: }

 71: /*TEST
 72:   test:
 73:     suffix: 2d
 74:     args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_box_faces 3,3 -u_loc_vec_view
 75:   test:
 76:     suffix: 3d
 77:     args: -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -u_loc_vec_view
 78:   test:
 79:     # The tensor degree is recovered from the dual space size with a floating-point cube root, which overshoots in
 80:     # single precision for order >= 4
 81:     suffix: 3d_q4
 82:     output_file: output/empty.out
 83:     args: -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 1,1,1 -order 4 -tensor_coords 0
 84: TEST*/