Actual source code: plexceed.c
1: #include <petsc/private/dmpleximpl.h>
3: PetscErrorCode DMGetPoints_Internal(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(PetscLogEventBegin(DMPLEX_GetLocalOffsets, dm, 0, 0, 0));
71: PetscCall(DMGetDimension(dm, &dim));
72: PetscCall(DMGetLocalSection(dm, §ion));
73: PetscCall(PetscSectionGetStorageSize(section, l_size));
74: {
75: IS field_is;
76: const PetscInt *fields;
77: PetscInt num_fields;
79: PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL));
80: PetscCheck(field_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Domain label does not have any fields associated with it");
81: // Translate dm_field to ds_field
82: PetscCall(ISGetIndices(field_is, &fields));
83: PetscCall(ISGetSize(field_is, &num_fields));
84: for (PetscInt i = 0; i < num_fields; i++) {
85: if (dm_field == fields[i]) {
86: ds_field = i;
87: break;
88: }
89: }
90: PetscCall(ISRestoreIndices(field_is, &fields));
91: }
92: PetscCheck(ds_field != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field);
94: PetscCall(DMGetPoints_Internal(dm, domain_label, label_value, height, &iter_is));
95: if (iter_is) {
96: PetscCall(ISGetLocalSize(iter_is, num_cells));
97: PetscCall(ISGetIndices(iter_is, &iter_indices));
98: } else {
99: *num_cells = 0;
100: iter_indices = NULL;
101: }
103: {
104: PetscDualSpace dual_space;
105: PetscInt num_dual_basis_vectors;
107: PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
108: PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
109: PetscCheck(fe, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Height %" PetscInt_FMT " is invalid for DG discretizations", height);
110: PetscCall(PetscFEGetDualSpace(fe, &dual_space));
111: PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
112: PetscCall(PetscDualSpaceGetNumComponents(dual_space, num_comp));
113: 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);
114: *cell_size = num_dual_basis_vectors / *num_comp;
115: }
116: PetscInt restr_size = (*num_cells) * (*cell_size);
117: PetscCall(PetscMalloc1(restr_size, &restr_indices));
118: PetscInt cell_offset = 0;
120: PetscInt P = dim - height ? (PetscInt)PetscPowReal(*cell_size, 1.0 / (dim - height)) : 0;
121: for (PetscInt p = 0; p < *num_cells; p++) {
122: PetscBool flip = PETSC_FALSE;
123: PetscInt c = iter_indices[p];
124: PetscInt num_indices, *indices;
125: PetscInt field_offsets[17]; // max number of fields plus 1
126: PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
127: if (height > 0) {
128: PetscInt num_cells_support, num_faces, start = -1;
129: const PetscInt *orients, *faces, *cells;
130: PetscCall(DMPlexGetSupport(dm, c, &cells));
131: PetscCall(DMPlexGetSupportSize(dm, c, &num_cells_support));
132: 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);
133: PetscCall(DMPlexGetCone(dm, cells[0], &faces));
134: PetscCall(DMPlexGetConeSize(dm, cells[0], &num_faces));
135: for (PetscInt i = 0; i < num_faces; i++) {
136: if (faces[i] == c) start = i;
137: }
138: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %" PetscInt_FMT " in cone of its support", c);
139: PetscCall(DMPlexGetConeOrientation(dm, cells[0], &orients));
140: if (orients[start] < 0) flip = PETSC_TRUE;
141: }
143: for (PetscInt i = 0; i < *cell_size; i++) {
144: PetscInt ii = i;
145: if (flip) {
146: if (*cell_size == P) ii = *cell_size - 1 - i;
147: else if (*cell_size == P * P) {
148: PetscInt row = i / P, col = i % P;
149: ii = row + col * P;
150: } 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);
151: }
152: // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
153: PetscInt loc = indices[field_offsets[dm_field] + ii * (*num_comp)];
154: loc = loc < 0 ? -(loc + 1) : loc;
155: restr_indices[cell_offset++] = loc;
156: PetscCheck(loc >= 0 && loc < *l_size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Location %" PetscInt_FMT " not in [0, %" PetscInt_FMT ") local vector", loc, *l_size);
157: }
158: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
159: }
160: 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);
161: if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
162: PetscCall(ISDestroy(&iter_is));
164: *offsets = restr_indices;
165: PetscCall(PetscLogEventEnd(DMPLEX_GetLocalOffsets, dm, 0, 0, 0));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: /*@C
170: DMPlexGetLocalOffsetsSupport - Allocate and populate arrays of local offsets for each face support.
172: Not collective
174: Input Parameters:
175: + dm - The `DMPLEX` object
176: . domain_label - label for `DMPLEX` domain, or NULL for whole domain
177: - label_value - Stratum value
179: Output Parameters:
180: + num_faces - Number of local, non-boundary faces
181: . num_comp - Number of components per dof
182: . l_size - Size of local vector
183: . offsetsNeg - Allocated offsets array for cells on the inward normal side of each face
184: - offsetsPos - Allocated offsets array for cells on the outward normal side of each face
186: Level: developer
188: Notes:
189: 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].
191: Caller is responsible for freeing the offsets array using `PetscFree()`.
193: .seealso: [](ch_unstructured), `DMPlexGetLocalOffsets()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
194: @*/
195: PetscErrorCode DMPlexGetLocalOffsetsSupport(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt *num_faces, PetscInt *num_comp, PetscInt *l_size, PetscInt **offsetsNeg, PetscInt **offsetsPos)
196: {
197: PetscDS ds = NULL;
198: PetscFV fv;
199: PetscSection section;
200: PetscInt dim, height = 1, dm_field = 0, ds_field = 0, Nf, NfInt = 0, Nc;
201: PetscInt *restr_indices_neg, *restr_indices_pos;
202: const PetscInt *iter_indices;
203: IS iter_is;
205: PetscFunctionBeginUser;
207: PetscCall(DMGetDimension(dm, &dim));
208: PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));
209: PetscCall(DMGetLocalSection(dm, §ion));
210: PetscCall(PetscSectionGetStorageSize(section, l_size));
212: PetscCall(DMGetPoints_Internal(dm, domain_label, label_value, height, &iter_is));
213: if (iter_is) {
214: PetscCall(ISGetIndices(iter_is, &iter_indices));
215: PetscCall(ISGetLocalSize(iter_is, &Nf));
216: for (PetscInt p = 0, Ns; p < Nf; ++p) {
217: PetscCall(DMPlexGetSupportSize(dm, iter_indices[p], &Ns));
218: if (Ns == 2) ++NfInt;
219: }
220: *num_faces = NfInt;
221: } else {
222: *num_faces = 0;
223: iter_indices = NULL;
224: }
226: PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fv));
227: PetscCall(PetscFVGetNumComponents(fv, &Nc));
228: PetscCall(PetscMalloc1(NfInt, &restr_indices_neg));
229: PetscCall(PetscMalloc1(NfInt, &restr_indices_pos));
230: PetscInt face_offset_neg = 0, face_offset_pos = 0;
232: for (PetscInt p = 0; p < Nf; ++p) {
233: const PetscInt face = iter_indices[p];
234: PetscInt num_indices, *indices;
235: PetscInt field_offsets[17]; // max number of fields plus 1
236: const PetscInt *supp;
237: PetscInt Ns, loc;
239: PetscCall(DMPlexGetSupport(dm, face, &supp));
240: PetscCall(DMPlexGetSupportSize(dm, face, &Ns));
241: // Ignore boundary faces
242: // TODO check for face on parallel boundary
243: if (Ns == 2) {
244: // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
245: PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
246: PetscCheck(num_indices == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of closure indices %" PetscInt_FMT " != %" PetscInt_FMT " number of FV components", num_indices, Nc);
247: loc = indices[field_offsets[dm_field]];
248: loc = loc < 0 ? -(loc + 1) : loc;
249: restr_indices_neg[face_offset_neg++] = loc;
250: PetscCheck(loc >= 0 && loc + Nc <= *l_size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Location %" PetscInt_FMT " + Nc not in [0, %" PetscInt_FMT ") local vector", loc, *l_size);
251: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
252: PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
253: PetscCheck(num_indices == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of closure indices %" PetscInt_FMT " != %" PetscInt_FMT " number of FV components", num_indices, Nc);
254: loc = indices[field_offsets[dm_field]];
255: loc = loc < 0 ? -(loc + 1) : loc;
256: restr_indices_pos[face_offset_pos++] = loc;
257: PetscCheck(loc >= 0 && loc + Nc <= *l_size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Location %" PetscInt_FMT " + Nc not in [0, %" PetscInt_FMT ") local vector", loc, *l_size);
258: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
259: }
260: }
261: PetscCheck(face_offset_neg == NfInt, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, neg offsets array of shape (%" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", NfInt, face_offset_neg);
262: PetscCheck(face_offset_pos == NfInt, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, pos offsets array of shape (%" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", NfInt, face_offset_pos);
263: if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
264: PetscCall(ISDestroy(&iter_is));
266: *num_comp = Nc;
267: *offsetsNeg = restr_indices_neg;
268: *offsetsPos = restr_indices_pos;
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: #if defined(PETSC_HAVE_LIBCEED)
273: #include <petscdmplexceed.h>
275: // Consumes the input petsc_indices and provides the output ceed_indices; no-copy when the sizes match.
276: static PetscErrorCode PetscIntArrayIntoCeedInt_Private(PetscInt length, PetscInt max_bound, const char *max_bound_name, PetscInt **petsc_indices, CeedInt **ceed_indices)
277: {
278: PetscFunctionBegin;
279: if (length) PetscAssertPointer(petsc_indices, 3);
280: PetscAssertPointer(ceed_indices, 4);
281: #if defined(PETSC_USE_64BIT_INDICES)
282: PetscCheck(max_bound <= PETSC_INT32_MAX, PETSC_COMM_SELF, PETSC_ERR_SUP, "%s %" PetscInt_FMT " does not fit in int32_t", max_bound_name, max_bound);
283: {
284: CeedInt *ceed_ind;
285: PetscCall(PetscMalloc1(length, &ceed_ind));
286: for (PetscInt i = 0; i < length; i++) ceed_ind[i] = (*petsc_indices)[i];
287: *ceed_indices = ceed_ind;
288: PetscCall(PetscFree(*petsc_indices));
289: }
290: #else
291: *ceed_indices = *petsc_indices;
292: *petsc_indices = NULL;
293: #endif
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@C
298: DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector)
300: Input Parameters:
301: + dm - The `DMPLEX` object
302: . domain_label - label for `DMPLEX` domain, or NULL for the whole domain
303: . label_value - Stratum value
304: . height - Height of target cells in `DMPLEX` topology
305: - dm_field - Index of `DMPLEX` field
307: Output Parameter:
308: . ERestrict - libCEED restriction from local vector to the cells
310: Level: developer
312: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetLocalOffsets()`
313: @*/
314: PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
315: {
316: PetscFunctionBeginUser;
318: PetscAssertPointer(ERestrict, 6);
319: if (!dm->ceedERestrict) {
320: PetscInt num_cells, cell_size, num_comp, lvec_size, *restr_indices;
321: CeedElemRestriction elem_restr;
322: Ceed ceed;
324: PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices));
325: PetscCall(DMGetCeed(dm, &ceed));
326: {
327: CeedInt *ceed_indices;
328: PetscCall(PetscIntArrayIntoCeedInt_Private(num_cells * cell_size, lvec_size, "lvec_size", &restr_indices, &ceed_indices));
329: PetscCallCEED(CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, ceed_indices, &elem_restr));
330: PetscCall(PetscFree(ceed_indices));
331: }
332: dm->ceedERestrict = elem_restr;
333: }
334: *ERestrict = dm->ceedERestrict;
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: PetscErrorCode DMPlexCreateCeedRestrictionFVM(DM dm, CeedElemRestriction *erL, CeedElemRestriction *erR)
339: {
340: Ceed ceed;
341: PetscInt *offL, *offR;
342: PetscInt num_faces, num_comp, lvec_size;
344: PetscFunctionBeginUser;
346: PetscAssertPointer(erL, 2);
347: PetscAssertPointer(erR, 3);
348: PetscCall(DMGetCeed(dm, &ceed));
349: PetscCall(DMPlexGetLocalOffsetsSupport(dm, NULL, 0, &num_faces, &num_comp, &lvec_size, &offL, &offR));
350: {
351: CeedInt *ceed_off;
352: PetscCall(PetscIntArrayIntoCeedInt_Private(num_faces * 1, lvec_size, "lvec_size", &offL, &ceed_off));
353: PetscCallCEED(CeedElemRestrictionCreate(ceed, num_faces, 1, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, ceed_off, erL));
354: PetscCall(PetscFree(ceed_off));
355: PetscCall(PetscIntArrayIntoCeedInt_Private(num_faces * 1, lvec_size, "lvec_size", &offR, &ceed_off));
356: PetscCallCEED(CeedElemRestrictionCreate(ceed, num_faces, 1, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, ceed_off, erR));
357: PetscCall(PetscFree(ceed_off));
358: }
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: // TODO DMPlexComputeGeometryFVM() also computes centroids and minimum radius
363: // TODO DMPlexComputeGeometryFVM() flips normal to match support orientation
364: // This function computes area-weights normals
365: PetscErrorCode DMPlexCeedComputeGeometryFVM(DM dm, CeedVector qd)
366: {
367: DMLabel domain_label = NULL;
368: PetscInt label_value = 0, height = 1, Nf, NfInt = 0, cdim;
369: const PetscInt *iter_indices;
370: IS iter_is;
371: CeedScalar *qdata;
373: PetscFunctionBegin;
374: PetscCall(DMGetCoordinateDim(dm, &cdim));
375: PetscCall(DMGetPoints_Internal(dm, domain_label, label_value, height, &iter_is));
376: if (iter_is) {
377: PetscCall(ISGetIndices(iter_is, &iter_indices));
378: PetscCall(ISGetLocalSize(iter_is, &Nf));
379: for (PetscInt p = 0, Ns; p < Nf; ++p) {
380: PetscCall(DMPlexGetSupportSize(dm, iter_indices[p], &Ns));
381: if (Ns == 2) ++NfInt;
382: }
383: } else {
384: iter_indices = NULL;
385: }
387: PetscCallCEED(CeedVectorSetValue(qd, 0.));
388: PetscCallCEED(CeedVectorGetArray(qd, CEED_MEM_HOST, &qdata));
389: for (PetscInt p = 0, off = 0; p < Nf; ++p) {
390: const PetscInt face = iter_indices[p];
391: const PetscInt *supp;
392: PetscInt suppSize;
394: PetscCall(DMPlexGetSupport(dm, face, &supp));
395: PetscCall(DMPlexGetSupportSize(dm, face, &suppSize));
396: // Ignore boundary faces
397: // TODO check for face on parallel boundary
398: if (suppSize == 2) {
399: DMPolytopeType ct;
400: PetscReal area, fcentroid[3], centroids[2][3];
402: PetscCall(DMPlexComputeCellGeometryFVM(dm, face, &area, fcentroid, &qdata[off]));
403: for (PetscInt d = 0; d < cdim; ++d) qdata[off + d] *= area;
404: off += cdim;
405: for (PetscInt s = 0; s < suppSize; ++s) {
406: PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
407: if (ct == DM_POLYTOPE_FV_GHOST) continue;
408: PetscCall(DMPlexComputeCellGeometryFVM(dm, supp[s], &qdata[off + s], centroids[s], NULL));
409: }
410: // Give FV ghosts the same volume as the opposite cell
411: for (PetscInt s = 0; s < suppSize; ++s) {
412: PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
413: if (ct != DM_POLYTOPE_FV_GHOST) continue;
414: qdata[off + s] = qdata[off + (1 - s)];
415: for (PetscInt d = 0; d < cdim; ++d) centroids[s][d] = fcentroid[d];
416: }
417: // Flip normal orientation if necessary to match ordering in support
418: {
419: CeedScalar *normal = &qdata[off - cdim];
420: PetscReal l[3], r[3], v[3];
422: PetscCall(DMLocalizeCoordinateReal_Internal(dm, cdim, fcentroid, centroids[0], l));
423: PetscCall(DMLocalizeCoordinateReal_Internal(dm, cdim, fcentroid, centroids[1], r));
424: DMPlex_WaxpyD_Internal(cdim, -1, l, r, v);
425: if (DMPlex_DotRealD_Internal(cdim, normal, v) < 0) {
426: for (PetscInt d = 0; d < cdim; ++d) normal[d] = -normal[d];
427: }
428: if (DMPlex_DotRealD_Internal(cdim, normal, v) <= 0) {
429: PetscCheck(cdim != 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g) v (%g,%g)", face, (double)normal[0], (double)normal[1], (double)v[0], (double)v[1]);
430: PetscCheck(cdim != 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", face, (double)normal[0], (double)normal[1], (double)normal[2], (double)v[0], (double)v[1], (double)v[2]);
431: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", face);
432: }
433: }
434: off += suppSize;
435: }
436: }
437: PetscCallCEED(CeedVectorRestoreArray(qd, &qdata));
438: if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
439: PetscCall(ISDestroy(&iter_is));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: #endif