Actual source code: fecomposite.c

  1: #include <petsc/private/petscfeimpl.h>
  2: #include <petsc/private/dtimpl.h>
  3: #include <petscblaslapack.h>
  4: #include <petscdmplextransform.h>

  6: static PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
  7: {
  8:   PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;

 10:   PetscFunctionBegin;
 11:   PetscCall(PetscFree(cmp->embedding));
 12:   PetscCall(PetscFree(cmp));
 13:   PetscFunctionReturn(PETSC_SUCCESS);
 14: }

 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;

 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, &section));
 44:   for (s = 0; s < cmp->numSubelements; ++s) {
 45:     PetscInt  sd = 0;
 46:     PetscInt  closureSize;
 47:     PetscInt *closure = NULL;

 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;

 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;

 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: }

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;

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;

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];

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;

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;

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;

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: }

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: }

211: /*MC
212:   PETSCFECOMPOSITE = "composite" - A `PetscFEType` that represents a composite element

214:   Level: intermediate

216: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
217: M*/
218: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
219: {
220:   PetscFE_Composite *cmp;

222:   PetscFunctionBegin;
224:   PetscCall(PetscNew(&cmp));
225:   fem->data = cmp;

227:   cmp->numSubelements = -1;
228:   cmp->v0             = NULL;
229:   cmp->jac            = NULL;

231:   PetscCall(PetscFEInitialize_Composite(fem));
232:   PetscFunctionReturn(PETSC_SUCCESS);
233: }

235: /*@C
236:   PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement

238:   Not Collective

240:   Input Parameter:
241: . fem - The `PetscFE` object

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

249:   Level: intermediate

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;

257:   PetscFunctionBegin;
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: }