Actual source code: plexceed.c

  1: #include <petsc/private/dmpleximpl.h>

  3: static PetscErrorCode DMGetPoints_Private(DM dm, DMLabel domainLabel, PetscInt labelVal, PetscInt height, IS *pointIS)
  4: {
  5:   PetscInt depth;
  6:   DMLabel  depthLabel;
  7:   IS       depthIS;

  9:   PetscFunctionBegin;
 10:   PetscCall(DMPlexGetDepth(dm, &depth));
 11:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
 12:   PetscCall(DMLabelGetStratumIS(depthLabel, depth - height, &depthIS));
 13:   if (domainLabel) {
 14:     IS domainIS;

 16:     PetscCall(DMLabelGetStratumIS(domainLabel, labelVal, &domainIS));
 17:     if (domainIS) { // domainIS is non-empty
 18:       PetscCall(ISIntersect(depthIS, domainIS, pointIS));
 19:       PetscCall(ISDestroy(&domainIS));
 20:     } else { // domainIS is NULL (empty)
 21:       *pointIS = NULL;
 22:     }
 23:     PetscCall(ISDestroy(&depthIS));
 24:   } else {
 25:     *pointIS = depthIS;
 26:   }
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: /*@C
 31:   DMPlexGetLocalOffsets - Allocate and populate array of local offsets for each cell closure.

 33:   Not collective

 35:   Input Parameters:
 36: +  dm - The `DMPLEX` object
 37: .  domain_label - label for `DMPLEX` domain, or NULL for whole domain
 38: .  label_value - Stratum value
 39: .  height - Height of target cells in `DMPLEX` topology
 40: -  dm_field - Index of `DMPLEX` field

 42:   Output Parameters:
 43: +  num_cells - Number of local cells
 44: .  cell_size - Size of each cell, given by cell_size * num_comp = num_dof
 45: .  num_comp - Number of components per dof
 46: .  l_size - Size of local vector
 47: -  offsets - Allocated offsets array for cells

 49:   Level: developer

 51:   Notes:
 52:   Allocate and populate array of shape [num_cells, cell_size] defining offsets for each value (cell, node) for local vector of the `DMPLEX` field. All offsets are in the range [0, l_size - 1].

 54:    Caller is responsible for freeing the offsets array using `PetscFree()`.

 56: .seealso: [](ch_unstructured), `DMPlexGetLocalOffsetsSupport()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
 57: @*/
 58: PetscErrorCode DMPlexGetLocalOffsets(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, PetscInt *num_cells, PetscInt *cell_size, PetscInt *num_comp, PetscInt *l_size, PetscInt **offsets)
 59: {
 60:   PetscDS         ds = NULL;
 61:   PetscFE         fe;
 62:   PetscSection    section;
 63:   PetscInt        dim, ds_field = -1;
 64:   PetscInt       *restr_indices;
 65:   const PetscInt *iter_indices;
 66:   IS              iter_is;

 68:   PetscFunctionBeginUser;
 70:   PetscCall(DMGetLocalSection(dm, &section));
 71:   PetscCall(DMGetDimension(dm, &dim));
 72:   {
 73:     IS              field_is;
 74:     const PetscInt *fields;
 75:     PetscInt        num_fields;

 77:     PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL));
 78:     // Translate dm_field to ds_field
 79:     PetscCall(ISGetIndices(field_is, &fields));
 80:     PetscCall(ISGetSize(field_is, &num_fields));
 81:     for (PetscInt i = 0; i < num_fields; i++) {
 82:       if (dm_field == fields[i]) {
 83:         ds_field = i;
 84:         break;
 85:       }
 86:     }
 87:     PetscCall(ISRestoreIndices(field_is, &fields));
 88:   }
 89:   PetscCheck(ds_field != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field);

 91:   PetscCall(DMGetPoints_Private(dm, domain_label, label_value, height, &iter_is));
 92:   if (iter_is) {
 93:     PetscCall(ISGetLocalSize(iter_is, num_cells));
 94:     PetscCall(ISGetIndices(iter_is, &iter_indices));
 95:   } else {
 96:     *num_cells   = 0;
 97:     iter_indices = NULL;
 98:   }

100:   {
101:     PetscDualSpace dual_space;
102:     PetscInt       num_dual_basis_vectors;

104:     PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
105:     PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
106:     PetscCheck(fe, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Height %" PetscInt_FMT " is invalid for DG coordinates", height);
107:     PetscCall(PetscFEGetDualSpace(fe, &dual_space));
108:     PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
109:     PetscCall(PetscDualSpaceGetNumComponents(dual_space, num_comp));
110:     PetscCheck(num_dual_basis_vectors % *num_comp == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for number of dual basis vectors %" PetscInt_FMT " not divisible by %" PetscInt_FMT " components", num_dual_basis_vectors, *num_comp);
111:     *cell_size = num_dual_basis_vectors / *num_comp;
112:   }
113:   PetscInt restr_size = (*num_cells) * (*cell_size);
114:   PetscCall(PetscMalloc1(restr_size, &restr_indices));
115:   PetscInt cell_offset = 0;

117:   PetscInt P = (PetscInt)PetscPowReal(*cell_size, 1.0 / (dim - height));
118:   for (PetscInt p = 0; p < *num_cells; p++) {
119:     PetscBool flip = PETSC_FALSE;
120:     PetscInt  c    = iter_indices[p];
121:     PetscInt  num_indices, *indices;
122:     PetscInt  field_offsets[17]; // max number of fields plus 1
123:     PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
124:     if (height > 0) {
125:       PetscInt        num_cells_support, num_faces, start = -1;
126:       const PetscInt *orients, *faces, *cells;
127:       PetscCall(DMPlexGetSupport(dm, c, &cells));
128:       PetscCall(DMPlexGetSupportSize(dm, c, &num_cells_support));
129:       PetscCheck(num_cells_support == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected one cell in support of exterior face, but got %" PetscInt_FMT " cells", num_cells_support);
130:       PetscCall(DMPlexGetCone(dm, cells[0], &faces));
131:       PetscCall(DMPlexGetConeSize(dm, cells[0], &num_faces));
132:       for (PetscInt i = 0; i < num_faces; i++) {
133:         if (faces[i] == c) start = i;
134:       }
135:       PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %" PetscInt_FMT " in cone of its support", c);
136:       PetscCall(DMPlexGetConeOrientation(dm, cells[0], &orients));
137:       if (orients[start] < 0) flip = PETSC_TRUE;
138:     }

140:     for (PetscInt i = 0; i < *cell_size; i++) {
141:       PetscInt ii = i;
142:       if (flip) {
143:         if (*cell_size == P) ii = *cell_size - 1 - i;
144:         else if (*cell_size == P * P) {
145:           PetscInt row = i / P, col = i % P;
146:           ii = row + col * P;
147:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for flipping point with cell size %" PetscInt_FMT " != P (%" PetscInt_FMT ") or P^2", *cell_size, P);
148:       }
149:       // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
150:       PetscInt loc                 = indices[field_offsets[dm_field] + ii * (*num_comp)];
151:       restr_indices[cell_offset++] = loc >= 0 ? loc : -(loc + 1);
152:     }
153:     PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
154:   }
155:   PetscCheck(cell_offset == restr_size, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, offsets array of shape (%" PetscInt_FMT ", %" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", *num_cells, (*cell_size), cell_offset);
156:   if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
157:   PetscCall(ISDestroy(&iter_is));

159:   *offsets = restr_indices;
160:   PetscCall(PetscSectionGetStorageSize(section, l_size));
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

164: /*@C
165:   DMPlexGetLocalOffsetsSupport - Allocate and populate arrays of local offsets for each face support.

167:   Not collective

169:   Input Parameters:
170: +  dm - The `DMPLEX` object
171: .  domain_label - label for `DMPLEX` domain, or NULL for whole domain
172: -  label_value - Stratum value

174:   Output Parameters:
175: +  num_faces - Number of local, non-boundary faces
176: .  num_comp - Number of components per dof
177: .  l_size - Size of local vector
178: .  offsetsNeg - Allocated offsets array for cells on the inward normal side of each face
179: -  offsetsPos - Allocated offsets array for cells on the outward normal side of each face

181:   Level: developer

183:   Notes:
184:   Allocate and populate array of shape [num_cells, num_comp] defining offsets for each cell for local vector of the `DMPLEX` field. All offsets are in the range [0, l_size - 1].

186:    Caller is responsible for freeing the offsets array using `PetscFree()`.

188: .seealso: [](ch_unstructured), `DMPlexGetLocalOffsets()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
189: @*/
190: PetscErrorCode DMPlexGetLocalOffsetsSupport(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt *num_faces, PetscInt *num_comp, PetscInt *l_size, PetscInt **offsetsNeg, PetscInt **offsetsPos)
191: {
192:   PetscDS         ds = NULL;
193:   PetscFV         fv;
194:   PetscSection    section;
195:   PetscInt        dim, height = 1, dm_field = 0, ds_field = 0, Nf, NfInt = 0, cell_size, restr_size;
196:   PetscInt       *restr_indices_neg, *restr_indices_pos;
197:   const PetscInt *iter_indices;
198:   IS              iter_is;

200:   PetscFunctionBeginUser;
202:   PetscCall(DMGetLocalSection(dm, &section));
203:   PetscCall(DMGetDimension(dm, &dim));
204:   PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));

206:   PetscCall(DMGetPoints_Private(dm, domain_label, label_value, height, &iter_is));
207:   if (iter_is) {
208:     PetscCall(ISGetIndices(iter_is, &iter_indices));
209:     PetscCall(ISGetLocalSize(iter_is, &Nf));
210:     for (PetscInt p = 0, Ns; p < Nf; ++p) {
211:       PetscCall(DMPlexGetSupportSize(dm, iter_indices[p], &Ns));
212:       if (Ns == 2) ++NfInt;
213:     }
214:     *num_faces = NfInt;
215:   } else {
216:     *num_faces   = 0;
217:     iter_indices = NULL;
218:   }

220:   PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fv));
221:   PetscCall(PetscFVGetNumComponents(fv, num_comp));
222:   cell_size  = *num_comp;
223:   restr_size = (*num_faces) * cell_size;
224:   PetscCall(PetscMalloc1(restr_size, &restr_indices_neg));
225:   PetscCall(PetscMalloc1(restr_size, &restr_indices_pos));
226:   PetscInt face_offset_neg = 0, face_offset_pos = 0;

228:   for (PetscInt p = 0; p < Nf; ++p) {
229:     const PetscInt  face = iter_indices[p];
230:     PetscInt        num_indices, *indices;
231:     PetscInt        field_offsets[17]; // max number of fields plus 1
232:     const PetscInt *supp;
233:     PetscInt        Ns;

235:     PetscCall(DMPlexGetSupport(dm, face, &supp));
236:     PetscCall(DMPlexGetSupportSize(dm, face, &Ns));
237:     // Ignore boundary faces
238:     //   TODO check for face on parallel boundary
239:     if (Ns == 2) {
240:       // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
241:       PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
242:       for (PetscInt i = 0; i < cell_size; i++) {
243:         const PetscInt loc                   = indices[field_offsets[dm_field] + i * (*num_comp)];
244:         restr_indices_neg[face_offset_neg++] = loc >= 0 ? loc : -(loc + 1);
245:       }
246:       PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
247:       PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
248:       for (PetscInt i = 0; i < cell_size; i++) {
249:         const PetscInt loc                   = indices[field_offsets[dm_field] + i * (*num_comp)];
250:         restr_indices_pos[face_offset_pos++] = loc >= 0 ? loc : -(loc + 1);
251:       }
252:       PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
253:     }
254:   }
255:   PetscCheck(face_offset_neg == restr_size, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, neg offsets array of shape (%" PetscInt_FMT ", %" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", *num_faces, cell_size, face_offset_neg);
256:   PetscCheck(face_offset_pos == restr_size, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, pos offsets array of shape (%" PetscInt_FMT ", %" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", *num_faces, cell_size, face_offset_pos);
257:   if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
258:   PetscCall(ISDestroy(&iter_is));

260:   *offsetsNeg = restr_indices_neg;
261:   *offsetsPos = restr_indices_pos;
262:   PetscCall(PetscSectionGetStorageSize(section, l_size));
263:   PetscFunctionReturn(PETSC_SUCCESS);
264: }

266: #if defined(PETSC_HAVE_LIBCEED)
267: #include <petscdmplexceed.h>

269: /*@C
270:   DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector)

272:   Input Parameters:
273: +  dm - The `DMPLEX` object
274: .  domain_label - label for `DMPLEX` domain, or NULL for the whole domain
275: .  label_value - Stratum value
276: .  height - Height of target cells in `DMPLEX` topology
277: -  dm_field - Index of `DMPLEX` field

279:   Output Parameter:
280: .  ERestrict - libCEED restriction from local vector to to the cells

282:   Level: developer

284: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetLocalOffsets()`
285: @*/
286: PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
287: {
288:   PetscFunctionBeginUser;
291:   if (!dm->ceedERestrict) {
292:     PetscInt            num_cells, cell_size, num_comp, lvec_size, *restr_indices;
293:     CeedElemRestriction elem_restr;
294:     Ceed                ceed;

296:     PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices));
297:     PetscCall(DMGetCeed(dm, &ceed));
298:     PetscCallCEED(CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices, &elem_restr));
299:     PetscCall(PetscFree(restr_indices));
300:     dm->ceedERestrict = elem_restr;
301:   }
302:   *ERestrict = dm->ceedERestrict;
303:   PetscFunctionReturn(PETSC_SUCCESS);
304: }

306: #endif