Actual source code: dualspacerefined.c

petsc-3.14.1 2020-11-03
Report Typos and Errors
  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, &section);
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: }