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;
9: PetscErrorCode ierr;
12: PetscFree(cmp->embedding);
13: PetscFree(cmp);
14: return(0);
15: }
17: static PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
18: {
19: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
20: DM K;
21: DMPolytopeType ct;
22: DMPlexTransform tr;
23: PetscReal *subpoint;
24: PetscBLASInt *pivots;
25: PetscBLASInt n, info;
26: PetscScalar *work, *invVscalar;
27: PetscInt dim, pdim, spdim, j, s;
28: PetscSection section;
29: PetscErrorCode ierr;
32: /* Get affine mapping from reference cell to each subcell */
33: PetscDualSpaceGetDM(fem->dualSpace, &K);
34: DMGetDimension(K, &dim);
35: DMPlexGetCellType(K, 0, &ct);
36: DMPlexTransformCreate(PETSC_COMM_SELF, &tr);
37: DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR);
38: DMPlexRefineRegularGetAffineTransforms(tr, ct, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
39: DMPlexTransformDestroy(&tr);
40: /* Determine dof embedding into subelements */
41: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
42: PetscSpaceGetDimension(fem->basisSpace, &spdim);
43: PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
44: DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);
45: PetscDualSpaceGetSection(fem->dualSpace, §ion);
46: for (s = 0; s < cmp->numSubelements; ++s) {
47: PetscInt sd = 0;
48: PetscInt closureSize;
49: PetscInt *closure = NULL;
51: DMPlexGetTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);
52: for (j = 0; j < closureSize; j++) {
53: PetscInt point = closure[2*j];
54: PetscInt dof, off, k;
56: PetscSectionGetDof(section, point, &dof);
57: PetscSectionGetOffset(section, point, &off);
58: for (k = 0; k < dof; k++) cmp->embedding[s*spdim+sd++] = off + k;
59: }
60: DMPlexRestoreTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure);
61: if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
62: }
63: DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint);
64: /* Construct the change of basis from prime basis to nodal basis for each subelement */
65: PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
66: PetscMalloc2(spdim,&pivots,spdim,&work);
67: #if defined(PETSC_USE_COMPLEX)
68: PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);
69: #else
70: invVscalar = fem->invV;
71: #endif
72: for (s = 0; s < cmp->numSubelements; ++s) {
73: for (j = 0; j < spdim; ++j) {
74: PetscReal *Bf;
75: PetscQuadrature f;
76: const PetscReal *points, *weights;
77: PetscInt Nc, Nq, q, k;
79: PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
80: PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
81: PetscMalloc1(f->numPoints*spdim*Nc,&Bf);
82: PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);
83: for (k = 0; k < spdim; ++k) {
84: /* n_j \cdot \phi_k */
85: invVscalar[(s*spdim + j)*spdim+k] = 0.0;
86: for (q = 0; q < Nq; ++q) {
87: invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*weights[q];
88: }
89: }
90: PetscFree(Bf);
91: }
92: n = spdim;
93: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info));
94: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info));
95: }
96: #if defined(PETSC_USE_COMPLEX)
97: for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
98: PetscFree(invVscalar);
99: #endif
100: PetscFree2(pivots,work);
101: return(0);
102: }
104: static PetscErrorCode PetscFECreateTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
105: {
106: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
107: DM dm;
108: DMPolytopeType ct;
109: PetscInt pdim; /* Dimension of FE space P */
110: PetscInt spdim; /* Dimension of subelement FE space P */
111: PetscInt dim; /* Spatial dimension */
112: PetscInt comp; /* Field components */
113: PetscInt *subpoints;
114: PetscReal *B = K >= 0 ? T->T[0] : NULL;
115: PetscReal *D = K >= 1 ? T->T[1] : NULL;
116: PetscReal *H = K >= 2 ? T->T[2] : NULL;
117: PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint;
118: PetscInt p, s, d, e, j, k;
119: PetscErrorCode ierr;
122: PetscDualSpaceGetDM(fem->dualSpace, &dm);
123: DMGetDimension(dm, &dim);
124: DMPlexGetCellType(dm, 0, &ct);
125: PetscSpaceGetDimension(fem->basisSpace, &spdim);
126: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
127: PetscFEGetNumComponents(fem, &comp);
128: /* Divide points into subelements */
129: DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);
130: DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);
131: for (p = 0; p < npoints; ++p) {
132: for (s = 0; s < cmp->numSubelements; ++s) {
133: PetscBool inside;
135: /* Apply transform, and check that point is inside cell */
136: for (d = 0; d < dim; ++d) {
137: subpoint[d] = -1.0;
138: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
139: }
140: DMPolytopeInCellTest(ct, subpoint, &inside);
141: if (inside) {subpoints[p] = s; break;}
142: }
143: if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
144: }
145: DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);
146: /* Evaluate the prime basis functions at all points */
147: if (K >= 0) {DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
148: if (K >= 1) {DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
149: if (K >= 2) {DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
150: PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
151: /* Translate to the nodal basis */
152: if (K >= 0) {PetscArrayzero(B, npoints*pdim*comp);}
153: if (K >= 1) {PetscArrayzero(D, npoints*pdim*comp*dim);}
154: if (K >= 2) {PetscArrayzero(H, npoints*pdim*comp*dim*dim);}
155: for (p = 0; p < npoints; ++p) {
156: const PetscInt s = subpoints[p];
158: if (B) {
159: /* Multiply by V^{-1} (spdim x spdim) */
160: for (j = 0; j < spdim; ++j) {
161: const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
163: B[i] = 0.0;
164: for (k = 0; k < spdim; ++k) {
165: B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
166: }
167: }
168: }
169: if (D) {
170: /* Multiply by V^{-1} (spdim x spdim) */
171: for (j = 0; j < spdim; ++j) {
172: for (d = 0; d < dim; ++d) {
173: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
175: D[i] = 0.0;
176: for (k = 0; k < spdim; ++k) {
177: D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
178: }
179: }
180: }
181: }
182: if (H) {
183: /* Multiply by V^{-1} (pdim x pdim) */
184: for (j = 0; j < spdim; ++j) {
185: for (d = 0; d < dim*dim; ++d) {
186: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
188: H[i] = 0.0;
189: for (k = 0; k < spdim; ++k) {
190: H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
191: }
192: }
193: }
194: }
195: }
196: DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);
197: if (K >= 0) {DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);}
198: if (K >= 1) {DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);}
199: if (K >= 2) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);}
200: return(0);
201: }
203: static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
204: {
206: fem->ops->setfromoptions = NULL;
207: fem->ops->setup = PetscFESetUp_Composite;
208: fem->ops->view = NULL;
209: fem->ops->destroy = PetscFEDestroy_Composite;
210: fem->ops->getdimension = PetscFEGetDimension_Basic;
211: fem->ops->createtabulation = PetscFECreateTabulation_Composite;
212: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
213: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
214: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
215: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
216: return(0);
217: }
219: /*MC
220: PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
222: Level: intermediate
224: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
225: M*/
226: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
227: {
228: PetscFE_Composite *cmp;
229: PetscErrorCode ierr;
233: PetscNewLog(fem, &cmp);
234: fem->data = cmp;
236: cmp->numSubelements = -1;
237: cmp->v0 = NULL;
238: cmp->jac = NULL;
240: PetscFEInitialize_Composite(fem);
241: return(0);
242: }
244: /*@C
245: PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
247: Not collective
249: Input Parameter:
250: . fem - The PetscFE object
252: Output Parameters:
253: + blockSize - The number of elements in a block
254: . numBlocks - The number of blocks in a batch
255: . batchSize - The number of elements in a batch
256: - numBatches - The number of batches in a chunk
258: Level: intermediate
260: .seealso: PetscFECreate()
261: @*/
262: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
263: {
264: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
272: return(0);
273: }