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, §ion);
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