Actual source code: dualspacerefined.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: typedef struct {
5: PetscInt dummy;
6: } PetscDualSpace_Refined;
8: /*@
9: PetscDualSpaceRefinedSetCellSpaces - Set the dual spaces for the closures of each of the cells
10: in the multicell `DM` of a `PetscDualSpace`
12: Collective
14: Input Parameters:
15: + sp - a `PetscDualSpace`
16: - cellSpaces - one `PetscDualSpace` for each of the cells. The reference count of each cell space will be incremented,
17: so the user is still responsible for these spaces afterwards
19: Level: intermediate
21: .seealso: `PETSCDUALSPACEREFINED`, `PetscDualSpace`, `PetscFERefine()`
22: @*/
23: PetscErrorCode PetscDualSpaceRefinedSetCellSpaces(PetscDualSpace sp, const PetscDualSpace cellSpaces[])
24: {
25: PetscFunctionBegin;
27: PetscAssertPointer(cellSpaces, 2);
28: PetscCheck(!sp->setupcalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change cell spaces after setup is called");
29: PetscTryMethod(sp, "PetscDualSpaceRefinedSetCellSpaces_C", (PetscDualSpace, const PetscDualSpace[]), (sp, cellSpaces));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: static PetscErrorCode PetscDualSpaceRefinedSetCellSpaces_Refined(PetscDualSpace sp, const PetscDualSpace cellSpaces[])
34: {
35: DM dm;
36: PetscInt pStart, pEnd;
37: PetscInt cStart, cEnd, c;
39: PetscFunctionBegin;
40: dm = sp->dm;
41: PetscCheck(dm, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "PetscDualSpace must have a DM (PetscDualSpaceSetDM()) before calling PetscDualSpaceRefinedSetCellSpaces");
42: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
43: if (!sp->pointSpaces) PetscCall(PetscCalloc1(pEnd - pStart, &sp->pointSpaces));
44: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
45: for (c = 0; c < cEnd - cStart; c++) {
46: PetscCall(PetscObjectReference((PetscObject)cellSpaces[c]));
47: PetscCall(PetscDualSpaceDestroy(&sp->pointSpaces[c + cStart - pStart]));
48: sp->pointSpaces[c + cStart - pStart] = cellSpaces[c];
49: }
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PetscDualSpaceDestroy_Refined(PetscDualSpace sp)
54: {
55: PetscDualSpace_Refined *ref = (PetscDualSpace_Refined *)sp->data;
57: PetscFunctionBegin;
58: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceRefinedSetCellSpaces_C", NULL));
59: PetscCall(PetscFree(ref));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode PetscDualSpaceSetUp_Refined(PetscDualSpace sp)
64: {
65: PetscInt pStart, pEnd, depth;
66: PetscInt cStart, cEnd, c, spdim;
67: PetscInt h;
68: DM dm;
69: PetscSection section;
71: PetscFunctionBegin;
72: PetscCall(PetscDualSpaceGetDM(sp, &dm));
73: PetscCall(DMPlexGetDepth(dm, &depth));
74: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
75: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
76: for (c = cStart; c < cEnd; c++) {
77: if (sp->pointSpaces[c - pStart]) {
78: PetscInt ccStart, ccEnd;
79: PetscCheck(sp->pointSpaces[c - pStart]->k == sp->k, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have the same form degree as the refined dual space");
80: PetscCheck(sp->pointSpaces[c - pStart]->Nc == sp->Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have the same number of components as the refined dual space");
81: PetscCall(DMPlexGetHeightStratum(sp->pointSpaces[c - pStart]->dm, 0, &ccStart, &ccEnd));
82: PetscCheck(ccEnd - ccStart == 1, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have a single cell themselves");
83: }
84: }
85: for (c = cStart; c < cEnd; c++) {
86: if (sp->pointSpaces[c - pStart]) {
87: PetscBool cUniform;
89: PetscCall(PetscDualSpaceGetUniform(sp->pointSpaces[c - pStart], &cUniform));
90: if (!cUniform) break;
91: }
92: if ((c > cStart) && sp->pointSpaces[c - pStart] != sp->pointSpaces[c - 1 - pStart]) break;
93: }
94: if (c < cEnd) sp->uniform = PETSC_FALSE;
95: for (h = 0; h < depth; h++) {
96: PetscInt hStart, hEnd;
98: PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
99: for (c = hStart; c < hEnd; c++) {
100: PetscInt coneSize, e;
101: PetscDualSpace cspace = sp->pointSpaces[c - pStart];
102: const PetscInt *cone;
103: const PetscInt *refCone;
105: if (!cspace) continue;
106: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
107: PetscCall(DMPlexGetCone(dm, c, &cone));
108: PetscCall(DMPlexGetCone(cspace->dm, 0, &refCone));
109: for (e = 0; e < coneSize; e++) {
110: PetscInt point = cone[e];
111: PetscInt refpoint = refCone[e];
112: PetscDualSpace espace;
114: PetscCall(PetscDualSpaceGetPointSubspace(cspace, refpoint, &espace));
115: if (sp->pointSpaces[point - pStart] == NULL) {
116: PetscCall(PetscObjectReference((PetscObject)espace));
117: sp->pointSpaces[point - pStart] = espace;
118: }
119: }
120: }
121: }
122: PetscCall(PetscDualSpaceGetSection(sp, §ion));
123: PetscCall(PetscDualSpaceGetDimension(sp, &spdim));
124: PetscCall(PetscMalloc1(spdim, &sp->functional));
125: PetscCall(PetscDualSpacePushForwardSubspaces_Internal(sp, pStart, pEnd));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode PetscDualSpaceRefinedView_Ascii(PetscDualSpace sp, PetscViewer viewer)
130: {
131: PetscFunctionBegin;
132: if (sp->dm && sp->pointSpaces) {
133: PetscInt pStart, pEnd;
134: PetscInt cStart, cEnd, c;
136: PetscCall(DMPlexGetChart(sp->dm, &pStart, &pEnd));
137: PetscCall(DMPlexGetHeightStratum(sp->dm, 0, &cStart, &cEnd));
138: PetscCall(PetscViewerASCIIPrintf(viewer, "Refined dual space:\n"));
139: PetscCall(PetscViewerASCIIPushTab(viewer));
140: for (c = cStart; c < cEnd; c++) {
141: if (!sp->pointSpaces[c - pStart]) {
142: PetscCall(PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT " not set yet\n", c));
143: } else {
144: PetscCall(PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT ":ot set yet\n", c));
145: PetscCall(PetscViewerASCIIPushTab(viewer));
146: PetscCall(PetscDualSpaceView(sp->pointSpaces[c - pStart], viewer));
147: PetscCall(PetscViewerASCIIPopTab(viewer));
148: }
149: }
150: PetscCall(PetscViewerASCIIPopTab(viewer));
151: } else PetscCall(PetscViewerASCIIPrintf(viewer, "Refined dual space: (cell spaces not set yet)\n"));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode PetscDualSpaceView_Refined(PetscDualSpace sp, PetscViewer viewer)
156: {
157: PetscBool iascii;
159: PetscFunctionBegin;
162: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
163: if (iascii) PetscCall(PetscDualSpaceRefinedView_Ascii(sp, viewer));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode PetscDualSpaceInitialize_Refined(PetscDualSpace sp)
168: {
169: PetscFunctionBegin;
170: sp->ops->destroy = PetscDualSpaceDestroy_Refined;
171: sp->ops->view = PetscDualSpaceView_Refined;
172: sp->ops->setfromoptions = NULL;
173: sp->ops->duplicate = NULL;
174: sp->ops->setup = PetscDualSpaceSetUp_Refined;
175: sp->ops->createheightsubspace = NULL;
176: sp->ops->createpointsubspace = NULL;
177: sp->ops->getsymmetries = NULL;
178: sp->ops->apply = PetscDualSpaceApplyDefault;
179: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
180: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
181: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
182: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*MC
187: PETSCDUALSPACEREFINED = "refined" - A `PetscDualSpaceType` that defines the joint dual space of a group of cells, usually refined from one larger cell
189: Level: intermediate
191: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceRefinedSetCellSpaces`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
192: M*/
193: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Refined(PetscDualSpace sp)
194: {
195: PetscDualSpace_Refined *ref;
197: PetscFunctionBegin;
199: PetscCall(PetscNew(&ref));
200: sp->data = ref;
202: PetscCall(PetscDualSpaceInitialize_Refined(sp));
203: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceRefinedSetCellSpaces_C", PetscDualSpaceRefinedSetCellSpaces_Refined));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }