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 on sp
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: `PetscDualSpace`, `PetscFERefine()`
22: @*/
23: PetscErrorCode PetscDualSpaceRefinedSetCellSpaces(PetscDualSpace sp, const PetscDualSpace cellSpaces[])
24: {
28: PetscTryMethod(sp, "PetscDualSpaceRefinedSetCellSpaces_C", (PetscDualSpace, const PetscDualSpace[]), (sp, cellSpaces));
29: return 0;
30: }
32: static PetscErrorCode PetscDualSpaceRefinedSetCellSpaces_Refined(PetscDualSpace sp, const PetscDualSpace cellSpaces[])
33: {
34: DM dm;
35: PetscInt pStart, pEnd;
36: PetscInt cStart, cEnd, c;
38: dm = sp->dm;
40: DMPlexGetChart(dm, &pStart, &pEnd);
41: if (!sp->pointSpaces) PetscCalloc1(pEnd - pStart, &(sp->pointSpaces));
42: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
43: for (c = 0; c < cEnd - cStart; c++) {
44: PetscObjectReference((PetscObject)cellSpaces[c]);
45: PetscDualSpaceDestroy(&(sp->pointSpaces[c + cStart - pStart]));
46: sp->pointSpaces[c + cStart - pStart] = cellSpaces[c];
47: }
48: return 0;
49: }
51: static PetscErrorCode PetscDualSpaceDestroy_Refined(PetscDualSpace sp)
52: {
53: PetscDualSpace_Refined *ref = (PetscDualSpace_Refined *)sp->data;
55: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceRefinedSetCellSpaces_C", NULL);
56: PetscFree(ref);
57: return 0;
58: }
60: static PetscErrorCode PetscDualSpaceSetUp_Refined(PetscDualSpace sp)
61: {
62: PetscInt pStart, pEnd, depth;
63: PetscInt cStart, cEnd, c, spdim;
64: PetscInt h;
65: DM dm;
66: PetscSection section;
68: PetscDualSpaceGetDM(sp, &dm);
69: DMPlexGetDepth(dm, &depth);
70: DMPlexGetChart(dm, &pStart, &pEnd);
71: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
72: for (c = cStart; c < cEnd; c++) {
73: if (sp->pointSpaces[c - pStart]) {
74: PetscInt ccStart, ccEnd;
77: DMPlexGetHeightStratum(sp->pointSpaces[c - pStart]->dm, 0, &ccStart, &ccEnd);
79: }
80: }
81: for (c = cStart; c < cEnd; c++) {
82: if (sp->pointSpaces[c - pStart]) {
83: PetscBool cUniform;
85: PetscDualSpaceGetUniform(sp->pointSpaces[c - pStart], &cUniform);
86: if (!cUniform) break;
87: }
88: if ((c > cStart) && sp->pointSpaces[c - pStart] != sp->pointSpaces[c - 1 - pStart]) break;
89: }
90: if (c < cEnd) sp->uniform = PETSC_FALSE;
91: for (h = 0; h < depth; h++) {
92: PetscInt hStart, hEnd;
94: DMPlexGetHeightStratum(dm, h, &hStart, &hEnd);
95: for (c = hStart; c < hEnd; c++) {
96: PetscInt coneSize, e;
97: PetscDualSpace cspace = sp->pointSpaces[c - pStart];
98: const PetscInt *cone;
99: const PetscInt *refCone;
101: if (!cspace) continue;
102: DMPlexGetConeSize(dm, c, &coneSize);
103: DMPlexGetCone(dm, c, &cone);
104: DMPlexGetCone(cspace->dm, 0, &refCone);
105: for (e = 0; e < coneSize; e++) {
106: PetscInt point = cone[e];
107: PetscInt refpoint = refCone[e];
108: PetscDualSpace espace;
110: PetscDualSpaceGetPointSubspace(cspace, refpoint, &espace);
111: if (sp->pointSpaces[point - pStart] == NULL) {
112: PetscObjectReference((PetscObject)espace);
113: sp->pointSpaces[point - pStart] = espace;
114: }
115: }
116: }
117: }
118: PetscDualSpaceGetSection(sp, §ion);
119: PetscDualSpaceGetDimension(sp, &spdim);
120: PetscMalloc1(spdim, &(sp->functional));
121: PetscDualSpacePushForwardSubspaces_Internal(sp, pStart, pEnd);
122: return 0;
123: }
125: static PetscErrorCode PetscDualSpaceRefinedView_Ascii(PetscDualSpace sp, PetscViewer viewer)
126: {
127: if (sp->dm && sp->pointSpaces) {
128: PetscInt pStart, pEnd;
129: PetscInt cStart, cEnd, c;
131: DMPlexGetChart(sp->dm, &pStart, &pEnd);
132: DMPlexGetHeightStratum(sp->dm, 0, &cStart, &cEnd);
133: PetscViewerASCIIPrintf(viewer, "Refined dual space:\n");
134: PetscViewerASCIIPushTab(viewer);
135: for (c = cStart; c < cEnd; c++) {
136: if (!sp->pointSpaces[c - pStart]) {
137: PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT " not set yet\n", c);
138: } else {
139: PetscViewerASCIIPrintf(viewer, "Cell space %" PetscInt_FMT ":ot set yet\n", c);
140: PetscViewerASCIIPushTab(viewer);
141: PetscDualSpaceView(sp->pointSpaces[c - pStart], viewer);
142: PetscViewerASCIIPopTab(viewer);
143: }
144: }
145: PetscViewerASCIIPopTab(viewer);
146: } else PetscViewerASCIIPrintf(viewer, "Refined dual space: (cell spaces not set yet)\n");
147: return 0;
148: }
150: static PetscErrorCode PetscDualSpaceView_Refined(PetscDualSpace sp, PetscViewer viewer)
151: {
152: PetscBool iascii;
156: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
157: if (iascii) PetscDualSpaceRefinedView_Ascii(sp, viewer);
158: return 0;
159: }
161: static PetscErrorCode PetscDualSpaceInitialize_Refined(PetscDualSpace sp)
162: {
163: sp->ops->destroy = PetscDualSpaceDestroy_Refined;
164: sp->ops->view = PetscDualSpaceView_Refined;
165: sp->ops->setfromoptions = NULL;
166: sp->ops->duplicate = NULL;
167: sp->ops->setup = PetscDualSpaceSetUp_Refined;
168: sp->ops->createheightsubspace = NULL;
169: sp->ops->createpointsubspace = NULL;
170: sp->ops->getsymmetries = NULL;
171: sp->ops->apply = PetscDualSpaceApplyDefault;
172: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
173: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
174: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
175: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
176: return 0;
177: }
179: /*MC
180: PETSCDUALSPACEREFINED = "refined" - A `PetscDualSpaceType` that defines the joint dual space of a group of cells, usually refined from one larger cell
182: Level: intermediate
184: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
185: M*/
186: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Refined(PetscDualSpace sp)
187: {
188: PetscDualSpace_Refined *ref;
191: PetscNew(&ref);
192: sp->data = ref;
194: PetscDualSpaceInitialize_Refined(sp);
195: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceRefinedSetCellSpaces_C", PetscDualSpaceRefinedSetCellSpaces_Refined);
196: return 0;
197: }