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