Actual source code: dualspacerefined.c
petsc-3.14.6 2021-03-30
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 Arguments:
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: {
30: if (sp->setupcalled) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change cell spaces after setup is called");
31: PetscTryMethod(sp, "PetscDualSpaceRefinedSetCellSpaces_C", (PetscDualSpace,const PetscDualSpace []),(sp,cellSpaces));
32: return(0);
33: }
35: static PetscErrorCode PetscDualSpaceRefinedSetCellSpaces_Refined(PetscDualSpace sp, const PetscDualSpace cellSpaces[])
36: {
37: DM dm;
38: PetscInt pStart, pEnd;
39: PetscInt cStart, cEnd, c;
43: dm = sp->dm;
44: if (!dm) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_WRONGSTATE, "PetscDualSpace must have a DM (PetscDualSpaceSetDM()) before calling PetscDualSpaceRefinedSetCellSpaces");
45: DMPlexGetChart(dm, &pStart, &pEnd);
46: if (!sp->pointSpaces) {
48: PetscCalloc1(pEnd-pStart,&(sp->pointSpaces));
49: }
50: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
51: for (c = 0; c < cEnd - cStart; c++) {
52: PetscObjectReference((PetscObject)cellSpaces[c]);
53: PetscDualSpaceDestroy(&(sp->pointSpaces[c + cStart - pStart]));
54: sp->pointSpaces[c+cStart-pStart] = cellSpaces[c];
55: }
56: return(0);
57: }
59: static PetscErrorCode PetscDualSpaceDestroy_Refined(PetscDualSpace sp)
60: {
61: PetscDualSpace_Refined *ref = (PetscDualSpace_Refined *) sp->data;
62: PetscErrorCode ierr;
65: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceRefinedSetCellSpaces_C", NULL);
66: PetscFree(ref);
67: return(0);
68: }
70: static PetscErrorCode PetscDualSpaceSetUp_Refined(PetscDualSpace sp)
71: {
72: PetscInt pStart, pEnd, depth;
73: PetscInt cStart, cEnd, c, spdim;
74: PetscInt h;
75: DM dm;
76: PetscSection section;
80: PetscDualSpaceGetDM(sp, &dm);
81: DMPlexGetDepth(dm, &depth);
82: DMPlexGetChart(dm, &pStart, &pEnd);
83: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
84: for (c = cStart; c < cEnd; c++) {
85: if (sp->pointSpaces[c-pStart]) {
86: PetscInt ccStart, ccEnd;
87: if (sp->pointSpaces[c-pStart]->k != sp->k) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have the same form degree as the refined dual space");
88: if (sp->pointSpaces[c-pStart]->Nc != sp->Nc) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have the same number of components as the refined dual space");
89: DMPlexGetHeightStratum(sp->pointSpaces[c-pStart]->dm, 0, &ccStart, &ccEnd);
90: if (ccEnd - ccStart != 1) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "All cell spaces must have a single cell themselves");
91: }
92: }
93: for (c = cStart; c < cEnd; c++) {
94: if (sp->pointSpaces[c-pStart]) {
95: PetscBool cUniform;
97: PetscDualSpaceGetUniform(sp->pointSpaces[c-pStart],&cUniform);
98: if (!cUniform) break;
99: }
100: if ((c > cStart) && sp->pointSpaces[c-pStart] != sp->pointSpaces[c-1-pStart]) break;
101: }
102: if (c < cEnd) sp->uniform = PETSC_FALSE;
103: for (h = 0; h < depth; h++) {
104: PetscInt hStart, hEnd;
106: DMPlexGetHeightStratum(dm, h, &hStart, &hEnd);
107: for (c = hStart; c < hEnd; c++) {
108: PetscInt coneSize, e;
109: PetscDualSpace cspace = sp->pointSpaces[c-pStart];
110: const PetscInt *cone;
111: const PetscInt *refCone;
113: if (!cspace) continue;
114: DMPlexGetConeSize(dm, c, &coneSize);
115: DMPlexGetCone(dm, c, &cone);
116: DMPlexGetCone(cspace->dm, 0, &refCone);
117: for (e = 0; e < coneSize; e++) {
118: PetscInt point = cone[e];
119: PetscInt refpoint = refCone[e];
120: PetscDualSpace espace;
122: PetscDualSpaceGetPointSubspace(cspace,refpoint,&espace);
123: if (sp->pointSpaces[point-pStart] == NULL) {
124: PetscObjectReference((PetscObject)espace);
125: sp->pointSpaces[point-pStart] = espace;
126: }
127: }
128: }
129: }
130: PetscDualSpaceGetSection(sp, §ion);
131: PetscDualSpaceGetDimension(sp, &spdim);
132: PetscMalloc1(spdim, &(sp->functional));
133: PetscDualSpacePushForwardSubspaces_Internal(sp, pStart, pEnd);
134: return(0);
135: }
137: static PetscErrorCode PetscDualSpaceRefinedView_Ascii(PetscDualSpace sp, PetscViewer viewer)
138: {
139: PetscErrorCode ierr;
142: if (sp->dm && sp->pointSpaces) {
143: PetscInt pStart, pEnd;
144: PetscInt cStart, cEnd, c;
146: DMPlexGetChart(sp->dm, &pStart, &pEnd);
147: DMPlexGetHeightStratum(sp->dm, 0, &cStart, &cEnd);
148: PetscViewerASCIIPrintf(viewer, "Refined dual space:\n");
149: PetscViewerASCIIPushTab(viewer);
150: for (c = cStart; c < cEnd; c++) {
151: if (!sp->pointSpaces[c-pStart]) {
152: PetscViewerASCIIPrintf(viewer, "Cell space %D not set yet\n", c);
153: } else {
154: PetscViewerASCIIPrintf(viewer, "Cell space %D:ot set yet\n", c);
155: PetscViewerASCIIPushTab(viewer);
156: PetscDualSpaceView(sp->pointSpaces[c-pStart],viewer);
157: PetscViewerASCIIPopTab(viewer);
158: }
159: }
160: PetscViewerASCIIPopTab(viewer);
161: } else {
162: PetscViewerASCIIPrintf(viewer, "Refined dual space: (cell spaces not set yet)\n");
163: }
164: return(0);
165: }
167: static PetscErrorCode PetscDualSpaceView_Refined(PetscDualSpace sp, PetscViewer viewer)
168: {
169: PetscBool iascii;
175: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
176: if (iascii) {PetscDualSpaceRefinedView_Ascii(sp, viewer);}
177: return(0);
178: }
180: static PetscErrorCode PetscDualSpaceInitialize_Refined(PetscDualSpace sp)
181: {
183: sp->ops->destroy = PetscDualSpaceDestroy_Refined;
184: sp->ops->view = PetscDualSpaceView_Refined;
185: sp->ops->setfromoptions = NULL;
186: sp->ops->duplicate = NULL;
187: sp->ops->setup = PetscDualSpaceSetUp_Refined;
188: sp->ops->createheightsubspace = NULL;
189: sp->ops->createpointsubspace = NULL;
190: sp->ops->getsymmetries = NULL;
191: sp->ops->apply = PetscDualSpaceApplyDefault;
192: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
193: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
194: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
195: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
196: return(0);
197: }
199: /*MC
200: PETSCDUALSPACEREFINED = "refined" - A PetscDualSpace object that defines the joint dual space of a group of cells, usually refined from one larger cell
202: Level: intermediate
204: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
205: M*/
206: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Refined(PetscDualSpace sp)
207: {
208: PetscDualSpace_Refined *ref;
209: PetscErrorCode ierr;
213: PetscNewLog(sp,&ref);
214: sp->data = ref;
216: PetscDualSpaceInitialize_Refined(sp);
217: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceRefinedSetCellSpaces_C", PetscDualSpaceRefinedSetCellSpaces_Refined);
218: return(0);
219: }