1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <../src/sys/utils/hash.h> 3 #include <petsc-private/isimpl.h> 4 #include <petscsf.h> 5 6 /** hierarchy routines */ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMPlexSetReferenceTree" 10 /*@ 11 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes. 12 13 Not collective 14 15 Input Parameters: 16 + dm - The DMPlex object 17 - ref - The reference tree DMPlex object 18 19 Level: beginner 20 21 .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree() 22 @*/ 23 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref) 24 { 25 DM_Plex *mesh = (DM_Plex *)dm->data; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 30 PetscValidHeaderSpecific(ref, DM_CLASSID, 2); 31 ierr = PetscObjectReference((PetscObject)ref);CHKERRQ(ierr); 32 ierr = DMDestroy(&mesh->referenceTree);CHKERRQ(ierr); 33 mesh->referenceTree = ref; 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "DMPlexGetReferenceTree" 39 /*@ 40 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes. 41 42 Not collective 43 44 Input Parameters: 45 . dm - The DMPlex object 46 47 Output Parameters 48 . ref - The reference tree DMPlex object 49 50 Level: beginner 51 52 .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree() 53 @*/ 54 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref) 55 { 56 DM_Plex *mesh = (DM_Plex *)dm->data; 57 58 PetscFunctionBegin; 59 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 60 PetscValidPointer(ref,2); 61 *ref = mesh->referenceTree; 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "DMPlexCreateDefaultReferenceTree" 67 /*@ 68 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 69 70 Collective on comm 71 72 Input Parameters: 73 + comm - the MPI communicator 74 . dim - the spatial dimension 75 - simplex - Flag for simplex, otherwise use a tensor-product cell 76 77 Output Parameters: 78 . ref - the reference tree DMPlex object 79 80 Level: beginner 81 82 .keywords: reference cell 83 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 84 @*/ 85 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 86 { 87 DM K, Kref; 88 PetscInt p, pStart, pEnd, pRefStart, pRefEnd, d, offset; 89 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 90 DMLabel identity, identityRef; 91 PetscSection unionSection, unionConeSection; 92 PetscScalar *unionCoords; 93 IS perm; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 /* create a reference element */ 98 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 99 ierr = DMPlexCreateLabel(K, "identity");CHKERRQ(ierr); 100 ierr = DMPlexGetLabel(K, "identity", &identity);CHKERRQ(ierr); 101 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 102 for (p = pStart; p < pEnd; p++) { 103 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 104 } 105 /* refine it */ 106 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 107 108 /* the reference tree is the union of these two, without duplicating 109 * points that appear in both */ 110 ierr = DMPlexGetLabel(Kref, "identity", &identityRef);CHKERRQ(ierr); 111 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 112 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 113 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 114 /* count points that will go in the union */ 115 for (p = pStart; p < pEnd; p++) { 116 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 117 } 118 for (p = pRefStart; p < pRefEnd; p++) { 119 PetscInt q, qSize; 120 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 121 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 122 if (qSize > 1) { 123 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 124 } 125 } 126 ierr = PetscMalloc1((pEnd - pStart) + (pRefEnd - pRefStart),&permvals);CHKERRQ(ierr); 127 offset = 0; 128 /* stratify points in the union by topological dimension */ 129 for (d = 0; d <= dim; d++) { 130 PetscInt cStart, cEnd, c; 131 132 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 133 for (c = cStart; c < cEnd; c++) { 134 permvals[offset++] = c; 135 } 136 137 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 138 for (c = cStart; c < cEnd; c++) { 139 permvals[offset++] = c + (pEnd - pStart); 140 } 141 } 142 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 143 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 144 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 145 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 146 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 147 /* count dimension points */ 148 for (d = 0; d <= dim; d++) { 149 PetscInt cStart, cOff, cOff2; 150 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 151 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 152 if (d < dim) { 153 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 154 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 155 } 156 else { 157 cOff2 = numUnionPoints; 158 } 159 numDimPoints[dim - d] = cOff2 - cOff; 160 } 161 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 162 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 163 /* count the cones in the union */ 164 for (p = pStart; p < pEnd; p++) { 165 PetscInt dof, uOff; 166 167 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 168 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 169 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 170 coneSizes[uOff] = dof; 171 } 172 for (p = pRefStart; p < pRefEnd; p++) { 173 PetscInt dof, uDof, uOff; 174 175 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 176 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 177 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 178 if (uDof) { 179 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 180 coneSizes[uOff] = dof; 181 } 182 } 183 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 184 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 185 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 186 /* write the cones in the union */ 187 for (p = pStart; p < pEnd; p++) { 188 PetscInt dof, uOff, c, cOff; 189 const PetscInt *cone, *orientation; 190 191 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 192 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 193 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 194 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 195 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 196 for (c = 0; c < dof; c++) { 197 PetscInt e, eOff; 198 e = cone[c]; 199 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 200 unionCones[cOff + c] = eOff; 201 unionOrientations[cOff + c] = orientation[c]; 202 } 203 } 204 for (p = pRefStart; p < pRefEnd; p++) { 205 PetscInt dof, uDof, uOff, c, cOff; 206 const PetscInt *cone, *orientation; 207 208 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 209 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 210 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 211 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 212 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 213 if (uDof) { 214 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 215 for (c = 0; c < dof; c++) { 216 PetscInt e, eOff, eDof; 217 218 e = cone[c]; 219 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 220 if (eDof) { 221 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 222 } 223 else { 224 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 225 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 226 } 227 unionCones[cOff + c] = eOff; 228 unionOrientations[cOff + c] = orientation[c]; 229 } 230 } 231 } 232 /* get the coordinates */ 233 { 234 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 235 PetscSection KcoordsSec, KrefCoordsSec; 236 Vec KcoordsVec, KrefCoordsVec; 237 PetscScalar *Kcoords; 238 239 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 240 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 241 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 242 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 243 244 numVerts = numDimPoints[0]; 245 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 246 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 247 248 offset = 0; 249 for (v = vStart; v < vEnd; v++) { 250 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 251 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 252 for (d = 0; d < dim; d++) { 253 unionCoords[offset * dim + d] = Kcoords[d]; 254 } 255 offset++; 256 } 257 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 258 for (v = vRefStart; v < vRefEnd; v++) { 259 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 260 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 261 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 262 if (vDof) { 263 for (d = 0; d < dim; d++) { 264 unionCoords[offset * dim + d] = Kcoords[d]; 265 } 266 offset++; 267 } 268 } 269 } 270 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 271 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 272 ierr = DMPlexSetDimension(*ref,dim);CHKERRQ(ierr); 273 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 274 /* clean up */ 275 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 276 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 277 ierr = ISDestroy(&perm);CHKERRQ(ierr); 278 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 279 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 280 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 281 ierr = DMDestroy(&K);CHKERRQ(ierr); 282 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286