1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 /*@C 4 DMPlexGetLocalOffsets - Allocate and populate array of local offsets. 5 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 12 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 19 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(). 21 22 Level: developer 23 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; 35 36 PetscFunctionBeginUser; 37 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 38 PetscCall(DMGetLocalSection(dm, §ion)); 39 PetscCall(DMGetDimension(dm, &dim)); 40 { 41 IS field_is; 42 const PetscInt *fields; 43 PetscInt num_fields; 44 45 PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds)); 46 // Translate dm_field to ds_field 47 PetscCall(ISGetIndices(field_is, &fields)); 48 PetscCall(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 PetscCall(ISRestoreIndices(field_is, &fields)); 56 } 57 PetscCheck(ds_field != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field); 58 59 { 60 PetscInt depth; 61 DMLabel depth_label; 62 IS depth_is; 63 64 PetscCall(DMPlexGetDepth(dm, &depth)); 65 PetscCall(DMPlexGetDepthLabel(dm, &depth_label)); 66 PetscCall(DMLabelGetStratumIS(depth_label, depth - height, &depth_is)); 67 if (domain_label) { 68 IS domain_is; 69 70 PetscCall(DMLabelGetStratumIS(domain_label, label_value, &domain_is)); 71 if (domain_is) { // domainIS is non-empty 72 PetscCall(ISIntersect(depth_is, domain_is, &iter_is)); 73 PetscCall(ISDestroy(&domain_is)); 74 } else { // domainIS is NULL (empty) 75 iter_is = NULL; 76 } 77 PetscCall(ISDestroy(&depth_is)); 78 } else { 79 iter_is = depth_is; 80 } 81 if (iter_is) { 82 PetscCall(ISGetLocalSize(iter_is, num_cells)); 83 PetscCall(ISGetIndices(iter_is, &iter_indices)); 84 } else { 85 *num_cells = 0; 86 iter_indices = NULL; 87 } 88 } 89 90 { 91 PetscDualSpace dual_space; 92 PetscInt num_dual_basis_vectors; 93 94 PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe)); 95 PetscCall(PetscFEGetHeightSubspace(fe, height, &fe)); 96 PetscCheck(fe, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Height %" PetscInt_FMT " is invalid for DG coordinates", height); 97 PetscCall(PetscFEGetDualSpace(fe, &dual_space)); 98 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors)); 99 PetscCall(PetscDualSpaceGetNumComponents(dual_space, num_comp)); 100 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); 101 *cell_size = num_dual_basis_vectors / *num_comp; 102 } 103 PetscInt restr_size = (*num_cells) * (*cell_size); 104 PetscCall(PetscMalloc1(restr_size, &restr_indices)); 105 PetscInt cell_offset = 0; 106 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 PetscCall(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 PetscCall(DMPlexGetSupport(dm, c, &cells)); 118 PetscCall(DMPlexGetSupportSize(dm, c, &num_cells_support)); 119 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); 120 PetscCall(DMPlexGetCone(dm, cells[0], &faces)); 121 PetscCall(DMPlexGetConeSize(dm, cells[0], &num_faces)); 122 for (PetscInt i = 0; i < num_faces; i++) { 123 if (faces[i] == c) start = i; 124 } 125 PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %" PetscInt_FMT " in cone of its support", c); 126 PetscCall(DMPlexGetConeOrientation(dm, cells[0], &orients)); 127 if (orients[start] < 0) flip = PETSC_TRUE; 128 } 129 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 PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL)); 144 } 145 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); 146 if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices)); 147 PetscCall(ISDestroy(&iter_is)); 148 149 *offsets = restr_indices; 150 PetscCall(PetscSectionGetStorageSize(section, l_size)); 151 PetscFunctionReturn(0); 152 } 153 154 #if defined(PETSC_HAVE_LIBCEED) 155 #include <petscdmplexceed.h> 156 157 /*@C 158 DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector) 159 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 166 167 Output Parameters: 168 ERestrict - libCEED restriction from local vector to to the cells 169 170 Level: developer 171 @*/ 172 PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict) 173 { 174 PetscFunctionBeginUser; 175 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 176 PetscValidPointer(ERestrict, 2); 177 if (!dm->ceedERestrict) { 178 PetscInt num_cells, cell_size, num_comp, lvec_size, *restr_indices; 179 CeedElemRestriction elem_restr; 180 Ceed ceed; 181 182 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices)); 183 PetscCall(DMGetCeed(dm, &ceed)); 184 PetscCallCEED(CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices, &elem_restr)); 185 PetscCall(PetscFree(restr_indices)); 186 dm->ceedERestrict = elem_restr; 187 } 188 *ERestrict = dm->ceedERestrict; 189 PetscFunctionReturn(0); 190 } 191 192 #endif 193