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:   Notes: 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]. Caller is responsible for freeing the offsets array using PetscFree().

 22:   Level: developer

 24: .seealso: `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
 25: @*/
 26: 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)
 27: {
 28:   PetscDS         ds = NULL;
 29:   PetscFE         fe;
 30:   PetscSection    section;
 31:   PetscInt        dim, ds_field = -1;
 32:   PetscInt       *restr_indices;
 33:   const PetscInt *iter_indices;
 34:   IS              iter_is;

 38:   DMGetLocalSection(dm, &section);
 39:   DMGetDimension(dm, &dim);
 40:   {
 41:     IS              field_is;
 42:     const PetscInt *fields;
 43:     PetscInt        num_fields;

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

 59:   {
 60:     PetscInt depth;
 61:     DMLabel  depth_label;
 62:     IS       depth_is;

 64:     DMPlexGetDepth(dm, &depth);
 65:     DMPlexGetDepthLabel(dm, &depth_label);
 66:     DMLabelGetStratumIS(depth_label, depth - height, &depth_is);
 67:     if (domain_label) {
 68:       IS domain_is;

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

 90:   {
 91:     PetscDualSpace dual_space;
 92:     PetscInt       num_dual_basis_vectors;

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

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

130:     for (PetscInt i = 0; i < *cell_size; i++) {
131:       PetscInt ii = i;
132:       if (flip) {
133:         if (*cell_size == P) ii = *cell_size - 1 - i;
134:         else if (*cell_size == P * P) {
135:           PetscInt row = i / P, col = i % P;
136:           ii = row + col * P;
137:         } 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);
138:       }
139:       // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
140:       PetscInt loc                 = indices[field_offsets[dm_field] + ii * (*num_comp)];
141:       restr_indices[cell_offset++] = loc >= 0 ? loc : -(loc + 1);
142:     }
143:     DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL);
144:   }
146:   if (iter_is) ISRestoreIndices(iter_is, &iter_indices);
147:   ISDestroy(&iter_is);

149:   *offsets = restr_indices;
150:   PetscSectionGetStorageSize(section, l_size);
151:   return 0;
152: }

154: #if defined(PETSC_HAVE_LIBCEED)
155: #include <petscdmplexceed.h>

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

160:   Input Parameters:
161:   dm - The DMPlex object
162:   domain_label - label for DMPlex domain, or NULL for the whole domain
163:   label_value - Stratum value
164:   height - Height of target cells in DMPlex topology
165:   dm_field - Index of DMPlex field

167:   Output Parameters:
168:   ERestrict - libCEED restriction from local vector to to the cells

170:   Level: developer
171: @*/
172: PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
173: {
177:   if (!dm->ceedERestrict) {
178:     PetscInt            num_cells, cell_size, num_comp, lvec_size, *restr_indices;
179:     CeedElemRestriction elem_restr;
180:     Ceed                ceed;

182:     DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices);
183:     DMGetCeed(dm, &ceed);
184:     CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices, &elem_restr);
185:     PetscFree(restr_indices);
186:     dm->ceedERestrict = elem_restr;
187:   }
188:   *ERestrict = dm->ceedERestrict;
189:   return 0;
190: }

192: #endif