Actual source code: plexceed.c
1: #include <petsc/private/dmpleximpl.h>
3: static PetscErrorCode DMGetPoints_Private(DM dm, DMLabel domainLabel, PetscInt labelVal, PetscInt height, IS *pointIS)
4: {
5: PetscInt depth;
6: DMLabel depthLabel;
7: IS depthIS;
9: PetscFunctionBegin;
10: PetscCall(DMPlexGetDepth(dm, &depth));
11: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
12: PetscCall(DMLabelGetStratumIS(depthLabel, depth - height, &depthIS));
13: if (domainLabel) {
14: IS domainIS;
16: PetscCall(DMLabelGetStratumIS(domainLabel, labelVal, &domainIS));
17: if (domainIS) { // domainIS is non-empty
18: PetscCall(ISIntersect(depthIS, domainIS, pointIS));
19: PetscCall(ISDestroy(&domainIS));
20: } else { // domainIS is NULL (empty)
21: *pointIS = NULL;
22: }
23: PetscCall(ISDestroy(&depthIS));
24: } else {
25: *pointIS = depthIS;
26: }
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: /*@C
31: DMPlexGetLocalOffsets - Allocate and populate array of local offsets for each cell closure.
33: Not collective
35: Input Parameters:
36: + dm - The `DMPLEX` object
37: . domain_label - label for `DMPLEX` domain, or NULL for whole domain
38: . label_value - Stratum value
39: . height - Height of target cells in `DMPLEX` topology
40: - dm_field - Index of `DMPLEX` field
42: Output Parameters:
43: + num_cells - Number of local cells
44: . cell_size - Size of each cell, given by cell_size * num_comp = num_dof
45: . num_comp - Number of components per dof
46: . l_size - Size of local vector
47: - offsets - Allocated offsets array for cells
49: Level: developer
51: Notes:
52: 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].
54: Caller is responsible for freeing the offsets array using `PetscFree()`.
56: .seealso: [](ch_unstructured), `DMPlexGetLocalOffsetsSupport()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
57: @*/
58: 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)
59: {
60: PetscDS ds = NULL;
61: PetscFE fe;
62: PetscSection section;
63: PetscInt dim, ds_field = -1;
64: PetscInt *restr_indices;
65: const PetscInt *iter_indices;
66: IS iter_is;
68: PetscFunctionBeginUser;
70: PetscCall(PetscLogEventBegin(DMPLEX_GetLocalOffsets, dm, 0, 0, 0));
71: PetscCall(DMGetLocalSection(dm, §ion));
72: PetscCall(DMGetDimension(dm, &dim));
73: {
74: IS field_is;
75: const PetscInt *fields;
76: PetscInt num_fields;
78: PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL));
79: // Translate dm_field to ds_field
80: PetscCall(ISGetIndices(field_is, &fields));
81: PetscCall(ISGetSize(field_is, &num_fields));
82: for (PetscInt i = 0; i < num_fields; i++) {
83: if (dm_field == fields[i]) {
84: ds_field = i;
85: break;
86: }
87: }
88: PetscCall(ISRestoreIndices(field_is, &fields));
89: }
90: PetscCheck(ds_field != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field);
92: PetscCall(DMGetPoints_Private(dm, domain_label, label_value, height, &iter_is));
93: if (iter_is) {
94: PetscCall(ISGetLocalSize(iter_is, num_cells));
95: PetscCall(ISGetIndices(iter_is, &iter_indices));
96: } else {
97: *num_cells = 0;
98: iter_indices = NULL;
99: }
101: {
102: PetscDualSpace dual_space;
103: PetscInt num_dual_basis_vectors;
105: PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
106: PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
107: PetscCheck(fe, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Height %" PetscInt_FMT " is invalid for DG coordinates", height);
108: PetscCall(PetscFEGetDualSpace(fe, &dual_space));
109: PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
110: PetscCall(PetscDualSpaceGetNumComponents(dual_space, num_comp));
111: PetscCheck(num_dual_basis_vectors % *num_comp == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for number of dual basis vectors %" PetscInt_FMT " not divisible by %" PetscInt_FMT " components", num_dual_basis_vectors, *num_comp);
112: *cell_size = num_dual_basis_vectors / *num_comp;
113: }
114: PetscInt restr_size = (*num_cells) * (*cell_size);
115: PetscCall(PetscMalloc1(restr_size, &restr_indices));
116: PetscInt cell_offset = 0;
118: PetscInt P = (PetscInt)PetscPowReal(*cell_size, 1.0 / (dim - height));
119: for (PetscInt p = 0; p < *num_cells; p++) {
120: PetscBool flip = PETSC_FALSE;
121: PetscInt c = iter_indices[p];
122: PetscInt num_indices, *indices;
123: PetscInt field_offsets[17]; // max number of fields plus 1
124: PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
125: if (height > 0) {
126: PetscInt num_cells_support, num_faces, start = -1;
127: const PetscInt *orients, *faces, *cells;
128: PetscCall(DMPlexGetSupport(dm, c, &cells));
129: PetscCall(DMPlexGetSupportSize(dm, c, &num_cells_support));
130: PetscCheck(num_cells_support == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected one cell in support of exterior face, but got %" PetscInt_FMT " cells", num_cells_support);
131: PetscCall(DMPlexGetCone(dm, cells[0], &faces));
132: PetscCall(DMPlexGetConeSize(dm, cells[0], &num_faces));
133: for (PetscInt i = 0; i < num_faces; i++) {
134: if (faces[i] == c) start = i;
135: }
136: PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %" PetscInt_FMT " in cone of its support", c);
137: PetscCall(DMPlexGetConeOrientation(dm, cells[0], &orients));
138: if (orients[start] < 0) flip = PETSC_TRUE;
139: }
141: for (PetscInt i = 0; i < *cell_size; i++) {
142: PetscInt ii = i;
143: if (flip) {
144: if (*cell_size == P) ii = *cell_size - 1 - i;
145: else if (*cell_size == P * P) {
146: PetscInt row = i / P, col = i % P;
147: ii = row + col * P;
148: } 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);
149: }
150: // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
151: PetscInt loc = indices[field_offsets[dm_field] + ii * (*num_comp)];
152: restr_indices[cell_offset++] = loc >= 0 ? loc : -(loc + 1);
153: }
154: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
155: }
156: PetscCheck(cell_offset == restr_size, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, offsets array of shape (%" PetscInt_FMT ", %" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", *num_cells, (*cell_size), cell_offset);
157: if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
158: PetscCall(ISDestroy(&iter_is));
160: *offsets = restr_indices;
161: PetscCall(PetscSectionGetStorageSize(section, l_size));
162: PetscCall(PetscLogEventEnd(DMPLEX_GetLocalOffsets, dm, 0, 0, 0));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /*@C
167: DMPlexGetLocalOffsetsSupport - Allocate and populate arrays of local offsets for each face support.
169: Not collective
171: Input Parameters:
172: + dm - The `DMPLEX` object
173: . domain_label - label for `DMPLEX` domain, or NULL for whole domain
174: - label_value - Stratum value
176: Output Parameters:
177: + num_faces - Number of local, non-boundary faces
178: . num_comp - Number of components per dof
179: . l_size - Size of local vector
180: . offsetsNeg - Allocated offsets array for cells on the inward normal side of each face
181: - offsetsPos - Allocated offsets array for cells on the outward normal side of each face
183: Level: developer
185: Notes:
186: Allocate and populate array of shape [num_cells, num_comp] defining offsets for each cell for local vector of the `DMPLEX` field. All offsets are in the range [0, l_size - 1].
188: Caller is responsible for freeing the offsets array using `PetscFree()`.
190: .seealso: [](ch_unstructured), `DMPlexGetLocalOffsets()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
191: @*/
192: PetscErrorCode DMPlexGetLocalOffsetsSupport(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt *num_faces, PetscInt *num_comp, PetscInt *l_size, PetscInt **offsetsNeg, PetscInt **offsetsPos)
193: {
194: PetscDS ds = NULL;
195: PetscFV fv;
196: PetscSection section;
197: PetscInt dim, height = 1, dm_field = 0, ds_field = 0, Nf, NfInt = 0, cell_size, restr_size;
198: PetscInt *restr_indices_neg, *restr_indices_pos;
199: const PetscInt *iter_indices;
200: IS iter_is;
202: PetscFunctionBeginUser;
204: PetscCall(DMGetLocalSection(dm, §ion));
205: PetscCall(DMGetDimension(dm, &dim));
206: PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));
208: PetscCall(DMGetPoints_Private(dm, domain_label, label_value, height, &iter_is));
209: if (iter_is) {
210: PetscCall(ISGetIndices(iter_is, &iter_indices));
211: PetscCall(ISGetLocalSize(iter_is, &Nf));
212: for (PetscInt p = 0, Ns; p < Nf; ++p) {
213: PetscCall(DMPlexGetSupportSize(dm, iter_indices[p], &Ns));
214: if (Ns == 2) ++NfInt;
215: }
216: *num_faces = NfInt;
217: } else {
218: *num_faces = 0;
219: iter_indices = NULL;
220: }
222: PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fv));
223: PetscCall(PetscFVGetNumComponents(fv, num_comp));
224: cell_size = *num_comp;
225: restr_size = (*num_faces) * cell_size;
226: PetscCall(PetscMalloc1(restr_size, &restr_indices_neg));
227: PetscCall(PetscMalloc1(restr_size, &restr_indices_pos));
228: PetscInt face_offset_neg = 0, face_offset_pos = 0;
230: for (PetscInt p = 0; p < Nf; ++p) {
231: const PetscInt face = iter_indices[p];
232: PetscInt num_indices, *indices;
233: PetscInt field_offsets[17]; // max number of fields plus 1
234: const PetscInt *supp;
235: PetscInt Ns;
237: PetscCall(DMPlexGetSupport(dm, face, &supp));
238: PetscCall(DMPlexGetSupportSize(dm, face, &Ns));
239: // Ignore boundary faces
240: // TODO check for face on parallel boundary
241: if (Ns == 2) {
242: // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
243: PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
244: for (PetscInt i = 0; i < cell_size; i++) {
245: const PetscInt loc = indices[field_offsets[dm_field] + i * (*num_comp)];
246: restr_indices_neg[face_offset_neg++] = loc >= 0 ? loc : -(loc + 1);
247: }
248: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
249: PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
250: for (PetscInt i = 0; i < cell_size; i++) {
251: const PetscInt loc = indices[field_offsets[dm_field] + i * (*num_comp)];
252: restr_indices_pos[face_offset_pos++] = loc >= 0 ? loc : -(loc + 1);
253: }
254: PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
255: }
256: }
257: PetscCheck(face_offset_neg == restr_size, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, neg offsets array of shape (%" PetscInt_FMT ", %" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", *num_faces, cell_size, face_offset_neg);
258: PetscCheck(face_offset_pos == restr_size, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, pos offsets array of shape (%" PetscInt_FMT ", %" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", *num_faces, cell_size, face_offset_pos);
259: if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
260: PetscCall(ISDestroy(&iter_is));
262: *offsetsNeg = restr_indices_neg;
263: *offsetsPos = restr_indices_pos;
264: PetscCall(PetscSectionGetStorageSize(section, l_size));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: #if defined(PETSC_HAVE_LIBCEED)
269: #include <petscdmplexceed.h>
271: /*@C
272: DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector)
274: Input Parameters:
275: + dm - The `DMPLEX` object
276: . domain_label - label for `DMPLEX` domain, or NULL for the whole domain
277: . label_value - Stratum value
278: . height - Height of target cells in `DMPLEX` topology
279: - dm_field - Index of `DMPLEX` field
281: Output Parameter:
282: . ERestrict - libCEED restriction from local vector to to the cells
284: Level: developer
286: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetLocalOffsets()`
287: @*/
288: PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
289: {
290: PetscFunctionBeginUser;
292: PetscAssertPointer(ERestrict, 6);
293: if (!dm->ceedERestrict) {
294: PetscInt num_cells, cell_size, num_comp, lvec_size, *restr_indices;
295: CeedElemRestriction elem_restr;
296: Ceed ceed;
298: PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices));
299: PetscCall(DMGetCeed(dm, &ceed));
300: PetscCallCEED(CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices, &elem_restr));
301: PetscCall(PetscFree(restr_indices));
302: dm->ceedERestrict = elem_restr;
303: }
304: *ERestrict = dm->ceedERestrict;
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: #endif