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