Actual source code: plexceed.c

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

  3: /*@C
  4:   DMPlexGetLocalOffsets - Allocate and populate array of local offsets.

  6:   Input Parameters:
  7: +  dm - The `DMPLEX` object
  8: .  domain_label - label for `DMPLEX` domain, or NULL for whole domain
  9: .  label_value - Stratum value
 10: .  height - Height of target cells in `DMPLEX` topology
 11: -  dm_field - Index of `DMPLEX` field

 13:   Output Parameters:
 14: +  num_cells - Number of local cells
 15: .  cell_size - Size of each cell, given by cell_size * num_comp = num_dof
 16: .  num_comp - Number of components per dof
 17: .  l_size - Size of local vector
 18: -  offsets - Allocated offsets array for cells

 20:   Level: developer

 22:   Notes:
 23:   Allocate and populate array of shape [num_cells, cell_size] defining offsets for each value (cell, node) for local vector of the `DMPLEX` field. All offsets are in the range [0, l_size - 1].

 25:    Caller is responsible for freeing the offsets array using `PetscFree()`.

 27: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
 28: @*/
 29: PetscErrorCode DMPlexGetLocalOffsets(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, PetscInt *num_cells, PetscInt *cell_size, PetscInt *num_comp, PetscInt *l_size, PetscInt **offsets)
 30: {
 31:   PetscDS         ds = NULL;
 32:   PetscFE         fe;
 33:   PetscSection    section;
 34:   PetscInt        dim, ds_field = -1;
 35:   PetscInt       *restr_indices;
 36:   const PetscInt *iter_indices;
 37:   IS              iter_is;

 41:   DMGetLocalSection(dm, &section);
 42:   DMGetDimension(dm, &dim);
 43:   {
 44:     IS              field_is;
 45:     const PetscInt *fields;
 46:     PetscInt        num_fields;

 48:     DMGetRegionDS(dm, domain_label, &field_is, &ds);
 49:     // Translate dm_field to ds_field
 50:     ISGetIndices(field_is, &fields);
 51:     ISGetSize(field_is, &num_fields);
 52:     for (PetscInt i = 0; i < num_fields; i++) {
 53:       if (dm_field == fields[i]) {
 54:         ds_field = i;
 55:         break;
 56:       }
 57:     }
 58:     ISRestoreIndices(field_is, &fields);
 59:   }

 62:   {
 63:     PetscInt depth;
 64:     DMLabel  depth_label;
 65:     IS       depth_is;

 67:     DMPlexGetDepth(dm, &depth);
 68:     DMPlexGetDepthLabel(dm, &depth_label);
 69:     DMLabelGetStratumIS(depth_label, depth - height, &depth_is);
 70:     if (domain_label) {
 71:       IS domain_is;

 73:       DMLabelGetStratumIS(domain_label, label_value, &domain_is);
 74:       if (domain_is) { // domainIS is non-empty
 75:         ISIntersect(depth_is, domain_is, &iter_is);
 76:         ISDestroy(&domain_is);
 77:       } else { // domainIS is NULL (empty)
 78:         iter_is = NULL;
 79:       }
 80:       ISDestroy(&depth_is);
 81:     } else {
 82:       iter_is = depth_is;
 83:     }
 84:     if (iter_is) {
 85:       ISGetLocalSize(iter_is, num_cells);
 86:       ISGetIndices(iter_is, &iter_indices);
 87:     } else {
 88:       *num_cells   = 0;
 89:       iter_indices = NULL;
 90:     }
 91:   }

 93:   {
 94:     PetscDualSpace dual_space;
 95:     PetscInt       num_dual_basis_vectors;

 97:     PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe);
 98:     PetscFEGetHeightSubspace(fe, height, &fe);
100:     PetscFEGetDualSpace(fe, &dual_space);
101:     PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors);
102:     PetscDualSpaceGetNumComponents(dual_space, num_comp);
104:     *cell_size = num_dual_basis_vectors / *num_comp;
105:   }
106:   PetscInt restr_size = (*num_cells) * (*cell_size);
107:   PetscMalloc1(restr_size, &restr_indices);
108:   PetscInt cell_offset = 0;

110:   PetscInt P = (PetscInt)pow(*cell_size, 1.0 / (dim - height));
111:   for (PetscInt p = 0; p < *num_cells; p++) {
112:     PetscBool flip = PETSC_FALSE;
113:     PetscInt  c    = iter_indices[p];
114:     PetscInt  num_indices, *indices;
115:     PetscInt  field_offsets[17]; // max number of fields plus 1
116:     DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL);
117:     if (height > 0) {
118:       PetscInt        num_cells_support, num_faces, start = -1;
119:       const PetscInt *orients, *faces, *cells;
120:       DMPlexGetSupport(dm, c, &cells);
121:       DMPlexGetSupportSize(dm, c, &num_cells_support);
123:       DMPlexGetCone(dm, cells[0], &faces);
124:       DMPlexGetConeSize(dm, cells[0], &num_faces);
125:       for (PetscInt i = 0; i < num_faces; i++) {
126:         if (faces[i] == c) start = i;
127:       }
129:       DMPlexGetConeOrientation(dm, cells[0], &orients);
130:       if (orients[start] < 0) flip = PETSC_TRUE;
131:     }

133:     for (PetscInt i = 0; i < *cell_size; i++) {
134:       PetscInt ii = i;
135:       if (flip) {
136:         if (*cell_size == P) ii = *cell_size - 1 - i;
137:         else if (*cell_size == P * P) {
138:           PetscInt row = i / P, col = i % P;
139:           ii = row + col * P;
140:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for flipping point with cell size %" PetscInt_FMT " != P (%" PetscInt_FMT ") or P^2", *cell_size, P);
141:       }
142:       // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
143:       PetscInt loc                 = indices[field_offsets[dm_field] + ii * (*num_comp)];
144:       restr_indices[cell_offset++] = loc >= 0 ? loc : -(loc + 1);
145:     }
146:     DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL);
147:   }
149:   if (iter_is) ISRestoreIndices(iter_is, &iter_indices);
150:   ISDestroy(&iter_is);

152:   *offsets = restr_indices;
153:   PetscSectionGetStorageSize(section, l_size);
154:   return 0;
155: }

157: #if defined(PETSC_HAVE_LIBCEED)
158: #include <petscdmplexceed.h>

160: /*@C
161:   DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector)

163:   Input Parameters:
164: +  dm - The `DMPLEX` object
165: .  domain_label - label for `DMPLEX` domain, or NULL for the whole domain
166: .  label_value - Stratum value
167: .  height - Height of target cells in `DMPLEX` topology
168: -  dm_field - Index of `DMPLEX` field

170:   Output Parameter:
171: .  ERestrict - libCEED restriction from local vector to to the cells

173:   Level: developer

175: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetLocalOffsets()`
176: @*/
177: PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
178: {
182:   if (!dm->ceedERestrict) {
183:     PetscInt            num_cells, cell_size, num_comp, lvec_size, *restr_indices;
184:     CeedElemRestriction elem_restr;
185:     Ceed                ceed;

187:     DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices);
188:     DMGetCeed(dm, &ceed);
189:     CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices, &elem_restr);
190:     PetscFree(restr_indices);
191:     dm->ceedERestrict = elem_restr;
192:   }
193:   *ERestrict = dm->ceedERestrict;
194:   return 0;
195: }

197: #endif