1 // Copyright (c) 2017-2022, 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 // Utility function to create local CEED restriction 18 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt label_value, PetscInt dm_field, 19 CeedElemRestriction *elem_restr) { 20 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets_petsc; 21 CeedInt *elem_restr_offsets_ceed; 22 23 PetscFunctionBeginUser; 24 PetscCall( 25 DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets_petsc)); 26 27 PetscCall(IntArrayP2C(num_elem * elem_size, &elem_restr_offsets_petsc, &elem_restr_offsets_ceed)); 28 PetscCallCeed(ceed, CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, 29 elem_restr_offsets_ceed, elem_restr)); 30 PetscCall(PetscFree(elem_restr_offsets_ceed)); 31 PetscFunctionReturn(PETSC_SUCCESS); 32 } 33 34 // Utility function to get Ceed Restriction for each domain 35 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt label_value, PetscInt dm_field, CeedInt Q, 36 CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, 37 CeedElemRestriction *elem_restr_qd_i) { 38 DM dm_coord; 39 CeedInt loc_num_elem; 40 PetscInt dim; 41 CeedElemRestriction elem_restr_tmp; 42 43 PetscFunctionBeginUser; 44 PetscCall(DMGetDimension(dm, &dim)); 45 dim -= height; 46 PetscCall(CreateRestrictionFromPlex(ceed, dm, height, domain_label, label_value, dm_field, &elem_restr_tmp)); 47 if (elem_restr_q) *elem_restr_q = elem_restr_tmp; 48 if (elem_restr_x) { 49 PetscCall(DMGetCellCoordinateDM(dm, &dm_coord)); 50 if (!dm_coord) { 51 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 52 } 53 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); 54 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, height, domain_label, label_value, dm_field, elem_restr_x)); 55 } 56 if (elem_restr_qd_i) { 57 PetscCallCeed(ceed, CeedElemRestrictionGetNumElements(elem_restr_tmp, &loc_num_elem)); 58 PetscCallCeed(ceed, CeedElemRestrictionCreateStrided(ceed, loc_num_elem, Q, q_data_size, q_data_size * loc_num_elem * Q, CEED_STRIDES_BACKEND, 59 elem_restr_qd_i)); 60 } 61 if (!elem_restr_q) PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_tmp)); 62 PetscFunctionReturn(PETSC_SUCCESS); 63 } 64 65 /** 66 @brief Convert `DM` field to `DS` field. 67 68 @param[in] dm `DM` holding mesh 69 @param[in] domain_label Label for `DM` domain 70 @param[in] dm_field Index of `DM` field 71 @param[out] ds_field Index of `DS` field 72 73 @return An error code: 0 - success, otherwise - failure 74 **/ 75 PetscErrorCode DMFieldToDSField(DM dm, DMLabel domain_label, PetscInt dm_field, PetscInt *ds_field) { 76 PetscDS ds; 77 IS field_is; 78 const PetscInt *fields; 79 PetscInt num_fields; 80 81 PetscFunctionBeginUser; 82 PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL)); 83 PetscCall(ISGetIndices(field_is, &fields)); 84 PetscCall(ISGetSize(field_is, &num_fields)); 85 for (PetscInt i = 0; i < num_fields; i++) { 86 if (dm_field == fields[i]) { 87 *ds_field = i; 88 break; 89 } 90 } 91 PetscCall(ISRestoreIndices(field_is, &fields)); 92 93 if (*ds_field == -1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 // ----------------------------------------------------------------------------- 98 // Utility function - convert from DMPolytopeType to CeedElemTopology 99 // ----------------------------------------------------------------------------- 100 static CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) { 101 switch (cell_type) { 102 case DM_POLYTOPE_TRIANGLE: 103 return CEED_TOPOLOGY_TRIANGLE; 104 case DM_POLYTOPE_QUADRILATERAL: 105 return CEED_TOPOLOGY_QUAD; 106 case DM_POLYTOPE_TETRAHEDRON: 107 return CEED_TOPOLOGY_TET; 108 case DM_POLYTOPE_HEXAHEDRON: 109 return CEED_TOPOLOGY_HEX; 110 default: 111 return 0; 112 } 113 } 114 115 // ----------------------------------------------------------------------------- 116 // Create libCEED Basis from PetscTabulation 117 // ----------------------------------------------------------------------------- 118 PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt face, PetscFE fe, 119 PetscTabulation basis_tabulation, PetscQuadrature quadrature, CeedBasis *basis) { 120 PetscInt first_point; 121 PetscInt ids[1] = {label_value}; 122 DMLabel depth_label; 123 DMPolytopeType cell_type; 124 CeedElemTopology elem_topo; 125 PetscScalar *q_points, *interp, *grad; 126 const PetscScalar *q_weights; 127 PetscDualSpace dual_space; 128 PetscInt num_dual_basis_vectors; 129 PetscInt dim, num_comp, P, Q; 130 131 PetscFunctionBeginUser; 132 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 133 PetscCall(PetscFEGetNumComponents(fe, &num_comp)); 134 PetscCall(PetscFEGetDualSpace(fe, &dual_space)); 135 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors)); 136 P = num_dual_basis_vectors / num_comp; 137 138 // Use depth label if no domain label present 139 if (!domain_label) { 140 PetscInt depth; 141 142 PetscCall(DMPlexGetDepth(dm, &depth)); 143 PetscCall(DMPlexGetDepthLabel(dm, &depth_label)); 144 ids[0] = depth - height; 145 } 146 147 // Get cell interp, grad, and quadrature data 148 PetscCall(DMGetFirstLabeledPoint(dm, dm, domain_label ? domain_label : depth_label, 1, ids, height, &first_point, NULL)); 149 PetscCall(DMPlexGetCellType(dm, first_point, &cell_type)); 150 elem_topo = ElemTopologyP2C(cell_type); 151 if (!elem_topo) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMPlex topology not supported"); 152 { 153 size_t q_points_size; 154 const PetscScalar *q_points_petsc; 155 PetscInt q_dim; 156 157 PetscCall(PetscQuadratureGetData(quadrature, &q_dim, NULL, &Q, &q_points_petsc, &q_weights)); 158 q_points_size = Q * dim * sizeof(CeedScalar); 159 PetscCall(PetscCalloc(q_points_size, &q_points)); 160 for (PetscInt q = 0; q < Q; q++) { 161 for (PetscInt d = 0; d < q_dim; d++) q_points[q * dim + d] = q_points_petsc[q * q_dim + d]; 162 } 163 } 164 165 { // Convert to libCEED orientation 166 PetscBool is_simplex = PETSC_FALSE; 167 IS permutation = NULL; 168 const PetscInt *permutation_indices; 169 170 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 171 if (!is_simplex) { 172 PetscSection section; 173 174 // -- Get permutation 175 PetscCall(DMGetLocalSection(dm, §ion)); 176 PetscCall(PetscSectionGetClosurePermutation(section, (PetscObject)dm, dim, num_comp * P, &permutation)); 177 PetscCall(ISGetIndices(permutation, &permutation_indices)); 178 } 179 180 // -- Copy interp, grad matrices 181 PetscCall(PetscCalloc(P * Q * sizeof(CeedScalar), &interp)); 182 PetscCall(PetscCalloc(P * Q * dim * sizeof(CeedScalar), &grad)); 183 const CeedInt c = 0; 184 for (CeedInt q = 0; q < Q; q++) { 185 for (CeedInt p_ceed = 0; p_ceed < P; p_ceed++) { 186 CeedInt p_petsc = is_simplex ? (p_ceed * num_comp) : permutation_indices[p_ceed * num_comp]; 187 188 interp[q * P + p_ceed] = basis_tabulation->T[0][((face * Q + q) * P * num_comp + p_petsc) * num_comp + c]; 189 for (CeedInt d = 0; d < dim; d++) { 190 grad[(d * Q + q) * P + p_ceed] = basis_tabulation->T[1][(((face * Q + q) * P * num_comp + p_petsc) * num_comp + c) * dim + d]; 191 } 192 } 193 } 194 195 // -- Cleanup 196 if (permutation) PetscCall(ISRestoreIndices(permutation, &permutation_indices)); 197 PetscCall(ISDestroy(&permutation)); 198 } 199 200 // Finally, create libCEED basis 201 PetscCallCeed(ceed, CeedBasisCreateH1(ceed, elem_topo, num_comp, P, Q, interp, grad, q_points, q_weights, basis)); 202 PetscCall(PetscFree(q_points)); 203 PetscCall(PetscFree(interp)); 204 PetscCall(PetscFree(grad)); 205 PetscFunctionReturn(PETSC_SUCCESS); 206 } 207 208 // ----------------------------------------------------------------------------- 209 // Get CEED Basis from DMPlex 210 // ----------------------------------------------------------------------------- 211 PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, CeedBasis *basis) { 212 PetscDS ds; 213 PetscFE fe; 214 PetscQuadrature quadrature; 215 PetscBool is_simplex = PETSC_TRUE; 216 PetscInt ds_field = -1; 217 218 PetscFunctionBeginUser; 219 // Get element information 220 PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL)); 221 PetscCall(DMFieldToDSField(dm, domain_label, dm_field, &ds_field)); 222 PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe)); 223 PetscCall(PetscFEGetHeightSubspace(fe, height, &fe)); 224 PetscCall(PetscFEGetQuadrature(fe, &quadrature)); 225 226 // Check if simplex or tensor-product mesh 227 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 228 229 // Build libCEED basis 230 if (is_simplex) { 231 PetscTabulation basis_tabulation; 232 PetscInt num_derivatives = 1, face = 0; 233 234 PetscCall(PetscFEGetCellTabulation(fe, num_derivatives, &basis_tabulation)); 235 PetscCall(BasisCreateFromTabulation(ceed, dm, domain_label, label_value, height, face, fe, basis_tabulation, quadrature, basis)); 236 } else { 237 PetscDualSpace dual_space; 238 PetscInt num_dual_basis_vectors; 239 PetscInt dim, num_comp, P, Q; 240 241 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 242 PetscCall(PetscFEGetNumComponents(fe, &num_comp)); 243 PetscCall(PetscFEGetDualSpace(fe, &dual_space)); 244 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors)); 245 P = num_dual_basis_vectors / num_comp; 246 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &Q, NULL, NULL)); 247 248 CeedInt P_1d = (CeedInt)round(pow(P, 1.0 / dim)); 249 CeedInt Q_1d = (CeedInt)round(pow(Q, 1.0 / dim)); 250 251 PetscCallCeed(ceed, CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, CEED_GAUSS, basis)); 252 } 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255 256 /** 257 @brief Setup `DM` with FE space of appropriate degree 258 259 Must be followed by `DMSetupByOrderEnd_FEM` 260 261 @param[in] setup_faces Flag to setup face geometry 262 @param[in] setup_coords Flag to setup coordinate spaces 263 @param[in] degree Polynomial orders of field 264 @param[in] coord_order Polynomial order of coordinate basis, or `RATEL_DECIDE` for default 265 @param[in] q_extra Additional quadrature order, or `RATEL_DECIDE` for default 266 @param[in] num_fields Number of fields in solution vector 267 @param[in] field_sizes Array of number of components for each field 268 @param[out] dm `DM` to setup 269 270 @return An error code: 0 - success, otherwise - failure 271 **/ 272 PetscErrorCode DMSetupByOrderBegin_FEM(PetscBool setup_faces, PetscBool setup_coords, PetscInt degree, PetscInt coord_order, PetscInt q_extra, 273 CeedInt num_fields, const CeedInt *field_sizes, DM dm) { 274 PetscInt dim, q_order = degree + q_extra; 275 PetscBool is_simplex = PETSC_TRUE; 276 PetscFE fe; 277 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 278 279 PetscFunctionBeginUser; 280 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 281 282 // Setup DM 283 PetscCall(DMGetDimension(dm, &dim)); 284 for (PetscInt i = 0; i < num_fields; i++) { 285 PetscFE fe_face; 286 PetscInt q_order = degree + q_extra; 287 288 PetscCall(PetscFECreateLagrange(comm, dim, field_sizes[i], is_simplex, degree, q_order, &fe)); 289 if (setup_faces) PetscCall(PetscFEGetHeightSubspace(fe, 1, &fe_face)); 290 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 291 PetscCall(PetscFEDestroy(&fe)); 292 } 293 PetscCall(DMCreateDS(dm)); 294 295 // Project coordinates to enrich quadrature space 296 if (setup_coords) { 297 DM dm_coord; 298 PetscDS ds_coord; 299 PetscFE fe_coord_current, fe_coord_new, fe_coord_face_new; 300 PetscDualSpace fe_coord_dual_space; 301 PetscInt fe_coord_order, num_comp_coord; 302 303 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 304 PetscCall(DMGetCoordinateDim(dm, &num_comp_coord)); 305 PetscCall(DMGetRegionDS(dm_coord, NULL, NULL, &ds_coord, NULL)); 306 PetscCall(PetscDSGetDiscretization(ds_coord, 0, (PetscObject *)&fe_coord_current)); 307 PetscCall(PetscFEGetDualSpace(fe_coord_current, &fe_coord_dual_space)); 308 PetscCall(PetscDualSpaceGetOrder(fe_coord_dual_space, &fe_coord_order)); 309 310 // Create FE for coordinates 311 PetscCheck(fe_coord_order == 1 || coord_order == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER_INPUT, 312 "Only linear mesh geometry supported. Recieved %d order", fe_coord_order); 313 PetscCall(PetscFECreateLagrange(comm, dim, num_comp_coord, is_simplex, fe_coord_order, q_order, &fe_coord_new)); 314 if (setup_faces) PetscCall(PetscFEGetHeightSubspace(fe_coord_new, 1, &fe_coord_face_new)); 315 PetscCall(DMProjectCoordinates(dm, fe_coord_new)); 316 PetscCall(PetscFEDestroy(&fe_coord_new)); 317 } 318 PetscFunctionReturn(PETSC_SUCCESS); 319 } 320 321 /** 322 @brief Finish setting up `DM` 323 324 Must be called after `DMSetupByOrderBegin_FEM` and all strong BCs have be declared via `DMAddBoundaries`. 325 326 @param[in] setup_coords Flag to setup coordinate spaces 327 @param[out] dm `DM` to setup 328 329 @return An error code: 0 - success, otherwise - failure 330 **/ 331 PetscErrorCode DMSetupByOrderEnd_FEM(PetscBool setup_coords, DM dm) { 332 PetscBool is_simplex; 333 334 PetscFunctionBeginUser; 335 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 336 // Set tensor permutation if needed 337 if (!is_simplex) { 338 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 339 if (setup_coords) { 340 DM dm_coord; 341 342 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 343 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); 344 } 345 } 346 PetscFunctionReturn(PETSC_SUCCESS); 347 } 348 349 /** 350 @brief Setup `DM` with FE space of appropriate degree with no boundary conditions 351 352 Calls `DMSetupByOrderBegin_FEM` and `DMSetupByOrderEnd_FEM` successively 353 354 @param[in] setup_faces Flag to setup face geometry 355 @param[in] setup_coords Flag to setup coordinate spaces 356 @param[in] degree Polynomial orders of field 357 @param[in] coord_order Polynomial order of coordinate basis, or `RATEL_DECIDE` for default 358 @param[in] q_extra Additional quadrature order, or `RATEL_DECIDE` for default 359 @param[in] num_fields Number of fields in solution vector 360 @param[in] field_sizes Array of number of components for each field 361 @param[out] dm `DM` to setup 362 363 @return An error code: 0 - success, otherwise - failure 364 **/ 365 PetscErrorCode DMSetupByOrder_FEM(PetscBool setup_faces, PetscBool setup_coords, PetscInt degree, PetscInt coord_order, PetscInt q_extra, 366 CeedInt num_fields, const CeedInt *field_sizes, DM dm) { 367 PetscFunctionBeginUser; 368 PetscCall(DMSetupByOrderBegin_FEM(setup_faces, setup_coords, degree, coord_order, q_extra, num_fields, field_sizes, dm)); 369 PetscCall(DMSetupByOrderEnd_FEM(setup_coords, dm)); 370 PetscFunctionReturn(PETSC_SUCCESS); 371 } 372