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 = 0, ncols, nnz = 0;
305: PetscInt Nc, num_subspaces;
307: PetscFunctionBegin;
308: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
309: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &num_subspaces));
310: PetscCall(PetscQuadratureGetData(fullquad, NULL, NULL, &npoints, NULL, NULL));
311: ncols = Nc * npoints;
312: for (PetscInt s = 0; s < num_subspaces; s++) {
313: // Get the COO data for each matrix, map the is and js, and append to growing COO data
314: PetscInt sNrows, sNcols;
315: Mat smat;
316: const PetscInt *si;
317: const PetscInt *sj;
318: PetscScalar *sv;
319: PetscMemType memtype;
320: PetscInt snz;
321: PetscInt snz_actual;
322: PetscInt *cooi;
323: PetscInt *cooj;
324: PetscScalar *coov;
325: PetscScalar *v_new;
326: PetscInt *i_new;
327: PetscInt *j_new;
329: if (!submats[s]) continue;
330: PetscCall(MatGetSize(submats[s], &sNrows, &sNcols));
331: nrows += sNrows;
332: PetscCall(MatConvert(submats[s], MATSEQAIJ, MAT_INITIAL_MATRIX, &smat));
333: PetscCall(MatBindToCPU(smat, PETSC_TRUE));
334: PetscCall(MatSeqAIJGetCSRAndMemType(smat, &si, &sj, &sv, &memtype));
335: PetscCheck(memtype == PETSC_MEMTYPE_HOST, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not convert matrix to host memory");
336: snz = si[sNrows];
338: PetscCall(PetscMalloc1(nnz + snz, &v_new));
339: PetscCall(PetscArraycpy(v_new, v, nnz));
340: PetscCall(PetscFree(v));
341: v = v_new;
343: PetscCall(PetscMalloc1(nnz + snz, &i_new));
344: PetscCall(PetscArraycpy(i_new, i, nnz));
345: PetscCall(PetscFree(i));
346: i = i_new;
348: PetscCall(PetscMalloc1(nnz + snz, &j_new));
349: PetscCall(PetscArraycpy(j_new, j, nnz));
350: PetscCall(PetscFree(j));
351: j = j_new;
353: PetscCall(PetscMalloc2(snz, &cooi, snz, &cooj));
355: coov = &v[nnz];
357: snz_actual = 0;
358: for (PetscInt row = 0; row < sNrows; row++) {
359: for (PetscInt k = si[row]; k < si[row + 1]; k++, snz_actual++) {
360: cooi[snz_actual] = row;
361: cooj[snz_actual] = sj[k];
362: coov[snz_actual] = sv[k];
363: }
364: }
365: PetscCall(MatDestroy(&smat));
367: if (snz_actual > 0) {
368: PetscCall(ISLocalToGlobalMappingApply(map_rows[s], snz_actual, cooi, &i[nnz]));
369: PetscCall(ISLocalToGlobalMappingApply(map_cols[s], snz_actual, cooj, &j[nnz]));
370: nnz += snz_actual;
371: }
372: PetscCall(PetscFree2(cooi, cooj));
373: }
374: PetscCall(MatCreate(PETSC_COMM_SELF, fullmat));
375: mat = *fullmat;
376: PetscCall(MatSetSizes(mat, nrows, ncols, nrows, ncols));
377: PetscCall(MatSetType(mat, MATSEQAIJ));
378: PetscCall(MatSetPreallocationCOO(mat, nnz, i, j));
379: PetscCall(MatSetValuesCOO(mat, v, INSERT_VALUES));
380: PetscCall(PetscFree(i));
381: PetscCall(PetscFree(j));
382: PetscCall(PetscFree(v));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode PetscDualSpaceSumCreateMappings(PetscDualSpace sp, PetscBool interior, PetscBool uniform_points, ISLocalToGlobalMapping map_row[], ISLocalToGlobalMapping map_col[])
387: {
388: PetscSection section;
389: PetscInt pStart, pEnd, Ns, Nc;
390: PetscInt num_points = 0;
391: PetscBool interleave_basis, interleave_components, concatenate;
392: PetscInt *roffset;
394: PetscFunctionBegin;
395: if (interior) {
396: PetscCall(PetscDualSpaceGetInteriorSection(sp, §ion));
397: } else {
398: PetscCall(PetscDualSpaceGetSection(sp, §ion));
399: }
400: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
401: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
402: PetscCall(PetscDualSpaceSumGetInterleave(sp, &interleave_basis, &interleave_components));
403: PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
404: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
405: for (PetscInt p = pStart; p < pEnd; p++) {
406: PetscInt dof;
408: PetscCall(PetscSectionGetDof(section, p, &dof));
409: num_points += (dof > 0);
410: }
411: PetscCall(PetscCalloc1(pEnd - pStart, &roffset));
412: for (PetscInt s = 0, coffset = 0; s < Ns; s++) {
413: PetscDualSpace subsp;
414: PetscSection subsection;
415: IS is_row, is_col;
416: PetscInt subNb;
417: PetscInt sNc, sNpoints, sNcols;
418: PetscQuadrature quad;
420: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
421: PetscCall(PetscDualSpaceGetNumComponents(subsp, &sNc));
422: if (interior) {
423: PetscCall(PetscDualSpaceGetInteriorSection(subsp, &subsection));
424: PetscCall(PetscDualSpaceGetInteriorData(subsp, &quad, NULL));
425: } else {
426: PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
427: PetscCall(PetscDualSpaceGetAllData(subsp, &quad, NULL));
428: }
429: PetscCall(PetscSectionGetStorageSize(subsection, &subNb));
430: if (num_points == 1) {
431: PetscInt rstride;
433: rstride = interleave_basis ? Ns : 1;
434: PetscCall(ISCreateStride(PETSC_COMM_SELF, subNb, roffset[0], rstride, &is_row));
435: roffset[0] += interleave_basis ? 1 : subNb;
436: } else {
437: PetscInt *rows;
439: PetscCall(PetscMalloc1(subNb, &rows));
440: for (PetscInt p = pStart; p < pEnd; p++) {
441: PetscInt subdof, dof, off, suboff, stride;
443: PetscCall(PetscSectionGetOffset(subsection, p, &suboff));
444: PetscCall(PetscSectionGetDof(subsection, p, &subdof));
445: PetscCall(PetscSectionGetOffset(section, p, &off));
446: PetscCall(PetscSectionGetDof(section, p, &dof));
447: PetscCheck(subdof * Ns == dof || !interleave_basis, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Basis cannot be interleaved");
448: stride = interleave_basis ? Ns : 1;
449: for (PetscInt k = 0; k < subdof; k++) { rows[suboff + k] = off + roffset[p - pStart] + k * stride; }
450: roffset[p - pStart] += interleave_basis ? 1 : subdof;
451: }
452: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, subNb, rows, PETSC_OWN_POINTER, &is_row));
453: }
455: sNpoints = 0;
456: if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &sNpoints, NULL, NULL));
457: sNcols = sNpoints * sNc;
459: if (!concatenate) {
460: PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, 1, &is_col));
461: coffset += uniform_points ? 0 : sNcols;
462: } else {
463: if (uniform_points && interleave_components) {
464: PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, Ns, &is_col));
465: coffset += 1;
466: } else {
467: PetscInt *cols;
469: PetscCall(PetscMalloc1(sNcols, &cols));
470: for (PetscInt p = 0, r = 0; p < sNpoints; p++) {
471: for (PetscInt c = 0; c < sNc; c++, r++) { cols[r] = coffset + p * Nc + c; }
472: }
473: coffset += uniform_points ? sNc : Nc * sNpoints + sNc;
474: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sNcols, cols, PETSC_OWN_POINTER, &is_col));
475: }
476: }
477: PetscCall(ISLocalToGlobalMappingCreateIS(is_row, &map_row[s]));
478: PetscCall(ISLocalToGlobalMappingCreateIS(is_col, &map_col[s]));
479: PetscCall(ISDestroy(&is_row));
480: PetscCall(ISDestroy(&is_col));
481: }
482: PetscCall(PetscFree(roffset));
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Sum(PetscDualSpace sp, PetscInt f, PetscDualSpace *bdsp)
487: {
488: PetscInt k, Nc, Nk, Nknew, Ncnew, Ns;
489: PetscInt dim, pointDim = -1;
490: PetscInt depth, Ncopies;
491: PetscBool any;
492: DM dm, K;
493: DMPolytopeType ct;
495: PetscFunctionBegin;
496: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
497: any = PETSC_FALSE;
498: for (PetscInt s = 0; s < Ns; s++) {
499: PetscDualSpace subsp, subpointsp;
501: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
502: PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
503: if (subpointsp) any = PETSC_TRUE;
504: }
505: if (!any) {
506: *bdsp = NULL;
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
509: PetscCall(PetscDualSpaceGetDM(sp, &dm));
510: PetscCall(DMGetDimension(dm, &dim));
511: PetscCall(DMPlexGetDepth(dm, &depth));
512: PetscCall(PetscDualSpaceDuplicate(sp, bdsp));
513: PetscCheck((depth == dim) || (depth == 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
514: pointDim = (depth == dim) ? (dim - 1) : 0;
515: PetscCall(DMPlexGetCellType(dm, f, &ct));
516: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
517: PetscCall(PetscDualSpaceSetDM(*bdsp, K));
518: PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
519: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
520: PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
521: Ncopies = Nc / Nk;
522: PetscCall(PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew));
523: Ncnew = Nknew * Ncopies;
524: PetscCall(PetscDualSpaceSetNumComponents(*bdsp, Ncnew));
525: for (PetscInt s = 0; s < Ns; s++) {
526: PetscDualSpace subsp, subpointsp;
528: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
529: PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
530: if (subpointsp) {
531: PetscCall(PetscObjectReference((PetscObject)subpointsp));
532: } else {
533: // make an empty dual space
534: PetscInt subNc, subNcopies, subpointNc;
536: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subsp), &subpointsp));
537: PetscCall(PetscDualSpaceGetNumComponents(subsp, &subNc));
538: subNcopies = subNc / Nk;
539: subpointNc = subNcopies * Nknew;
540: PetscCall(PetscDualSpaceSetType(subpointsp, PETSCDUALSPACESIMPLE));
541: PetscCall(PetscDualSpaceSimpleSetDimension(subpointsp, 0));
542: PetscCall(PetscDualSpaceSetFormDegree(subpointsp, k));
543: PetscCall(PetscDualSpaceSetNumComponents(subpointsp, subpointNc));
544: }
545: // K should be equal to subpointsp->dm, but we want it to be exactly the same
546: PetscCall(PetscObjectReference((PetscObject)K));
547: PetscCall(DMDestroy(&subpointsp->dm));
548: subpointsp->dm = K;
549: PetscCall(PetscDualSpaceSetUp(subpointsp));
550: PetscCall(PetscDualSpaceSumSetSubspace(*bdsp, s, subpointsp));
551: PetscCall(PetscDualSpaceDestroy(&subpointsp));
552: }
553: PetscCall(DMDestroy(&K));
554: PetscCall(PetscDualSpaceSetUp(*bdsp));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: static PetscErrorCode PetscDualSpaceSumIsUniform(PetscDualSpace sp, PetscBool *is_uniform)
559: {
560: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
561: PetscBool uniform = PETSC_TRUE;
563: PetscFunctionBegin;
564: for (PetscInt s = 1; s < sum->numSumSpaces; s++) {
565: if (sum->sumspaces[s] != sum->sumspaces[0]) {
566: uniform = PETSC_FALSE;
567: break;
568: }
569: }
570: *is_uniform = uniform;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: static PetscErrorCode PetscDualSpaceGetSymmetries_Sum(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
575: {
576: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
578: PetscFunctionBegin;
579: if (!sum->symComputed) {
580: PetscInt Ns;
581: PetscBool any_perms = PETSC_FALSE;
582: PetscBool any_flips = PETSC_FALSE;
583: PetscInt ***symperms = NULL;
584: PetscScalar ***symflips = NULL;
586: sum->symComputed = PETSC_TRUE;
587: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
588: for (PetscInt s = 0; s < Ns; s++) {
589: PetscDualSpace subsp;
590: const PetscInt ***sub_perms;
591: const PetscScalar ***sub_flips;
593: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
594: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
595: if (sub_perms) any_perms = PETSC_TRUE;
596: if (sub_flips) any_flips = PETSC_TRUE;
597: }
598: if (any_perms || any_flips) {
599: DM K;
600: PetscInt pStart, pEnd, numPoints;
601: PetscInt spintdim;
603: PetscCall(PetscDualSpaceGetDM(sp, &K));
604: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
605: numPoints = pEnd - pStart;
606: PetscCall(PetscCalloc1(numPoints, &symperms));
607: PetscCall(PetscCalloc1(numPoints, &symflips));
608: PetscCall(PetscDualSpaceGetBoundarySymmetries_Internal(sp, symperms, symflips));
609: // get interior symmetries
610: PetscCall(PetscDualSpaceGetInteriorDimension(sp, &spintdim));
611: if (spintdim) {
612: PetscInt groupSize;
613: PetscInt **cellPerms;
614: PetscScalar **cellFlips;
615: DMPolytopeType ct;
617: PetscCall(DMPlexGetCellType(K, 0, &ct));
618: groupSize = DMPolytopeTypeGetNumArrangements(ct);
619: sum->numSelfSym = groupSize;
620: sum->selfSymOff = groupSize / 2;
621: PetscCall(PetscCalloc1(groupSize, &cellPerms));
622: PetscCall(PetscCalloc1(groupSize, &cellFlips));
623: symperms[0] = &cellPerms[groupSize / 2];
624: symflips[0] = &cellFlips[groupSize / 2];
625: for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
626: PetscBool any_o_perms = PETSC_FALSE;
627: PetscBool any_o_flips = PETSC_FALSE;
629: for (PetscInt s = 0; s < Ns; s++) {
630: PetscDualSpace subsp;
631: const PetscInt ***sub_perms;
632: const PetscScalar ***sub_flips;
634: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
635: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
636: if (sub_perms && sub_perms[0] && sub_perms[0][o]) any_o_perms = PETSC_TRUE;
637: if (sub_flips && sub_flips[0] && sub_flips[0][o]) any_o_flips = PETSC_TRUE;
638: }
639: if (any_o_perms) {
640: PetscInt *o_perm;
642: PetscCall(PetscMalloc1(spintdim, &o_perm));
643: for (PetscInt i = 0; i < spintdim; i++) o_perm[i] = i;
644: for (PetscInt s = 0; s < Ns; s++) {
645: PetscDualSpace subsp;
646: const PetscInt ***sub_perms;
647: const PetscScalar ***sub_flips;
649: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
650: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
651: if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
652: PetscInt subspdim;
653: PetscInt *range, *domain;
654: PetscInt *range_mapped, *domain_mapped;
656: PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
657: PetscCall(PetscMalloc4(subspdim, &range, subspdim, &range_mapped, subspdim, &domain, subspdim, &domain_mapped));
658: for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
659: PetscCall(PetscArraycpy(range, sub_perms[0][o], subspdim));
660: PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
661: PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, range, range_mapped));
662: for (PetscInt i = 0; i < subspdim; i++) o_perm[domain_mapped[i]] = range_mapped[i];
663: PetscCall(PetscFree4(range, range_mapped, domain, domain_mapped));
664: }
665: }
666: symperms[0][o] = o_perm;
667: }
668: if (any_o_flips) {
669: PetscScalar *o_flip;
671: PetscCall(PetscMalloc1(spintdim, &o_flip));
672: for (PetscInt i = 0; i < spintdim; i++) o_flip[i] = 1.0;
673: for (PetscInt s = 0; s < Ns; s++) {
674: PetscDualSpace subsp;
675: const PetscInt ***sub_perms;
676: const PetscScalar ***sub_flips;
678: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
679: PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
680: if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
681: PetscInt subspdim;
682: PetscInt *domain;
683: PetscInt *domain_mapped;
685: PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
686: PetscCall(PetscMalloc2(subspdim, &domain, subspdim, &domain_mapped));
687: for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
688: PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
689: for (PetscInt i = 0; i < subspdim; i++) o_flip[domain_mapped[i]] = sub_perms[0][o][i];
690: PetscCall(PetscFree2(domain, domain_mapped));
691: }
692: }
693: symflips[0][o] = o_flip;
694: }
695: }
696: {
697: PetscBool any_perms = PETSC_FALSE;
698: PetscBool any_flips = PETSC_FALSE;
699: for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
700: if (symperms[0][o]) any_perms = PETSC_TRUE;
701: if (symflips[0][o]) any_flips = PETSC_TRUE;
702: }
703: if (!any_perms) {
704: PetscCall(PetscFree(cellPerms));
705: symperms[0] = NULL;
706: }
707: if (!any_flips) {
708: PetscCall(PetscFree(cellFlips));
709: symflips[0] = NULL;
710: }
711: }
712: }
713: if (!any_perms) {
714: PetscCall(PetscFree(symperms));
715: symperms = NULL;
716: }
717: if (!any_flips) {
718: PetscCall(PetscFree(symflips));
719: symflips = NULL;
720: }
721: }
722: sum->symperms = symperms;
723: sum->symflips = symflips;
724: }
725: if (perms) *perms = (const PetscInt ***)sum->symperms;
726: if (flips) *flips = (const PetscScalar ***)sum->symflips;
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: static PetscErrorCode PetscDualSpaceSetUp_Sum(PetscDualSpace sp)
731: {
732: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
733: PetscBool concatenate = PETSC_TRUE;
734: PetscBool uniform;
735: PetscInt Ns, Nc, i, sum_Nc = 0;
736: PetscInt minNc, maxNc;
737: PetscInt minForm, maxForm, cdim, depth;
738: DM K;
739: PetscQuadrature *all_quads = NULL;
740: PetscQuadrature *int_quads = NULL;
741: Mat *all_mats = NULL;
742: Mat *int_mats = NULL;
744: PetscFunctionBegin;
745: if (sum->setupCalled) PetscFunctionReturn(PETSC_SUCCESS);
746: sum->setupCalled = PETSC_TRUE;
748: PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
749: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have %" PetscInt_FMT " subspaces", Ns);
751: // step 1: make sure they share a DM
752: PetscCall(PetscDualSpaceGetDM(sp, &K));
753: PetscCall(DMGetCoordinateDim(K, &cdim));
754: PetscCall(DMPlexGetDepth(K, &depth));
755: PetscCall(PetscDualSpaceSumIsUniform(sp, &sp->uniform));
756: uniform = sp->uniform;
757: {
758: for (PetscInt s = 0; s < Ns; s++) {
759: PetscDualSpace subsp;
760: DM sub_K;
762: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
763: PetscCall(PetscDualSpaceSetUp(subsp));
764: PetscCall(PetscDualSpaceGetDM(subsp, &sub_K));
765: if (s == 0 && K == NULL) {
766: PetscCall(PetscDualSpaceSetDM(sp, sub_K));
767: K = sub_K;
768: }
769: 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);
770: }
771: }
773: // step 2: count components
774: PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
775: PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
776: minNc = Nc;
777: maxNc = Nc;
778: minForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MAX : sp->k;
779: maxForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MIN : sp->k;
780: for (i = 0; i < Ns; ++i) {
781: PetscInt sNc, formDegree;
782: PetscDualSpace si;
784: PetscCall(PetscDualSpaceSumGetSubspace(sp, i, &si));
785: PetscCall(PetscDualSpaceSetUp(si));
786: PetscCall(PetscDualSpaceGetNumComponents(si, &sNc));
787: 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");
788: minNc = PetscMin(minNc, sNc);
789: maxNc = PetscMax(maxNc, sNc);
790: sum_Nc += sNc;
791: PetscCall(PetscDualSpaceGetFormDegree(si, &formDegree));
792: minForm = PetscMin(minForm, formDegree);
793: maxForm = PetscMax(maxForm, formDegree);
794: }
795: sp->k = (minForm != maxForm) ? PETSC_FORM_DEGREE_UNDEFINED : minForm;
797: 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);
798: else PetscCheck(minNc == Nc && maxNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspaces must have same number of components as the target space.");
800: /* A PetscDualSpace should have a fixed number of components, but
801: if the spaces we are combining have different form degrees, they will not
802: have the same number of components on subcomponents of the boundary,
803: so we do not try to create boundary dual spaces in this case */
804: if (sp->k != PETSC_FORM_DEGREE_UNDEFINED && depth > 0) {
805: PetscInt pStart, pEnd;
806: PetscInt *pStratStart, *pStratEnd;
808: PetscCall(DMPlexGetDepth(K, &depth));
809: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
810: PetscCall(PetscCalloc1(pEnd, &sp->pointSpaces));
811: PetscCall(PetscMalloc2(depth + 1, &pStratStart, depth + 1, &pStratEnd));
812: for (PetscInt d = 0; d <= depth; ++d) PetscCall(DMPlexGetDepthStratum(K, d, &pStratStart[d], &pStratEnd[d]));
814: for (PetscInt p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
815: PetscReal v0[3];
816: DMPolytopeType ptype;
817: PetscReal J[9], detJ;
818: PetscInt q;
820: PetscCall(DMPlexComputeCellGeometryAffineFEM(K, p, v0, J, NULL, &detJ));
821: PetscCall(DMPlexGetCellType(K, p, &ptype));
823: /* compare to previous facets: if computed, reference that dualspace */
824: for (q = pStratStart[depth - 1]; q < p; q++) {
825: DMPolytopeType qtype;
827: PetscCall(DMPlexGetCellType(K, q, &qtype));
828: if (qtype == ptype) break;
829: }
830: if (q < p) { /* this facet has the same dual space as that one */
831: PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[q]));
832: sp->pointSpaces[p] = sp->pointSpaces[q];
833: continue;
834: }
835: /* if not, recursively compute this dual space */
836: PetscCall(PetscDualSpaceCreateFacetSubspace_Sum(sp, p, &sp->pointSpaces[p]));
837: }
838: for (PetscInt h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
839: PetscInt hd = depth - h;
841: for (PetscInt p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
842: PetscInt suppSize;
843: const PetscInt *supp;
844: DM qdm;
845: PetscDualSpace qsp, psp;
846: PetscInt c, coneSize, q;
847: const PetscInt *cone;
848: const PetscInt *refCone;
850: PetscCall(DMPlexGetSupportSize(K, p, &suppSize));
851: PetscCall(DMPlexGetSupport(K, p, &supp));
852: q = supp[0];
853: qsp = sp->pointSpaces[q];
854: PetscCall(DMPlexGetConeSize(K, q, &coneSize));
855: PetscCall(DMPlexGetCone(K, q, &cone));
856: for (c = 0; c < coneSize; c++)
857: if (cone[c] == p) break;
858: PetscCheck(c != coneSize, PetscObjectComm((PetscObject)K), PETSC_ERR_PLIB, "cone/support mismatch");
859: if (!qsp) {
860: sp->pointSpaces[p] = NULL;
861: continue;
862: }
863: PetscCall(PetscDualSpaceGetDM(qsp, &qdm));
864: PetscCall(DMPlexGetCone(qdm, 0, &refCone));
865: /* get the equivalent dual space from the support dual space */
866: PetscCall(PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp));
867: PetscCall(PetscObjectReference((PetscObject)psp));
868: sp->pointSpaces[p] = psp;
869: }
870: }
871: PetscCall(PetscFree2(pStratStart, pStratEnd));
872: }
874: sum->uniform = uniform;
875: PetscCall(PetscCalloc1(Ns, &sum->all_rows));
876: PetscCall(PetscCalloc1(Ns, &sum->all_cols));
877: PetscCall(PetscCalloc1(Ns, &sum->int_rows));
878: PetscCall(PetscCalloc1(Ns, &sum->int_cols));
879: PetscCall(PetscMalloc4(Ns, &all_quads, Ns, &all_mats, Ns, &int_quads, Ns, &int_mats));
880: {
881: // test for uniform all points and uniform interior points
882: PetscBool uniform_all = PETSC_TRUE;
883: PetscBool uniform_interior = PETSC_TRUE;
884: PetscQuadrature quad_all_first = NULL;
885: PetscQuadrature quad_interior_first = NULL;
886: PetscInt pStart, pEnd;
888: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
889: PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &sp->pointSection));
891: for (PetscInt p = pStart; p < pEnd; p++) {
892: PetscInt full_dof = 0;
894: for (PetscInt s = 0; s < Ns; s++) {
895: PetscDualSpace subsp;
896: PetscSection subsection;
897: PetscInt dof;
899: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
900: PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
901: PetscCall(PetscSectionGetDof(subsection, p, &dof));
902: full_dof += dof;
903: }
904: PetscCall(PetscSectionSetDof(sp->pointSection, p, full_dof));
905: }
906: PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection));
908: for (PetscInt s = 0; s < Ns; s++) {
909: PetscDualSpace subsp;
910: PetscQuadrature subquad_all;
911: PetscQuadrature subquad_interior;
912: Mat submat_all;
913: Mat submat_interior;
915: PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
916: PetscCall(PetscDualSpaceGetAllData(subsp, &subquad_all, &submat_all));
917: PetscCall(PetscDualSpaceGetInteriorData(subsp, &subquad_interior, &submat_interior));
918: if (!s) {
919: quad_all_first = subquad_all;
920: quad_interior_first = subquad_interior;
921: } else {
922: if (subquad_all != quad_all_first) uniform_all = PETSC_FALSE;
923: if (subquad_interior != quad_interior_first) uniform_interior = PETSC_FALSE;
924: }
925: all_quads[s] = subquad_all;
926: int_quads[s] = subquad_interior;
927: all_mats[s] = submat_all;
928: int_mats[s] = submat_interior;
929: }
930: sum->uniform_all_points = uniform_all;
931: sum->uniform_interior_points = uniform_interior;
933: PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_TRUE, uniform_interior, sum->int_rows, sum->int_cols));
934: PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_interior, int_quads, &sp->intNodes));
935: if (sp->intNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, int_mats, sum->int_rows, sum->int_cols, sp->intNodes, &sp->intMat));
937: PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_FALSE, uniform_all, sum->all_rows, sum->all_cols));
938: PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_all, all_quads, &sp->allNodes));
939: if (sp->allNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, all_mats, sum->all_rows, sum->all_cols, sp->allNodes, &sp->allMat));
940: }
941: PetscCall(PetscFree4(all_quads, all_mats, int_quads, int_mats));
942: PetscCall(PetscDualSpaceComputeFunctionalsFromAllData(sp));
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: static PetscErrorCode PetscDualSpaceSumView_Ascii(PetscDualSpace sp, PetscViewer v)
947: {
948: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
949: PetscBool concatenate = sum->concatenate;
950: PetscInt i, Ns = sum->numSumSpaces;
952: PetscFunctionBegin;
953: if (concatenate) PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
954: else PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
955: for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
956: PetscCall(PetscViewerASCIIPushTab(v));
957: PetscCall(PetscDualSpaceView(sum->sumspaces[i], v));
958: PetscCall(PetscViewerASCIIPopTab(v));
959: }
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: static PetscErrorCode PetscDualSpaceView_Sum(PetscDualSpace sp, PetscViewer viewer)
964: {
965: PetscBool iascii;
967: PetscFunctionBegin;
968: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
969: if (iascii) PetscCall(PetscDualSpaceSumView_Ascii(sp, viewer));
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
973: static PetscErrorCode PetscDualSpaceDestroy_Sum(PetscDualSpace sp)
974: {
975: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
976: PetscInt i, Ns = sum->numSumSpaces;
978: PetscFunctionBegin;
979: if (sum->symperms) {
980: PetscInt **selfSyms = sum->symperms[0];
982: if (selfSyms) {
983: PetscInt i, **allocated = &selfSyms[-sum->selfSymOff];
985: for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
986: PetscCall(PetscFree(allocated));
987: }
988: PetscCall(PetscFree(sum->symperms));
989: }
990: if (sum->symflips) {
991: PetscScalar **selfSyms = sum->symflips[0];
993: if (selfSyms) {
994: PetscInt i;
995: PetscScalar **allocated = &selfSyms[-sum->selfSymOff];
997: for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
998: PetscCall(PetscFree(allocated));
999: }
1000: PetscCall(PetscFree(sum->symflips));
1001: }
1002: for (i = 0; i < Ns; ++i) {
1003: PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[i]));
1004: if (sum->all_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_rows[i]));
1005: if (sum->all_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_cols[i]));
1006: if (sum->int_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_rows[i]));
1007: if (sum->int_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_cols[i]));
1008: }
1009: PetscCall(PetscFree(sum->sumspaces));
1010: PetscCall(PetscFree(sum->all_rows));
1011: PetscCall(PetscFree(sum->all_cols));
1012: PetscCall(PetscFree(sum->int_rows));
1013: PetscCall(PetscFree(sum->int_cols));
1014: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", NULL));
1015: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", NULL));
1016: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", NULL));
1017: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", NULL));
1018: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", NULL));
1019: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", NULL));
1020: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", NULL));
1021: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", NULL));
1022: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL));
1023: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL));
1024: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL));
1025: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL));
1026: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL));
1027: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL));
1028: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", NULL));
1029: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", NULL));
1030: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL));
1031: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL));
1032: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL));
1033: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL));
1034: PetscCall(PetscFree(sum));
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: /*@
1039: PetscDualSpaceSumSetInterleave - Set whether the basis functions and components of a uniform sum are interleaved
1041: Logically collective
1043: Input Parameters:
1044: + sp - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1045: . interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1046: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1047: interleave the concatenated components
1049: Level: developer
1051: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumGetInterleave()`
1052: @*/
1053: PetscErrorCode PetscDualSpaceSumSetInterleave(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1054: {
1055: PetscFunctionBegin;
1057: PetscTryMethod(sp, "PetscDualSpaceSumSetInterleave_C", (PetscDualSpace, PetscBool, PetscBool), (sp, interleave_basis, interleave_components));
1058: PetscFunctionReturn(PETSC_SUCCESS);
1059: }
1061: static PetscErrorCode PetscDualSpaceSumSetInterleave_Sum(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1062: {
1063: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
1065: PetscFunctionBegin;
1066: sum->interleave_basis = interleave_basis;
1067: sum->interleave_components = interleave_components;
1068: PetscFunctionReturn(PETSC_SUCCESS);
1069: }
1071: /*@
1072: PetscDualSpaceSumGetInterleave - Get whether the basis functions and components of a uniform sum are interleaved
1074: Logically collective
1076: Input Parameter:
1077: . sp - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1079: Output Parameters:
1080: + interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1081: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1082: interleave the concatenated components
1084: Level: developer
1086: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumSetInterleave()`
1087: @*/
1088: PetscErrorCode PetscDualSpaceSumGetInterleave(PetscDualSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
1089: {
1090: PetscFunctionBegin;
1092: if (interleave_basis) PetscAssertPointer(interleave_basis, 2);
1093: if (interleave_components) PetscAssertPointer(interleave_components, 3);
1094: PetscTryMethod(sp, "PetscDualSpaceSumGetInterleave_C", (PetscDualSpace, PetscBool *, PetscBool *), (sp, interleave_basis, interleave_components));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: static PetscErrorCode PetscDualSpaceSumGetInterleave_Sum(PetscDualSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
1099: {
1100: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
1102: PetscFunctionBegin;
1103: if (interleave_basis) *interleave_basis = sum->interleave_basis;
1104: if (interleave_components) *interleave_components = sum->interleave_components;
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: #define PetscDualSpaceSumPassthrough(sp, func, ...) \
1109: do { \
1110: PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data; \
1111: PetscBool is_uniform; \
1112: PetscCall(PetscDualSpaceSumIsUniform(sp, &is_uniform)); \
1113: if (is_uniform && sum->numSumSpaces > 0) { \
1114: PetscDualSpace subsp; \
1115: PetscCall(PetscDualSpaceSumGetSubspace(sp, 0, &subsp)); \
1116: PetscCall(func(subsp, __VA_ARGS__)); \
1117: } \
1118: } while (0)
1120: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Sum(PetscDualSpace sp, PetscBool *value)
1121: {
1122: PetscFunctionBegin;
1123: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetContinuity, value);
1124: PetscFunctionReturn(PETSC_SUCCESS);
1125: }
1127: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Sum(PetscDualSpace sp, PetscBool value)
1128: {
1129: PetscFunctionBegin;
1130: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetContinuity, value);
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: }
1134: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Sum(PetscDualSpace sp, PetscBool *value)
1135: {
1136: PetscFunctionBegin;
1137: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTensor, value);
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Sum(PetscDualSpace sp, PetscBool value)
1142: {
1143: PetscFunctionBegin;
1144: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTensor, value);
1145: PetscFunctionReturn(PETSC_SUCCESS);
1146: }
1148: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Sum(PetscDualSpace sp, PetscBool *value)
1149: {
1150: PetscFunctionBegin;
1151: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTrimmed, value);
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Sum(PetscDualSpace sp, PetscBool value)
1156: {
1157: PetscFunctionBegin;
1158: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTrimmed, value);
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Sum(PetscDualSpace sp, PetscBool *value)
1163: {
1164: PetscFunctionBegin;
1165: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetUseMoments, value);
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Sum(PetscDualSpace sp, PetscBool value)
1170: {
1171: PetscFunctionBegin;
1172: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetUseMoments, value);
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Sum(PetscDualSpace sp, PetscInt *value)
1177: {
1178: PetscFunctionBegin;
1179: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetMomentOrder, value);
1180: PetscFunctionReturn(PETSC_SUCCESS);
1181: }
1183: static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Sum(PetscDualSpace sp, PetscInt value)
1184: {
1185: PetscFunctionBegin;
1186: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetMomentOrder, value);
1187: PetscFunctionReturn(PETSC_SUCCESS);
1188: }
1190: static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType *node_type, PetscBool *include_endpoints, PetscReal *exponent)
1191: {
1192: PetscFunctionBegin;
1193: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetNodeType, node_type, include_endpoints, exponent);
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType node_type, PetscBool include_endpoints, PetscReal exponent)
1198: {
1199: PetscFunctionBegin;
1200: PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetNodeType, node_type, include_endpoints, exponent);
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: static PetscErrorCode PetscDualSpaceInitialize_Sum(PetscDualSpace sp)
1205: {
1206: PetscFunctionBegin;
1207: sp->ops->destroy = PetscDualSpaceDestroy_Sum;
1208: sp->ops->view = PetscDualSpaceView_Sum;
1209: sp->ops->setfromoptions = NULL;
1210: sp->ops->duplicate = PetscDualSpaceDuplicate_Sum;
1211: sp->ops->setup = PetscDualSpaceSetUp_Sum;
1212: sp->ops->createheightsubspace = NULL;
1213: sp->ops->createpointsubspace = NULL;
1214: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Sum;
1215: sp->ops->apply = PetscDualSpaceApplyDefault;
1216: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
1217: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
1218: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
1219: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
1221: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", PetscDualSpaceSumGetNumSubspaces_Sum));
1222: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", PetscDualSpaceSumSetNumSubspaces_Sum));
1223: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", PetscDualSpaceSumGetSubspace_Sum));
1224: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", PetscDualSpaceSumSetSubspace_Sum));
1225: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", PetscDualSpaceSumGetConcatenate_Sum));
1226: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", PetscDualSpaceSumSetConcatenate_Sum));
1227: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", PetscDualSpaceSumGetInterleave_Sum));
1228: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", PetscDualSpaceSumSetInterleave_Sum));
1229: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Sum));
1230: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Sum));
1231: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Sum));
1232: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Sum));
1233: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Sum));
1234: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Sum));
1235: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Sum));
1236: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Sum));
1237: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Sum));
1238: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Sum));
1239: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Sum));
1240: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Sum));
1241: PetscFunctionReturn(PETSC_SUCCESS);
1242: }
1244: /*MC
1245: PETSCDUALSPACESUM = "sum" - A `PetscDualSpace` object that encapsulates a sum of subspaces.
1247: Level: intermediate
1249: Note:
1250: That sum can either be direct or a concatenation. For example if A and B are spaces each with 2 components,
1251: 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
1252: same reference element.
1254: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`, `PetscDualSpaceSumGetNumSubspaces()`, `PetscDualSpaceSumSetNumSubspaces()`,
1255: `PetscDualSpaceSumGetConcatenate()`, `PetscDualSpaceSumSetConcatenate()`, `PetscDualSpaceSumSetInterleave()`, `PetscDualSpaceSumGetInterleave()`
1256: M*/
1257: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Sum(PetscDualSpace sp)
1258: {
1259: PetscDualSpace_Sum *sum;
1261: PetscFunctionBegin;
1263: PetscCall(PetscNew(&sum));
1264: sum->numSumSpaces = PETSC_DEFAULT;
1265: sp->data = sum;
1266: sp->k = PETSC_FORM_DEGREE_UNDEFINED;
1267: PetscCall(PetscDualSpaceInitialize_Sum(sp));
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: }
1271: /*@
1272: PetscDualSpaceCreateSum - Create a finite element dual basis that is the sum of other dual bases
1274: Collective
1276: Input Parameters:
1277: + numSubspaces - the number of spaces that will be added together
1278: . subspaces - an array of length `numSubspaces` of spaces
1279: - 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
1281: Output Parameter:
1282: . sumSpace - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1284: Level: advanced
1286: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCSPACESUM`
1287: @*/
1288: PetscErrorCode PetscDualSpaceCreateSum(PetscInt numSubspaces, const PetscDualSpace subspaces[], PetscBool concatenate, PetscDualSpace *sumSpace)
1289: {
1290: PetscInt i, Nc = 0;
1292: PetscFunctionBegin;
1293: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace));
1294: PetscCall(PetscDualSpaceSetType(*sumSpace, PETSCDUALSPACESUM));
1295: PetscCall(PetscDualSpaceSumSetNumSubspaces(*sumSpace, numSubspaces));
1296: PetscCall(PetscDualSpaceSumSetConcatenate(*sumSpace, concatenate));
1297: for (i = 0; i < numSubspaces; ++i) {
1298: PetscInt sNc;
1300: PetscCall(PetscDualSpaceSumSetSubspace(*sumSpace, i, subspaces[i]));
1301: PetscCall(PetscDualSpaceGetNumComponents(subspaces[i], &sNc));
1302: if (concatenate) Nc += sNc;
1303: else Nc = sNc;
1305: if (i == 0) {
1306: DM dm;
1308: PetscCall(PetscDualSpaceGetDM(subspaces[i], &dm));
1309: PetscCall(PetscDualSpaceSetDM(*sumSpace, dm));
1310: }
1311: }
1312: PetscCall(PetscDualSpaceSetNumComponents(*sumSpace, Nc));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }