xref: /libCEED/examples/fluids/src/dm_utils.c (revision 5ebd836c59d60a2e5e1cb67f6731404c7da26f85)
1 // Copyright (c) 2017-2024, 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 **/
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 **/
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(
71       DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_elem, &elem_size, &num_comp, &num_dof, &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 **/
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 **/
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 **/
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 **/
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 // -----------------------------------------------------------------------------
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 // -----------------------------------------------------------------------------
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, &section));
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 // -----------------------------------------------------------------------------
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 **/
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(DMLocalizeCoordinates(dm));  // Update CellCoordinateDM with projected coordinates
461     PetscCall(PetscFEDestroy(&fe_coord_new));
462   }
463   PetscFunctionReturn(PETSC_SUCCESS);
464 }
465 
466 /**
467   @brief Finish setting up `DM`
468 
469   Must be called after `DMSetupByOrderBegin_FEM` and all strong BCs have be declared via `DMAddBoundaries`.
470 
471   @param[in]   setup_coords    Flag to setup coordinate spaces
472   @param[out]  dm              `DM` to setup
473 
474   @return An error code: 0 - success, otherwise - failure
475 **/
476 PetscErrorCode DMSetupByOrderEnd_FEM(PetscBool setup_coords, DM dm) {
477   PetscBool is_simplex;
478 
479   PetscFunctionBeginUser;
480   PetscCall(DMPlexIsSimplex(dm, &is_simplex));
481   // Set tensor permutation if needed
482   if (!is_simplex) {
483     PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
484     if (setup_coords) {
485       DM dm_coord;
486 
487       PetscCall(DMGetCoordinateDM(dm, &dm_coord));
488       PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL));
489       PetscCall(DMGetCellCoordinateDM(dm, &dm_coord));
490       if (dm_coord) PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL));
491     }
492   }
493   PetscFunctionReturn(PETSC_SUCCESS);
494 }
495 
496 /**
497   @brief Setup `DM` with FE space of appropriate degree with no boundary conditions
498 
499   Calls `DMSetupByOrderBegin_FEM` and `DMSetupByOrderEnd_FEM` successively
500 
501   @param[in]   setup_faces     Flag to setup face geometry
502   @param[in]   setup_coords    Flag to setup coordinate spaces
503   @param[in]   degree          Polynomial orders of field
504   @param[in]   coord_order     Polynomial order of coordinate basis, or `PETSC_DECIDE` for default
505   @param[in]   q_extra         Additional quadrature order
506   @param[in]   num_fields      Number of fields in solution vector
507   @param[in]   field_sizes     Array of number of components for each field
508   @param[out]  dm              `DM` to setup
509 
510   @return An error code: 0 - success, otherwise - failure
511 **/
512 PetscErrorCode DMSetupByOrder_FEM(PetscBool setup_faces, PetscBool setup_coords, PetscInt degree, PetscInt coord_order, PetscInt q_extra,
513                                   PetscInt num_fields, const PetscInt *field_sizes, DM dm) {
514   PetscFunctionBeginUser;
515   PetscCall(DMSetupByOrderBegin_FEM(setup_faces, setup_coords, degree, coord_order, q_extra, num_fields, field_sizes, dm));
516   PetscCall(DMSetupByOrderEnd_FEM(setup_coords, dm));
517   PetscFunctionReturn(PETSC_SUCCESS);
518 }
519