1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 /// @file
9 /// Utilities for setting up DM and PetscFE
10
11 #include <ceed.h>
12 #include <petscdmplex.h>
13 #include <petscds.h>
14
15 #include "../navierstokes.h"
16
17 /**
18 @brief Convert `DM` field to `DS` field.
19
20 @param[in] dm `DM` holding mesh
21 @param[in] domain_label Label for `DM` domain
22 @param[in] dm_field Index of `DM` field
23 @param[out] ds_field Index of `DS` field
24
25 @return An error code: 0 - success, otherwise - failure
26 **/
DMFieldToDSField(DM dm,DMLabel domain_label,PetscInt dm_field,PetscInt * ds_field)27 PetscErrorCode DMFieldToDSField(DM dm, DMLabel domain_label, PetscInt dm_field, PetscInt *ds_field) {
28 PetscDS ds;
29 IS field_is;
30 const PetscInt *fields;
31 PetscInt num_fields;
32
33 PetscFunctionBeginUser;
34 PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL));
35 PetscCall(ISGetIndices(field_is, &fields));
36 PetscCall(ISGetSize(field_is, &num_fields));
37 for (PetscInt i = 0; i < num_fields; i++) {
38 if (dm_field == fields[i]) {
39 *ds_field = i;
40 break;
41 }
42 }
43 PetscCall(ISRestoreIndices(field_is, &fields));
44
45 PetscCheck(*ds_field != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field);
46 PetscFunctionReturn(PETSC_SUCCESS);
47 }
48
49 /**
50 @brief Create `CeedElemRestriction` from `DMPlex`.
51
52 Not collective across MPI processes.
53
54 @param[in] ceed `Ceed` context
55 @param[in] dm `DMPlex` holding mesh
56 @param[in] domain_label `DMLabel` for `DMPlex` domain
57 @param[in] label_value Stratum value
58 @param[in] height Height of `DMPlex` topology
59 @param[in] dm_field Index of `DMPlex` field
60 @param[out] restriction `CeedElemRestriction` for `DMPlex`
61
62 @return An error code: 0 - success, otherwise - failure
63 **/
DMPlexCeedElemRestrictionCreate(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt dm_field,CeedElemRestriction * restriction)64 PetscErrorCode DMPlexCeedElemRestrictionCreate(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field,
65 CeedElemRestriction *restriction) {
66 PetscInt num_elem, elem_size, num_dof, num_comp, *restriction_offsets_petsc;
67 CeedInt *restriction_offsets_ceed = NULL;
68
69 PetscFunctionBeginUser;
70 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_elem, &elem_size, &num_comp, &num_dof,
71 &restriction_offsets_petsc));
72 PetscCall(IntArrayPetscToCeed(num_elem * elem_size, &restriction_offsets_petsc, &restriction_offsets_ceed));
73 PetscCallCeed(ceed, CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES,
74 restriction_offsets_ceed, restriction));
75 PetscCall(PetscFree(restriction_offsets_ceed));
76
77 PetscFunctionReturn(PETSC_SUCCESS);
78 }
79
80 /**
81 @brief Create `CeedElemRestriction` from `DMPlex` domain for mesh coordinates.
82
83 Not collective across MPI processes.
84
85 @param[in] ceed `Ceed` context
86 @param[in] dm `DMPlex` holding mesh
87 @param[in] domain_label Label for `DMPlex` domain
88 @param[in] label_value Stratum value
89 @param[in] height Height of `DMPlex` topology
90 @param[out] restriction `CeedElemRestriction` for mesh
91
92 @return An error code: 0 - success, otherwise - failure
93 **/
DMPlexCeedElemRestrictionCoordinateCreate(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,CeedElemRestriction * restriction)94 PetscErrorCode DMPlexCeedElemRestrictionCoordinateCreate(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height,
95 CeedElemRestriction *restriction) {
96 DM dm_coord;
97
98 PetscFunctionBeginUser;
99 PetscCall(DMGetCellCoordinateDM(dm, &dm_coord));
100 if (!dm_coord) {
101 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
102 }
103 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm_coord, domain_label, label_value, height, 0, restriction));
104
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
108 /**
109 @brief Create `CeedElemRestriction` from `DMPlex` domain for auxilury `QFunction` data.
110
111 Not collective across MPI processes.
112
113 @param[in] ceed `Ceed` context
114 @param[in] dm `DMPlex` holding mesh
115 @param[in] domain_label Label for `DMPlex` domain
116 @param[in] label_value Stratum value
117 @param[in] height Height of `DMPlex` topology
118 @param[in] q_data_size Number of components for `QFunction` data
119 @param[in] is_collocated Boolean flag indicating if the data is collocated on the nodes (`PETSC_TRUE`) or on quadrature points (`PETSC_FALSE`)
120 @param[out] restriction Strided `CeedElemRestriction` for `QFunction` data
121
122 @return An error code: 0 - success, otherwise - failure
123 **/
DMPlexCeedElemRestrictionStridedCreate(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt q_data_size,PetscBool is_collocated,CeedElemRestriction * restriction)124 static PetscErrorCode DMPlexCeedElemRestrictionStridedCreate(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height,
125 PetscInt q_data_size, PetscBool is_collocated, CeedElemRestriction *restriction) {
126 PetscInt num_elem, num_qpts, dm_field = 0;
127
128 PetscFunctionBeginUser;
129 { // Get number of elements
130 PetscInt depth;
131 DMLabel depth_label;
132 IS point_is, depth_is;
133
134 PetscCall(DMPlexGetDepth(dm, &depth));
135 PetscCall(DMPlexGetDepthLabel(dm, &depth_label));
136 PetscCall(DMLabelGetStratumIS(depth_label, depth - height, &depth_is));
137 if (domain_label) {
138 IS domain_is;
139
140 PetscCall(DMLabelGetStratumIS(domain_label, label_value, &domain_is));
141 if (domain_is) {
142 PetscCall(ISIntersect(depth_is, domain_is, &point_is));
143 PetscCall(ISDestroy(&domain_is));
144 } else {
145 point_is = NULL;
146 }
147 PetscCall(ISDestroy(&depth_is));
148 } else {
149 point_is = depth_is;
150 }
151 if (point_is) {
152 PetscCall(ISGetLocalSize(point_is, &num_elem));
153 } else {
154 num_elem = 0;
155 }
156 PetscCall(ISDestroy(&point_is));
157 }
158
159 { // Get number of quadrature points
160 PetscDS ds;
161 PetscFE fe;
162 PetscInt ds_field = -1;
163
164 PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));
165 PetscCall(DMFieldToDSField(dm, domain_label, dm_field, &ds_field));
166 PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
167 PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
168 if (is_collocated) {
169 PetscDualSpace dual_space;
170 PetscInt num_dual_basis_vectors, dim, num_comp;
171
172 PetscCall(PetscFEGetSpatialDimension(fe, &dim));
173 PetscCall(PetscFEGetNumComponents(fe, &num_comp));
174 PetscCall(PetscFEGetDualSpace(fe, &dual_space));
175 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
176 num_qpts = num_dual_basis_vectors / num_comp;
177 } else {
178 PetscQuadrature quadrature;
179
180 PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));
181 PetscCall(DMFieldToDSField(dm, domain_label, dm_field, &ds_field));
182 PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
183 PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
184 PetscCall(PetscFEGetQuadrature(fe, &quadrature));
185 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &num_qpts, NULL, NULL));
186 }
187 }
188
189 // Create the restriction
190 PetscCallCeed(ceed, CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, num_elem * num_qpts * q_data_size, CEED_STRIDES_BACKEND,
191 restriction));
192
193 PetscFunctionReturn(PETSC_SUCCESS);
194 }
195
196 /**
197 @brief Create `CeedElemRestriction` from `DMPlex` domain for auxilury `QFunction` data.
198
199 Not collective across MPI processes.
200
201 @param[in] ceed `Ceed` context
202 @param[in] dm `DMPlex` holding mesh
203 @param[in] domain_label Label for `DMPlex` domain
204 @param[in] label_value Stratum value
205 @param[in] height Height of `DMPlex` topology
206 @param[in] q_data_size Number of components for `QFunction` data
207 @param[out] restriction Strided `CeedElemRestriction` for `QFunction` data
208
209 @return An error code: 0 - success, otherwise - failure
210 **/
DMPlexCeedElemRestrictionQDataCreate(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt q_data_size,CeedElemRestriction * restriction)211 PetscErrorCode DMPlexCeedElemRestrictionQDataCreate(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height,
212 PetscInt q_data_size, CeedElemRestriction *restriction) {
213 PetscFunctionBeginUser;
214 PetscCall(DMPlexCeedElemRestrictionStridedCreate(ceed, dm, domain_label, label_value, height, q_data_size, PETSC_FALSE, restriction));
215
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
219 /**
220 @brief Create `CeedElemRestriction` from `DMPlex` domain for nodally collocated auxilury `QFunction` data.
221
222 Not collective across MPI processes.
223
224 @param[in] ceed `Ceed` context
225 @param[in] dm `DMPlex` holding mesh
226 @param[in] domain_label Label for `DMPlex` domain
227 @param[in] label_value Stratum value
228 @param[in] height Height of `DMPlex` topology
229 @param[in] q_data_size Number of components for `QFunction` data
230 @param[out] restriction Strided `CeedElemRestriction` for `QFunction` data
231
232 @return An error code: 0 - success, otherwise - failure
233 **/
DMPlexCeedElemRestrictionCollocatedCreate(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt q_data_size,CeedElemRestriction * restriction)234 PetscErrorCode DMPlexCeedElemRestrictionCollocatedCreate(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height,
235 PetscInt q_data_size, CeedElemRestriction *restriction) {
236 PetscFunctionBeginUser;
237 PetscCall(DMPlexCeedElemRestrictionStridedCreate(ceed, dm, domain_label, label_value, height, q_data_size, PETSC_TRUE, restriction));
238
239 PetscFunctionReturn(PETSC_SUCCESS);
240 }
241
242 // -----------------------------------------------------------------------------
243 // Utility function - convert from DMPolytopeType to CeedElemTopology
244 // -----------------------------------------------------------------------------
ElemTopologyP2C(DMPolytopeType cell_type)245 static CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) {
246 switch (cell_type) {
247 case DM_POLYTOPE_TRIANGLE:
248 return CEED_TOPOLOGY_TRIANGLE;
249 case DM_POLYTOPE_QUADRILATERAL:
250 return CEED_TOPOLOGY_QUAD;
251 case DM_POLYTOPE_TETRAHEDRON:
252 return CEED_TOPOLOGY_TET;
253 case DM_POLYTOPE_HEXAHEDRON:
254 return CEED_TOPOLOGY_HEX;
255 default:
256 return 0;
257 }
258 }
259
260 // -----------------------------------------------------------------------------
261 // Create libCEED Basis from PetscTabulation
262 // -----------------------------------------------------------------------------
BasisCreateFromTabulation(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt face,PetscFE fe,PetscTabulation basis_tabulation,PetscQuadrature quadrature,CeedBasis * basis)263 PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt face, PetscFE fe,
264 PetscTabulation basis_tabulation, PetscQuadrature quadrature, CeedBasis *basis) {
265 PetscInt first_point;
266 PetscInt ids[1] = {label_value};
267 DMLabel depth_label;
268 DMPolytopeType cell_type;
269 CeedElemTopology elem_topo;
270 PetscScalar *q_points, *interp, *grad;
271 const PetscScalar *q_weights;
272 PetscDualSpace dual_space;
273 PetscInt num_dual_basis_vectors;
274 PetscInt dim, num_comp, P, Q;
275
276 PetscFunctionBeginUser;
277 PetscCall(PetscFEGetSpatialDimension(fe, &dim));
278 PetscCall(PetscFEGetNumComponents(fe, &num_comp));
279 PetscCall(PetscFEGetDualSpace(fe, &dual_space));
280 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
281 P = num_dual_basis_vectors / num_comp;
282
283 // Use depth label if no domain label present
284 if (!domain_label) {
285 PetscInt depth;
286
287 PetscCall(DMPlexGetDepth(dm, &depth));
288 PetscCall(DMPlexGetDepthLabel(dm, &depth_label));
289 ids[0] = depth - height;
290 }
291
292 // Get cell interp, grad, and quadrature data
293 PetscCall(DMGetFirstLabeledPoint(dm, dm, domain_label ? domain_label : depth_label, 1, ids, height, &first_point, NULL));
294 PetscCall(DMPlexGetCellType(dm, first_point, &cell_type));
295 elem_topo = ElemTopologyP2C(cell_type);
296 PetscCheck(elem_topo, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMPlex topology not supported");
297 {
298 size_t q_points_size;
299 const PetscScalar *q_points_petsc;
300 PetscInt q_dim;
301
302 PetscCall(PetscQuadratureGetData(quadrature, &q_dim, NULL, &Q, &q_points_petsc, &q_weights));
303 q_points_size = Q * dim * sizeof(CeedScalar);
304 PetscCall(PetscCalloc(q_points_size, &q_points));
305 for (PetscInt q = 0; q < Q; q++) {
306 for (PetscInt d = 0; d < q_dim; d++) q_points[q * dim + d] = q_points_petsc[q * q_dim + d];
307 }
308 }
309
310 { // Convert to libCEED orientation
311 PetscBool is_simplex = PETSC_FALSE;
312 IS permutation = NULL;
313 const PetscInt *permutation_indices;
314
315 PetscCall(DMPlexIsSimplex(dm, &is_simplex));
316 if (!is_simplex) {
317 PetscSection section;
318
319 // -- Get permutation
320 PetscCall(DMGetLocalSection(dm, §ion));
321 PetscCall(PetscSectionGetClosurePermutation(section, (PetscObject)dm, dim, num_comp * P, &permutation));
322 PetscCall(ISGetIndices(permutation, &permutation_indices));
323 }
324
325 // -- Copy interp, grad matrices
326 PetscCall(PetscCalloc(P * Q * sizeof(CeedScalar), &interp));
327 PetscCall(PetscCalloc(P * Q * dim * sizeof(CeedScalar), &grad));
328 const CeedInt c = 0;
329 for (CeedInt q = 0; q < Q; q++) {
330 for (CeedInt p_ceed = 0; p_ceed < P; p_ceed++) {
331 CeedInt p_petsc = is_simplex ? (p_ceed * num_comp) : permutation_indices[p_ceed * num_comp];
332
333 interp[q * P + p_ceed] = basis_tabulation->T[0][((face * Q + q) * P * num_comp + p_petsc) * num_comp + c];
334 for (CeedInt d = 0; d < dim; d++) {
335 grad[(d * Q + q) * P + p_ceed] = basis_tabulation->T[1][(((face * Q + q) * P * num_comp + p_petsc) * num_comp + c) * dim + d];
336 }
337 }
338 }
339
340 // -- Cleanup
341 if (permutation) PetscCall(ISRestoreIndices(permutation, &permutation_indices));
342 PetscCall(ISDestroy(&permutation));
343 }
344
345 // Finally, create libCEED basis
346 PetscCallCeed(ceed, CeedBasisCreateH1(ceed, elem_topo, num_comp, P, Q, interp, grad, q_points, q_weights, basis));
347 PetscCall(PetscFree(q_points));
348 PetscCall(PetscFree(interp));
349 PetscCall(PetscFree(grad));
350 PetscFunctionReturn(PETSC_SUCCESS);
351 }
352
353 // -----------------------------------------------------------------------------
354 // Get CEED Basis from DMPlex
355 // -----------------------------------------------------------------------------
CreateBasisFromPlex(Ceed ceed,DM dm,DMLabel domain_label,CeedInt label_value,CeedInt height,CeedInt dm_field,CeedBasis * basis)356 PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, CeedBasis *basis) {
357 PetscDS ds;
358 PetscFE fe;
359 PetscQuadrature quadrature;
360 PetscBool is_simplex = PETSC_TRUE;
361 PetscInt ds_field = -1;
362
363 PetscFunctionBeginUser;
364 // Get element information
365 PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));
366 PetscCall(DMFieldToDSField(dm, domain_label, dm_field, &ds_field));
367 PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
368 PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
369 PetscCall(PetscFEGetQuadrature(fe, &quadrature));
370
371 // Check if simplex or tensor-product mesh
372 PetscCall(DMPlexIsSimplex(dm, &is_simplex));
373
374 // Build libCEED basis
375 if (is_simplex) {
376 PetscTabulation basis_tabulation;
377 PetscInt num_derivatives = 1, face = 0;
378
379 PetscCall(PetscFEGetCellTabulation(fe, num_derivatives, &basis_tabulation));
380 PetscCall(BasisCreateFromTabulation(ceed, dm, domain_label, label_value, height, face, fe, basis_tabulation, quadrature, basis));
381 } else {
382 PetscDualSpace dual_space;
383 PetscInt num_dual_basis_vectors;
384 PetscInt dim, num_comp, P, Q;
385
386 PetscCall(PetscFEGetSpatialDimension(fe, &dim));
387 PetscCall(PetscFEGetNumComponents(fe, &num_comp));
388 PetscCall(PetscFEGetDualSpace(fe, &dual_space));
389 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
390 P = num_dual_basis_vectors / num_comp;
391 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &Q, NULL, NULL));
392
393 CeedInt P_1d = (CeedInt)round(pow(P, 1.0 / dim));
394 CeedInt Q_1d = (CeedInt)round(pow(Q, 1.0 / dim));
395
396 PetscCallCeed(ceed, CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, CEED_GAUSS, basis));
397 }
398 PetscFunctionReturn(PETSC_SUCCESS);
399 }
400
401 /**
402 @brief Setup `DM` with FE space of appropriate degree
403
404 Must be followed by `DMSetupByOrderEnd_FEM`
405
406 @param[in] setup_faces Flag to setup face geometry
407 @param[in] setup_coords Flag to setup coordinate spaces
408 @param[in] degree Polynomial orders of field
409 @param[in] coord_order Polynomial order of coordinate basis, or `PETSC_DECIDE` for default
410 @param[in] q_extra Additional quadrature order
411 @param[in] num_fields Number of fields in solution vector
412 @param[in] field_sizes Array of number of components for each field
413 @param[out] dm `DM` to setup
414
415 @return An error code: 0 - success, otherwise - failure
416 **/
DMSetupByOrderBegin_FEM(PetscBool setup_faces,PetscBool setup_coords,PetscInt degree,PetscInt coord_order,PetscInt q_extra,PetscInt num_fields,const PetscInt * field_sizes,DM dm)417 PetscErrorCode DMSetupByOrderBegin_FEM(PetscBool setup_faces, PetscBool setup_coords, PetscInt degree, PetscInt coord_order, PetscInt q_extra,
418 PetscInt num_fields, const PetscInt *field_sizes, DM dm) {
419 PetscInt dim, q_order = degree + q_extra;
420 PetscBool is_simplex = PETSC_TRUE;
421 PetscFE fe;
422 MPI_Comm comm = PetscObjectComm((PetscObject)dm);
423
424 PetscFunctionBeginUser;
425 PetscCall(DMPlexIsSimplex(dm, &is_simplex));
426
427 // Setup DM
428 PetscCall(DMGetDimension(dm, &dim));
429 for (PetscInt i = 0; i < num_fields; i++) {
430 PetscFE fe_face;
431 PetscInt q_order = degree + q_extra;
432
433 PetscCall(PetscFECreateLagrange(comm, dim, field_sizes[i], is_simplex, degree, q_order, &fe));
434 if (setup_faces) PetscCall(PetscFEGetHeightSubspace(fe, 1, &fe_face));
435 PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
436 PetscCall(PetscFEDestroy(&fe));
437 }
438 PetscCall(DMCreateDS(dm));
439
440 // Project coordinates to enrich quadrature space
441 if (setup_coords) {
442 DM dm_coord;
443 PetscDS ds_coord;
444 PetscFE fe_coord_current, fe_coord_new, fe_coord_face_new;
445 PetscDualSpace fe_coord_dual_space;
446 PetscInt fe_coord_order, num_comp_coord;
447
448 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
449 PetscCall(DMGetCoordinateDim(dm, &num_comp_coord));
450 PetscCall(DMGetRegionDS(dm_coord, NULL, NULL, &ds_coord, NULL));
451 PetscCall(PetscDSGetDiscretization(ds_coord, 0, (PetscObject *)&fe_coord_current));
452 PetscCall(PetscFEGetDualSpace(fe_coord_current, &fe_coord_dual_space));
453 PetscCall(PetscDualSpaceGetOrder(fe_coord_dual_space, &fe_coord_order));
454
455 // Create FE for coordinates
456 if (coord_order != PETSC_DECIDE) fe_coord_order = coord_order;
457 PetscCall(PetscFECreateLagrange(comm, dim, num_comp_coord, is_simplex, fe_coord_order, q_order, &fe_coord_new));
458 if (setup_faces) PetscCall(PetscFEGetHeightSubspace(fe_coord_new, 1, &fe_coord_face_new));
459 PetscCall(DMSetCoordinateDisc(dm, fe_coord_new, PETSC_TRUE));
460 PetscCall(PetscFEDestroy(&fe_coord_new));
461 }
462 PetscFunctionReturn(PETSC_SUCCESS);
463 }
464
465 /**
466 @brief Finish setting up `DM`
467
468 Must be called after `DMSetupByOrderBegin_FEM` and all strong BCs have be declared via `DMAddBoundaries`.
469
470 @param[in] setup_coords Flag to setup coordinate spaces
471 @param[out] dm `DM` to setup
472
473 @return An error code: 0 - success, otherwise - failure
474 **/
DMSetupByOrderEnd_FEM(PetscBool setup_coords,DM dm)475 PetscErrorCode DMSetupByOrderEnd_FEM(PetscBool setup_coords, DM dm) {
476 PetscBool is_simplex;
477
478 PetscFunctionBeginUser;
479 PetscCall(DMPlexIsSimplex(dm, &is_simplex));
480 // Set tensor permutation if needed
481 if (!is_simplex) {
482 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
483 if (setup_coords) {
484 DM dm_coord;
485
486 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
487 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL));
488 }
489 }
490 PetscCall(DMLocalizeCoordinates(dm)); // Must localize after tensor closure setting
491 PetscFunctionReturn(PETSC_SUCCESS);
492 }
493
494 /**
495 @brief Setup `DM` with FE space of appropriate degree with no boundary conditions
496
497 Calls `DMSetupByOrderBegin_FEM` and `DMSetupByOrderEnd_FEM` successively
498
499 @param[in] setup_faces Flag to setup face geometry
500 @param[in] setup_coords Flag to setup coordinate spaces
501 @param[in] degree Polynomial orders of field
502 @param[in] coord_order Polynomial order of coordinate basis, or `PETSC_DECIDE` for default
503 @param[in] q_extra Additional quadrature order
504 @param[in] num_fields Number of fields in solution vector
505 @param[in] field_sizes Array of number of components for each field
506 @param[out] dm `DM` to setup
507
508 @return An error code: 0 - success, otherwise - failure
509 **/
DMSetupByOrder_FEM(PetscBool setup_faces,PetscBool setup_coords,PetscInt degree,PetscInt coord_order,PetscInt q_extra,PetscInt num_fields,const PetscInt * field_sizes,DM dm)510 PetscErrorCode DMSetupByOrder_FEM(PetscBool setup_faces, PetscBool setup_coords, PetscInt degree, PetscInt coord_order, PetscInt q_extra,
511 PetscInt num_fields, const PetscInt *field_sizes, DM dm) {
512 PetscFunctionBeginUser;
513 PetscCall(DMSetupByOrderBegin_FEM(setup_faces, setup_coords, degree, coord_order, q_extra, num_fields, field_sizes, dm));
514 PetscCall(DMSetupByOrderEnd_FEM(setup_coords, dm));
515 PetscFunctionReturn(PETSC_SUCCESS);
516 }
517