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 PetscErrorCode ierr; 29 30 PetscFunctionBeginUser; 31 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 32 PetscDS ds = NULL; 33 PetscFE fe; 34 PetscSection section; 35 PetscInt dim; 36 PetscInt *restr_indices; 37 const PetscInt *iter_indices; 38 IS iter_is; 39 40 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 41 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 42 if (domain_label) { 43 ierr = DMGetFirstLabelEntry_Internal(dm, dm, domain_label, 1, &label_value, dim, NULL, &ds);CHKERRQ(ierr); 44 } 45 46 // Translate dm_field to ds_field 47 PetscInt ds_field = -1; 48 for (PetscInt i=0; i<dm->Nds; i++) { 49 if (!domain_label || domain_label == dm->probs[i].label) { 50 ds = dm->probs[i].ds; 51 } 52 if (ds == dm->probs[i].ds) { 53 const PetscInt *arr; 54 PetscInt nf; 55 IS is = dm->probs[i].fields; 56 ierr = ISGetIndices(is, &arr);CHKERRQ(ierr); 57 ierr = ISGetSize(is, &nf);CHKERRQ(ierr); 58 for (PetscInt j=0; j<nf; j++) { 59 if (dm_field == arr[j]) { 60 ds_field = j; 61 break; 62 } 63 } 64 ierr = ISRestoreIndices(is, &arr);CHKERRQ(ierr); 65 } 66 } 67 if (ds_field == -1) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Could not find dm_field %D in DS", dm_field); 68 69 { 70 PetscInt depth; 71 DMLabel depth_label; 72 IS depth_is; 73 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 74 ierr = DMPlexGetDepthLabel(dm, &depth_label);CHKERRQ(ierr); 75 ierr = DMLabelGetStratumIS(depth_label, depth - height, &depth_is);CHKERRQ(ierr); 76 if (domain_label) { 77 IS domain_is; 78 ierr = DMLabelGetStratumIS(domain_label, label_value, &domain_is);CHKERRQ(ierr); 79 if (domain_is) { // domainIS is non-empty 80 ierr = ISIntersect(depth_is, domain_is, &iter_is);CHKERRQ(ierr); 81 ierr = ISDestroy(&domain_is);CHKERRQ(ierr); 82 } else { // domainIS is NULL (empty) 83 iter_is = NULL; 84 } 85 ierr = ISDestroy(&depth_is);CHKERRQ(ierr); 86 } else { 87 iter_is = depth_is; 88 } 89 if (iter_is) { 90 ierr = ISGetLocalSize(iter_is, num_cells);CHKERRQ(ierr); 91 ierr = ISGetIndices(iter_is, &iter_indices);CHKERRQ(ierr); 92 } else { 93 *num_cells = 0; 94 iter_indices = NULL; 95 } 96 } 97 98 { 99 PetscDualSpace dual_space; 100 PetscInt num_dual_basis_vectors; 101 ierr = PetscDSGetDiscretization(ds, ds_field, (PetscObject*)&fe);CHKERRQ(ierr); 102 ierr = PetscFEGetHeightSubspace(fe, height, &fe);CHKERRQ(ierr); 103 ierr = PetscFEGetDualSpace(fe, &dual_space);CHKERRQ(ierr); 104 ierr = PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors);CHKERRQ(ierr); 105 ierr = PetscDualSpaceGetNumComponents(dual_space, num_comp);CHKERRQ(ierr); 106 if (num_dual_basis_vectors % *num_comp != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for number of dual basis vectors %D not divisible by %D components", num_dual_basis_vectors, *num_comp); 107 *cell_size = num_dual_basis_vectors / *num_comp; 108 } 109 PetscInt restr_size = (*num_cells)*(*cell_size); 110 ierr = PetscMalloc1(restr_size, &restr_indices);CHKERRQ(ierr); 111 PetscInt cell_offset = 0; 112 113 PetscInt P = (PetscInt) pow(*cell_size, 1.0 / (dim - height)); 114 for (PetscInt p = 0; p < *num_cells; p++) { 115 PetscBool flip = PETSC_FALSE; 116 PetscInt c = iter_indices[p]; 117 PetscInt num_indices, *indices; 118 PetscInt field_offsets[17]; // max number of fields plus 1 119 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL);CHKERRQ(ierr); 120 if (height > 0) { 121 PetscInt num_cells_support, num_faces, start = -1; 122 const PetscInt *orients, *faces, *cells; 123 ierr = DMPlexGetSupport(dm, c, &cells);CHKERRQ(ierr); 124 ierr = DMPlexGetSupportSize(dm, c, &num_cells_support);CHKERRQ(ierr); 125 if (num_cells_support != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected one cell in support of exterior face, but got %D cells", num_cells_support); 126 ierr = DMPlexGetCone(dm, cells[0], &faces);CHKERRQ(ierr); 127 ierr = DMPlexGetConeSize(dm, cells[0], &num_faces);CHKERRQ(ierr); 128 for (PetscInt i=0; i<num_faces; i++) {if (faces[i] == c) start = i;} 129 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %D in cone of its support", c); 130 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients);CHKERRQ(ierr); 131 if (orients[start] < 0) flip = PETSC_TRUE; 132 } 133 134 for (PetscInt i = 0; i < *cell_size; i++) { 135 PetscInt ii = i; 136 if (flip) { 137 if (*cell_size == P) ii = *cell_size - 1 - i; 138 else if (*cell_size == P*P) { 139 PetscInt row = i / P, col = i % P; 140 ii = row + col * P; 141 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for flipping point with cell size %D != P (%D) or P^2", *cell_size, P); 142 } 143 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 144 PetscInt loc = indices[field_offsets[dm_field] + ii*(*num_comp)]; 145 restr_indices[cell_offset++] = loc >= 0 ? loc : -(loc + 1); 146 } 147 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL);CHKERRQ(ierr); 148 } 149 if (cell_offset != restr_size) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, offsets array of shape (%D, %D) initialized for %D nodes", *num_cells, (*cell_size), cell_offset); 150 if (iter_is) { ierr = ISRestoreIndices(iter_is, &iter_indices);CHKERRQ(ierr); } 151 ierr = ISDestroy(&iter_is);CHKERRQ(ierr); 152 153 *offsets = restr_indices; 154 ierr = PetscSectionGetStorageSize(section, l_size);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #if defined(PETSC_HAVE_LIBCEED) 159 #include <petscdmplexceed.h> 160 161 /*@C 162 DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector) 163 164 Input Parameters: 165 dm - The DMPlex object 166 domain_label - label for DMPlex domain, or NULL for the whole domain 167 label_value - Stratum value 168 height - Height of target cells in DMPlex topology 169 dm_field - Index of DMPlex field 170 171 Output Parameters: 172 ERestrict - libCEED restriction from local vector to to the cells 173 174 Level: developer 175 @*/ 176 PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict) 177 { 178 PetscErrorCode ierr; 179 180 PetscFunctionBeginUser; 181 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 182 PetscValidPointer(ERestrict, 2); 183 if (!dm->ceedERestrict) { 184 PetscInt num_cells, cell_size, num_comp, lvec_size, *restr_indices; 185 CeedElemRestriction elem_restr; 186 Ceed ceed; 187 188 ierr = DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices);CHKERRQ(ierr); 189 190 ierr = DMGetCeed(dm, &ceed);CHKERRQ(ierr); 191 ierr = CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices, &elem_restr);CHKERRQ_CEED(ierr); 192 ierr = PetscFree(restr_indices);CHKERRQ(ierr); 193 dm->ceedERestrict = elem_restr; 194 } 195 *ERestrict = dm->ceedERestrict; 196 PetscFunctionReturn(0); 197 } 198 199 #endif 200