1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/petscfeimpl.h> /* For PetscFEInterpolate_Static() */ 3 4 #include <petscdmplextransform.h> 5 #include <petscsf.h> 6 7 /*@ 8 DMPlexCreateProcessSF - Create an SF which just has process connectivity 9 10 Collective on dm 11 12 Input Parameters: 13 + dm - The DM 14 - sfPoint - The PetscSF which encodes point connectivity 15 16 Output Parameters: 17 + processRanks - A list of process neighbors, or NULL 18 - sfProcess - An SF encoding the process connectivity, or NULL 19 20 Level: developer 21 22 .seealso: PetscSFCreate(), DMPlexCreateTwoSidedProcessSF() 23 @*/ 24 PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess) 25 { 26 PetscInt numRoots, numLeaves, l; 27 const PetscInt *localPoints; 28 const PetscSFNode *remotePoints; 29 PetscInt *localPointsNew; 30 PetscSFNode *remotePointsNew; 31 PetscInt *ranks, *ranksNew; 32 PetscMPIInt size; 33 34 PetscFunctionBegin; 35 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 36 PetscValidHeaderSpecific(sfPoint, PETSCSF_CLASSID, 2); 37 if (processRanks) {PetscValidPointer(processRanks, 3);} 38 if (sfProcess) {PetscValidPointer(sfProcess, 4);} 39 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 40 CHKERRQ(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints)); 41 CHKERRQ(PetscMalloc1(numLeaves, &ranks)); 42 for (l = 0; l < numLeaves; ++l) { 43 ranks[l] = remotePoints[l].rank; 44 } 45 CHKERRQ(PetscSortRemoveDupsInt(&numLeaves, ranks)); 46 CHKERRQ(PetscMalloc1(numLeaves, &ranksNew)); 47 CHKERRQ(PetscMalloc1(numLeaves, &localPointsNew)); 48 CHKERRQ(PetscMalloc1(numLeaves, &remotePointsNew)); 49 for (l = 0; l < numLeaves; ++l) { 50 ranksNew[l] = ranks[l]; 51 localPointsNew[l] = l; 52 remotePointsNew[l].index = 0; 53 remotePointsNew[l].rank = ranksNew[l]; 54 } 55 CHKERRQ(PetscFree(ranks)); 56 if (processRanks) CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks)); 57 else CHKERRQ(PetscFree(ranksNew)); 58 if (sfProcess) { 59 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess)); 60 CHKERRQ(PetscObjectSetName((PetscObject) *sfProcess, "Process SF")); 61 CHKERRQ(PetscSFSetFromOptions(*sfProcess)); 62 CHKERRQ(PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 63 } 64 PetscFunctionReturn(0); 65 } 66 67 /*@ 68 DMPlexCreateCoarsePointIS - Creates an IS covering the coarse DM chart with the fine points as data 69 70 Collective on dm 71 72 Input Parameter: 73 . dm - The coarse DM 74 75 Output Parameter: 76 . fpointIS - The IS of all the fine points which exist in the original coarse mesh 77 78 Level: developer 79 80 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetSubpointIS() 81 @*/ 82 PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS) 83 { 84 DMPlexTransform tr; 85 PetscInt *fpoints; 86 PetscInt pStart, pEnd, p, vStart, vEnd, v; 87 88 PetscFunctionBegin; 89 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 90 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 91 CHKERRQ(DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr)); 92 CHKERRQ(DMPlexTransformSetUp(tr)); 93 CHKERRQ(PetscMalloc1(pEnd-pStart, &fpoints)); 94 for (p = 0; p < pEnd-pStart; ++p) fpoints[p] = -1; 95 for (v = vStart; v < vEnd; ++v) { 96 PetscInt vNew = -1; /* quiet overzealous may be used uninitialized check */ 97 98 CHKERRQ(DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew)); 99 fpoints[v-pStart] = vNew; 100 } 101 CHKERRQ(DMPlexTransformDestroy(&tr)); 102 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, pEnd-pStart, fpoints, PETSC_OWN_POINTER, fpointIS)); 103 PetscFunctionReturn(0); 104 } 105 106 /*@C 107 DMPlexSetTransformType - Set the transform type for uniform refinement 108 109 Input Parameters: 110 + dm - The DM 111 - type - The transform type for uniform refinement 112 113 Level: developer 114 115 .seealso: DMPlexTransformType, DMRefine(), DMPlexGetTransformType(), DMPlexSetRefinementUniform() 116 @*/ 117 PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type) 118 { 119 DM_Plex *mesh = (DM_Plex*) dm->data; 120 121 PetscFunctionBegin; 122 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMPLEX); 123 if (type) PetscValidCharPointer(type, 2); 124 CHKERRQ(PetscFree(mesh->transformType)); 125 CHKERRQ(PetscStrallocpy(type, &mesh->transformType)); 126 PetscFunctionReturn(0); 127 } 128 129 /*@C 130 DMPlexGetTransformType - Retrieve the transform type for uniform refinement 131 132 Input Parameter: 133 . dm - The DM 134 135 Output Parameter: 136 . type - The transform type for uniform refinement 137 138 Level: developer 139 140 .seealso: DMPlexTransformType, DMRefine(), DMPlexSetTransformType(), DMPlexGetRefinementUniform() 141 @*/ 142 PetscErrorCode DMPlexGetTransformType(DM dm, DMPlexTransformType *type) 143 { 144 DM_Plex *mesh = (DM_Plex*) dm->data; 145 146 PetscFunctionBegin; 147 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMPLEX); 148 PetscValidPointer(type, 2); 149 *type = mesh->transformType; 150 PetscFunctionReturn(0); 151 } 152 153 /*@ 154 DMPlexSetRefinementUniform - Set the flag for uniform refinement 155 156 Input Parameters: 157 + dm - The DM 158 - refinementUniform - The flag for uniform refinement 159 160 Level: developer 161 162 .seealso: DMRefine(), DMPlexGetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 163 @*/ 164 PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform) 165 { 166 DM_Plex *mesh = (DM_Plex*) dm->data; 167 168 PetscFunctionBegin; 169 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMPLEX); 170 mesh->refinementUniform = refinementUniform; 171 PetscFunctionReturn(0); 172 } 173 174 /*@ 175 DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement 176 177 Input Parameter: 178 . dm - The DM 179 180 Output Parameter: 181 . refinementUniform - The flag for uniform refinement 182 183 Level: developer 184 185 .seealso: DMRefine(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 186 @*/ 187 PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform) 188 { 189 DM_Plex *mesh = (DM_Plex*) dm->data; 190 191 PetscFunctionBegin; 192 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMPLEX); 193 PetscValidPointer(refinementUniform, 2); 194 *refinementUniform = mesh->refinementUniform; 195 PetscFunctionReturn(0); 196 } 197 198 /*@ 199 DMPlexSetRefinementLimit - Set the maximum cell volume for refinement 200 201 Input Parameters: 202 + dm - The DM 203 - refinementLimit - The maximum cell volume in the refined mesh 204 205 Level: developer 206 207 .seealso: DMRefine(), DMPlexGetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform() 208 @*/ 209 PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit) 210 { 211 DM_Plex *mesh = (DM_Plex*) dm->data; 212 213 PetscFunctionBegin; 214 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 215 mesh->refinementLimit = refinementLimit; 216 PetscFunctionReturn(0); 217 } 218 219 /*@ 220 DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement 221 222 Input Parameter: 223 . dm - The DM 224 225 Output Parameter: 226 . refinementLimit - The maximum cell volume in the refined mesh 227 228 Level: developer 229 230 .seealso: DMRefine(), DMPlexSetRefinementLimit(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform() 231 @*/ 232 PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit) 233 { 234 DM_Plex *mesh = (DM_Plex*) dm->data; 235 236 PetscFunctionBegin; 237 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 238 PetscValidPointer(refinementLimit, 2); 239 /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */ 240 *refinementLimit = mesh->refinementLimit; 241 PetscFunctionReturn(0); 242 } 243 244 /*@ 245 DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement 246 247 Input Parameters: 248 + dm - The DM 249 - refinementFunc - Function giving the maximum cell volume in the refined mesh 250 251 Note: The calling sequence is refinementFunc(coords, limit) 252 $ coords - Coordinates of the current point, usually a cell centroid 253 $ limit - The maximum cell volume for a cell containing this point 254 255 Level: developer 256 257 .seealso: DMRefine(), DMPlexGetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 258 @*/ 259 PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *)) 260 { 261 DM_Plex *mesh = (DM_Plex*) dm->data; 262 263 PetscFunctionBegin; 264 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 265 mesh->refinementFunc = refinementFunc; 266 PetscFunctionReturn(0); 267 } 268 269 /*@ 270 DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement 271 272 Input Parameter: 273 . dm - The DM 274 275 Output Parameter: 276 . refinementFunc - Function giving the maximum cell volume in the refined mesh 277 278 Note: The calling sequence is refinementFunc(coords, limit) 279 $ coords - Coordinates of the current point, usually a cell centroid 280 $ limit - The maximum cell volume for a cell containing this point 281 282 Level: developer 283 284 .seealso: DMRefine(), DMPlexSetRefinementFunction(), DMPlexGetRefinementUniform(), DMPlexSetRefinementUniform(), DMPlexGetRefinementLimit(), DMPlexSetRefinementLimit() 285 @*/ 286 PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal [], PetscReal *)) 287 { 288 DM_Plex *mesh = (DM_Plex*) dm->data; 289 290 PetscFunctionBegin; 291 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 292 PetscValidPointer(refinementFunc, 2); 293 *refinementFunc = mesh->refinementFunc; 294 PetscFunctionReturn(0); 295 } 296 297 PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm) 298 { 299 PetscBool isUniform; 300 301 PetscFunctionBegin; 302 CHKERRQ(DMPlexGetRefinementUniform(dm, &isUniform)); 303 CHKERRQ(DMViewFromOptions(dm, NULL, "-initref_dm_view")); 304 if (isUniform) { 305 DMPlexTransform tr; 306 DM cdm, rcdm; 307 DMPlexTransformType trType; 308 const char *prefix; 309 PetscOptions options; 310 311 CHKERRQ(DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr)); 312 CHKERRQ(DMPlexTransformSetDM(tr, dm)); 313 CHKERRQ(DMPlexGetTransformType(dm, &trType)); 314 if (trType) CHKERRQ(DMPlexTransformSetType(tr, trType)); 315 CHKERRQ(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 316 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) tr, prefix)); 317 CHKERRQ(PetscObjectGetOptions((PetscObject) dm, &options)); 318 CHKERRQ(PetscObjectSetOptions((PetscObject) tr, options)); 319 CHKERRQ(DMPlexTransformSetFromOptions(tr)); 320 CHKERRQ(PetscObjectSetOptions((PetscObject) tr, NULL)); 321 CHKERRQ(DMPlexTransformSetUp(tr)); 322 CHKERRQ(PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view")); 323 CHKERRQ(DMPlexTransformApply(tr, dm, rdm)); 324 CHKERRQ(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE)); 325 CHKERRQ(DMCopyDisc(dm, *rdm)); 326 CHKERRQ(DMGetCoordinateDM(dm, &cdm)); 327 CHKERRQ(DMGetCoordinateDM(*rdm, &rcdm)); 328 CHKERRQ(DMCopyDisc(cdm, rcdm)); 329 CHKERRQ(DMPlexTransformCreateDiscLabels(tr, *rdm)); 330 CHKERRQ(DMPlexTransformDestroy(&tr)); 331 } else { 332 CHKERRQ(DMPlexRefine_Internal(dm, NULL, NULL, NULL, rdm)); 333 } 334 if (*rdm) { 335 ((DM_Plex *) (*rdm)->data)->printFEM = ((DM_Plex *) dm->data)->printFEM; 336 ((DM_Plex *) (*rdm)->data)->printL2 = ((DM_Plex *) dm->data)->printL2; 337 } 338 CHKERRQ(DMViewFromOptions(dm, NULL, "-postref_dm_view")); 339 PetscFunctionReturn(0); 340 } 341 342 PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[]) 343 { 344 DM cdm = dm; 345 PetscInt r; 346 PetscBool isUniform, localized; 347 348 PetscFunctionBegin; 349 CHKERRQ(DMPlexGetRefinementUniform(dm, &isUniform)); 350 CHKERRQ(DMGetCoordinatesLocalized(dm, &localized)); 351 if (isUniform) { 352 for (r = 0; r < nlevels; ++r) { 353 DMPlexTransform tr; 354 DM codm, rcodm; 355 const char *prefix; 356 357 CHKERRQ(DMPlexTransformCreate(PetscObjectComm((PetscObject) cdm), &tr)); 358 CHKERRQ(PetscObjectGetOptionsPrefix((PetscObject) cdm, &prefix)); 359 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) tr, prefix)); 360 CHKERRQ(DMPlexTransformSetDM(tr, cdm)); 361 CHKERRQ(DMPlexTransformSetFromOptions(tr)); 362 CHKERRQ(DMPlexTransformSetUp(tr)); 363 CHKERRQ(DMPlexTransformApply(tr, cdm, &rdm[r])); 364 CHKERRQ(DMSetCoarsenLevel(rdm[r], cdm->leveldown)); 365 CHKERRQ(DMSetRefineLevel(rdm[r], cdm->levelup+1)); 366 CHKERRQ(DMCopyDisc(cdm, rdm[r])); 367 CHKERRQ(DMGetCoordinateDM(dm, &codm)); 368 CHKERRQ(DMGetCoordinateDM(rdm[r], &rcodm)); 369 CHKERRQ(DMCopyDisc(codm, rcodm)); 370 CHKERRQ(DMPlexTransformCreateDiscLabels(tr, rdm[r])); 371 CHKERRQ(DMSetCoarseDM(rdm[r], cdm)); 372 CHKERRQ(DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE)); 373 if (rdm[r]) { 374 ((DM_Plex *) (rdm[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM; 375 ((DM_Plex *) (rdm[r])->data)->printL2 = ((DM_Plex *) dm->data)->printL2; 376 } 377 cdm = rdm[r]; 378 CHKERRQ(DMPlexTransformDestroy(&tr)); 379 } 380 } else { 381 for (r = 0; r < nlevels; ++r) { 382 CHKERRQ(DMRefine(cdm, PetscObjectComm((PetscObject) dm), &rdm[r])); 383 CHKERRQ(DMCopyDisc(cdm, rdm[r])); 384 if (localized) CHKERRQ(DMLocalizeCoordinates(rdm[r])); 385 CHKERRQ(DMSetCoarseDM(rdm[r], cdm)); 386 if (rdm[r]) { 387 ((DM_Plex *) (rdm[r])->data)->printFEM = ((DM_Plex *) dm->data)->printFEM; 388 ((DM_Plex *) (rdm[r])->data)->printL2 = ((DM_Plex *) dm->data)->printL2; 389 } 390 cdm = rdm[r]; 391 } 392 } 393 PetscFunctionReturn(0); 394 } 395