Actual source code: ex49.c

  1: static char help[] = "Tests dof numberings for external integrators such as LibCEED.\n\n";

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

  6: typedef struct {
  7:   PetscInt dummy;
  8: } AppCtx;

 10: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 11: {
 13:   options->dummy = 1;
 14:   PetscOptionsBegin(comm, "", "Dof Ordering Options", "DMPLEX");
 15:   //PetscOptionsInt("-num_refine", "Refine cycles", "ex46.c", options->Nr, &options->Nr, NULL);
 16:   PetscOptionsEnd();
 17:   return 0;
 18: }

 20: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 21: {
 23:   DMCreate(comm, dm);
 24:   DMSetType(*dm, DMPLEX);
 25:   DMSetFromOptions(*dm);
 26:   DMSetApplicationContext(*dm, user);
 27:   DMViewFromOptions(*dm, NULL, "-dm_view");
 28:   return 0;
 29: }

 31: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
 32: {
 33:   DM       cdm = dm;
 34:   PetscFE  fe;
 35:   PetscInt dim;

 38:   DMGetDimension(dm, &dim);
 39:   PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe);
 40:   PetscObjectSetName((PetscObject)fe, "scalar");
 41:   DMSetField(dm, 0, NULL, (PetscObject)fe);
 42:   DMSetField(dm, 1, NULL, (PetscObject)fe);
 43:   PetscFEDestroy(&fe);
 44:   DMCreateDS(dm);
 45:   while (cdm) {
 46:     DMCopyDisc(dm, cdm);
 47:     DMGetCoarseDM(cdm, &cdm);
 48:   }
 49:   return 0;
 50: }

 52: static PetscErrorCode CheckOffsets(DM dm, const char *domain_name, PetscInt label_value, PetscInt height)
 53: {
 54:   const char *height_name[] = {"cells", "faces"};
 55:   DMLabel     domain_label  = NULL;
 56:   DM          cdm;
 57:   IS          offIS;
 58:   PetscInt   *offsets, Ncell, Ncl, Nc, n;
 59:   PetscInt    Nf, f;

 62:   if (domain_name) DMGetLabel(dm, domain_name, &domain_label);
 63:   PetscPrintf(PETSC_COMM_WORLD, "## %s: '%s' {%" PetscInt_FMT "}\n", height_name[height], domain_name ? domain_name : "default", label_value);
 64:   // Offsets for cell closures
 65:   DMGetNumFields(dm, &Nf);
 66:   for (f = 0; f < Nf; ++f) {
 67:     char name[PETSC_MAX_PATH_LEN];

 69:     DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets);
 70:     ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS);
 71:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f);
 72:     PetscObjectSetName((PetscObject)offIS, name);
 73:     ISViewFromOptions(offIS, NULL, "-offsets_view");
 74:     ISDestroy(&offIS);
 75:   }
 76:   // Offsets for coordinates
 77:   {
 78:     Vec                X;
 79:     PetscSection       s;
 80:     const PetscScalar *x;
 81:     const char        *cname;
 82:     PetscInt           cdim;
 83:     PetscBool          isDG = PETSC_FALSE;

 85:     DMGetCellCoordinateDM(dm, &cdm);
 86:     if (!cdm) {
 87:       DMGetCoordinateDM(dm, &cdm);
 88:       cname = "Coordinates";
 89:       DMGetCoordinatesLocal(dm, &X);
 90:     } else {
 91:       isDG  = PETSC_TRUE;
 92:       cname = "DG Coordinates";
 93:       DMGetCellCoordinatesLocal(dm, &X);
 94:     }
 95:     if (isDG && height) return 0;
 96:     if (domain_name) DMGetLabel(cdm, domain_name, &domain_label);
 97:     DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets);
 98:     DMGetCoordinateDim(dm, &cdim);
100:     DMGetLocalSection(cdm, &s);
101:     VecGetArrayRead(X, &x);
102:     PetscPrintf(PETSC_COMM_WORLD, "%s by element\n", cname);
103:     for (PetscInt c = 0; c < Ncell; ++c) {
104:       for (PetscInt v = 0; v < Ncl; ++v) {
105:         PetscInt           off = offsets[c * Ncl + v], dgdof;
106:         const PetscScalar *vx  = &x[off];

108:         if (isDG) {
109:           PetscSectionGetDof(s, c, &dgdof);
111:         }
112:         switch (cdim) {
113:         case 2:
114:           PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1]));
115:           break;
116:         case 3:
117:           PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1]), (double)PetscRealPart(vx[2]));
118:           break;
119:         }
120:       }
121:     }
122:     VecRestoreArrayRead(X, &x);
123:     PetscFree(offsets);
124:   }
125:   return 0;
126: }

128: int main(int argc, char **argv)
129: {
130:   DM     dm;
131:   AppCtx user;

134:   PetscInitialize(&argc, &argv, NULL, help);
135:   ProcessOptions(PETSC_COMM_WORLD, &user);
136:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
137:   SetupDiscretization(dm, &user);
138:   CheckOffsets(dm, NULL, 0, 0);
139:   CheckOffsets(dm, "Face Sets", 1, 1);
140:   DMDestroy(&dm);
141:   PetscFinalize();
142:   return 0;
143: }

145: /*TEST

147:   test:
148:     suffix: 0
149:     requires: triangle
150:     args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view

152:   test:
153:     suffix: 1
154:     args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \
155:           -dm_view -offsets_view

157:   test:
158:     suffix: cg_2d
159:     args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \
160:           -dm_view -offsets_view
161: TEST*/