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