1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
DMGetPoints_Internal(DM dm,DMLabel domainLabel,PetscInt labelVal,PetscInt height,IS * pointIS)3 PetscErrorCode DMGetPoints_Internal(DM dm, DMLabel domainLabel, PetscInt labelVal, PetscInt height, IS *pointIS)
4 {
5 PetscInt depth;
6 DMLabel depthLabel;
7 IS depthIS;
8
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;
15
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 }
29
30 /*@C
31 DMPlexGetLocalOffsets - Allocate and populate array of local offsets for each cell closure.
32
33 Not collective
34
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
41
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
48
49 Level: developer
50
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].
53
54 Caller is responsible for freeing the offsets array using `PetscFree()`.
55
56 .seealso: [](ch_unstructured), `DMPlexGetLocalOffsetsSupport()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
57 @*/
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[])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;
67
68 PetscFunctionBeginUser;
69 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
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;
78
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);
93
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 }
102
103 {
104 PetscDualSpace dual_space;
105 PetscInt num_dual_basis_vectors;
106
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;
119
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 }
142
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));
163
164 *offsets = restr_indices;
165 PetscCall(PetscLogEventEnd(DMPLEX_GetLocalOffsets, dm, 0, 0, 0));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
169 /*@C
170 DMPlexGetLocalOffsetsSupport - Allocate and populate arrays of local offsets for each face support.
171
172 Not collective
173
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
178
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
185
186 Level: developer
187
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].
190
191 Caller is responsible for freeing the offsets array using `PetscFree()`.
192
193 .seealso: [](ch_unstructured), `DMPlexGetLocalOffsets()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
194 @*/
DMPlexGetLocalOffsetsSupport(DM dm,DMLabel domain_label,PetscInt label_value,PetscInt * num_faces,PetscInt * num_comp,PetscInt * l_size,PetscInt ** offsetsNeg,PetscInt ** offsetsPos)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;
204
205 PetscFunctionBeginUser;
206 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
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));
211
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 }
225
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;
231
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;
238
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));
265
266 *num_comp = Nc;
267 *offsetsNeg = restr_indices_neg;
268 *offsetsPos = restr_indices_pos;
269 PetscFunctionReturn(PETSC_SUCCESS);
270 }
271
272 #if defined(PETSC_HAVE_LIBCEED)
273 #include <petscdmplexceed.h>
274
275 // Consumes the input petsc_indices and provides the output ceed_indices; no-copy when the sizes match.
PetscIntArrayIntoCeedInt_Private(PetscInt length,PetscInt max_bound,const char * max_bound_name,PetscInt ** petsc_indices,CeedInt ** ceed_indices)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 }
296
297 /*@C
298 DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector)
299
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
306
307 Output Parameter:
308 . ERestrict - libCEED restriction from local vector to the cells
309
310 Level: developer
311
312 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetLocalOffsets()`
313 @*/
DMPlexGetCeedRestriction(DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt dm_field,CeedElemRestriction * ERestrict)314 PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
315 {
316 PetscFunctionBeginUser;
317 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
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;
323
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 }
337
DMPlexCreateCeedRestrictionFVM(DM dm,CeedElemRestriction * erL,CeedElemRestriction * erR)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;
343
344 PetscFunctionBeginUser;
345 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
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 }
361
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
DMPlexCeedComputeGeometryFVM(DM dm,CeedVector qd)365 PetscErrorCode DMPlexCeedComputeGeometryFVM(DM dm, CeedVector qd)
366 {
367 DMLabel domain_label = NULL;
368 PetscInt label_value = 0, height = 1, Nf, cdim;
369 const PetscInt *iter_indices;
370 IS iter_is;
371 CeedScalar *qdata;
372
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) PetscCall(DMPlexGetSupportSize(dm, iter_indices[p], &Ns));
380 } else {
381 iter_indices = NULL;
382 }
383
384 PetscCallCEED(CeedVectorSetValue(qd, 0.));
385 PetscCallCEED(CeedVectorGetArray(qd, CEED_MEM_HOST, &qdata));
386 for (PetscInt p = 0, off = 0; p < Nf; ++p) {
387 const PetscInt face = iter_indices[p];
388 const PetscInt *supp;
389 PetscInt suppSize;
390
391 PetscCall(DMPlexGetSupport(dm, face, &supp));
392 PetscCall(DMPlexGetSupportSize(dm, face, &suppSize));
393 // Ignore boundary faces
394 // TODO check for face on parallel boundary
395 if (suppSize == 2) {
396 DMPolytopeType ct;
397 PetscReal area, fcentroid[3], centroids[2][3];
398
399 PetscCall(DMPlexComputeCellGeometryFVM(dm, face, &area, fcentroid, &qdata[off]));
400 for (PetscInt d = 0; d < cdim; ++d) qdata[off + d] *= area;
401 off += cdim;
402 for (PetscInt s = 0; s < suppSize; ++s) {
403 PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
404 if (ct == DM_POLYTOPE_FV_GHOST) continue;
405 PetscCall(DMPlexComputeCellGeometryFVM(dm, supp[s], &qdata[off + s], centroids[s], NULL));
406 }
407 // Give FV ghosts the same volume as the opposite cell
408 for (PetscInt s = 0; s < suppSize; ++s) {
409 PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
410 if (ct != DM_POLYTOPE_FV_GHOST) continue;
411 qdata[off + s] = qdata[off + (1 - s)];
412 for (PetscInt d = 0; d < cdim; ++d) centroids[s][d] = fcentroid[d];
413 }
414 // Flip normal orientation if necessary to match ordering in support
415 {
416 CeedScalar *normal = &qdata[off - cdim];
417 PetscReal l[3], r[3], v[3];
418
419 PetscCall(DMLocalizeCoordinateReal_Internal(dm, cdim, fcentroid, centroids[0], l));
420 PetscCall(DMLocalizeCoordinateReal_Internal(dm, cdim, fcentroid, centroids[1], r));
421 DMPlex_WaxpyD_Internal(cdim, -1, l, r, v);
422 if (DMPlex_DotRealD_Internal(cdim, normal, v) < 0) {
423 for (PetscInt d = 0; d < cdim; ++d) normal[d] = -normal[d];
424 }
425 if (DMPlex_DotRealD_Internal(cdim, normal, v) <= 0) {
426 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]);
427 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]);
428 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", face);
429 }
430 }
431 off += suppSize;
432 }
433 }
434 PetscCallCEED(CeedVectorRestoreArray(qd, &qdata));
435 if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
436 PetscCall(ISDestroy(&iter_is));
437 PetscFunctionReturn(PETSC_SUCCESS);
438 }
439
440 #endif
441