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: DM dm;
68: PetscSection section;
70: PetscFunctionBegin;
71: PetscCall(PetscDualSpaceGetDM(sp, &dm));
72: PetscCall(DMPlexGetDepth(dm, &depth));
73: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
74: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
75: for (c = cStart; c < cEnd; c++) {
76: if (sp->pointSpaces[c - pStart]) {
77: PetscInt ccStart, ccEnd;
78: 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");
79: 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");
80: PetscCall(DMPlexGetHeightStratum(sp->pointSpaces[c - pStart]->dm, 0, &ccStart, &ccEnd));
81: PetscCheck(ccEnd - ccStart == 1, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have a single cell themselves");
82: }
83: }
84: for (c = cStart; c < cEnd; c++) {
85: if (sp->pointSpaces[c - pStart]) {
86: PetscBool cUniform;
88: PetscCall(PetscDualSpaceGetUniform(sp->pointSpaces[c - pStart], &cUniform));
89: if (!cUniform) break;
90: }
91: if ((c > cStart) && sp->pointSpaces[c - pStart] != sp->pointSpaces[c - 1 - pStart]) break;
92: }
93: if (c < cEnd) sp->uniform = PETSC_FALSE;
94: for (PetscInt h = 0; h < depth; h++) {
95: PetscInt hStart, hEnd;
97: PetscCall(DMPlexGetHeightStratum(dm, h, &hStart, &hEnd));
98: for (c = hStart; c < hEnd; c++) {
99: PetscInt coneSize;
100: PetscDualSpace cspace = sp->pointSpaces[c - pStart];
101: const PetscInt *cone;
102: const PetscInt *refCone;
104: if (!cspace) continue;
105: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
106: PetscCall(DMPlexGetCone(dm, c, &cone));
107: PetscCall(DMPlexGetCone(cspace->dm, 0, &refCone));
108: for (PetscInt e = 0; e < coneSize; e++) {
109: PetscInt point = cone[e];
110: PetscInt refpoint = refCone[e];
111: PetscDualSpace espace;
113: PetscCall(PetscDualSpaceGetPointSubspace(cspace, refpoint, &espace));
114: if (sp->pointSpaces[point - pStart] == NULL) {
115: PetscCall(PetscObjectReference((PetscObject)espace));
116: sp->pointSpaces[point - pStart] = espace;
117: }
118: }
119: }
120: }
121: PetscCall(PetscDualSpaceGetSection(sp, §ion));
122: PetscCall(PetscDualSpaceGetDimension(sp, &spdim));
123: PetscCall(PetscMalloc1(spdim, &sp->functional));
124: PetscCall(PetscDualSpacePushForwardSubspaces_Internal(sp, pStart, pEnd));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: static PetscErrorCode PetscDualSpaceRefinedView_Ascii(PetscDualSpace sp, PetscViewer viewer)
129: {
130: PetscFunctionBegin;
131: if (sp->dm && sp->pointSpaces) {
132: PetscInt pStart, pEnd;
133: PetscInt cStart, cEnd, c;
135: PetscCall(DMPlexGetChart(sp->dm, &pStart, &pEnd));
136: PetscCall(DMPlexGetHeightStratum(sp->dm, 0, &cStart, &cEnd));
137: PetscCall(PetscViewerASCIIPrintf(viewer, "Refined dual space:\n"));
138: PetscCall(PetscViewerASCIIPushTab(viewer));
139: for (c = cStart; c < cEnd; c++) {
140: if (!sp->pointSpaces[c - pStart]) {
141: PetscCall(PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT " not set yet\n", c));
142: } else {
143: PetscCall(PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT ":ot set yet\n", c));
144: PetscCall(PetscViewerASCIIPushTab(viewer));
145: PetscCall(PetscDualSpaceView(sp->pointSpaces[c - pStart], viewer));
146: PetscCall(PetscViewerASCIIPopTab(viewer));
147: }
148: }
149: PetscCall(PetscViewerASCIIPopTab(viewer));
150: } else PetscCall(PetscViewerASCIIPrintf(viewer, "Refined dual space: (cell spaces not set yet)\n"));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode PetscDualSpaceView_Refined(PetscDualSpace sp, PetscViewer viewer)
155: {
156: PetscBool isascii;
158: PetscFunctionBegin;
161: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
162: if (isascii) PetscCall(PetscDualSpaceRefinedView_Ascii(sp, viewer));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode PetscDualSpaceInitialize_Refined(PetscDualSpace sp)
167: {
168: PetscFunctionBegin;
169: sp->ops->destroy = PetscDualSpaceDestroy_Refined;
170: sp->ops->view = PetscDualSpaceView_Refined;
171: sp->ops->setfromoptions = NULL;
172: sp->ops->duplicate = NULL;
173: sp->ops->setup = PetscDualSpaceSetUp_Refined;
174: sp->ops->createheightsubspace = NULL;
175: sp->ops->createpointsubspace = NULL;
176: sp->ops->getsymmetries = NULL;
177: sp->ops->apply = PetscDualSpaceApplyDefault;
178: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
179: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
180: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
181: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*MC
186: PETSCDUALSPACEREFINED = "refined" - A `PetscDualSpaceType` that defines the joint dual space of a group of cells, usually refined from one larger cell
188: Level: intermediate
190: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceRefinedSetCellSpaces`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
191: M*/
192: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Refined(PetscDualSpace sp)
193: {
194: PetscDualSpace_Refined *ref;
196: PetscFunctionBegin;
198: PetscCall(PetscNew(&ref));
199: sp->data = ref;
201: PetscCall(PetscDualSpaceInitialize_Refined(sp));
202: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceRefinedSetCellSpaces_C", PetscDualSpaceRefinedSetCellSpaces_Refined));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }