Actual source code: dualspacesum.c
1: #include <petsc/private/petscfeimpl.h>
3: /*@
4: PetscDualSpaceSumGetNumSubspaces - Get the number of spaces in the sum space
6: Input Parameter:
7: . sp - the dual space object
9: Output Parameter:
10: . numSumSpaces - the number of spaces
12: Level: intermediate
14: Note:
15: The name NumSubspaces is slightly misleading because it is actually getting the number of defining spaces of the sum, not a number of Subspaces of it
17: .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetNumSubspaces()`
18: @*/
19: PetscErrorCode PetscDualSpaceSumGetNumSubspaces(PetscDualSpace sp, PetscInt *numSumSpaces)
20: {
21: PetscFunctionBegin;
23: PetscAssertPointer(numSumSpaces, 2);
24: PetscTryMethod(sp, "PetscDualSpaceSumGetNumSubspaces_C", (PetscDualSpace, PetscInt *), (sp, numSumSpaces));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /*@
29: PetscDualSpaceSumSetNumSubspaces - Set the number of spaces in the sum space
31: Input Parameters:
32: + sp - the dual space object
33: - numSumSpaces - the number of spaces
35: Level: intermediate
37: Note:
38: The name NumSubspaces is slightly misleading because it is actually setting the number of defining spaces of the sum, not a number of Subspaces of it
40: .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetNumSubspaces()`
41: @*/
42: PetscErrorCode PetscDualSpaceSumSetNumSubspaces(PetscDualSpace sp, PetscInt numSumSpaces)
43: {
44: PetscFunctionBegin;
46: PetscTryMethod(sp, "PetscDualSpaceSumSetNumSubspaces_C", (PetscDualSpace, PetscInt), (sp, numSumSpaces));
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: /*@
51: PetscDualSpaceSumGetConcatenate - Get the concatenate flag for this space.
53: Input Parameter:
54: . sp - the dual space object
56: Output Parameter:
57: . concatenate - flag indicating whether subspaces are concatenated.
59: Level: intermediate
61: Note:
62: A concatenated sum space will have the number of components equal to the sum of the number of
63: components of all subspaces. A non-concatenated, or direct sum space will have the same
64: number of components as its subspaces.
66: .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetConcatenate()`
67: @*/
68: PetscErrorCode PetscDualSpaceSumGetConcatenate(PetscDualSpace sp, PetscBool *concatenate)
69: {
70: PetscFunctionBegin;
72: PetscTryMethod(sp, "PetscDualSpaceSumGetConcatenate_C", (PetscDualSpace, PetscBool *), (sp, concatenate));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@
77: PetscDualSpaceSumSetConcatenate - Sets the concatenate flag for this space.
79: Input Parameters:
80: + sp - the dual space object
81: - concatenate - are subspaces concatenated components (true) or direct summands (false)
83: Level: intermediate
85: Notes:
86: A concatenated sum space will have the number of components equal to the sum of the number of
87: components of all subspaces. A non-concatenated, or direct sum space will have the same
88: number of components as its subspaces .
90: .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetConcatenate()`
91: @*/
92: PetscErrorCode PetscDualSpaceSumSetConcatenate(PetscDualSpace sp, PetscBool concatenate)
93: {
94: PetscFunctionBegin;
96: PetscTryMethod(sp, "PetscDualSpaceSumSetConcatenate_C", (PetscDualSpace, PetscBool), (sp, concatenate));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: /*@
101: PetscDualSpaceSumGetSubspace - Get a space in the sum space
103: Input Parameters:
104: + sp - the dual space object
105: - s - The space number
107: Output Parameter:
108: . subsp - the `PetscDualSpace`
110: Level: intermediate
112: Note:
113: The name GetSubspace is slightly misleading because it is actually getting one of the defining spaces of the sum, not a Subspace of it
115: .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetSubspace()`
116: @*/
117: PetscErrorCode PetscDualSpaceSumGetSubspace(PetscDualSpace sp, PetscInt s, PetscDualSpace *subsp)
118: {
119: PetscFunctionBegin;
121: PetscAssertPointer(subsp, 3);
122: PetscTryMethod(sp, "PetscDualSpaceSumGetSubspace_C", (PetscDualSpace, PetscInt, PetscDualSpace *), (sp, s, subsp));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /*@
127: PetscDualSpaceSumSetSubspace - Set a space in the sum space
129: Input Parameters:
130: + sp - the dual space object
131: . s - The space number
132: - subsp - the number of spaces
134: Level: intermediate
136: Note:
137: The name SetSubspace is slightly misleading because it is actually setting one of the defining spaces of the sum, not a Subspace of it
139: .seealso: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetSubspace()`
140: @*/
141: PetscErrorCode PetscDualSpaceSumSetSubspace(PetscDualSpace sp, PetscInt s, PetscDualSpace subsp)
142: {
143: PetscFunctionBegin;
146: PetscTryMethod(sp, "PetscDualSpaceSumSetSubspace_C", (PetscDualSpace, PetscInt, PetscDualSpace), (sp, s, subsp));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode PetscDualSpaceSumGetNumSubspaces_Sum(PetscDualSpace space, PetscInt *numSumSpaces)
151: {
152: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
154: PetscFunctionBegin;
155: *numSumSpaces = sum->numSumSpaces;
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode PetscDualSpaceSumSetNumSubspaces_Sum(PetscDualSpace space, PetscInt numSumSpaces)
160: {
161: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
162: PetscInt Ns = sum->numSumSpaces;
164: PetscFunctionBegin;
165: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
166: if (numSumSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
167: for (PetscInt s = 0; s < Ns; ++s) PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[s]));
168: PetscCall(PetscFree(sum->sumspaces));
170: Ns = sum->numSumSpaces = numSumSpaces;
171: PetscCall(PetscCalloc1(Ns, &sum->sumspaces));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: static PetscErrorCode PetscDualSpaceSumGetConcatenate_Sum(PetscDualSpace sp, PetscBool *concatenate)
176: {
177: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
179: PetscFunctionBegin;
180: *concatenate = sum->concatenate;
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode PetscDualSpaceSumSetConcatenate_Sum(PetscDualSpace sp, PetscBool concatenate)
185: {
186: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
188: PetscFunctionBegin;
189: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change space concatenation after setup called.");
191: sum->concatenate = concatenate;
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: static PetscErrorCode PetscDualSpaceSumGetSubspace_Sum(PetscDualSpace space, PetscInt s, PetscDualSpace *subspace)
196: {
197: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
198: PetscInt Ns = sum->numSumSpaces;
200: PetscFunctionBegin;
201: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscDualSpaceSumSetNumSubspaces() first");
202: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
204: *subspace = sum->sumspaces[s];
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode PetscDualSpaceSumSetSubspace_Sum(PetscDualSpace space, PetscInt s, PetscDualSpace subspace)
209: {
210: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
211: PetscInt Ns = sum->numSumSpaces;
213: PetscFunctionBegin;
214: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
215: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscDualSpaceSumSetNumSubspaces() first");
216: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
218: PetscCall(PetscObjectReference((PetscObject)subspace));
219: PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[s]));
220: sum->sumspaces[s] = subspace;
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: static PetscErrorCode PetscDualSpaceDuplicate_Sum(PetscDualSpace sp, PetscDualSpace spNew)
225: {
226: PetscInt num_subspaces, Nc;
227: PetscBool concatenate, interleave_basis, interleave_components;
228: PetscDualSpace subsp_first = NULL;
229: PetscDualSpace subsp_dup_first = NULL;
230: DM K;
232: PetscFunctionBegin;
233: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &num_subspaces));
234: PetscCall(PetscDualSpaceSumSetNumSubspaces(spNew, num_subspaces));
235: PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
236: PetscCall(PetscDualSpaceSumSetConcatenate(spNew, concatenate));
237: PetscCall(PetscDualSpaceSumGetInterleave(sp, &interleave_basis, &interleave_components));
238: PetscCall(PetscDualSpaceSumSetInterleave(spNew, interleave_basis, interleave_components));
239: PetscCall(PetscDualSpaceGetDM(sp, &K));
240: PetscCall(PetscDualSpaceSetDM(spNew, K));
241: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
242: PetscCall(PetscDualSpaceSetNumComponents(spNew, Nc));
243: for (PetscInt s = 0; s < num_subspaces; s++) {
244: PetscDualSpace subsp, subspNew;
246: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
247: if (s == 0) {
248: subsp_first = subsp;
249: PetscCall(PetscDualSpaceDuplicate(subsp, &subsp_dup_first));
250: subspNew = subsp_dup_first;
251: } else if (subsp == subsp_first) {
252: PetscCall(PetscObjectReference((PetscObject)subsp_dup_first));
253: subspNew = subsp_dup_first;
254: } else {
255: PetscCall(PetscDualSpaceDuplicate(subsp, &subspNew));
256: }
257: PetscCall(PetscDualSpaceSumSetSubspace(spNew, s, subspNew));
258: PetscCall(PetscDualSpaceDestroy(&subspNew));
259: }
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode PetscDualSpaceSumCreateQuadrature(PetscInt Ns, PetscInt cdim, PetscBool uniform_points, PetscQuadrature subquads[], PetscQuadrature *fullquad)
264: {
265: PetscReal *points;
266: PetscInt Npoints;
268: PetscFunctionBegin;
269: if (uniform_points) {
270: PetscCall(PetscObjectReference((PetscObject)subquads[0]));
271: *fullquad = subquads[0];
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
274: Npoints = 0;
275: for (PetscInt s = 0; s < Ns; s++) {
276: PetscInt subNpoints;
278: if (!subquads[s]) continue;
279: PetscCall(PetscQuadratureGetData(subquads[s], NULL, NULL, &subNpoints, NULL, NULL));
280: Npoints += subNpoints;
281: }
282: *fullquad = NULL;
283: if (!Npoints) PetscFunctionReturn(PETSC_SUCCESS);
284: PetscCall(PetscMalloc1(Npoints * cdim, &points));
285: for (PetscInt s = 0, offset = 0; s < Ns; s++) {
286: PetscInt subNpoints;
287: const PetscReal *subpoints;
289: PetscCall(PetscQuadratureGetData(subquads[s], NULL, NULL, &subNpoints, &subpoints, NULL));
290: PetscCall(PetscArraycpy(&points[offset], subpoints, cdim * subNpoints));
291: offset += cdim * subNpoints;
292: }
293: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, fullquad));
294: PetscCall(PetscQuadratureSetData(*fullquad, cdim, 0, Npoints, points, NULL));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode PetscDualSpaceSumCreateMatrix(PetscDualSpace sp, Mat submats[], ISLocalToGlobalMapping map_rows[], ISLocalToGlobalMapping map_cols[], PetscQuadrature fullquad, Mat *fullmat)
299: {
300: Mat mat;
301: PetscInt *i = NULL, *j = NULL;
302: PetscScalar *v = NULL;
303: PetscInt npoints;
304: PetscInt nrows, ncols, nnz;
305: PetscInt Nc, num_subspaces;
307: PetscFunctionBegin;
308: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
309: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &num_subspaces));
310: nrows = 0;
311: ncols = 0;
312: nnz = 0;
313: PetscCall(PetscQuadratureGetData(fullquad, NULL, NULL, &npoints, NULL, NULL));
314: ncols = Nc * npoints;
315: for (PetscInt s = 0; s < num_subspaces; s++) {
316: // Get the COO data for each matrix, map the is and js, and append to growing COO data
317: PetscInt sNrows, sNcols;
318: Mat smat;
319: const PetscInt *si;
320: const PetscInt *sj;
321: PetscScalar *sv;
322: PetscMemType memtype;
323: PetscInt snz;
324: PetscInt snz_actual;
325: PetscInt *cooi;
326: PetscInt *cooj;
327: PetscScalar *coov;
328: PetscScalar *v_new;
329: PetscInt *i_new;
330: PetscInt *j_new;
332: if (!submats[s]) continue;
333: PetscCall(MatGetSize(submats[s], &sNrows, &sNcols));
334: nrows += sNrows;
335: PetscCall(MatConvert(submats[s], MATSEQAIJ, MAT_INITIAL_MATRIX, &smat));
336: PetscCall(MatBindToCPU(smat, PETSC_TRUE));
337: PetscCall(MatSeqAIJGetCSRAndMemType(smat, &si, &sj, &sv, &memtype));
338: PetscCheck(memtype == PETSC_MEMTYPE_HOST, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not convert matrix to host memory");
339: snz = si[sNrows];
341: PetscCall(PetscMalloc1(nnz + snz, &v_new));
342: PetscCall(PetscArraycpy(v_new, v, nnz));
343: PetscCall(PetscFree(v));
344: v = v_new;
346: PetscCall(PetscMalloc1(nnz + snz, &i_new));
347: PetscCall(PetscArraycpy(i_new, i, nnz));
348: PetscCall(PetscFree(i));
349: i = i_new;
351: PetscCall(PetscMalloc1(nnz + snz, &j_new));
352: PetscCall(PetscArraycpy(j_new, j, nnz));
353: PetscCall(PetscFree(j));
354: j = j_new;
356: PetscCall(PetscMalloc2(snz, &cooi, snz, &cooj));
358: coov = &v[nnz];
360: snz_actual = 0;
361: for (PetscInt row = 0; row < sNrows; row++) {
362: for (PetscInt k = si[row]; k < si[row + 1]; k++, snz_actual++) {
363: cooi[snz_actual] = row;
364: cooj[snz_actual] = sj[k];
365: coov[snz_actual] = sv[k];
366: }
367: }
368: PetscCall(MatDestroy(&smat));
370: if (snz_actual > 0) {
371: PetscCall(ISLocalToGlobalMappingApply(map_rows[s], snz_actual, cooi, &i[nnz]));
372: PetscCall(ISLocalToGlobalMappingApply(map_cols[s], snz_actual, cooj, &j[nnz]));
373: nnz += snz_actual;
374: }
375: PetscCall(PetscFree2(cooi, cooj));
376: }
377: PetscCall(MatCreate(PETSC_COMM_SELF, fullmat));
378: mat = *fullmat;
379: PetscCall(MatSetSizes(mat, nrows, ncols, nrows, ncols));
380: PetscCall(MatSetType(mat, MATSEQAIJ));
381: PetscCall(MatSetPreallocationCOO(mat, nnz, i, j));
382: PetscCall(MatSetValuesCOO(mat, v, INSERT_VALUES));
383: PetscCall(PetscFree(i));
384: PetscCall(PetscFree(j));
385: PetscCall(PetscFree(v));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode PetscDualSpaceSumCreateMappings(PetscDualSpace sp, PetscBool interior, PetscBool uniform_points, ISLocalToGlobalMapping map_row[], ISLocalToGlobalMapping map_col[])
390: {
391: PetscSection section;
392: PetscInt pStart, pEnd, Ns, Nc;
393: PetscInt num_points = 0;
394: PetscBool interleave_basis, interleave_components, concatenate;
395: PetscInt *roffset;
397: PetscFunctionBegin;
398: if (interior) {
399: PetscCall(PetscDualSpaceGetInteriorSection(sp, §ion));
400: } else {
401: PetscCall(PetscDualSpaceGetSection(sp, §ion));
402: }
403: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
404: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
405: PetscCall(PetscDualSpaceSumGetInterleave(sp, &interleave_basis, &interleave_components));
406: PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
407: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
408: for (PetscInt p = pStart; p < pEnd; p++) {
409: PetscInt dof;
411: PetscCall(PetscSectionGetDof(section, p, &dof));
412: num_points += (dof > 0);
413: }
414: PetscCall(PetscCalloc1(pEnd - pStart, &roffset));
415: for (PetscInt s = 0, coffset = 0; s < Ns; s++) {
416: PetscDualSpace subsp;
417: PetscSection subsection;
418: IS is_row, is_col;
419: PetscInt subNb;
420: PetscInt sNc, sNpoints, sNcols;
421: PetscQuadrature quad;
423: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
424: PetscCall(PetscDualSpaceGetNumComponents(subsp, &sNc));
425: if (interior) {
426: PetscCall(PetscDualSpaceGetInteriorSection(subsp, &subsection));
427: PetscCall(PetscDualSpaceGetInteriorData(subsp, &quad, NULL));
428: } else {
429: PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
430: PetscCall(PetscDualSpaceGetAllData(subsp, &quad, NULL));
431: }
432: PetscCall(PetscSectionGetStorageSize(subsection, &subNb));
433: if (num_points == 1) {
434: PetscInt rstride;
436: rstride = interleave_basis ? Ns : 1;
437: PetscCall(ISCreateStride(PETSC_COMM_SELF, subNb, roffset[0], rstride, &is_row));
438: roffset[0] += interleave_basis ? 1 : subNb;
439: } else {
440: PetscInt *rows;
442: PetscCall(PetscMalloc1(subNb, &rows));
443: for (PetscInt p = pStart; p < pEnd; p++) {
444: PetscInt subdof, dof, off, suboff, stride;
446: PetscCall(PetscSectionGetOffset(subsection, p, &suboff));
447: PetscCall(PetscSectionGetDof(subsection, p, &subdof));
448: PetscCall(PetscSectionGetOffset(section, p, &off));
449: PetscCall(PetscSectionGetDof(section, p, &dof));
450: PetscCheck(subdof * Ns == dof || !interleave_basis, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Basis cannot be interleaved\n");
451: stride = interleave_basis ? Ns : 1;
452: for (PetscInt k = 0; k < subdof; k++) { rows[suboff + k] = off + roffset[p - pStart] + k * stride; }
453: roffset[p - pStart] += interleave_basis ? 1 : subdof;
454: }
455: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, subNb, rows, PETSC_OWN_POINTER, &is_row));
456: }
458: sNpoints = 0;
459: if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &sNpoints, NULL, NULL));
460: sNcols = sNpoints * sNc;
462: if (!concatenate) {
463: PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, 1, &is_col));
464: coffset += uniform_points ? 0 : sNcols;
465: } else {
466: if (uniform_points && interleave_components) {
467: PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, Ns, &is_col));
468: coffset += 1;
469: } else {
470: PetscInt *cols;
472: PetscCall(PetscMalloc1(sNcols, &cols));
473: for (PetscInt p = 0, r = 0; p < sNpoints; p++) {
474: for (PetscInt c = 0; c < sNc; c++, r++) { cols[r] = coffset + p * Nc + c; }
475: }
476: coffset += uniform_points ? sNc : Nc * sNpoints + sNc;
477: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sNcols, cols, PETSC_OWN_POINTER, &is_col));
478: }
479: }
480: PetscCall(ISLocalToGlobalMappingCreateIS(is_row, &map_row[s]));
481: PetscCall(ISLocalToGlobalMappingCreateIS(is_col, &map_col[s]));
482: PetscCall(ISDestroy(&is_row));
483: PetscCall(ISDestroy(&is_col));
484: }
485: PetscCall(PetscFree(roffset));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Sum(PetscDualSpace sp, PetscInt f, PetscDualSpace *bdsp)
490: {
491: PetscInt k, Nc, Nk, Nknew, Ncnew, Ns;
492: PetscInt dim, pointDim = -1;
493: PetscInt depth, Ncopies;
494: PetscBool any;
495: DM dm, K;
496: DMPolytopeType ct;
498: PetscFunctionBegin;
499: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
500: any = PETSC_FALSE;
501: for (PetscInt s = 0; s < Ns; s++) {
502: PetscDualSpace subsp, subpointsp;
504: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
505: PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
506: if (subpointsp) any = PETSC_TRUE;
507: }
508: if (!any) {
509: *bdsp = NULL;
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
512: PetscCall(PetscDualSpaceGetDM(sp, &dm));
513: PetscCall(DMGetDimension(dm, &dim));
514: PetscCall(DMPlexGetDepth(dm, &depth));
515: PetscCall(PetscDualSpaceDuplicate(sp, bdsp));
516: PetscCheck((depth == dim) || (depth == 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
517: pointDim = (depth == dim) ? (dim - 1) : 0;
518: PetscCall(DMPlexGetCellType(dm, f, &ct));
519: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
520: PetscCall(PetscDualSpaceSetDM(*bdsp, K));
521: PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
522: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
523: PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
524: Ncopies = Nc / Nk;
525: PetscCall(PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew));
526: Ncnew = Nknew * Ncopies;
527: PetscCall(PetscDualSpaceSetNumComponents(*bdsp, Ncnew));
528: for (PetscInt s = 0; s < Ns; s++) {
529: PetscDualSpace subsp, subpointsp;
531: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
532: PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
533: if (subpointsp) {
534: PetscCall(PetscObjectReference((PetscObject)subpointsp));
535: } else {
536: // make an empty dual space
537: PetscInt subNc, subNcopies, subpointNc;
539: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subsp), &subpointsp));
540: PetscCall(PetscDualSpaceGetNumComponents(subsp, &subNc));
541: subNcopies = subNc / Nk;
542: subpointNc = subNcopies * Nknew;
543: PetscCall(PetscDualSpaceSetType(subpointsp, PETSCDUALSPACESIMPLE));
544: PetscCall(PetscDualSpaceSimpleSetDimension(subpointsp, 0));
545: PetscCall(PetscDualSpaceSetFormDegree(subpointsp, k));
546: PetscCall(PetscDualSpaceSetNumComponents(subpointsp, subpointNc));
547: }
548: // K should be equal to subpointsp->dm, but we want it to be exactly the same
549: PetscCall(PetscObjectReference((PetscObject)K));
550: PetscCall(DMDestroy(&subpointsp->dm));
551: subpointsp->dm = K;
552: PetscCall(PetscDualSpaceSetUp(subpointsp));
553: PetscCall(PetscDualSpaceSumSetSubspace(*bdsp, s, subpointsp));
554: PetscCall(PetscDualSpaceDestroy(&subpointsp));
555: }
556: PetscCall(DMDestroy(&K));
557: PetscCall(PetscDualSpaceSetUp(*bdsp));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode PetscDualSpaceSumIsUniform(PetscDualSpace sp, PetscBool *is_uniform)
562: {
563: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
564: PetscBool uniform = PETSC_TRUE;
566: PetscFunctionBegin;
567: for (PetscInt s = 1; s < sum->numSumSpaces; s++) {
568: if (sum->sumspaces[s] != sum->sumspaces[0]) {
569: uniform = PETSC_FALSE;
570: break;
571: }
572: }
573: *is_uniform = uniform;
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode PetscDualSpaceGetSymmetries_Sum(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
578: {
579: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
581: PetscFunctionBegin;
582: if (!sum->symComputed) {
583: PetscInt Ns;
584: PetscBool any_perms = PETSC_FALSE;
585: PetscBool any_flips = PETSC_FALSE;
586: PetscInt ***symperms = NULL;
587: PetscScalar ***symflips = NULL;
589: sum->symComputed = PETSC_TRUE;
590: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
591: for (PetscInt s = 0; s < Ns; s++) {
592: PetscDualSpace subsp;
593: const PetscInt ***sub_perms;
594: const PetscScalar ***sub_flips;
596: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
597: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
598: if (sub_perms) any_perms = PETSC_TRUE;
599: if (sub_flips) any_flips = PETSC_TRUE;
600: }
601: if (any_perms || any_flips) {
602: DM K;
603: PetscInt pStart, pEnd, numPoints;
604: PetscInt spintdim;
606: PetscCall(PetscDualSpaceGetDM(sp, &K));
607: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
608: numPoints = pEnd - pStart;
609: PetscCall(PetscCalloc1(numPoints, &symperms));
610: PetscCall(PetscCalloc1(numPoints, &symflips));
611: PetscCall(PetscDualSpaceGetBoundarySymmetries_Internal(sp, symperms, symflips));
612: // get interior symmetries
613: PetscCall(PetscDualSpaceGetInteriorDimension(sp, &spintdim));
614: if (spintdim) {
615: PetscInt groupSize;
616: PetscInt **cellPerms;
617: PetscScalar **cellFlips;
618: DMPolytopeType ct;
620: PetscCall(DMPlexGetCellType(K, 0, &ct));
621: groupSize = DMPolytopeTypeGetNumArrangments(ct);
622: sum->numSelfSym = groupSize;
623: sum->selfSymOff = groupSize / 2;
624: PetscCall(PetscCalloc1(groupSize, &cellPerms));
625: PetscCall(PetscCalloc1(groupSize, &cellFlips));
626: symperms[0] = &cellPerms[groupSize / 2];
627: symflips[0] = &cellFlips[groupSize / 2];
628: for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
629: PetscBool any_o_perms = PETSC_FALSE;
630: PetscBool any_o_flips = PETSC_FALSE;
632: for (PetscInt s = 0; s < Ns; s++) {
633: PetscDualSpace subsp;
634: const PetscInt ***sub_perms;
635: const PetscScalar ***sub_flips;
637: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
638: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
639: if (sub_perms && sub_perms[0] && sub_perms[0][o]) any_o_perms = PETSC_TRUE;
640: if (sub_flips && sub_flips[0] && sub_flips[0][o]) any_o_flips = PETSC_TRUE;
641: }
642: if (any_o_perms) {
643: PetscInt *o_perm;
645: PetscCall(PetscMalloc1(spintdim, &o_perm));
646: for (PetscInt i = 0; i < spintdim; i++) o_perm[i] = i;
647: for (PetscInt s = 0; s < Ns; s++) {
648: PetscDualSpace subsp;
649: const PetscInt ***sub_perms;
650: const PetscScalar ***sub_flips;
652: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
653: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
654: if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
655: PetscInt subspdim;
656: PetscInt *range, *domain;
657: PetscInt *range_mapped, *domain_mapped;
659: PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
660: PetscCall(PetscMalloc4(subspdim, &range, subspdim, &range_mapped, subspdim, &domain, subspdim, &domain_mapped));
661: for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
662: PetscCall(PetscArraycpy(range, sub_perms[0][o], subspdim));
663: PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
664: PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, range, range_mapped));
665: for (PetscInt i = 0; i < subspdim; i++) o_perm[domain_mapped[i]] = range_mapped[i];
666: PetscCall(PetscFree4(range, range_mapped, domain, domain_mapped));
667: }
668: }
669: symperms[0][o] = o_perm;
670: }
671: if (any_o_flips) {
672: PetscScalar *o_flip;
674: PetscCall(PetscMalloc1(spintdim, &o_flip));
675: for (PetscInt i = 0; i < spintdim; i++) o_flip[i] = 1.0;
676: for (PetscInt s = 0; s < Ns; s++) {
677: PetscDualSpace subsp;
678: const PetscInt ***sub_perms;
679: const PetscScalar ***sub_flips;
681: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
682: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
683: if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
684: PetscInt subspdim;
685: PetscInt *domain;
686: PetscInt *domain_mapped;
688: PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
689: PetscCall(PetscMalloc2(subspdim, &domain, subspdim, &domain_mapped));
690: for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
691: PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
692: for (PetscInt i = 0; i < subspdim; i++) o_flip[domain_mapped[i]] = sub_perms[0][o][i];
693: PetscCall(PetscFree2(domain, domain_mapped));
694: }
695: }
696: symflips[0][o] = o_flip;
697: }
698: }
699: {
700: PetscBool any_perms = PETSC_FALSE;
701: PetscBool any_flips = PETSC_FALSE;
702: for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
703: if (symperms[0][o]) any_perms = PETSC_TRUE;
704: if (symflips[0][o]) any_flips = PETSC_TRUE;
705: }
706: if (!any_perms) {
707: PetscCall(PetscFree(cellPerms));
708: symperms[0] = NULL;
709: }
710: if (!any_flips) {
711: PetscCall(PetscFree(cellFlips));
712: symflips[0] = NULL;
713: }
714: }
715: }
716: if (!any_perms) {
717: PetscCall(PetscFree(symperms));
718: symperms = NULL;
719: }
720: if (!any_flips) {
721: PetscCall(PetscFree(symflips));
722: symflips = NULL;
723: }
724: }
725: sum->symperms = symperms;
726: sum->symflips = symflips;
727: }
728: if (perms) *perms = (const PetscInt ***)sum->symperms;
729: if (flips) *flips = (const PetscScalar ***)sum->symflips;
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: static PetscErrorCode PetscDualSpaceSetUp_Sum(PetscDualSpace sp)
734: {
735: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
736: PetscBool concatenate = PETSC_TRUE;
737: PetscBool uniform;
738: PetscInt Ns, Nc, i, sum_Nc = 0;
739: PetscInt minNc, maxNc;
740: PetscInt minForm, maxForm, cdim, depth;
741: DM K;
742: PetscQuadrature *all_quads = NULL;
743: PetscQuadrature *int_quads = NULL;
744: Mat *all_mats = NULL;
745: Mat *int_mats = NULL;
747: PetscFunctionBegin;
748: if (sum->setupCalled) PetscFunctionReturn(PETSC_SUCCESS);
749: sum->setupCalled = PETSC_TRUE;
751: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
752: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have %" PetscInt_FMT " subspaces", Ns);
754: // step 1: make sure they share a DM
755: PetscCall(PetscDualSpaceGetDM(sp, &K));
756: PetscCall(DMGetCoordinateDim(K, &cdim));
757: PetscCall(DMPlexGetDepth(K, &depth));
758: PetscCall(PetscDualSpaceSumIsUniform(sp, &sp->uniform));
759: uniform = sp->uniform;
760: {
761: for (PetscInt s = 0; s < Ns; s++) {
762: PetscDualSpace subsp;
763: DM sub_K;
765: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
766: PetscCall(PetscDualSpaceSetUp(subsp));
767: PetscCall(PetscDualSpaceGetDM(subsp, &sub_K));
768: if (s == 0 && K == NULL) {
769: PetscCall(PetscDualSpaceSetDM(sp, sub_K));
770: K = sub_K;
771: }
772: PetscCheck(sub_K == K, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %d does not have the same DM as the sum space", (int)s);
773: }
774: }
776: // step 2: count components
777: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
778: PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
779: minNc = Nc;
780: maxNc = Nc;
781: minForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MAX : sp->k;
782: maxForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MIN : sp->k;
783: for (i = 0; i < Ns; ++i) {
784: PetscInt sNc, formDegree;
785: PetscDualSpace si;
787: PetscCall(PetscDualSpaceSumGetSubspace(sp, i, &si));
788: PetscCall(PetscDualSpaceSetUp(si));
789: PetscCall(PetscDualSpaceGetNumComponents(si, &sNc));
790: if (sNc != Nc) PetscCheck(concatenate, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace as a different number of components but space does not concatenate components");
791: minNc = PetscMin(minNc, sNc);
792: maxNc = PetscMax(maxNc, sNc);
793: sum_Nc += sNc;
794: PetscCall(PetscDualSpaceGetFormDegree(si, &formDegree));
795: minForm = PetscMin(minForm, formDegree);
796: maxForm = PetscMax(maxForm, formDegree);
797: }
798: sp->k = (minForm != maxForm) ? PETSC_FORM_DEGREE_UNDEFINED : minForm;
800: if (concatenate) PetscCheck(sum_Nc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Total number of subspace components (%" PetscInt_FMT ") does not match number of target space components (%" PetscInt_FMT ").", sum_Nc, Nc);
801: else PetscCheck(minNc == Nc && maxNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspaces must have same number of components as the target space.");
803: /* A PetscDualSpace should have a fixed number of components, but
804: if the spaces we are combining have different form degrees, they will not
805: have the same number of components on subcomponents of the boundary,
806: so we do not try to create boundary dual spaces in this case */
807: if (sp->k != PETSC_FORM_DEGREE_UNDEFINED && depth > 0) {
808: PetscInt pStart, pEnd;
809: PetscInt *pStratStart, *pStratEnd;
811: PetscCall(DMPlexGetDepth(K, &depth));
812: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
813: PetscCall(PetscCalloc1(pEnd, &(sp->pointSpaces)));
814: PetscCall(PetscMalloc2(depth + 1, &pStratStart, depth + 1, &pStratEnd));
815: for (PetscInt d = 0; d <= depth; ++d) PetscCall(DMPlexGetDepthStratum(K, d, &pStratStart[d], &pStratEnd[d]));
817: for (PetscInt p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
818: PetscReal v0[3];
819: DMPolytopeType ptype;
820: PetscReal J[9], detJ;
821: PetscInt q;
823: PetscCall(DMPlexComputeCellGeometryAffineFEM(K, p, v0, J, NULL, &detJ));
824: PetscCall(DMPlexGetCellType(K, p, &ptype));
826: /* compare to previous facets: if computed, reference that dualspace */
827: for (q = pStratStart[depth - 1]; q < p; q++) {
828: DMPolytopeType qtype;
830: PetscCall(DMPlexGetCellType(K, q, &qtype));
831: if (qtype == ptype) break;
832: }
833: if (q < p) { /* this facet has the same dual space as that one */
834: PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[q]));
835: sp->pointSpaces[p] = sp->pointSpaces[q];
836: continue;
837: }
838: /* if not, recursively compute this dual space */
839: PetscCall(PetscDualSpaceCreateFacetSubspace_Sum(sp, p, &sp->pointSpaces[p]));
840: }
841: for (PetscInt h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
842: PetscInt hd = depth - h;
844: for (PetscInt p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
845: PetscInt suppSize;
846: const PetscInt *supp;
847: DM qdm;
848: PetscDualSpace qsp, psp;
849: PetscInt c, coneSize, q;
850: const PetscInt *cone;
851: const PetscInt *refCone;
853: PetscCall(DMPlexGetSupportSize(K, p, &suppSize));
854: PetscCall(DMPlexGetSupport(K, p, &supp));
855: q = supp[0];
856: qsp = sp->pointSpaces[q];
857: PetscCall(DMPlexGetConeSize(K, q, &coneSize));
858: PetscCall(DMPlexGetCone(K, q, &cone));
859: for (c = 0; c < coneSize; c++)
860: if (cone[c] == p) break;
861: PetscCheck(c != coneSize, PetscObjectComm((PetscObject)K), PETSC_ERR_PLIB, "cone/support mismatch");
862: if (!qsp) {
863: sp->pointSpaces[p] = NULL;
864: continue;
865: }
866: PetscCall(PetscDualSpaceGetDM(qsp, &qdm));
867: PetscCall(DMPlexGetCone(qdm, 0, &refCone));
868: /* get the equivalent dual space from the support dual space */
869: PetscCall(PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp));
870: PetscCall(PetscObjectReference((PetscObject)psp));
871: sp->pointSpaces[p] = psp;
872: }
873: }
874: PetscCall(PetscFree2(pStratStart, pStratEnd));
875: }
877: sum->uniform = uniform;
878: PetscCall(PetscCalloc1(Ns, &sum->all_rows));
879: PetscCall(PetscCalloc1(Ns, &sum->all_cols));
880: PetscCall(PetscCalloc1(Ns, &sum->int_rows));
881: PetscCall(PetscCalloc1(Ns, &sum->int_cols));
882: PetscCall(PetscMalloc4(Ns, &all_quads, Ns, &all_mats, Ns, &int_quads, Ns, &int_mats));
883: {
884: // test for uniform all points and uniform interior points
885: PetscBool uniform_all = PETSC_TRUE;
886: PetscBool uniform_interior = PETSC_TRUE;
887: PetscQuadrature quad_all_first = NULL;
888: PetscQuadrature quad_interior_first = NULL;
889: PetscInt pStart, pEnd;
891: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
892: PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &sp->pointSection));
894: for (PetscInt p = pStart; p < pEnd; p++) {
895: PetscInt full_dof = 0;
897: for (PetscInt s = 0; s < Ns; s++) {
898: PetscDualSpace subsp;
899: PetscSection subsection;
900: PetscInt dof;
902: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
903: PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
904: PetscCall(PetscSectionGetDof(subsection, p, &dof));
905: full_dof += dof;
906: }
907: PetscCall(PetscSectionSetDof(sp->pointSection, p, full_dof));
908: }
909: PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection));
911: for (PetscInt s = 0; s < Ns; s++) {
912: PetscDualSpace subsp;
913: PetscQuadrature subquad_all;
914: PetscQuadrature subquad_interior;
915: Mat submat_all;
916: Mat submat_interior;
918: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
919: PetscCall(PetscDualSpaceGetAllData(subsp, &subquad_all, &submat_all));
920: PetscCall(PetscDualSpaceGetInteriorData(subsp, &subquad_interior, &submat_interior));
921: if (!s) {
922: quad_all_first = subquad_all;
923: quad_interior_first = subquad_interior;
924: } else {
925: if (subquad_all != quad_all_first) uniform_all = PETSC_FALSE;
926: if (subquad_interior != quad_interior_first) uniform_interior = PETSC_FALSE;
927: }
928: all_quads[s] = subquad_all;
929: int_quads[s] = subquad_interior;
930: all_mats[s] = submat_all;
931: int_mats[s] = submat_interior;
932: }
933: sum->uniform_all_points = uniform_all;
934: sum->uniform_interior_points = uniform_interior;
936: PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_TRUE, uniform_interior, sum->int_rows, sum->int_cols));
937: PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_interior, int_quads, &sp->intNodes));
938: if (sp->intNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, int_mats, sum->int_rows, sum->int_cols, sp->intNodes, &sp->intMat));
940: PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_FALSE, uniform_all, sum->all_rows, sum->all_cols));
941: PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_all, all_quads, &sp->allNodes));
942: if (sp->allNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, all_mats, sum->all_rows, sum->all_cols, sp->allNodes, &sp->allMat));
943: }
944: PetscCall(PetscFree4(all_quads, all_mats, int_quads, int_mats));
945: PetscCall(PetscDualSpaceComputeFunctionalsFromAllData(sp));
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: static PetscErrorCode PetscDualSpaceSumView_Ascii(PetscDualSpace sp, PetscViewer v)
950: {
951: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
952: PetscBool concatenate = sum->concatenate;
953: PetscInt i, Ns = sum->numSumSpaces;
955: PetscFunctionBegin;
956: if (concatenate) PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
957: else PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
958: for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
959: PetscCall(PetscViewerASCIIPushTab(v));
960: PetscCall(PetscDualSpaceView(sum->sumspaces[i], v));
961: PetscCall(PetscViewerASCIIPopTab(v));
962: }
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: static PetscErrorCode PetscDualSpaceView_Sum(PetscDualSpace sp, PetscViewer viewer)
967: {
968: PetscBool iascii;
970: PetscFunctionBegin;
971: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
972: if (iascii) PetscCall(PetscDualSpaceSumView_Ascii(sp, viewer));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: static PetscErrorCode PetscDualSpaceDestroy_Sum(PetscDualSpace sp)
977: {
978: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
979: PetscInt i, Ns = sum->numSumSpaces;
981: PetscFunctionBegin;
982: if (sum->symperms) {
983: PetscInt **selfSyms = sum->symperms[0];
985: if (selfSyms) {
986: PetscInt i, **allocated = &selfSyms[-sum->selfSymOff];
988: for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
989: PetscCall(PetscFree(allocated));
990: }
991: PetscCall(PetscFree(sum->symperms));
992: }
993: if (sum->symflips) {
994: PetscScalar **selfSyms = sum->symflips[0];
996: if (selfSyms) {
997: PetscInt i;
998: PetscScalar **allocated = &selfSyms[-sum->selfSymOff];
1000: for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
1001: PetscCall(PetscFree(allocated));
1002: }
1003: PetscCall(PetscFree(sum->symflips));
1004: }
1005: for (i = 0; i < Ns; ++i) {
1006: PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[i]));
1007: if (sum->all_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_rows[i]));
1008: if (sum->all_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_cols[i]));
1009: if (sum->int_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_rows[i]));
1010: if (sum->int_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_cols[i]));
1011: }
1012: PetscCall(PetscFree(sum->sumspaces));
1013: PetscCall(PetscFree(sum->all_rows));
1014: PetscCall(PetscFree(sum->all_cols));
1015: PetscCall(PetscFree(sum->int_rows));
1016: PetscCall(PetscFree(sum->int_cols));
1017: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", NULL));
1018: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", NULL));
1019: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", NULL));
1020: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", NULL));
1021: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", NULL));
1022: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", NULL));
1023: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", NULL));
1024: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", NULL));
1025: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL));
1026: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL));
1027: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL));
1028: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL));
1029: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL));
1030: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL));
1031: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", NULL));
1032: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", NULL));
1033: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL));
1034: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL));
1035: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL));
1036: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL));
1037: PetscCall(PetscFree(sum));
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: /*@
1042: PetscDualSpaceSumSetInterleave - Set whether the basis functions and components of a uniform sum are interleaved
1044: Logically collective
1046: Input Parameters:
1047: + sp - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1048: . interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1049: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1050: interleave the concatenated components
1052: Level: developer
1054: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumGetInterleave()`
1055: @*/
1056: PetscErrorCode PetscDualSpaceSumSetInterleave(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1057: {
1058: PetscFunctionBegin;
1060: PetscTryMethod(sp, "PetscDualSpaceSumSetInterleave_C", (PetscDualSpace, PetscBool, PetscBool), (sp, interleave_basis, interleave_components));
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: static PetscErrorCode PetscDualSpaceSumSetInterleave_Sum(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1065: {
1066: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
1067: PetscFunctionBegin;
1068: sum->interleave_basis = interleave_basis;
1069: sum->interleave_components = interleave_components;
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1073: /*@
1074: PetscDualSpaceSumGetInterleave - Get whether the basis functions and components of a uniform sum are interleaved
1076: Logically collective
1078: Input Parameter:
1079: . sp - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1081: Output Parameters:
1082: + interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1083: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1084: interleave the concatenated components
1086: Level: developer
1088: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumSetInterleave()`
1089: @*/
1090: PetscErrorCode PetscDualSpaceSumGetInterleave(PetscDualSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
1091: {
1092: PetscFunctionBegin;
1094: if (interleave_basis) PetscAssertPointer(interleave_basis, 2);
1095: if (interleave_components) PetscAssertPointer(interleave_components, 3);
1096: PetscTryMethod(sp, "PetscDualSpaceSumGetInterleave_C", (PetscDualSpace, PetscBool *, PetscBool *), (sp, interleave_basis, interleave_components));
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: static PetscErrorCode PetscDualSpaceSumGetInterleave_Sum(PetscDualSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
1101: {
1102: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
1103: PetscFunctionBegin;
1104: if (interleave_basis) *interleave_basis = sum->interleave_basis;
1105: if (interleave_components) *interleave_components = sum->interleave_components;
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }
1109: #define PetscDualSpaceSumPassthrough(sp, func, ...) \
1110: do { \
1111: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data; \
1112: PetscBool is_uniform; \
1113: PetscCall(PetscDualSpaceSumIsUniform(sp, &is_uniform)); \
1114: if (is_uniform && sum->numSumSpaces > 0) { \
1115: PetscDualSpace subsp; \
1116: PetscCall(PetscDualSpaceSumGetSubspace(sp, 0, &subsp)); \
1117: PetscCall(func(subsp, __VA_ARGS__)); \
1118: } \
1119: } while (0)
1121: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Sum(PetscDualSpace sp, PetscBool *value)
1122: {
1123: PetscFunctionBegin;
1124: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetContinuity, value);
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Sum(PetscDualSpace sp, PetscBool value)
1129: {
1130: PetscFunctionBegin;
1131: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetContinuity, value);
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Sum(PetscDualSpace sp, PetscBool *value)
1136: {
1137: PetscFunctionBegin;
1138: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTensor, value);
1139: PetscFunctionReturn(PETSC_SUCCESS);
1140: }
1142: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Sum(PetscDualSpace sp, PetscBool value)
1143: {
1144: PetscFunctionBegin;
1145: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTensor, value);
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Sum(PetscDualSpace sp, PetscBool *value)
1150: {
1151: PetscFunctionBegin;
1152: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTrimmed, value);
1153: PetscFunctionReturn(PETSC_SUCCESS);
1154: }
1156: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Sum(PetscDualSpace sp, PetscBool value)
1157: {
1158: PetscFunctionBegin;
1159: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTrimmed, value);
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Sum(PetscDualSpace sp, PetscBool *value)
1164: {
1165: PetscFunctionBegin;
1166: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetUseMoments, value);
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Sum(PetscDualSpace sp, PetscBool value)
1171: {
1172: PetscFunctionBegin;
1173: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetUseMoments, value);
1174: PetscFunctionReturn(PETSC_SUCCESS);
1175: }
1177: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Sum(PetscDualSpace sp, PetscInt *value)
1178: {
1179: PetscFunctionBegin;
1180: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetMomentOrder, value);
1181: PetscFunctionReturn(PETSC_SUCCESS);
1182: }
1184: static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Sum(PetscDualSpace sp, PetscInt value)
1185: {
1186: PetscFunctionBegin;
1187: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetMomentOrder, value);
1188: PetscFunctionReturn(PETSC_SUCCESS);
1189: }
1191: static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType *node_type, PetscBool *include_endpoints, PetscReal *exponent)
1192: {
1193: PetscFunctionBegin;
1194: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetNodeType, node_type, include_endpoints, exponent);
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType node_type, PetscBool include_endpoints, PetscReal exponent)
1199: {
1200: PetscFunctionBegin;
1201: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetNodeType, node_type, include_endpoints, exponent);
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: static PetscErrorCode PetscDualSpaceInitialize_Sum(PetscDualSpace sp)
1206: {
1207: PetscFunctionBegin;
1208: sp->ops->destroy = PetscDualSpaceDestroy_Sum;
1209: sp->ops->view = PetscDualSpaceView_Sum;
1210: sp->ops->setfromoptions = NULL;
1211: sp->ops->duplicate = PetscDualSpaceDuplicate_Sum;
1212: sp->ops->setup = PetscDualSpaceSetUp_Sum;
1213: sp->ops->createheightsubspace = NULL;
1214: sp->ops->createpointsubspace = NULL;
1215: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Sum;
1216: sp->ops->apply = PetscDualSpaceApplyDefault;
1217: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
1218: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
1219: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
1220: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
1222: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", PetscDualSpaceSumGetNumSubspaces_Sum));
1223: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", PetscDualSpaceSumSetNumSubspaces_Sum));
1224: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", PetscDualSpaceSumGetSubspace_Sum));
1225: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", PetscDualSpaceSumSetSubspace_Sum));
1226: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", PetscDualSpaceSumGetConcatenate_Sum));
1227: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", PetscDualSpaceSumSetConcatenate_Sum));
1228: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", PetscDualSpaceSumGetInterleave_Sum));
1229: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", PetscDualSpaceSumSetInterleave_Sum));
1230: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Sum));
1231: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Sum));
1232: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Sum));
1233: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Sum));
1234: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Sum));
1235: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Sum));
1236: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Sum));
1237: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Sum));
1238: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Sum));
1239: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Sum));
1240: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Sum));
1241: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Sum));
1242: PetscFunctionReturn(PETSC_SUCCESS);
1243: }
1245: /*MC
1246: PETSCDUALSPACESUM = "sum" - A `PetscDualSpace` object that encapsulates a sum of subspaces.
1248: Level: intermediate
1250: Note:
1251: That sum can either be direct or a concatenation. For example if A and B are spaces each with 2 components,
1252: the direct sum of A and B will also have 2 components while the concatenated sum will have 4 components. In both cases A and B must be defined over the
1253: same reference element.
1255: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`, `PetscDualSpaceSumGetNumSubspaces()`, `PetscDualSpaceSumSetNumSubspaces()`,
1256: `PetscDualSpaceSumGetConcatenate()`, `PetscDualSpaceSumSetConcatenate()`, `PetscDualSpaceSumSetInterleave()`, `PetscDualSpaceSumGetInterleave()`
1257: M*/
1258: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Sum(PetscDualSpace sp)
1259: {
1260: PetscDualSpace_Sum *sum;
1262: PetscFunctionBegin;
1264: PetscCall(PetscNew(&sum));
1265: sum->numSumSpaces = PETSC_DEFAULT;
1266: sp->data = sum;
1267: sp->k = PETSC_FORM_DEGREE_UNDEFINED;
1268: PetscCall(PetscDualSpaceInitialize_Sum(sp));
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: /*@
1273: PetscDualSpaceCreateSum - Create a finite element dual basis that is the sum of other dual bases
1275: Collective
1277: Input Parameters:
1278: + numSubspaces - the number of spaces that will be added together
1279: . subspaces - an array of length `numSubspaces` of spaces
1280: - concatenate - if `PETSC_FALSE`, the sum-space has the same components as the individual dual spaces (`PetscDualSpaceGetNumComponents()`); if `PETSC_TRUE`, the individual components are concatenated to create a dual space with more components
1282: Output Parameter:
1283: . sumSpace - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1285: Level: advanced
1287: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCSPACESUM`
1288: @*/
1289: PetscErrorCode PetscDualSpaceCreateSum(PetscInt numSubspaces, const PetscDualSpace subspaces[], PetscBool concatenate, PetscDualSpace *sumSpace)
1290: {
1291: PetscInt i, Nc = 0;
1293: PetscFunctionBegin;
1294: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace));
1295: PetscCall(PetscDualSpaceSetType(*sumSpace, PETSCDUALSPACESUM));
1296: PetscCall(PetscDualSpaceSumSetNumSubspaces(*sumSpace, numSubspaces));
1297: PetscCall(PetscDualSpaceSumSetConcatenate(*sumSpace, concatenate));
1298: for (i = 0; i < numSubspaces; ++i) {
1299: PetscInt sNc;
1301: PetscCall(PetscDualSpaceSumSetSubspace(*sumSpace, i, subspaces[i]));
1302: PetscCall(PetscDualSpaceGetNumComponents(subspaces[i], &sNc));
1303: if (concatenate) Nc += sNc;
1304: else Nc = sNc;
1306: if (i == 0) {
1307: DM dm;
1309: PetscCall(PetscDualSpaceGetDM(subspaces[i], &dm));
1310: PetscCall(PetscDualSpaceSetDM(*sumSpace, dm));
1311: }
1312: }
1313: PetscCall(PetscDualSpaceSetNumComponents(*sumSpace, Nc));
1314: PetscFunctionReturn(PETSC_SUCCESS);
1315: }