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);
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