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 PetscDualSpace
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: 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) {
43: PetscCalloc1(pEnd-pStart,&(sp->pointSpaces));
44: }
45: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
46: for (c = 0; c < cEnd - cStart; c++) {
47: PetscObjectReference((PetscObject)cellSpaces[c]);
48: PetscDualSpaceDestroy(&(sp->pointSpaces[c + cStart - pStart]));
49: sp->pointSpaces[c+cStart-pStart] = cellSpaces[c];
50: }
51: return 0;
52: }
54: static PetscErrorCode PetscDualSpaceDestroy_Refined(PetscDualSpace sp)
55: {
56: PetscDualSpace_Refined *ref = (PetscDualSpace_Refined *) sp->data;
58: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceRefinedSetCellSpaces_C", NULL);
59: PetscFree(ref);
60: return 0;
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: PetscDualSpaceGetDM(sp, &dm);
72: DMPlexGetDepth(dm, &depth);
73: DMPlexGetChart(dm, &pStart, &pEnd);
74: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
75: for (c = cStart; c < cEnd; c++) {
76: if (sp->pointSpaces[c-pStart]) {
77: PetscInt ccStart, ccEnd;
80: DMPlexGetHeightStratum(sp->pointSpaces[c-pStart]->dm, 0, &ccStart, &ccEnd);
82: }
83: }
84: for (c = cStart; c < cEnd; c++) {
85: if (sp->pointSpaces[c-pStart]) {
86: PetscBool cUniform;
88: 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 (h = 0; h < depth; h++) {
95: PetscInt hStart, hEnd;
97: DMPlexGetHeightStratum(dm, h, &hStart, &hEnd);
98: for (c = hStart; c < hEnd; c++) {
99: PetscInt coneSize, e;
100: PetscDualSpace cspace = sp->pointSpaces[c-pStart];
101: const PetscInt *cone;
102: const PetscInt *refCone;
104: if (!cspace) continue;
105: DMPlexGetConeSize(dm, c, &coneSize);
106: DMPlexGetCone(dm, c, &cone);
107: DMPlexGetCone(cspace->dm, 0, &refCone);
108: for (e = 0; e < coneSize; e++) {
109: PetscInt point = cone[e];
110: PetscInt refpoint = refCone[e];
111: PetscDualSpace espace;
113: PetscDualSpaceGetPointSubspace(cspace,refpoint,&espace);
114: if (sp->pointSpaces[point-pStart] == NULL) {
115: PetscObjectReference((PetscObject)espace);
116: sp->pointSpaces[point-pStart] = espace;
117: }
118: }
119: }
120: }
121: PetscDualSpaceGetSection(sp, §ion);
122: PetscDualSpaceGetDimension(sp, &spdim);
123: PetscMalloc1(spdim, &(sp->functional));
124: PetscDualSpacePushForwardSubspaces_Internal(sp, pStart, pEnd);
125: return 0;
126: }
128: static PetscErrorCode PetscDualSpaceRefinedView_Ascii(PetscDualSpace sp, PetscViewer viewer)
129: {
130: if (sp->dm && sp->pointSpaces) {
131: PetscInt pStart, pEnd;
132: PetscInt cStart, cEnd, c;
134: DMPlexGetChart(sp->dm, &pStart, &pEnd);
135: DMPlexGetHeightStratum(sp->dm, 0, &cStart, &cEnd);
136: PetscViewerASCIIPrintf(viewer, "Refined dual space:\n");
137: PetscViewerASCIIPushTab(viewer);
138: for (c = cStart; c < cEnd; c++) {
139: if (!sp->pointSpaces[c-pStart]) {
140: PetscViewerASCIIPrintf(viewer, "Cell space %D not set yet\n", c);
141: } else {
142: PetscViewerASCIIPrintf(viewer, "Cell space %D:ot set yet\n", c);
143: PetscViewerASCIIPushTab(viewer);
144: PetscDualSpaceView(sp->pointSpaces[c-pStart],viewer);
145: PetscViewerASCIIPopTab(viewer);
146: }
147: }
148: PetscViewerASCIIPopTab(viewer);
149: } else {
150: PetscViewerASCIIPrintf(viewer, "Refined dual space: (cell spaces not set yet)\n");
151: }
152: return 0;
153: }
155: static PetscErrorCode PetscDualSpaceView_Refined(PetscDualSpace sp, PetscViewer viewer)
156: {
157: PetscBool iascii;
161: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
162: if (iascii) PetscDualSpaceRefinedView_Ascii(sp, viewer);
163: return 0;
164: }
166: static PetscErrorCode PetscDualSpaceInitialize_Refined(PetscDualSpace sp)
167: {
168: sp->ops->destroy = PetscDualSpaceDestroy_Refined;
169: sp->ops->view = PetscDualSpaceView_Refined;
170: sp->ops->setfromoptions = NULL;
171: sp->ops->duplicate = NULL;
172: sp->ops->setup = PetscDualSpaceSetUp_Refined;
173: sp->ops->createheightsubspace = NULL;
174: sp->ops->createpointsubspace = NULL;
175: sp->ops->getsymmetries = NULL;
176: sp->ops->apply = PetscDualSpaceApplyDefault;
177: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
178: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
179: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
180: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
181: return 0;
182: }
184: /*MC
185: PETSCDUALSPACEREFINED = "refined" - A PetscDualSpace object that defines the joint dual space of a group of cells, usually refined from one larger cell
187: Level: intermediate
189: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
190: M*/
191: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Refined(PetscDualSpace sp)
192: {
193: PetscDualSpace_Refined *ref;
196: PetscNewLog(sp,&ref);
197: sp->data = ref;
199: PetscDualSpaceInitialize_Refined(sp);
200: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceRefinedSetCellSpaces_C", PetscDualSpaceRefinedSetCellSpaces_Refined);
201: return 0;
202: }