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: DMLabelSetValue(trType, p, 0); break;
 23:       case DM_POLYTOPE_SEGMENT:
 24:       case DM_POLYTOPE_POINT_PRISM_TENSOR:
 25:         DMLabelGetValue(active, p, &val);
 26:         if (val == 1) DMLabelSetValue(trType, p, val);
 27:         else          DMLabelSetValue(trType, p, 2);
 28:         break;
 29:       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle points of type %s", DMPolytopeTypes[ct]);
 30:     }
 31:   }
 32:   return 0;
 33: }

 35: static PetscErrorCode DMPlexTransformGetSubcellOrientation_1D(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
 36: {
 37:   PetscInt       rt;

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

 51: static PetscErrorCode DMPlexTransformCellTransform_1D(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
 52: {
 53:   DMLabel        trType = tr->trType;
 54:   PetscInt       val;

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

 74: static PetscErrorCode DMPlexTransformSetFromOptions_1D(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr)
 75: {
 76:   PetscInt       cells[256], n = 256, i;
 77:   PetscBool      flg;

 80:   PetscOptionsHead(PetscOptionsObject,"DMPlex Options");
 81:   PetscOptionsIntArray("-dm_plex_transform_1d_ref_cell", "Mark cells for refinement", "", cells, &n, &flg);
 82:   if (flg) {
 83:     DMLabel active;

 85:     DMLabelCreate(PETSC_COMM_SELF, "Adaptation Label", &active);
 86:     for (i = 0; i < n; ++i) DMLabelSetValue(active, cells[i], DM_ADAPT_REFINE);
 87:     DMPlexTransformSetActive(tr, active);
 88:     DMLabelDestroy(&active);
 89:   }
 90:   PetscOptionsTail();
 91:   return 0;
 92: }

 94: static PetscErrorCode DMPlexTransformView_1D(DMPlexTransform tr, PetscViewer viewer)
 95: {
 96:   PetscBool      isascii;

100:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
101:   if (isascii) {
102:     PetscViewerFormat format;
103:     const char       *name;

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

117: static PetscErrorCode DMPlexTransformDestroy_1D(DMPlexTransform tr)
118: {
119:   PetscFree(tr->data);
120:   return 0;
121: }

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

135: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform tr)
136: {
137:   DMPlexRefine_1D *f;

140:   PetscNewLog(tr, &f);
141:   tr->data = f;

143:   DMPlexTransformInitialize_1D(tr);
144:   return 0;
145: }