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
PetscFEDestroy_Composite(PetscFE fem)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
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 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 PetscCall(PetscBLASIntCast(spdim, &n));
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
PetscFEComputeTabulation_Composite(PetscFE fem,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation T)100 static PetscErrorCode PetscFEComputeTabulation_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
PetscFEInitialize_Composite(PetscFE fem)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->computetabulation = PetscFEComputeTabulation_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*/
PetscFECreate_Composite(PetscFE fem)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, an array of length $dim * Nc$. Pass `NULL` to ignore.
246 . jac - The Jacobian for each element, an array of length $dim^2 * Nc$. Pass `NULL` to ignore.
247 - invjac - The inverse of the Jacobian, an array of length $dim^2 * Nc$. Pass `NULL` to ignore.
248
249 Level: intermediate
250
251 Note:
252 Do not free the output arrays.
253
254 .seealso: `PetscFE`, `PetscFECreate()`
255 @*/
PetscFECompositeGetMapping(PetscFE fem,PeOp PetscInt * numSubelements,PeOp const PetscReal * v0[],PeOp const PetscReal * jac[],PeOp const PetscReal * invjac[])256 PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PeOp PetscInt *numSubelements, PeOp const PetscReal *v0[], PeOp const PetscReal *jac[], PeOp const PetscReal *invjac[])
257 {
258 PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
259
260 PetscFunctionBegin;
261 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
262 if (numSubelements) {
263 PetscAssertPointer(numSubelements, 2);
264 *numSubelements = cmp->numSubelements;
265 }
266 if (v0) {
267 PetscAssertPointer(v0, 3);
268 *v0 = cmp->v0;
269 }
270 if (jac) {
271 PetscAssertPointer(jac, 4);
272 *jac = cmp->jac;
273 }
274 if (invjac) {
275 PetscAssertPointer(invjac, 5);
276 *invjac = cmp->invjac;
277 }
278 PetscFunctionReturn(PETSC_SUCCESS);
279 }
280