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(PETSC_SUCCESS); 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) invVscalar[(s * spdim + j) * spdim + k] += Bf[q * spdim + k] * weights[q]; 85 } 86 PetscCall(PetscFree(Bf)); 87 } 88 n = spdim; 89 PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s * spdim * spdim], &n, pivots, &info)); 90 PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s * spdim * spdim], &n, pivots, work, &n, &info)); 91 } 92 #if defined(PETSC_USE_COMPLEX) 93 for (s = 0; s < cmp->numSubelements * spdim * spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]); 94 PetscCall(PetscFree(invVscalar)); 95 #endif 96 PetscCall(PetscFree2(pivots, work)); 97 PetscFunctionReturn(PETSC_SUCCESS); 98 } 99 100 static PetscErrorCode PetscFECreateTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) 101 { 102 PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data; 103 DM dm; 104 DMPolytopeType ct; 105 PetscInt pdim; /* Dimension of FE space P */ 106 PetscInt spdim; /* Dimension of subelement FE space P */ 107 PetscInt dim; /* Spatial dimension */ 108 PetscInt comp; /* Field components */ 109 PetscInt *subpoints; 110 PetscReal *B = K >= 0 ? T->T[0] : NULL; 111 PetscReal *D = K >= 1 ? T->T[1] : NULL; 112 PetscReal *H = K >= 2 ? T->T[2] : NULL; 113 PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint; 114 PetscInt p, s, d, e, j, k; 115 116 PetscFunctionBegin; 117 PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm)); 118 PetscCall(DMGetDimension(dm, &dim)); 119 PetscCall(DMPlexGetCellType(dm, 0, &ct)); 120 PetscCall(PetscSpaceGetDimension(fem->basisSpace, &spdim)); 121 PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim)); 122 PetscCall(PetscFEGetNumComponents(fem, &comp)); 123 /* Divide points into subelements */ 124 PetscCall(DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints)); 125 PetscCall(DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint)); 126 for (p = 0; p < npoints; ++p) { 127 for (s = 0; s < cmp->numSubelements; ++s) { 128 PetscBool inside; 129 130 /* Apply transform, and check that point is inside cell */ 131 for (d = 0; d < dim; ++d) { 132 subpoint[d] = -1.0; 133 for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s * dim + d) * dim + e] * (points[p * dim + e] - cmp->v0[s * dim + e]); 134 } 135 PetscCall(DMPolytopeInCellTest(ct, subpoint, &inside)); 136 if (inside) { 137 subpoints[p] = s; 138 break; 139 } 140 } 141 PetscCheck(s < cmp->numSubelements, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " was not found in any subelement", p); 142 } 143 PetscCall(DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint)); 144 /* Evaluate the prime basis functions at all points */ 145 if (K >= 0) PetscCall(DMGetWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB)); 146 if (K >= 1) PetscCall(DMGetWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD)); 147 if (K >= 2) PetscCall(DMGetWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH)); 148 PetscCall(PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH)); 149 /* Translate to the nodal basis */ 150 if (K >= 0) PetscCall(PetscArrayzero(B, npoints * pdim * comp)); 151 if (K >= 1) PetscCall(PetscArrayzero(D, npoints * pdim * comp * dim)); 152 if (K >= 2) PetscCall(PetscArrayzero(H, npoints * pdim * comp * dim * dim)); 153 for (p = 0; p < npoints; ++p) { 154 const PetscInt s = subpoints[p]; 155 156 if (B) { 157 /* Multiply by V^{-1} (spdim x spdim) */ 158 for (j = 0; j < spdim; ++j) { 159 const PetscInt i = (p * pdim + cmp->embedding[s * spdim + j]) * comp; 160 161 B[i] = 0.0; 162 for (k = 0; k < spdim; ++k) B[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpB[p * spdim + k]; 163 } 164 } 165 if (D) { 166 /* Multiply by V^{-1} (spdim x spdim) */ 167 for (j = 0; j < spdim; ++j) { 168 for (d = 0; d < dim; ++d) { 169 const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim + d; 170 171 D[i] = 0.0; 172 for (k = 0; k < spdim; ++k) D[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpD[(p * spdim + k) * dim + d]; 173 } 174 } 175 } 176 if (H) { 177 /* Multiply by V^{-1} (pdim x pdim) */ 178 for (j = 0; j < spdim; ++j) { 179 for (d = 0; d < dim * dim; ++d) { 180 const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim * dim + d; 181 182 H[i] = 0.0; 183 for (k = 0; k < spdim; ++k) H[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpH[(p * spdim + k) * dim * dim + d]; 184 } 185 } 186 } 187 } 188 PetscCall(DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints)); 189 if (K >= 0) PetscCall(DMRestoreWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB)); 190 if (K >= 1) PetscCall(DMRestoreWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD)); 191 if (K >= 2) PetscCall(DMRestoreWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH)); 192 PetscFunctionReturn(PETSC_SUCCESS); 193 } 194 195 static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem) 196 { 197 PetscFunctionBegin; 198 fem->ops->setfromoptions = NULL; 199 fem->ops->setup = PetscFESetUp_Composite; 200 fem->ops->view = NULL; 201 fem->ops->destroy = PetscFEDestroy_Composite; 202 fem->ops->getdimension = PetscFEGetDimension_Basic; 203 fem->ops->createtabulation = PetscFECreateTabulation_Composite; 204 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 205 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 206 fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */; 207 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 208 PetscFunctionReturn(PETSC_SUCCESS); 209 } 210 211 /*MC 212 PETSCFECOMPOSITE = "composite" - A `PetscFEType` that represents a composite element 213 214 Level: intermediate 215 216 .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()` 217 M*/ 218 PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem) 219 { 220 PetscFE_Composite *cmp; 221 222 PetscFunctionBegin; 223 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 224 PetscCall(PetscNew(&cmp)); 225 fem->data = cmp; 226 227 cmp->numSubelements = -1; 228 cmp->v0 = NULL; 229 cmp->jac = NULL; 230 231 PetscCall(PetscFEInitialize_Composite(fem)); 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 /*@C 236 PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement 237 238 Not Collective 239 240 Input Parameter: 241 . fem - The `PetscFE` object 242 243 Output Parameters: 244 + numSubelements - The number of sub elements 245 . v0 - The affine transformation for each element 246 . jac - The Jacobian for each element 247 - invjac - The inverse of the Jacobian 248 249 Level: intermediate 250 251 .seealso: `PetscFE`, `PetscFECreate()` 252 @*/ 253 PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[]) 254 { 255 PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 259 if (numSubelements) { 260 PetscAssertPointer(numSubelements, 2); 261 *numSubelements = cmp->numSubelements; 262 } 263 if (v0) { 264 PetscAssertPointer(v0, 3); 265 *v0 = cmp->v0; 266 } 267 if (jac) { 268 PetscAssertPointer(jac, 4); 269 *jac = cmp->jac; 270 } 271 if (invjac) { 272 PetscAssertPointer(invjac, 5); 273 *invjac = cmp->invjac; 274 } 275 PetscFunctionReturn(PETSC_SUCCESS); 276 } 277