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 
PetscFEDestroy_Composite(PetscFE fem)5 static PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
6 {
7   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
8   PetscErrorCode     ierr;
9 
10   PetscFunctionBegin;
11   ierr = PetscFree(cmp->embedding);CHKERRQ(ierr);
12   ierr = PetscFree(cmp);CHKERRQ(ierr);
13   PetscFunctionReturn(0);
14 }
15 
PetscFESetUp_Composite(PetscFE fem)16 static PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
17 {
18   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
19   DM                 K;
20   DMPolytopeType     ct;
21   DMPlexCellRefiner  cr;
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   PetscErrorCode     ierr;
29 
30   PetscFunctionBegin;
31   /* Get affine mapping from reference cell to each subcell */
32   ierr = PetscDualSpaceGetDM(fem->dualSpace, &K);CHKERRQ(ierr);
33   ierr = DMGetDimension(K, &dim);CHKERRQ(ierr);
34   ierr = DMPlexGetCellType(K, 0, &ct);CHKERRQ(ierr);
35   ierr = DMPlexCellRefinerCreate(K, &cr);CHKERRQ(ierr);
36   ierr = DMPlexCellRefinerGetAffineTransforms(cr, ct, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);CHKERRQ(ierr);
37   ierr = DMPlexCellRefinerDestroy(&cr);CHKERRQ(ierr);
38   /* Determine dof embedding into subelements */
39   ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr);
40   ierr = PetscSpaceGetDimension(fem->basisSpace, &spdim);CHKERRQ(ierr);
41   ierr = PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);CHKERRQ(ierr);
42   ierr = DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr);
43   ierr = PetscDualSpaceGetSection(fem->dualSpace, &section);CHKERRQ(ierr);
44   for (s = 0; s < cmp->numSubelements; ++s) {
45     PetscInt sd = 0;
46     PetscInt closureSize;
47     PetscInt *closure = NULL;
48 
49     ierr = DMPlexGetTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
50     for (j = 0; j < closureSize; j++) {
51       PetscInt point = closure[2*j];
52       PetscInt dof, off, k;
53 
54       ierr = PetscSectionGetDof(section, point, &dof);CHKERRQ(ierr);
55       ierr = PetscSectionGetOffset(section, point, &off);CHKERRQ(ierr);
56       for (k = 0; k < dof; k++) cmp->embedding[s*spdim+sd++] = off + k;
57     }
58     ierr = DMPlexRestoreTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
59     if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
60   }
61   ierr = DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr);
62   /* Construct the change of basis from prime basis to nodal basis for each subelement */
63   ierr = PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);CHKERRQ(ierr);
64   ierr = PetscMalloc2(spdim,&pivots,spdim,&work);CHKERRQ(ierr);
65 #if defined(PETSC_USE_COMPLEX)
66   ierr = PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);CHKERRQ(ierr);
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       ierr = PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);CHKERRQ(ierr);
78       ierr = PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);CHKERRQ(ierr);
79       ierr = PetscMalloc1(f->numPoints*spdim*Nc,&Bf);CHKERRQ(ierr);
80       ierr = PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);CHKERRQ(ierr);
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       ierr = PetscFree(Bf);CHKERRQ(ierr);
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   ierr = PetscFree(invVscalar);CHKERRQ(ierr);
97 #endif
98   ierr = PetscFree2(pivots,work);CHKERRQ(ierr);
99   PetscFunctionReturn(0);
100 }
101 
PetscFECreateTabulation_Composite(PetscFE fem,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation T)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   PetscErrorCode     ierr;
118 
119   PetscFunctionBegin;
120   ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr);
121   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
122   ierr = DMPlexGetCellType(dm, 0, &ct);CHKERRQ(ierr);
123   ierr = PetscSpaceGetDimension(fem->basisSpace, &spdim);CHKERRQ(ierr);
124   ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr);
125   ierr = PetscFEGetNumComponents(fem, &comp);CHKERRQ(ierr);
126   /* Divide points into subelements */
127   ierr = DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);CHKERRQ(ierr);
128   ierr = DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr);
129   for (p = 0; p < npoints; ++p) {
130     for (s = 0; s < cmp->numSubelements; ++s) {
131       PetscBool inside;
132 
133       /* Apply transform, and check that point is inside cell */
134       for (d = 0; d < dim; ++d) {
135         subpoint[d] = -1.0;
136         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
137       }
138       ierr = CellRefinerInCellTest_Internal(ct, subpoint, &inside);CHKERRQ(ierr);
139       if (inside) {subpoints[p] = s; break;}
140     }
141     if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
142   }
143   ierr = DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr);
144   /* Evaluate the prime basis functions at all points */
145   if (K >= 0) {ierr = DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
146   if (K >= 1) {ierr = DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
147   if (K >= 2) {ierr = DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
148   ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);CHKERRQ(ierr);
149   /* Translate to the nodal basis */
150   if (K >= 0) {ierr = PetscArrayzero(B, npoints*pdim*comp);CHKERRQ(ierr);}
151   if (K >= 1) {ierr = PetscArrayzero(D, npoints*pdim*comp*dim);CHKERRQ(ierr);}
152   if (K >= 2) {ierr = PetscArrayzero(H, npoints*pdim*comp*dim*dim);CHKERRQ(ierr);}
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) {
163           B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
164         }
165       }
166     }
167     if (D) {
168       /* Multiply by V^{-1} (spdim x spdim) */
169       for (j = 0; j < spdim; ++j) {
170         for (d = 0; d < dim; ++d) {
171           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
172 
173           D[i] = 0.0;
174           for (k = 0; k < spdim; ++k) {
175             D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
176           }
177         }
178       }
179     }
180     if (H) {
181       /* Multiply by V^{-1} (pdim x pdim) */
182       for (j = 0; j < spdim; ++j) {
183         for (d = 0; d < dim*dim; ++d) {
184           const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
185 
186           H[i] = 0.0;
187           for (k = 0; k < spdim; ++k) {
188             H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
189           }
190         }
191       }
192     }
193   }
194   ierr = DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);CHKERRQ(ierr);
195   if (K >= 0) {ierr = DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
196   if (K >= 1) {ierr = DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
197   if (K >= 2) {ierr = DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
198   PetscFunctionReturn(0);
199 }
200 
PetscFEInitialize_Composite(PetscFE fem)201 static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
202 {
203   PetscFunctionBegin;
204   fem->ops->setfromoptions          = NULL;
205   fem->ops->setup                   = PetscFESetUp_Composite;
206   fem->ops->view                    = NULL;
207   fem->ops->destroy                 = PetscFEDestroy_Composite;
208   fem->ops->getdimension            = PetscFEGetDimension_Basic;
209   fem->ops->createtabulation        = PetscFECreateTabulation_Composite;
210   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
211   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
212   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
213   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
214   PetscFunctionReturn(0);
215 }
216 
217 /*MC
218   PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
219 
220   Level: intermediate
221 
222 .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
223 M*/
PetscFECreate_Composite(PetscFE fem)224 PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
225 {
226   PetscFE_Composite *cmp;
227   PetscErrorCode     ierr;
228 
229   PetscFunctionBegin;
230   PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
231   ierr      = PetscNewLog(fem, &cmp);CHKERRQ(ierr);
232   fem->data = cmp;
233 
234   cmp->numSubelements = -1;
235   cmp->v0             = NULL;
236   cmp->jac            = NULL;
237 
238   ierr = PetscFEInitialize_Composite(fem);CHKERRQ(ierr);
239   PetscFunctionReturn(0);
240 }
241 
242 /*@C
243   PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
244 
245   Not collective
246 
247   Input Parameter:
248 . fem - The PetscFE object
249 
250   Output Parameters:
251 + blockSize - The number of elements in a block
252 . numBlocks - The number of blocks in a batch
253 . batchSize - The number of elements in a batch
254 - numBatches - The number of batches in a chunk
255 
256   Level: intermediate
257 
258 .seealso: PetscFECreate()
259 @*/
PetscFECompositeGetMapping(PetscFE fem,PetscInt * numSubelements,const PetscReal * v0[],const PetscReal * jac[],const PetscReal * invjac[])260 PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
261 {
262   PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
263 
264   PetscFunctionBegin;
265   PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
266   if (numSubelements) {PetscValidPointer(numSubelements, 2); *numSubelements = cmp->numSubelements;}
267   if (v0)             {PetscValidPointer(v0, 3);             *v0             = cmp->v0;}
268   if (jac)            {PetscValidPointer(jac, 4);            *jac            = cmp->jac;}
269   if (invjac)         {PetscValidPointer(invjac, 5);         *invjac         = cmp->invjac;}
270   PetscFunctionReturn(0);
271 }
272