Actual source code: fevector.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petsc/private/petscimpl.h>
4: typedef struct _n_PetscFE_Vec {
5: PetscFE scalar_fe;
6: PetscInt num_copies;
7: PetscBool interleave_basis;
8: PetscBool interleave_components;
9: } PetscFE_Vec;
11: static PetscErrorCode PetscFEDestroy_Vector(PetscFE fe)
12: {
13: PetscFE_Vec *v;
15: PetscFunctionBegin;
16: v = (PetscFE_Vec *)fe->data;
17: PetscCall(PetscFEDestroy(&v->scalar_fe));
18: PetscCall(PetscFree(v));
19: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: static PetscErrorCode PetscFEView_Vector_Ascii(PetscFE fe, PetscViewer v)
23: {
24: PetscInt dim, Nc, scalar_Nc;
25: PetscSpace basis = NULL;
26: PetscDualSpace dual = NULL;
27: PetscQuadrature quad = NULL;
28: PetscFE_Vec *vec;
29: PetscViewerFormat fmt;
31: PetscFunctionBegin;
32: vec = (PetscFE_Vec *)fe->data;
33: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
34: PetscCall(PetscFEGetNumComponents(fe, &Nc));
35: PetscCall(PetscFEGetNumComponents(vec->scalar_fe, &scalar_Nc));
36: PetscCall(PetscFEGetBasisSpace(fe, &basis));
37: PetscCall(PetscFEGetDualSpace(fe, &dual));
38: PetscCall(PetscFEGetQuadrature(fe, &quad));
39: PetscCall(PetscViewerGetFormat(v, &fmt));
40: PetscCall(PetscViewerASCIIPushTab(v));
41: if (scalar_Nc == 1) {
42: PetscCall(PetscViewerASCIIPrintf(v, "Vector Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components\n", dim, Nc));
43: } else {
44: PetscCall(PetscViewerASCIIPrintf(v, "Vector Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components (%" PetscInt_FMT " copies of finite element with %" PetscInt_FMT " components)\n", dim, Nc, vec->num_copies, scalar_Nc));
45: }
46: if (fmt == PETSC_VIEWER_ASCII_INFO_DETAIL) {
47: PetscCall(PetscViewerASCIIPrintf(v, "Original finite element:\n"));
48: PetscCall(PetscViewerASCIIPushTab(v));
49: PetscCall(PetscFEView(vec->scalar_fe, v));
50: PetscCall(PetscViewerASCIIPopTab(v));
51: }
52: if (basis) PetscCall(PetscSpaceView(basis, v));
53: if (dual) PetscCall(PetscDualSpaceView(dual, v));
54: if (quad) PetscCall(PetscQuadratureView(quad, v));
55: PetscCall(PetscViewerASCIIPopTab(v));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode PetscFEView_Vector(PetscFE fe, PetscViewer v)
60: {
61: PetscBool iascii;
63: PetscFunctionBegin;
64: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
65: if (iascii) PetscCall(PetscFEView_Vector_Ascii(fe, v));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode PetscFESetUp_Vector(PetscFE fe)
70: {
71: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
72: PetscDualSpace dsp;
73: PetscInt n, Ncopies = v->num_copies;
74: PetscInt scalar_n;
75: PetscInt *d, *d_mapped;
76: PetscDualSpace_Sum *sum;
77: PetscBool is_sum;
79: PetscFunctionBegin;
80: PetscCall(PetscFESetUp(v->scalar_fe));
81: PetscCall(PetscFEGetDimension(v->scalar_fe, &scalar_n));
82: PetscCall(PetscFEGetDualSpace(fe, &dsp));
83: PetscCall(PetscObjectTypeCompare((PetscObject)dsp, PETSCDUALSPACESUM, &is_sum));
84: PetscCheck(is_sum, PetscObjectComm((PetscObject)fe), PETSC_ERR_ARG_INCOMP, "Expected PETSCDUALSPACESUM dual space");
85: sum = (PetscDualSpace_Sum *)dsp->data;
86: n = Ncopies * scalar_n;
87: PetscCall(PetscCalloc1(n * n, &fe->invV));
88: PetscCall(PetscMalloc2(scalar_n, &d, scalar_n, &d_mapped));
89: for (PetscInt i = 0; i < scalar_n; i++) d[i] = i;
90: for (PetscInt c = 0; c < Ncopies; c++) {
91: PetscCall(ISLocalToGlobalMappingApply(sum->all_rows[c], scalar_n, d, d_mapped));
92: for (PetscInt i = 0; i < scalar_n; i++) {
93: PetscInt iw = d_mapped[i];
94: PetscReal *row_w = &fe->invV[iw * n];
95: const PetscReal *row_r = &v->scalar_fe->invV[i * scalar_n];
96: PetscInt j0 = v->interleave_basis ? c : c * scalar_n;
97: PetscInt jstride = v->interleave_basis ? Ncopies : 1;
99: for (PetscInt j = 0; j < scalar_n; j++) row_w[j0 + j * jstride] = row_r[j];
100: }
101: }
102: PetscCall(PetscFree2(d, d_mapped));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: static PetscErrorCode PetscFEGetDimension_Vector(PetscFE fe, PetscInt *dim)
107: {
108: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
110: PetscFunctionBegin;
111: PetscCall(PetscFEGetDimension(v->scalar_fe, dim));
112: *dim *= v->num_copies;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode PetscFEVectorInsertTabulation(PetscFE fe, PetscInt npoints, const PetscReal points[], PetscInt k, PetscInt scalar_Nb, PetscInt scalar_point_stride, const PetscReal scalar_Tk[], PetscReal Tk[])
117: {
118: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
119: PetscInt scalar_Nc;
120: PetscInt Nc;
121: PetscInt cdim;
122: PetscInt dblock;
123: PetscInt block;
124: PetscInt scalar_block;
126: PetscFunctionBegin;
127: PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
128: Nc = scalar_Nc * v->num_copies;
129: PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
130: dblock = PetscPowInt(cdim, k);
131: block = Nc * dblock;
132: scalar_block = scalar_Nc * dblock;
133: for (PetscInt p = 0; p < npoints; p++) {
134: const PetscReal *scalar_Tp = &scalar_Tk[p * scalar_point_stride];
135: PetscReal *Tp = &Tk[p * scalar_point_stride * v->num_copies * v->num_copies];
137: for (PetscInt j = 0; j < v->num_copies; j++) {
138: for (PetscInt scalar_i = 0; scalar_i < scalar_Nb; scalar_i++) {
139: PetscInt i = v->interleave_basis ? (scalar_i * v->num_copies + j) : (j * scalar_Nb + scalar_i);
140: const PetscReal *scalar_Ti = &scalar_Tp[scalar_i * scalar_block];
141: PetscReal *Ti = &Tp[i * block];
143: for (PetscInt scalar_c = 0; scalar_c < scalar_Nc; scalar_c++) {
144: PetscInt c = v->interleave_components ? (scalar_c * v->num_copies + j) : (j * scalar_Nc + scalar_c);
145: const PetscReal *scalar_Tc = &scalar_Ti[scalar_c * dblock];
146: PetscReal *Tc = &Ti[c * dblock];
148: PetscCall(PetscArraycpy(Tc, scalar_Tc, dblock));
149: }
150: }
151: }
152: }
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode PetscFECreateTabulation_Vector(PetscFE fe, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
157: {
158: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
159: PetscInt scalar_Nc;
160: PetscInt scalar_Nb;
161: PetscInt cdim;
162: PetscTabulation scalar_T;
164: PetscFunctionBegin;
165: PetscAssert(npoints == T->Nr * T->Np, PetscObjectComm((PetscObject)fe), PETSC_ERR_PLIB, "Expected to be able decode PetscFECreateTabulation() from PetscTabulation fields");
166: PetscCall(PetscFECreateTabulation(v->scalar_fe, T->Nr, T->Np, points, K, &scalar_T));
167: PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
168: PetscCall(PetscFEGetDimension(v->scalar_fe, &scalar_Nb));
169: PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
170: for (PetscInt k = 0; k <= T->K; k++) {
171: PetscReal *Tk = T->T[k];
172: const PetscReal *scalar_Tk = scalar_T->T[k];
173: PetscInt dblock = PetscPowInt(cdim, k);
174: PetscInt scalar_block = scalar_Nc * dblock;
175: PetscInt scalar_point_stride = scalar_Nb * scalar_block;
177: if (v->interleave_basis) {
178: PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, scalar_Nb, scalar_point_stride, scalar_Tk, Tk));
179: } else {
180: PetscDualSpace scalar_dsp;
181: PetscSection ref_section;
182: PetscInt pStart, pEnd;
184: PetscCall(PetscFEGetDualSpace(v->scalar_fe, &scalar_dsp));
185: PetscCall(PetscDualSpaceGetSection(scalar_dsp, &ref_section));
186: PetscCall(PetscSectionGetChart(ref_section, &pStart, &pEnd));
187: for (PetscInt p = pStart; p < pEnd; p++) {
188: PetscInt dof, off;
189: PetscReal *Tp;
190: const PetscReal *scalar_Tp;
192: PetscCall(PetscSectionGetDof(ref_section, p, &dof));
193: PetscCall(PetscSectionGetOffset(ref_section, p, &off));
194: scalar_Tp = &scalar_Tk[off * scalar_block];
195: Tp = &Tk[off * scalar_block * v->num_copies * v->num_copies];
196: PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, dof, scalar_point_stride, scalar_Tp, Tp));
197: }
198: }
199: }
200: PetscCall(PetscTabulationDestroy(&scalar_T));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode PetscFECreatePointTrace_Vector(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
205: {
206: PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
207: PetscFE scalar_trFE;
208: const char *name;
210: PetscFunctionBegin;
211: PetscCall(PetscFECreatePointTrace(v->scalar_fe, refPoint, &scalar_trFE));
212: PetscCall(PetscFECreateVector(scalar_trFE, v->num_copies, v->interleave_basis, v->interleave_components, trFE));
213: PetscCall(PetscFEDestroy(&scalar_trFE));
214: PetscCall(PetscObjectGetName((PetscObject)fe, &name));
215: if (name) PetscCall(PetscFESetName(*trFE, name));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: PETSC_INTERN PetscErrorCode PetscFEIntegrate_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
220: PETSC_INTERN PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
221: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS, PetscDS, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
222: PETSC_INTERN PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS, PetscWeakForm, PetscFEJacobianType, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
223: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS, PetscDS, PetscFEJacobianType, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
225: static PetscErrorCode PetscFEInitialize_Vector(PetscFE fe)
226: {
227: PetscFunctionBegin;
228: fe->ops->setfromoptions = NULL;
229: fe->ops->setup = PetscFESetUp_Vector;
230: fe->ops->view = PetscFEView_Vector;
231: fe->ops->destroy = PetscFEDestroy_Vector;
232: fe->ops->getdimension = PetscFEGetDimension_Vector;
233: fe->ops->createpointtrace = PetscFECreatePointTrace_Vector;
234: fe->ops->createtabulation = PetscFECreateTabulation_Vector;
235: fe->ops->integrate = PetscFEIntegrate_Basic;
236: fe->ops->integratebd = PetscFEIntegrateBd_Basic;
237: fe->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
238: fe->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
239: fe->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
240: fe->ops->integratejacobianaction = NULL;
241: fe->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
242: fe->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
243: fe->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*MC
248: PETSCFEVECTOR = "vector" - A vector-valued `PetscFE` object that is repeated copies
249: of the same underlying finite element.
251: Level: intermediate
253: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PetscFECreateVector()`
254: M*/
255: PETSC_EXTERN PetscErrorCode PetscFECreate_Vector(PetscFE fe)
256: {
257: PetscFE_Vec *v;
259: PetscFunctionBegin;
261: PetscCall(PetscNew(&v));
262: fe->data = v;
264: PetscCall(PetscFEInitialize_Vector(fe));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: PetscFECreateVector - Create a vector-valued `PetscFE` from multiple copies of an underlying
270: `PetscFE`.
272: Collective
274: Input Parameters:
275: + scalar_fe - a `PetscFE` finite element
276: . num_copies - a positive integer
277: . interleave_basis - if `PETSC_TRUE`, the first `num_copies` basis vectors
278: of the output finite element will be copies of the
279: first basis vector of `scalar_fe`, and so on for the
280: other basis vectors; otherwise all of the first-copy
281: basis vectors will come first, followed by all of the
282: second-copy, and so on.
283: - interleave_components - if `PETSC_TRUE`, the first `num_copies` components
284: of the output finite element will be copies of the
285: first component of `scalar_fe`, and so on for the
286: other components; otherwise all of the first-copy
287: components will come first, followed by all of the
288: second-copy, and so on.
290: Output Parameter:
291: . vector_fe - a `PetscFE` of type `PETSCFEVECTOR` that represent a discretization space with `num_copies` copies of `scalar_fe`
293: Level: intermediate
295: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PETSCFEVECTOR`
296: @*/
297: PetscErrorCode PetscFECreateVector(PetscFE scalar_fe, PetscInt num_copies, PetscBool interleave_basis, PetscBool interleave_components, PetscFE *vector_fe)
298: {
299: MPI_Comm comm;
300: PetscFE fe_vec;
301: PetscFE_Vec *v;
302: PetscInt scalar_Nc;
303: PetscQuadrature quad;
305: PetscFunctionBegin;
307: PetscCall(PetscObjectGetComm((PetscObject)scalar_fe, &comm));
308: PetscCheck(num_copies >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Expected positive number of copies, got %" PetscInt_FMT, num_copies);
309: PetscCall(PetscFECreate(comm, vector_fe));
310: fe_vec = *vector_fe;
311: PetscCall(PetscFESetType(fe_vec, PETSCFEVECTOR));
312: v = (PetscFE_Vec *)fe_vec->data;
313: PetscCall(PetscObjectReference((PetscObject)scalar_fe));
314: v->scalar_fe = scalar_fe;
315: v->num_copies = num_copies;
316: v->interleave_basis = interleave_basis;
317: v->interleave_components = interleave_components;
318: PetscCall(PetscFEGetNumComponents(scalar_fe, &scalar_Nc));
319: PetscCall(PetscFESetNumComponents(fe_vec, scalar_Nc * num_copies));
320: PetscCall(PetscFEGetQuadrature(scalar_fe, &quad));
321: PetscCall(PetscFESetQuadrature(fe_vec, quad));
322: PetscCall(PetscFEGetFaceQuadrature(scalar_fe, &quad));
323: PetscCall(PetscFESetFaceQuadrature(fe_vec, quad));
324: {
325: PetscSpace scalar_sp;
326: PetscSpace *copies;
327: PetscSpace sp;
329: PetscCall(PetscFEGetBasisSpace(scalar_fe, &scalar_sp));
330: PetscCall(PetscMalloc1(num_copies, &copies));
331: for (PetscInt i = 0; i < num_copies; i++) {
332: PetscCall(PetscObjectReference((PetscObject)scalar_sp));
333: copies[i] = scalar_sp;
334: }
335: PetscCall(PetscSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
336: PetscCall(PetscSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
337: PetscCall(PetscSpaceSetUp(sp));
338: PetscCall(PetscFESetBasisSpace(fe_vec, sp));
339: PetscCall(PetscSpaceDestroy(&sp));
340: for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscSpaceDestroy(&copies[i]));
341: PetscCall(PetscFree(copies));
342: }
343: {
344: PetscDualSpace scalar_sp;
345: PetscDualSpace *copies;
346: PetscDualSpace sp;
348: PetscCall(PetscFEGetDualSpace(scalar_fe, &scalar_sp));
349: PetscCall(PetscMalloc1(num_copies, &copies));
350: for (PetscInt i = 0; i < num_copies; i++) {
351: PetscCall(PetscObjectReference((PetscObject)scalar_sp));
352: copies[i] = scalar_sp;
353: }
354: PetscCall(PetscDualSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
355: PetscCall(PetscDualSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
356: PetscCall(PetscDualSpaceSetUp(sp));
357: PetscCall(PetscFESetDualSpace(fe_vec, sp));
358: PetscCall(PetscDualSpaceDestroy(&sp));
359: for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscDualSpaceDestroy(&copies[i]));
360: PetscCall(PetscFree(copies));
361: }
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }