1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petsc/private/dtimpl.h> /*I "petscdt.h" I*/ 3 #include <petscblaslapack.h> 4 #include <petscdmplextransform.h> 5 6 static PetscErrorCode PetscFEDestroy_Composite(PetscFE fem) 7 { 8 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 9 10 PetscFunctionBegin; 11 PetscCall(PetscFree(cmp->embedding)); 12 PetscCall(PetscFree(cmp)); 13 PetscFunctionReturn(0); 14 } 15 16 static PetscErrorCode PetscFESetUp_Composite(PetscFE fem) 17 { 18 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 19 DM K; 20 DMPolytopeType ct; 21 DMPlexTransform tr; 22 PetscReal *subpoint; 23 PetscBLASInt *pivots; 24 PetscBLASInt n, info; 25 PetscScalar *work, *invVscalar; 26 PetscInt dim, pdim, spdim, j, s; 27 PetscSection section; 28 29 PetscFunctionBegin; 30 /* Get affine mapping from reference cell to each subcell */ 31 PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &K)); 32 PetscCall(DMGetDimension(K, &dim)); 33 PetscCall(DMPlexGetCellType(K, 0, &ct)); 34 PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr)); 35 PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR)); 36 PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac)); 37 PetscCall(DMPlexTransformDestroy(&tr)); 38 /* Determine dof embedding into subelements */ 39 PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim)); 40 PetscCall(PetscSpaceGetDimension(fem->basisSpace, &spdim)); 41 PetscCall(PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding)); 42 PetscCall(DMGetWorkArray(K, dim, MPIU_REAL, &subpoint)); 43 PetscCall(PetscDualSpaceGetSection(fem->dualSpace, §ion)); 44 for (s = 0; s < cmp->numSubelements; ++s) { 45 PetscInt sd = 0; 46 PetscInt closureSize; 47 PetscInt *closure = NULL; 48 49 PetscCall(DMPlexGetTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure)); 50 for (j = 0; j < closureSize; j++) { 51 PetscInt point = closure[2*j]; 52 PetscInt dof, off, k; 53 54 PetscCall(PetscSectionGetDof(section, point, &dof)); 55 PetscCall(PetscSectionGetOffset(section, point, &off)); 56 for (k = 0; k < dof; k++) cmp->embedding[s*spdim+sd++] = off + k; 57 } 58 PetscCall(DMPlexRestoreTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure)); 59 PetscCheck(sd == spdim,PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %" PetscInt_FMT " has %" PetscInt_FMT " dual basis vectors != %" PetscInt_FMT, s, sd, spdim); 60 } 61 PetscCall(DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint)); 62 /* Construct the change of basis from prime basis to nodal basis for each subelement */ 63 PetscCall(PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV)); 64 PetscCall(PetscMalloc2(spdim,&pivots,spdim,&work)); 65 #if defined(PETSC_USE_COMPLEX) 66 PetscCall(PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar)); 67 #else 68 invVscalar = fem->invV; 69 #endif 70 for (s = 0; s < cmp->numSubelements; ++s) { 71 for (j = 0; j < spdim; ++j) { 72 PetscReal *Bf; 73 PetscQuadrature f; 74 const PetscReal *points, *weights; 75 PetscInt Nc, Nq, q, k; 76 77 PetscCall(PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f)); 78 PetscCall(PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights)); 79 PetscCall(PetscMalloc1(f->numPoints*spdim*Nc,&Bf)); 80 PetscCall(PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL)); 81 for (k = 0; k < spdim; ++k) { 82 /* n_j \cdot \phi_k */ 83 invVscalar[(s*spdim + j)*spdim+k] = 0.0; 84 for (q = 0; q < Nq; ++q) { 85 invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*weights[q]; 86 } 87 } 88 PetscCall(PetscFree(Bf)); 89 } 90 n = spdim; 91 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info)); 92 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info)); 93 } 94 #if defined(PETSC_USE_COMPLEX) 95 for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]); 96 PetscCall(PetscFree(invVscalar)); 97 #endif 98 PetscCall(PetscFree2(pivots,work)); 99 PetscFunctionReturn(0); 100 } 101 102 static PetscErrorCode PetscFECreateTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) 103 { 104 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 105 DM dm; 106 DMPolytopeType ct; 107 PetscInt pdim; /* Dimension of FE space P */ 108 PetscInt spdim; /* Dimension of subelement FE space P */ 109 PetscInt dim; /* Spatial dimension */ 110 PetscInt comp; /* Field components */ 111 PetscInt *subpoints; 112 PetscReal *B = K >= 0 ? T->T[0] : NULL; 113 PetscReal *D = K >= 1 ? T->T[1] : NULL; 114 PetscReal *H = K >= 2 ? T->T[2] : NULL; 115 PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint; 116 PetscInt p, s, d, e, j, k; 117 118 PetscFunctionBegin; 119 PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm)); 120 PetscCall(DMGetDimension(dm, &dim)); 121 PetscCall(DMPlexGetCellType(dm, 0, &ct)); 122 PetscCall(PetscSpaceGetDimension(fem->basisSpace, &spdim)); 123 PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim)); 124 PetscCall(PetscFEGetNumComponents(fem, &comp)); 125 /* Divide points into subelements */ 126 PetscCall(DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints)); 127 PetscCall(DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint)); 128 for (p = 0; p < npoints; ++p) { 129 for (s = 0; s < cmp->numSubelements; ++s) { 130 PetscBool inside; 131 132 /* Apply transform, and check that point is inside cell */ 133 for (d = 0; d < dim; ++d) { 134 subpoint[d] = -1.0; 135 for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]); 136 } 137 PetscCall(DMPolytopeInCellTest(ct, subpoint, &inside)); 138 if (inside) {subpoints[p] = s; break;} 139 } 140 PetscCheck(s < cmp->numSubelements,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " was not found in any subelement", p); 141 } 142 PetscCall(DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint)); 143 /* Evaluate the prime basis functions at all points */ 144 if (K >= 0) PetscCall(DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB)); 145 if (K >= 1) PetscCall(DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD)); 146 if (K >= 2) PetscCall(DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH)); 147 PetscCall(PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH)); 148 /* Translate to the nodal basis */ 149 if (K >= 0) PetscCall(PetscArrayzero(B, npoints*pdim*comp)); 150 if (K >= 1) PetscCall(PetscArrayzero(D, npoints*pdim*comp*dim)); 151 if (K >= 2) PetscCall(PetscArrayzero(H, npoints*pdim*comp*dim*dim)); 152 for (p = 0; p < npoints; ++p) { 153 const PetscInt s = subpoints[p]; 154 155 if (B) { 156 /* Multiply by V^{-1} (spdim x spdim) */ 157 for (j = 0; j < spdim; ++j) { 158 const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp; 159 160 B[i] = 0.0; 161 for (k = 0; k < spdim; ++k) { 162 B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k]; 163 } 164 } 165 } 166 if (D) { 167 /* Multiply by V^{-1} (spdim x spdim) */ 168 for (j = 0; j < spdim; ++j) { 169 for (d = 0; d < dim; ++d) { 170 const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d; 171 172 D[i] = 0.0; 173 for (k = 0; k < spdim; ++k) { 174 D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d]; 175 } 176 } 177 } 178 } 179 if (H) { 180 /* Multiply by V^{-1} (pdim x pdim) */ 181 for (j = 0; j < spdim; ++j) { 182 for (d = 0; d < dim*dim; ++d) { 183 const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d; 184 185 H[i] = 0.0; 186 for (k = 0; k < spdim; ++k) { 187 H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d]; 188 } 189 } 190 } 191 } 192 } 193 PetscCall(DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints)); 194 if (K >= 0) PetscCall(DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB)); 195 if (K >= 1) PetscCall(DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD)); 196 if (K >= 2) PetscCall(DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH)); 197 PetscFunctionReturn(0); 198 } 199 200 static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem) 201 { 202 PetscFunctionBegin; 203 fem->ops->setfromoptions = NULL; 204 fem->ops->setup = PetscFESetUp_Composite; 205 fem->ops->view = NULL; 206 fem->ops->destroy = PetscFEDestroy_Composite; 207 fem->ops->getdimension = PetscFEGetDimension_Basic; 208 fem->ops->createtabulation = PetscFECreateTabulation_Composite; 209 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 210 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 211 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 212 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 213 PetscFunctionReturn(0); 214 } 215 216 /*MC 217 PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element 218 219 Level: intermediate 220 221 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 222 M*/ 223 PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem) 224 { 225 PetscFE_Composite *cmp; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 229 PetscCall(PetscNewLog(fem, &cmp)); 230 fem->data = cmp; 231 232 cmp->numSubelements = -1; 233 cmp->v0 = NULL; 234 cmp->jac = NULL; 235 236 PetscCall(PetscFEInitialize_Composite(fem)); 237 PetscFunctionReturn(0); 238 } 239 240 /*@C 241 PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement 242 243 Not collective 244 245 Input Parameter: 246 . fem - The PetscFE object 247 248 Output Parameters: 249 + blockSize - The number of elements in a block 250 . numBlocks - The number of blocks in a batch 251 . batchSize - The number of elements in a batch 252 - numBatches - The number of batches in a chunk 253 254 Level: intermediate 255 256 .seealso: PetscFECreate() 257 @*/ 258 PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[]) 259 { 260 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 261 262 PetscFunctionBegin; 263 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 264 if (numSubelements) {PetscValidIntPointer(numSubelements, 2); *numSubelements = cmp->numSubelements;} 265 if (v0) {PetscValidPointer(v0, 3); *v0 = cmp->v0;} 266 if (jac) {PetscValidPointer(jac, 4); *jac = cmp->jac;} 267 if (invjac) {PetscValidPointer(invjac, 5); *invjac = cmp->invjac;} 268 PetscFunctionReturn(0); 269 } 270