1 #include "../include/petscutils.h" 2 3 // ----------------------------------------------------------------------------- 4 // Convert PETSc MemType to libCEED MemType 5 // ----------------------------------------------------------------------------- 6 CeedMemType MemTypeP2C(PetscMemType mtype) { 7 return PetscMemTypeDevice(mtype) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 8 } 9 10 // ----------------------------------------------------------------------------- 11 // Utility function taken from petsc/src/dm/impls/plex/examples/tutorials/ex7.c 12 // ----------------------------------------------------------------------------- 13 PetscErrorCode ProjectToUnitSphere(DM dm) { 14 Vec coordinates; 15 PetscScalar *coords; 16 PetscInt Nv, v, dim, d; 17 PetscErrorCode ierr; 18 19 PetscFunctionBeginUser; 20 ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr); 21 ierr = VecGetLocalSize(coordinates, &Nv); CHKERRQ(ierr); 22 ierr = VecGetBlockSize(coordinates, &dim); CHKERRQ(ierr); 23 Nv /= dim; 24 ierr = VecGetArray(coordinates, &coords); CHKERRQ(ierr); 25 for (v = 0; v < Nv; ++v) { 26 PetscReal r = 0.0; 27 28 for (d = 0; d < dim; ++d) r += PetscSqr(PetscRealPart(coords[v*dim+d])); 29 r = PetscSqrtReal(r); 30 for (d = 0; d < dim; ++d) coords[v*dim+d] /= r; 31 } 32 ierr = VecRestoreArray(coordinates, &coords); CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 }; 35 36 // ----------------------------------------------------------------------------- 37 // PETSc FE Boilerplate 38 // ----------------------------------------------------------------------------- 39 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 40 PetscBool isSimplex, const char prefix[], 41 PetscInt order, PetscFE *fem) { 42 PetscQuadrature q, fq; 43 DM K; 44 PetscSpace P; 45 PetscDualSpace Q; 46 PetscInt quadPointsPerEdge; 47 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 48 PetscErrorCode ierr; 49 50 PetscFunctionBeginUser; 51 /* Create space */ 52 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P); CHKERRQ(ierr); 53 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix); CHKERRQ(ierr); 54 ierr = PetscSpacePolynomialSetTensor(P, tensor); CHKERRQ(ierr); 55 ierr = PetscSpaceSetFromOptions(P); CHKERRQ(ierr); 56 ierr = PetscSpaceSetNumComponents(P, Nc); CHKERRQ(ierr); 57 ierr = PetscSpaceSetNumVariables(P, dim); CHKERRQ(ierr); 58 ierr = PetscSpaceSetDegree(P, order, order); CHKERRQ(ierr); 59 ierr = PetscSpaceSetUp(P); CHKERRQ(ierr); 60 ierr = PetscSpacePolynomialGetTensor(P, &tensor); CHKERRQ(ierr); 61 /* Create dual space */ 62 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q); 63 CHKERRQ(ierr); 64 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE); CHKERRQ(ierr); 65 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix); CHKERRQ(ierr); 66 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K); CHKERRQ(ierr); 67 ierr = PetscDualSpaceSetDM(Q, K); CHKERRQ(ierr); 68 ierr = DMDestroy(&K); CHKERRQ(ierr); 69 ierr = PetscDualSpaceSetNumComponents(Q, Nc); CHKERRQ(ierr); 70 ierr = PetscDualSpaceSetOrder(Q, order); CHKERRQ(ierr); 71 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor); CHKERRQ(ierr); 72 ierr = PetscDualSpaceSetFromOptions(Q); CHKERRQ(ierr); 73 ierr = PetscDualSpaceSetUp(Q); CHKERRQ(ierr); 74 /* Create element */ 75 ierr = PetscFECreate(PetscObjectComm((PetscObject) dm), fem); CHKERRQ(ierr); 76 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix); CHKERRQ(ierr); 77 ierr = PetscFESetFromOptions(*fem); CHKERRQ(ierr); 78 ierr = PetscFESetBasisSpace(*fem, P); CHKERRQ(ierr); 79 ierr = PetscFESetDualSpace(*fem, Q); CHKERRQ(ierr); 80 ierr = PetscFESetNumComponents(*fem, Nc); CHKERRQ(ierr); 81 ierr = PetscFESetUp(*fem); CHKERRQ(ierr); 82 ierr = PetscSpaceDestroy(&P); CHKERRQ(ierr); 83 ierr = PetscDualSpaceDestroy(&Q); CHKERRQ(ierr); 84 /* Create quadrature */ 85 quadPointsPerEdge = PetscMax(order + 1,1); 86 if (isSimplex) { 87 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, 88 &q); CHKERRQ(ierr); 89 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, 90 &fq); CHKERRQ(ierr); 91 } else { 92 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, 93 &q); CHKERRQ(ierr); 94 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, 95 &fq); CHKERRQ(ierr); 96 } 97 ierr = PetscFESetQuadrature(*fem, q); CHKERRQ(ierr); 98 ierr = PetscFESetFaceQuadrature(*fem, fq); CHKERRQ(ierr); 99 ierr = PetscQuadratureDestroy(&q); CHKERRQ(ierr); 100 ierr = PetscQuadratureDestroy(&fq); CHKERRQ(ierr); 101 102 PetscFunctionReturn(0); 103 }; 104 105 // ----------------------------------------------------------------------------- 106 // Create BC label 107 // ----------------------------------------------------------------------------- 108 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) { 109 int ierr; 110 DMLabel label; 111 112 PetscFunctionBeginUser; 113 114 ierr = DMCreateLabel(dm, name); CHKERRQ(ierr); 115 ierr = DMGetLabel(dm, name, &label); CHKERRQ(ierr); 116 ierr = DMPlexMarkBoundaryFaces(dm, 1, label); CHKERRQ(ierr); 117 ierr = DMPlexLabelComplete(dm, label); CHKERRQ(ierr); 118 119 PetscFunctionReturn(0); 120 }; 121 122 // ----------------------------------------------------------------------------- 123 // This function sets up a DM for a given degree 124 // ----------------------------------------------------------------------------- 125 PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncompu, 126 PetscInt dim, bool enforcebc, BCFunction bcsfunc) { 127 PetscInt ierr, marker_ids[1] = {1}; 128 PetscFE fe; 129 130 PetscFunctionBeginUser; 131 132 // Setup FE 133 ierr = PetscFECreateByDegree(dm, dim, ncompu, PETSC_FALSE, NULL, degree, &fe); 134 CHKERRQ(ierr); 135 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 136 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 137 138 // Setup DM 139 ierr = DMCreateDS(dm); CHKERRQ(ierr); 140 if (enforcebc) { 141 PetscBool hasLabel; 142 DMHasLabel(dm, "marker", &hasLabel); 143 if (!hasLabel) {CreateBCLabel(dm, "marker");} 144 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, 145 (void(*)(void))bcsfunc, NULL, 1, marker_ids, NULL); 146 CHKERRQ(ierr); 147 } 148 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 149 CHKERRQ(ierr); 150 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 151 152 PetscFunctionReturn(0); 153 }; 154 155 // ----------------------------------------------------------------------------- 156 // Utility function - essential BC dofs are encoded in closure indices as -(i+1) 157 // ----------------------------------------------------------------------------- 158 PetscInt Involute(PetscInt i) { 159 return i >= 0 ? i : -(i + 1); 160 }; 161 162 // ----------------------------------------------------------------------------- 163 // Get CEED restriction data from DMPlex 164 // ----------------------------------------------------------------------------- 165 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 166 CeedInt topodim, CeedInt height, DMLabel domainLabel, CeedInt value, 167 CeedElemRestriction *Erestrict) { 168 PetscSection section; 169 PetscInt p, Nelem, Ndof, *erestrict, eoffset, nfields, dim, depth; 170 DMLabel depthLabel; 171 IS depthIS, iterIS; 172 Vec Uloc; 173 const PetscInt *iterIndices; 174 PetscErrorCode ierr; 175 176 PetscFunctionBeginUser; 177 178 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 179 dim -= height; 180 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 181 ierr = PetscSectionGetNumFields(section, &nfields); CHKERRQ(ierr); 182 PetscInt ncomp[nfields], fieldoff[nfields+1]; 183 fieldoff[0] = 0; 184 for (PetscInt f = 0; f < nfields; f++) { 185 ierr = PetscSectionGetFieldComponents(section, f, &ncomp[f]); CHKERRQ(ierr); 186 fieldoff[f+1] = fieldoff[f] + ncomp[f]; 187 } 188 189 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 190 ierr = DMPlexGetDepthLabel(dm, &depthLabel); CHKERRQ(ierr); 191 ierr = DMLabelGetStratumIS(depthLabel, depth - height, &depthIS); CHKERRQ(ierr); 192 if (domainLabel) { 193 IS domainIS; 194 ierr = DMLabelGetStratumIS(domainLabel, value, &domainIS); CHKERRQ(ierr); 195 if (domainIS) { // domainIS is non-empty 196 ierr = ISIntersect(depthIS, domainIS, &iterIS); CHKERRQ(ierr); 197 ierr = ISDestroy(&domainIS); CHKERRQ(ierr); 198 } else { // domainIS is NULL (empty) 199 iterIS = NULL; 200 } 201 ierr = ISDestroy(&depthIS); CHKERRQ(ierr); 202 } else { 203 iterIS = depthIS; 204 } 205 if (iterIS) { 206 ierr = ISGetLocalSize(iterIS, &Nelem); CHKERRQ(ierr); 207 ierr = ISGetIndices(iterIS, &iterIndices); CHKERRQ(ierr); 208 } else { 209 Nelem = 0; 210 iterIndices = NULL; 211 } 212 ierr = PetscMalloc1(Nelem*PetscPowInt(P, topodim), &erestrict); CHKERRQ(ierr); 213 for (p = 0, eoffset = 0; p < Nelem; p++) { 214 PetscInt c = iterIndices[p]; 215 PetscInt numindices, *indices, nnodes; 216 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, 217 &numindices, &indices, NULL, NULL); 218 CHKERRQ(ierr); 219 bool flip = false; 220 if (height > 0) { 221 PetscInt numCells, numFaces, start = -1; 222 const PetscInt *orients, *faces, *cells; 223 ierr = DMPlexGetSupport(dm, c, &cells); CHKERRQ(ierr); 224 ierr = DMPlexGetSupportSize(dm, c, &numCells); CHKERRQ(ierr); 225 if (numCells != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 226 "Expected one cell in support of exterior face, but got %D cells", 227 numCells); 228 ierr = DMPlexGetCone(dm, cells[0], &faces); CHKERRQ(ierr); 229 ierr = DMPlexGetConeSize(dm, cells[0], &numFaces); CHKERRQ(ierr); 230 for (PetscInt i=0; i<numFaces; i++) {if (faces[i] == c) start = i;} 231 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, 232 "Could not find face %D in cone of its support", 233 c); 234 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients); CHKERRQ(ierr); 235 if (orients[start] < 0) flip = true; 236 } 237 if (numindices % fieldoff[nfields]) SETERRQ1(PETSC_COMM_SELF, 238 PETSC_ERR_ARG_INCOMP, "Number of closure indices not compatible with Cell %D", 239 c); 240 nnodes = numindices / fieldoff[nfields]; 241 for (PetscInt i = 0; i < nnodes; i++) { 242 PetscInt ii = i; 243 if (flip) { 244 if (P == nnodes) ii = nnodes - 1 - i; 245 else if (P*P == nnodes) { 246 PetscInt row = i / P, col = i % P; 247 ii = row + col * P; 248 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, 249 "No support for flipping point with %D nodes != P (%D) or P^2", 250 nnodes, P); 251 } 252 // Check that indices are blocked by node and thus can be coalesced as a single field with 253 // fieldoff[nfields] = sum(ncomp) components. 254 for (PetscInt f = 0; f < nfields; f++) { 255 for (PetscInt j = 0; j < ncomp[f]; j++) { 256 if (Involute(indices[fieldoff[f]*nnodes + ii*ncomp[f] + j]) 257 != Involute(indices[ii*ncomp[0]]) + fieldoff[f] + j) 258 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 259 "Cell %D closure indices not interlaced for node %D field %D component %D", 260 c, ii, f, j); 261 } 262 } 263 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 264 PetscInt loc = Involute(indices[ii*ncomp[0]]); 265 erestrict[eoffset++] = loc; 266 } 267 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, 268 &numindices, &indices, NULL, NULL); 269 CHKERRQ(ierr); 270 } 271 if (eoffset != Nelem*PetscPowInt(P, topodim)) 272 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 273 "ElemRestriction of size (%D,%D) initialized %D nodes", Nelem, 274 PetscPowInt(P, topodim),eoffset); 275 if (iterIS) { 276 ierr = ISRestoreIndices(iterIS, &iterIndices); CHKERRQ(ierr); 277 } 278 ierr = ISDestroy(&iterIS); CHKERRQ(ierr); 279 280 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 281 ierr = VecGetLocalSize(Uloc, &Ndof); CHKERRQ(ierr); 282 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 283 CeedElemRestrictionCreate(ceed, Nelem, PetscPowInt(P, topodim), 284 fieldoff[nfields], 285 1, Ndof, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 286 Erestrict); 287 ierr = PetscFree(erestrict); CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 }; 290 291 // ----------------------------------------------------------------------------- 292