Actual source code: fecomposite.c

petsc-master 2019-11-16
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>
  2:  #include <petsc/private/dtimpl.h>
  3:  #include <petsc/private/dmpleximpl.h>
  4:  #include <petscblaslapack.h>

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

 12:   CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
 13:   PetscFree(cmp->embedding);
 14:   PetscFree(cmp);
 15:   return(0);
 16: }

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

 30:   /* Get affine mapping from reference cell to each subcell */
 31:   PetscDualSpaceGetDM(fem->dualSpace, &K);
 32:   DMGetDimension(K, &dim);
 33:   DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
 34:   CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
 35:   /* Determine dof embedding into subelements */
 36:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
 37:   PetscSpaceGetDimension(fem->basisSpace, &spdim);
 38:   PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
 39:   DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);
 40:   for (s = 0; s < cmp->numSubelements; ++s) {
 41:     PetscInt sd = 0;

 43:     for (j = 0; j < pdim; ++j) {
 44:       PetscBool       inside;
 45:       PetscQuadrature f;
 46:       PetscInt        d, e;

 48:       PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
 49:       /* Apply transform to first point, and check that point is inside subcell */
 50:       for (d = 0; d < dim; ++d) {
 51:         subpoint[d] = -1.0;
 52:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
 53:       }
 54:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
 55:       if (inside) {cmp->embedding[s*spdim+sd++] = j;}
 56:     }
 57:     if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
 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) {
 83:           invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*weights[q];
 84:         }
 85:       }
 86:       PetscFree(Bf);
 87:     }
 88:     n = spdim;
 89:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info));
 90:     PetscStackCallBLAS("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:   PetscFree(invVscalar);
 95: #endif
 96:   PetscFree2(pivots,work);
 97:   return(0);
 98: }

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

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

126:       /* Apply transform, and check that point is inside cell */
127:       for (d = 0; d < dim; ++d) {
128:         subpoint[d] = -1.0;
129:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
130:       }
131:       CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
132:       if (inside) {subpoints[p] = s; break;}
133:     }
134:     if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
135:   }
136:   DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);
137:   /* Evaluate the prime basis functions at all points */
138:   if (B) {DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
139:   if (D) {DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
140:   if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
141:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
142:   /* Translate to the nodal basis */
143:   if (B) {PetscArrayzero(B, npoints*pdim*comp);}
144:   if (D) {PetscArrayzero(D, npoints*pdim*comp*dim);}
145:   if (H) {PetscArrayzero(H, npoints*pdim*comp*dim*dim);}
146:   for (p = 0; p < npoints; ++p) {
147:     const PetscInt s = subpoints[p];

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

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

166:           D[i] = 0.0;
167:           for (k = 0; k < spdim; ++k) {
168:             D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
169:           }
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) {
181:             H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
182:           }
183:         }
184:       }
185:     }
186:   }
187:   DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);
188:   if (B) {DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
189:   if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
190:   if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
191:   return(0);
192: }

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

210: /*MC
211:   PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element

213:   Level: intermediate

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

224:   PetscNewLog(fem, &cmp);
225:   fem->data = cmp;

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

232:   PetscFEInitialize_Composite(fem);
233:   return(0);
234: }

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

239:   Not collective

241:   Input Parameter:
242: . fem - The PetscFE object

244:   Output Parameters:
245: + blockSize - The number of elements in a block
246: . numBlocks - The number of blocks in a batch
247: . batchSize - The number of elements in a batch
248: - numBatches - The number of batches in a chunk

250:   Level: intermediate

252: .seealso: PetscFECreate()
253: @*/
254: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
255: {
256:   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;

264:   return(0);
265: }