Actual source code: ex8.c

  1: static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n";

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

  6: static PetscErrorCode ViewOffsets(DM dm, Vec X)
  7: {
  8:   PetscInt           num_elem, elem_size, num_comp, num_dof;
  9:   PetscInt          *elem_restr_offsets;
 10:   const PetscScalar *x = NULL;
 11:   const char        *name;

 13:   PetscFunctionBegin;
 14:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
 15:   PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets));
 16:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "DM %s offsets: num_elem %" PetscInt_FMT ", size %" PetscInt_FMT ", comp %" PetscInt_FMT ", dof %" PetscInt_FMT "\n", name, num_elem, elem_size, num_comp, num_dof));
 17:   if (X) PetscCall(VecGetArrayRead(X, &x));
 18:   for (PetscInt c = 0; c < num_elem; c++) {
 19:     PetscCall(PetscIntView(elem_size, &elem_restr_offsets[c * elem_size], PETSC_VIEWER_STDOUT_SELF));
 20:     if (x) {
 21:       for (PetscInt i = 0; i < elem_size; i++) PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c * elem_size + i]], PETSC_VIEWER_STDERR_SELF));
 22:     }
 23:   }
 24:   if (X) PetscCall(VecRestoreArrayRead(X, &x));
 25:   PetscCall(PetscFree(elem_restr_offsets));
 26:   PetscFunctionReturn(PETSC_SUCCESS);
 27: }

 29: int main(int argc, char **argv)
 30: {
 31:   DM           dm;
 32:   PetscSection section;
 33:   PetscFE      fe;
 34:   PetscInt     dim, Nf = 1, c, cStart, cEnd;
 35:   PetscBool    view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE;

 37:   PetscFunctionBeginUser;
 38:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 39:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX");
 40:   PetscCall(PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL));
 41:   PetscCall(PetscOptionsBool("-project_coordinates", "Call DMSetCoordinateDisc() explicitly", "ex8.c", project, &project, NULL));
 42:   PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL));
 43:   PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields to use", "ex8.c", Nf, &Nf, NULL, 1));
 44:   PetscOptionsEnd();

 46:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 47:   PetscCall(DMSetType(dm, DMPLEX));
 48:   PetscCall(DMSetFromOptions(dm));
 49:   if (project) {
 50:     PetscFE  fe_coords;
 51:     PetscInt cdim;
 52:     PetscCall(DMGetCoordinateDim(dm, &cdim));
 53:     PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords));
 54:     PetscCall(DMSetCoordinateDisc(dm, fe_coords, PETSC_TRUE));
 55:     PetscCall(PetscFEDestroy(&fe_coords));
 56:   }
 57:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
 58:   PetscCall(DMGetDimension(dm, &dim));

 60:   if (Nf == 1) {
 61:     PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe));
 62:     PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
 63:     PetscCall(PetscFEDestroy(&fe));
 64:   } else {
 65:     for (PetscInt f = 0; f < Nf; ++f) {
 66:       char prefix[16];
 67:       PetscCall(PetscSNPrintf(prefix, 16, "f%" PetscInt_FMT "_", f));
 68:       PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, prefix, PETSC_DETERMINE, &fe));
 69:       PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
 70:       PetscCall(PetscFEDestroy(&fe));
 71:     }
 72:   }
 73:   PetscCall(DMCreateDS(dm));
 74:   if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
 75:   PetscCall(DMGetLocalSection(dm, &section));
 76:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 77:   for (c = cStart; c < cEnd; c++) {
 78:     PetscInt numindices, *indices;
 79:     PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
 80:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "\n", c - cStart));
 81:     PetscCall(PetscIntView(numindices, indices, PETSC_VIEWER_STDOUT_SELF));
 82:     PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL));
 83:   }
 84:   {
 85:     DMLabel         domain_label;
 86:     IS              valueIS;
 87:     const PetscInt *values;
 88:     PetscInt        Nv;

 90:     PetscCall(DMGetLabel(dm, "Face Sets", &domain_label));
 91:     PetscCall(DMLabelGetNumValues(domain_label, &Nv));
 92:     // Check that discretization bas any values on faces
 93:     {
 94:       PetscDS         ds;
 95:       PetscFE         fe;
 96:       IS              fieldIS;
 97:       const PetscInt *fields;
 98:       PetscInt        Nf, dmf = 0, dsf = -1;

100:       PetscCall(DMGetRegionDS(dm, domain_label, &fieldIS, &ds, NULL));
101:       // Translate DM field number to DS field number
102:       PetscCall(ISGetIndices(fieldIS, &fields));
103:       PetscCall(ISGetSize(fieldIS, &Nf));
104:       for (PetscInt f = 0; f < Nf; ++f) {
105:         if (dmf == fields[f]) {
106:           dsf = f;
107:           break;
108:         }
109:       }
110:       PetscCall(ISRestoreIndices(fieldIS, &fields));
111:       PetscCall(PetscDSGetDiscretization(ds, dsf, (PetscObject *)&fe));
112:       PetscCall(PetscFEGetHeightSubspace(fe, 1, &fe));
113:       if (!fe) Nv = 0;
114:     }
115:     PetscCall(DMLabelGetValueIS(domain_label, &valueIS));
116:     PetscCall(ISGetIndices(valueIS, &values));
117:     for (PetscInt v = 0; v < Nv; ++v) {
118:       PetscInt *elem_restr_offsets_face;
119:       PetscInt  num_elem, elem_size, num_comp, num_dof;

121:       PetscCall(DMPlexGetLocalOffsets(dm, domain_label, values[v], 1, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets_face));
122:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "========= Face restriction; marker: %" PetscInt_FMT " ========\n", values[v]));
123:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "num_elem: %" PetscInt_FMT ", elem_size: %" PetscInt_FMT ", num_dof:  %" PetscInt_FMT "\n", num_elem, elem_size, num_dof));
124:       for (PetscInt c = 0; c < num_elem; c++) PetscCall(PetscIntView(elem_size, &elem_restr_offsets_face[c * elem_size], PETSC_VIEWER_STDOUT_SELF));
125:       PetscCall(PetscFree(elem_restr_offsets_face));
126:     }
127:     PetscCall(ISRestoreIndices(valueIS, &values));
128:     PetscCall(ISDestroy(&valueIS));
129:   }
130:   if (view_coord) {
131:     DM       cdm;
132:     Vec      X;
133:     PetscInt cdim;

135:     PetscCall(DMGetCoordinatesLocalSetUp(dm));
136:     PetscCall(DMGetCoordinateDim(dm, &cdim));
137:     PetscCall(DMGetCoordinateDM(dm, &cdm));
138:     PetscCall(PetscObjectSetName((PetscObject)cdm, "coords"));
139:     if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
140:     for (c = cStart; c < cEnd; ++c) {
141:       const PetscScalar *array;
142:       PetscScalar       *x = NULL;
143:       PetscInt           ndof;
144:       PetscBool          isDG;

146:       PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
147:       PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim");
148:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart));
149:       for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF));
150:       PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
151:     }
152:     PetscCall(ViewOffsets(dm, NULL));
153:     PetscCall(DMGetCoordinatesLocal(dm, &X));
154:     PetscCall(ViewOffsets(cdm, X));
155:   }
156:   PetscCall(DMDestroy(&dm));
157:   PetscCall(PetscFinalize());
158:   return 0;
159: }

161: /*TEST

163:   test:
164:     suffix: 1d_q2
165:     args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2
166:   test:
167:     suffix: 2d_q1
168:     args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
169:   test:
170:     suffix: 2d_q1d
171:     args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscdualspace_lagrange_continuity 0
172:   test:
173:     suffix: 2d_q1_q1d
174:     args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 1 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0
175:   test:
176:     suffix: 2d_q2
177:     args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2
178:   test:
179:     suffix: 2d_q2_q1
180:     args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1
181:   test:
182:     suffix: 2d_q2_p0
183:     args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 0 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0
184:   test:
185:     suffix: 2d_q2_q1d
186:     args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
187:             -f0_petscspace_degree 2 \
188:             -f1_petscspace_degree 1 -f1_petscdualspace_lagrange_continuity 0
189:   test:
190:     suffix: 2d_q2_p1d
191:     args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \
192:             -f0_petscspace_degree 2 \
193:             -f1_petscspace_degree 1 -f1_petscspace_poly_tensor 0 -f1_petscdualspace_lagrange_continuity 0
194:   test:
195:     suffix: 2d_q3
196:     args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1

198:   test:
199:     suffix: 3d_q1
200:     args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
201:   test:
202:     suffix: 1d_q1_periodic
203:     requires: !complex
204:     args: -dm_plex_dim 1 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3 -dm_plex_box_bd periodic -dm_view -view_coord
205:   test:
206:     suffix: 2d_q1_periodic
207:     requires: !complex
208:     args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2 -dm_plex_box_bd periodic,none -dm_view -view_coord
209:   test:
210:     suffix: 3d_q1_periodic
211:     requires: !complex
212:     args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2,1 -dm_plex_box_bd periodic,none,none -dm_view -view_coord
213:   test:
214:     suffix: 3d_q1_periodic_project
215:     requires: !complex
216:     args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,3 -dm_plex_box_bd none,none,periodic -dm_view -view_coord -project_coordinates

218:   test:
219:     suffix: 3d_q2_periodic # not actually periodic because only 2 cells
220:     args: -dm_plex_dim 3 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2 -dm_plex_box_bd periodic,none,periodic -dm_view

222: TEST*/