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, §ion));
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: PetscCall(PetscBLASIntCast(spdim, &n));
89: // The next two lines are likely long-standing incorrect code. Do not check info as it was not previously checked but is nonzero
90: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s * spdim * spdim], &n, pivots, &info));
91: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s * spdim * spdim], &n, pivots, work, &n, &info));
92: }
93: #if defined(PETSC_USE_COMPLEX)
94: for (s = 0; s < cmp->numSubelements * spdim * spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
95: PetscCall(PetscFree(invVscalar));
96: #endif
97: PetscCall(PetscFree2(pivots, work));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode PetscFEComputeTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
102: {
103: PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
104: DM dm;
105: DMPolytopeType ct;
106: PetscInt pdim; /* Dimension of FE space P */
107: PetscInt spdim; /* Dimension of subelement FE space P */
108: PetscInt dim; /* Spatial dimension */
109: PetscInt comp; /* Field components */
110: PetscInt *subpoints;
111: PetscReal *B = K >= 0 ? T->T[0] : NULL;
112: PetscReal *D = K >= 1 ? T->T[1] : NULL;
113: PetscReal *H = K >= 2 ? T->T[2] : NULL;
114: PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint;
115: PetscInt p, s, d, e, j, k;
117: PetscFunctionBegin;
118: PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm));
119: PetscCall(DMGetDimension(dm, &dim));
120: PetscCall(DMPlexGetCellType(dm, 0, &ct));
121: PetscCall(PetscSpaceGetDimension(fem->basisSpace, &spdim));
122: PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim));
123: PetscCall(PetscFEGetNumComponents(fem, &comp));
124: /* Divide points into subelements */
125: PetscCall(DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints));
126: PetscCall(DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint));
127: for (p = 0; p < npoints; ++p) {
128: for (s = 0; s < cmp->numSubelements; ++s) {
129: PetscBool inside;
131: /* Apply transform, and check that point is inside cell */
132: for (d = 0; d < dim; ++d) {
133: subpoint[d] = -1.0;
134: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s * dim + d) * dim + e] * (points[p * dim + e] - cmp->v0[s * dim + e]);
135: }
136: PetscCall(DMPolytopeInCellTest(ct, subpoint, &inside));
137: if (inside) {
138: subpoints[p] = s;
139: break;
140: }
141: }
142: PetscCheck(s < cmp->numSubelements, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " was not found in any subelement", p);
143: }
144: PetscCall(DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint));
145: /* Evaluate the prime basis functions at all points */
146: if (K >= 0) PetscCall(DMGetWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB));
147: if (K >= 1) PetscCall(DMGetWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD));
148: if (K >= 2) PetscCall(DMGetWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH));
149: PetscCall(PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH));
150: /* Translate to the nodal basis */
151: if (K >= 0) PetscCall(PetscArrayzero(B, npoints * pdim * comp));
152: if (K >= 1) PetscCall(PetscArrayzero(D, npoints * pdim * comp * dim));
153: if (K >= 2) PetscCall(PetscArrayzero(H, npoints * pdim * comp * dim * dim));
154: for (p = 0; p < npoints; ++p) {
155: const PetscInt s = subpoints[p];
157: if (B) {
158: /* Multiply by V^{-1} (spdim x spdim) */
159: for (j = 0; j < spdim; ++j) {
160: const PetscInt i = (p * pdim + cmp->embedding[s * spdim + j]) * comp;
162: B[i] = 0.0;
163: for (k = 0; k < spdim; ++k) B[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpB[p * spdim + k];
164: }
165: }
166: if (D) {
167: /* Multiply by V^{-1} (spdim x spdim) */
168: for (j = 0; j < spdim; ++j) {
169: for (d = 0; d < dim; ++d) {
170: const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim + d;
172: D[i] = 0.0;
173: for (k = 0; k < spdim; ++k) D[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpD[(p * spdim + k) * dim + d];
174: }
175: }
176: }
177: if (H) {
178: /* Multiply by V^{-1} (pdim x pdim) */
179: for (j = 0; j < spdim; ++j) {
180: for (d = 0; d < dim * dim; ++d) {
181: const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim * dim + d;
183: H[i] = 0.0;
184: for (k = 0; k < spdim; ++k) H[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpH[(p * spdim + k) * dim * dim + d];
185: }
186: }
187: }
188: }
189: PetscCall(DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints));
190: if (K >= 0) PetscCall(DMRestoreWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB));
191: if (K >= 1) PetscCall(DMRestoreWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD));
192: if (K >= 2) PetscCall(DMRestoreWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
197: {
198: PetscFunctionBegin;
199: fem->ops->setfromoptions = NULL;
200: fem->ops->setup = PetscFESetUp_Composite;
201: fem->ops->view = NULL;
202: fem->ops->destroy = PetscFEDestroy_Composite;
203: fem->ops->getdimension = PetscFEGetDimension_Basic;
204: fem->ops->computetabulation = PetscFEComputeTabulation_Composite;
205: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
206: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
207: fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
208: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*MC
213: PETSCFECOMPOSITE = "composite" - A `PetscFEType` that represents a composite element
215: Level: intermediate
217: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
218: M*/
219: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
220: {
221: PetscFE_Composite *cmp;
223: PetscFunctionBegin;
225: PetscCall(PetscNew(&cmp));
226: fem->data = cmp;
228: cmp->numSubelements = -1;
229: cmp->v0 = NULL;
230: cmp->jac = NULL;
232: PetscCall(PetscFEInitialize_Composite(fem));
233: PetscFunctionReturn(PETSC_SUCCESS);
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: + numSubelements - The number of sub elements
246: . v0 - The affine transformation for each element, an array of length $dim * Nc$. Pass `NULL` to ignore.
247: . jac - The Jacobian for each element, an array of length $dim^2 * Nc$. Pass `NULL` to ignore.
248: - invjac - The inverse of the Jacobian, an array of length $dim^2 * Nc$. Pass `NULL` to ignore.
250: Level: intermediate
252: Note:
253: Do not free the output arrays.
255: .seealso: `PetscFE`, `PetscFECreate()`
256: @*/
257: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PeOp PetscInt *numSubelements, PeOp const PetscReal *v0[], PeOp const PetscReal *jac[], PeOp const PetscReal *invjac[])
258: {
259: PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
261: PetscFunctionBegin;
263: if (numSubelements) {
264: PetscAssertPointer(numSubelements, 2);
265: *numSubelements = cmp->numSubelements;
266: }
267: if (v0) {
268: PetscAssertPointer(v0, 3);
269: *v0 = cmp->v0;
270: }
271: if (jac) {
272: PetscAssertPointer(jac, 4);
273: *jac = cmp->jac;
274: }
275: if (invjac) {
276: PetscAssertPointer(invjac, 5);
277: *invjac = cmp->invjac;
278: }
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }