Actual source code: fecomposite.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/petscfeimpl.h>
2: #include <petsc/private/dtimpl.h>
3: #include <petscblaslapack.h>
5: static PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
6: {
7: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
8: PetscErrorCode ierr;
11: PetscFree(cmp->embedding);
12: PetscFree(cmp);
13: return(0);
14: }
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;
31: /* Get affine mapping from reference cell to each subcell */
32: PetscDualSpaceGetDM(fem->dualSpace, &K);
33: DMGetDimension(K, &dim);
34: DMPlexGetCellType(K, 0, &ct);
35: DMPlexCellRefinerCreate(K, &cr);
36: DMPlexCellRefinerGetAffineTransforms(cr, ct, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
37: DMPlexCellRefinerDestroy(&cr);
38: /* Determine dof embedding into subelements */
39: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
40: PetscSpaceGetDimension(fem->basisSpace, &spdim);
41: PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
42: DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);
43: PetscDualSpaceGetSection(fem->dualSpace, §ion);
44: for (s = 0; s < cmp->numSubelements; ++s) {
45: PetscInt sd = 0;
46: PetscInt closureSize;
47: PetscInt *closure = NULL;
49: 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: PetscSectionGetDof(section, point, &dof);
55: PetscSectionGetOffset(section, point, &off);
56: for (k = 0; k < dof; k++) cmp->embedding[s*spdim+sd++] = off + k;
57: }
58: DMPlexRestoreTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);
59: if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
60: }
61: DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint);
62: /* Construct the change of basis from prime basis to nodal basis for each subelement */
63: PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
64: PetscMalloc2(spdim,&pivots,spdim,&work);
65: #if defined(PETSC_USE_COMPLEX)
66: 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: PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
78: PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
79: PetscMalloc1(f->numPoints*spdim*Nc,&Bf);
80: 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) {
85: invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*weights[q];
86: }
87: }
88: PetscFree(Bf);
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: PetscFree(invVscalar);
97: #endif
98: PetscFree2(pivots,work);
99: return(0);
100: }
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;
120: PetscDualSpaceGetDM(fem->dualSpace, &dm);
121: DMGetDimension(dm, &dim);
122: DMPlexGetCellType(dm, 0, &ct);
123: PetscSpaceGetDimension(fem->basisSpace, &spdim);
124: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
125: PetscFEGetNumComponents(fem, &comp);
126: /* Divide points into subelements */
127: DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);
128: DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);
129: for (p = 0; p < npoints; ++p) {
130: for (s = 0; s < cmp->numSubelements; ++s) {
131: PetscBool inside;
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: CellRefinerInCellTest_Internal(ct, subpoint, &inside);
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: DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);
144: /* Evaluate the prime basis functions at all points */
145: if (K >= 0) {DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
146: if (K >= 1) {DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
147: if (K >= 2) {DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
148: PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
149: /* Translate to the nodal basis */
150: if (K >= 0) {PetscArrayzero(B, npoints*pdim*comp);}
151: if (K >= 1) {PetscArrayzero(D, npoints*pdim*comp*dim);}
152: if (K >= 2) {PetscArrayzero(H, npoints*pdim*comp*dim*dim);}
153: for (p = 0; p < npoints; ++p) {
154: const PetscInt s = subpoints[p];
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;
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;
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;
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: DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);
195: if (K >= 0) {DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
196: if (K >= 1) {DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
197: if (K >= 2) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
198: return(0);
199: }
201: static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
202: {
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: return(0);
215: }
217: /*MC
218: PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
220: Level: intermediate
222: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
223: M*/
224: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
225: {
226: PetscFE_Composite *cmp;
227: PetscErrorCode ierr;
231: PetscNewLog(fem, &cmp);
232: fem->data = cmp;
234: cmp->numSubelements = -1;
235: cmp->v0 = NULL;
236: cmp->jac = NULL;
238: PetscFEInitialize_Composite(fem);
239: return(0);
240: }
242: /*@C
243: PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
245: Not collective
247: Input Parameter:
248: . fem - The PetscFE object
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
256: Level: intermediate
258: .seealso: PetscFECreate()
259: @*/
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;
270: return(0);
271: }