Actual source code: plexref1d.c

  1: #include <petsc/private/dmplextransformimpl.h>

  3: static PetscErrorCode DMPlexTransformSetUp_1D(DMPlexTransform tr)
  4: {
  5:   DM       dm;
  6:   DMLabel  active;
  7:   PetscInt pStart, pEnd, p;

  9:   DMPlexTransformGetDM(tr, &dm);
 10:   DMPlexTransformGetActive(tr, &active);
 12:   /* Calculate refineType for each cell */
 13:   DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType);
 14:   DMPlexGetChart(dm, &pStart, &pEnd);
 15:   for (p = pStart; p < pEnd; ++p) {
 16:     DMLabel        trType = tr->trType;
 17:     DMPolytopeType ct;
 18:     PetscInt       val;

 20:     DMPlexGetCellType(dm, p, &ct);
 21:     switch (ct) {
 22:     case DM_POLYTOPE_POINT:
 23:       DMLabelSetValue(trType, p, 0);
 24:       break;
 25:     case DM_POLYTOPE_SEGMENT:
 26:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
 27:       DMLabelGetValue(active, p, &val);
 28:       if (val == 1) DMLabelSetValue(trType, p, val);
 29:       else DMLabelSetValue(trType, p, 2);
 30:       break;
 31:     default:
 32:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
 33:     }
 34:   }
 35:   return 0;
 36: }

 38: static PetscErrorCode DMPlexTransformGetSubcellOrientation_1D(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
 39: {
 40:   PetscInt rt;

 43:   DMLabelGetValue(tr->trType, sp, &rt);
 44:   *rnew = r;
 45:   *onew = o;
 46:   switch (rt) {
 47:   case 1:
 48:     DMPlexTransformGetSubcellOrientation_Regular(tr, sct, sp, so, tct, r, o, rnew, onew);
 49:     break;
 50:   default:
 51:     DMPlexTransformGetSubcellOrientationIdentity(tr, sct, sp, so, tct, r, o, rnew, onew);
 52:   }
 53:   return 0;
 54: }

 56: static PetscErrorCode DMPlexTransformCellTransform_1D(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
 57: {
 58:   DMLabel  trType = tr->trType;
 59:   PetscInt val;

 63:   DMLabelGetValue(trType, p, &val);
 64:   if (rt) *rt = val;
 65:   switch (source) {
 66:   case DM_POLYTOPE_POINT:
 67:     DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt);
 68:     break;
 69:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
 70:   case DM_POLYTOPE_SEGMENT:
 71:     if (val == 1) DMPlexTransformCellRefine_Regular(tr, source, p, NULL, Nt, target, size, cone, ornt);
 72:     else DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt);
 73:     break;
 74:   default:
 75:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
 76:   }
 77:   return 0;
 78: }

 80: static PetscErrorCode DMPlexTransformSetFromOptions_1D(DMPlexTransform tr, PetscOptionItems *PetscOptionsObject)
 81: {
 82:   PetscInt  cells[256], n = 256, i;
 83:   PetscBool flg;

 85:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Options");
 86:   PetscOptionsIntArray("-dm_plex_transform_1d_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
 87:   if (flg) {
 88:     DMLabel active;

 90:     DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
 91:     for (i = 0; i < n; ++i) DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);
 92:     DMPlexTransformSetActive(tr, active);
 93:     DMLabelDestroy(&active);
 94:   }
 95:   PetscOptionsHeadEnd();
 96:   return 0;
 97: }

 99: static PetscErrorCode DMPlexTransformView_1D(DMPlexTransform tr, PetscViewer viewer)
100: {
101:   PetscBool isascii;

105:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
106:   if (isascii) {
107:     PetscViewerFormat format;
108:     const char       *name;

110:     PetscObjectGetName((PetscObject)tr, &name);
111:     PetscViewerASCIIPrintf(viewer, "1D refinement %s\n", name ? name : "");
112:     PetscViewerGetFormat(viewer, &format);
113:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) DMLabelView(tr->trType, viewer);
114:   } else {
115:     SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name);
116:   }
117:   return 0;
118: }

120: static PetscErrorCode DMPlexTransformDestroy_1D(DMPlexTransform tr)
121: {
122:   PetscFree(tr->data);
123:   return 0;
124: }

126: static PetscErrorCode DMPlexTransformInitialize_1D(DMPlexTransform tr)
127: {
128:   tr->ops->view                  = DMPlexTransformView_1D;
129:   tr->ops->setfromoptions        = DMPlexTransformSetFromOptions_1D;
130:   tr->ops->setup                 = DMPlexTransformSetUp_1D;
131:   tr->ops->destroy               = DMPlexTransformDestroy_1D;
132:   tr->ops->celltransform         = DMPlexTransformCellTransform_1D;
133:   tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_1D;
134:   tr->ops->mapcoordinates        = DMPlexTransformMapCoordinatesBarycenter_Internal;
135:   return 0;
136: }

138: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform tr)
139: {
140:   DMPlexRefine_1D *f;

143:   PetscNew(&f);
144:   tr->data = f;

146:   DMPlexTransformInitialize_1D(tr);
147:   return 0;
148: }