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, &section);
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: }