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:   PetscFree(cmp->embedding);
 11:   PetscFree(cmp);
 12:   return 0;
 13: }

 15: static PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
 16: {
 17:   PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
 18:   DM                 K;
 19:   DMPolytopeType     ct;
 20:   DMPlexTransform    tr;
 21:   PetscReal         *subpoint;
 22:   PetscBLASInt      *pivots;
 23:   PetscBLASInt       n, info;
 24:   PetscScalar       *work, *invVscalar;
 25:   PetscInt           dim, pdim, spdim, j, s;
 26:   PetscSection       section;

 28:   /* Get affine mapping from reference cell to each subcell */
 29:   PetscDualSpaceGetDM(fem->dualSpace, &K);
 30:   DMGetDimension(K, &dim);
 31:   DMPlexGetCellType(K, 0, &ct);
 32:   DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
 33:   DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
 34:   DMPlexRefineRegularGetAffineTransforms(tr, ct, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
 35:   DMPlexTransformDestroy(&tr);
 36:   /* Determine dof embedding into subelements */
 37:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
 38:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
 39:   PetscMalloc1(cmp->numSubelements * spdim, &cmp->embedding);
 40:   DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);
 41:   PetscDualSpaceGetSection(fem->dualSpace, &section);
 42:   for (s = 0; s < cmp->numSubelements; ++s) {
 43:     PetscInt  sd = 0;
 44:     PetscInt  closureSize;
 45:     PetscInt *closure = NULL;

 47:     DMPlexGetTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);
 48:     for (j = 0; j < closureSize; j++) {
 49:       PetscInt point = closure[2 * j];
 50:       PetscInt dof, off, k;

 52:       PetscSectionGetDof(section, point, &dof);
 53:       PetscSectionGetOffset(section, point, &off);
 54:       for (k = 0; k < dof; k++) cmp->embedding[s * spdim + sd++] = off + k;
 55:     }
 56:     DMPlexRestoreTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);
 58:   }
 59:   DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint);
 60:   /* Construct the change of basis from prime basis to nodal basis for each subelement */
 61:   PetscMalloc1(cmp->numSubelements * spdim * spdim, &fem->invV);
 62:   PetscMalloc2(spdim, &pivots, spdim, &work);
 63: #if defined(PETSC_USE_COMPLEX)
 64:   PetscMalloc1(cmp->numSubelements * spdim * spdim, &invVscalar);
 65: #else
 66:   invVscalar = fem->invV;
 67: #endif
 68:   for (s = 0; s < cmp->numSubelements; ++s) {
 69:     for (j = 0; j < spdim; ++j) {
 70:       PetscReal       *Bf;
 71:       PetscQuadrature  f;
 72:       const PetscReal *points, *weights;
 73:       PetscInt         Nc, Nq, q, k;

 75:       PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s * spdim + j], &f);
 76:       PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
 77:       PetscMalloc1(f->numPoints * spdim * Nc, &Bf);
 78:       PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);
 79:       for (k = 0; k < spdim; ++k) {
 80:         /* n_j \cdot \phi_k */
 81:         invVscalar[(s * spdim + j) * spdim + k] = 0.0;
 82:         for (q = 0; q < Nq; ++q) invVscalar[(s * spdim + j) * spdim + k] += Bf[q * spdim + k] * weights[q];
 83:       }
 84:       PetscFree(Bf);
 85:     }
 86:     n = spdim;
 87:     PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s * spdim * spdim], &n, pivots, &info));
 88:     PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s * spdim * spdim], &n, pivots, work, &n, &info));
 89:   }
 90: #if defined(PETSC_USE_COMPLEX)
 91:   for (s = 0; s < cmp->numSubelements * spdim * spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
 92:   PetscFree(invVscalar);
 93: #endif
 94:   PetscFree2(pivots, work);
 95:   return 0;
 96: }

 98: static PetscErrorCode PetscFECreateTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
 99: {
100:   PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
101:   DM                 dm;
102:   DMPolytopeType     ct;
103:   PetscInt           pdim;  /* Dimension of FE space P */
104:   PetscInt           spdim; /* Dimension of subelement FE space P */
105:   PetscInt           dim;   /* Spatial dimension */
106:   PetscInt           comp;  /* Field components */
107:   PetscInt          *subpoints;
108:   PetscReal         *B    = K >= 0 ? T->T[0] : NULL;
109:   PetscReal         *D    = K >= 1 ? T->T[1] : NULL;
110:   PetscReal         *H    = K >= 2 ? T->T[2] : NULL;
111:   PetscReal         *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint;
112:   PetscInt           p, s, d, e, j, k;

114:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
115:   DMGetDimension(dm, &dim);
116:   DMPlexGetCellType(dm, 0, &ct);
117:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
118:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
119:   PetscFEGetNumComponents(fem, &comp);
120:   /* Divide points into subelements */
121:   DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);
122:   DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);
123:   for (p = 0; p < npoints; ++p) {
124:     for (s = 0; s < cmp->numSubelements; ++s) {
125:       PetscBool inside;

127:       /* Apply transform, and check that point is inside cell */
128:       for (d = 0; d < dim; ++d) {
129:         subpoint[d] = -1.0;
130:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s * dim + d) * dim + e] * (points[p * dim + e] - cmp->v0[s * dim + e]);
131:       }
132:       DMPolytopeInCellTest(ct, subpoint, &inside);
133:       if (inside) {
134:         subpoints[p] = s;
135:         break;
136:       }
137:     }
139:   }
140:   DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);
141:   /* Evaluate the prime basis functions at all points */
142:   if (K >= 0) DMGetWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB);
143:   if (K >= 1) DMGetWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD);
144:   if (K >= 2) DMGetWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH);
145:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
146:   /* Translate to the nodal basis */
147:   if (K >= 0) PetscArrayzero(B, npoints * pdim * comp);
148:   if (K >= 1) PetscArrayzero(D, npoints * pdim * comp * dim);
149:   if (K >= 2) PetscArrayzero(H, npoints * pdim * comp * dim * dim);
150:   for (p = 0; p < npoints; ++p) {
151:     const PetscInt s = subpoints[p];

153:     if (B) {
154:       /* Multiply by V^{-1} (spdim x spdim) */
155:       for (j = 0; j < spdim; ++j) {
156:         const PetscInt i = (p * pdim + cmp->embedding[s * spdim + j]) * comp;

158:         B[i] = 0.0;
159:         for (k = 0; k < spdim; ++k) B[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpB[p * spdim + k];
160:       }
161:     }
162:     if (D) {
163:       /* Multiply by V^{-1} (spdim x spdim) */
164:       for (j = 0; j < spdim; ++j) {
165:         for (d = 0; d < dim; ++d) {
166:           const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim + d;

168:           D[i] = 0.0;
169:           for (k = 0; k < spdim; ++k) D[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpD[(p * spdim + k) * dim + d];
170:         }
171:       }
172:     }
173:     if (H) {
174:       /* Multiply by V^{-1} (pdim x pdim) */
175:       for (j = 0; j < spdim; ++j) {
176:         for (d = 0; d < dim * dim; ++d) {
177:           const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim * dim + d;

179:           H[i] = 0.0;
180:           for (k = 0; k < spdim; ++k) H[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpH[(p * spdim + k) * dim * dim + d];
181:         }
182:       }
183:     }
184:   }
185:   DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);
186:   if (K >= 0) DMRestoreWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB);
187:   if (K >= 1) DMRestoreWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD);
188:   if (K >= 2) DMRestoreWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH);
189:   return 0;
190: }

192: static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
193: {
194:   fem->ops->setfromoptions          = NULL;
195:   fem->ops->setup                   = PetscFESetUp_Composite;
196:   fem->ops->view                    = NULL;
197:   fem->ops->destroy                 = PetscFEDestroy_Composite;
198:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
199:   fem->ops->createtabulation        = PetscFECreateTabulation_Composite;
200:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
201:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
202:   fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
203:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
204:   return 0;
205: }

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

210:   Level: intermediate

212: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
213: M*/
214: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
215: {
216:   PetscFE_Composite *cmp;

219:   PetscNew(&cmp);
220:   fem->data = cmp;

222:   cmp->numSubelements = -1;
223:   cmp->v0             = NULL;
224:   cmp->jac            = NULL;

226:   PetscFEInitialize_Composite(fem);
227:   return 0;
228: }

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

233:   Not collective

235:   Input Parameter:
236: . fem - The `PetscFE` object

238:   Output Parameters:
239: + blockSize - The number of elements in a block
240: . numBlocks - The number of blocks in a batch
241: . batchSize - The number of elements in a batch
242: - numBatches - The number of batches in a chunk

244:   Level: intermediate

246: .seealso: `PetscFE`, `PetscFECreate()`
247: @*/
248: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
249: {
250:   PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;

253:   if (numSubelements) {
255:     *numSubelements = cmp->numSubelements;
256:   }
257:   if (v0) {
259:     *v0 = cmp->v0;
260:   }
261:   if (jac) {
263:     *jac = cmp->jac;
264:   }
265:   if (invjac) {
267:     *invjac = cmp->invjac;
268:   }
269:   return 0;
270: }