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);
 96:     PetscFEGetDualSpace(fe, &dual_space);
 97:     PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors);
 98:     PetscDualSpaceGetNumComponents(dual_space, num_comp);
100:     *cell_size = num_dual_basis_vectors / *num_comp;
101:   }
102:   PetscInt restr_size = (*num_cells)*(*cell_size);
103:   PetscMalloc1(restr_size, &restr_indices);
104:   PetscInt cell_offset = 0;

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

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

146:   *offsets = restr_indices;
147:   PetscSectionGetStorageSize(section, l_size);
148:   return 0;
149: }

151: #if defined(PETSC_HAVE_LIBCEED)
152: #include <petscdmplexceed.h>

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

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

164:   Output Parameters:
165:   ERestrict - libCEED restriction from local vector to to the cells

167:   Level: developer
168: @*/
169: PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
170: {
174:   if (!dm->ceedERestrict) {
175:     PetscInt     num_cells, cell_size, num_comp, lvec_size, *restr_indices;
176:     CeedElemRestriction elem_restr;
177:     Ceed         ceed;

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

189: #endif