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