1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/ 2 #include <petscds.h> 3 4 // Greatest common divisor of two nonnegative integers 5 PetscInt PetscGCD(PetscInt a, PetscInt b) { 6 while (b != 0) { 7 PetscInt tmp = b; 8 b = a % b; 9 a = tmp; 10 } 11 return a; 12 } 13 14 PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm,Vec *vec) 15 { 16 PetscSection gSection; 17 PetscInt localSize, bs, blockSize = -1, pStart, pEnd, p; 18 PetscInt in[2],out[2]; 19 20 PetscFunctionBegin; 21 CHKERRQ(DMGetGlobalSection(dm, &gSection)); 22 CHKERRQ(PetscSectionGetChart(gSection, &pStart, &pEnd)); 23 for (p = pStart; p < pEnd; ++p) { 24 PetscInt dof, cdof; 25 26 CHKERRQ(PetscSectionGetDof(gSection, p, &dof)); 27 CHKERRQ(PetscSectionGetConstraintDof(gSection, p, &cdof)); 28 29 if (dof - cdof > 0) { 30 if (blockSize < 0) { 31 /* set blockSize */ 32 blockSize = dof-cdof; 33 } else { 34 blockSize = PetscGCD(dof - cdof, blockSize); 35 } 36 } 37 } 38 39 in[0] = blockSize < 0 ? PETSC_MIN_INT : -blockSize; 40 in[1] = blockSize; 41 CHKERRMPI(MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm))); 42 /* -out[0] = min(blockSize), out[1] = max(blockSize) */ 43 if (-out[0] == out[1]) { 44 bs = out[1]; 45 } else bs = 1; 46 47 if (bs < 0) { /* Everyone was empty */ 48 blockSize = 1; 49 bs = 1; 50 } 51 52 CHKERRQ(PetscSectionGetConstrainedStorageSize(gSection, &localSize)); 53 PetscCheck(localSize%blockSize == 0,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %" PetscInt_FMT " and local storage size %" PetscInt_FMT, blockSize, localSize); 54 CHKERRQ(VecCreate(PetscObjectComm((PetscObject)dm), vec)); 55 CHKERRQ(VecSetSizes(*vec, localSize, PETSC_DETERMINE)); 56 CHKERRQ(VecSetBlockSize(*vec, bs)); 57 CHKERRQ(VecSetType(*vec,dm->vectype)); 58 CHKERRQ(VecSetDM(*vec, dm)); 59 /* CHKERRQ(VecSetLocalToGlobalMapping(*vec, dm->ltogmap)); */ 60 PetscFunctionReturn(0); 61 } 62 63 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec) 64 { 65 PetscSection section; 66 PetscInt localSize, blockSize = -1, pStart, pEnd, p; 67 68 PetscFunctionBegin; 69 CHKERRQ(DMGetLocalSection(dm, §ion)); 70 CHKERRQ(PetscSectionGetChart(section, &pStart, &pEnd)); 71 for (p = pStart; p < pEnd; ++p) { 72 PetscInt dof; 73 74 CHKERRQ(PetscSectionGetDof(section, p, &dof)); 75 if ((blockSize < 0) && (dof > 0)) blockSize = dof; 76 if (dof > 0) blockSize = PetscGCD(dof, blockSize); 77 } 78 CHKERRQ(PetscSectionGetStorageSize(section, &localSize)); 79 CHKERRQ(VecCreate(PETSC_COMM_SELF, vec)); 80 CHKERRQ(VecSetSizes(*vec, localSize, localSize)); 81 CHKERRQ(VecSetBlockSize(*vec, blockSize)); 82 CHKERRQ(VecSetType(*vec,dm->vectype)); 83 CHKERRQ(VecSetDM(*vec, dm)); 84 PetscFunctionReturn(0); 85 } 86 87 /*@C 88 DMCreateSectionSubDM - Returns an IS and subDM+subSection encapsulating a subproblem defined by the fields in a PetscSection in the DM. 89 90 Not collective 91 92 Input Parameters: 93 + dm - The DM object 94 . numFields - The number of fields in this subproblem 95 - fields - The field numbers of the selected fields 96 97 Output Parameters: 98 + is - The global indices for the subproblem 99 - subdm - The DM for the subproblem, which must already have be cloned from dm 100 101 Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes, 102 such as Plex and Forest. 103 104 Level: intermediate 105 106 .seealso DMCreateSubDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView() 107 @*/ 108 PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm) 109 { 110 PetscSection section, sectionGlobal; 111 PetscInt *subIndices; 112 PetscInt subSize = 0, subOff = 0, Nf, f, pStart, pEnd, p; 113 114 PetscFunctionBegin; 115 if (!numFields) PetscFunctionReturn(0); 116 CHKERRQ(DMGetLocalSection(dm, §ion)); 117 CHKERRQ(DMGetGlobalSection(dm, §ionGlobal)); 118 PetscCheck(section,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 119 PetscCheck(sectionGlobal,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 120 CHKERRQ(PetscSectionGetNumFields(section, &Nf)); 121 PetscCheck(numFields <= Nf,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of DM fields %" PetscInt_FMT, numFields, Nf); 122 if (is) { 123 PetscInt bs, bsLocal[2], bsMinMax[2]; 124 125 for (f = 0, bs = 0; f < numFields; ++f) { 126 PetscInt Nc; 127 128 CHKERRQ(PetscSectionGetFieldComponents(section, fields[f], &Nc)); 129 bs += Nc; 130 } 131 CHKERRQ(PetscSectionGetChart(sectionGlobal, &pStart, &pEnd)); 132 for (p = pStart; p < pEnd; ++p) { 133 PetscInt gdof, pSubSize = 0; 134 135 CHKERRQ(PetscSectionGetDof(sectionGlobal, p, &gdof)); 136 if (gdof > 0) { 137 for (f = 0; f < numFields; ++f) { 138 PetscInt fdof, fcdof; 139 140 CHKERRQ(PetscSectionGetFieldDof(section, p, fields[f], &fdof)); 141 CHKERRQ(PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof)); 142 pSubSize += fdof-fcdof; 143 } 144 subSize += pSubSize; 145 if (pSubSize && bs != pSubSize) { 146 /* Layout does not admit a pointwise block size */ 147 bs = 1; 148 } 149 } 150 } 151 /* Must have same blocksize on all procs (some might have no points) */ 152 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 153 CHKERRQ(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax)); 154 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 155 else {bs = bsMinMax[0];} 156 CHKERRQ(PetscMalloc1(subSize, &subIndices)); 157 for (p = pStart; p < pEnd; ++p) { 158 PetscInt gdof, goff; 159 160 CHKERRQ(PetscSectionGetDof(sectionGlobal, p, &gdof)); 161 if (gdof > 0) { 162 CHKERRQ(PetscSectionGetOffset(sectionGlobal, p, &goff)); 163 for (f = 0; f < numFields; ++f) { 164 PetscInt fdof, fcdof, fc, f2, poff = 0; 165 166 /* Can get rid of this loop by storing field information in the global section */ 167 for (f2 = 0; f2 < fields[f]; ++f2) { 168 CHKERRQ(PetscSectionGetFieldDof(section, p, f2, &fdof)); 169 CHKERRQ(PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof)); 170 poff += fdof-fcdof; 171 } 172 CHKERRQ(PetscSectionGetFieldDof(section, p, fields[f], &fdof)); 173 CHKERRQ(PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof)); 174 for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) { 175 subIndices[subOff] = goff+poff+fc; 176 } 177 } 178 } 179 } 180 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is)); 181 if (bs > 1) { 182 /* We need to check that the block size does not come from non-contiguous fields */ 183 PetscInt i, j, set = 1; 184 for (i = 0; i < subSize; i += bs) { 185 for (j = 0; j < bs; ++j) { 186 if (subIndices[i+j] != subIndices[i]+j) {set = 0; break;} 187 } 188 } 189 if (set) CHKERRQ(ISSetBlockSize(*is, bs)); 190 } 191 } 192 if (subdm) { 193 PetscSection subsection; 194 PetscBool haveNull = PETSC_FALSE; 195 PetscInt f, nf = 0, of = 0; 196 197 CHKERRQ(PetscSectionCreateSubsection(section, numFields, fields, &subsection)); 198 CHKERRQ(DMSetLocalSection(*subdm, subsection)); 199 CHKERRQ(PetscSectionDestroy(&subsection)); 200 for (f = 0; f < numFields; ++f) { 201 (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]]; 202 if ((*subdm)->nullspaceConstructors[f]) { 203 haveNull = PETSC_TRUE; 204 nf = f; 205 of = fields[f]; 206 } 207 } 208 if (dm->probs) { 209 CHKERRQ(DMSetNumFields(*subdm, numFields)); 210 for (f = 0; f < numFields; ++f) { 211 PetscObject disc; 212 213 CHKERRQ(DMGetField(dm, fields[f], NULL, &disc)); 214 CHKERRQ(DMSetField(*subdm, f, NULL, disc)); 215 } 216 CHKERRQ(DMCreateDS(*subdm)); 217 if (numFields == 1 && is) { 218 PetscObject disc, space, pmat; 219 220 CHKERRQ(DMGetField(*subdm, 0, NULL, &disc)); 221 CHKERRQ(PetscObjectQuery(disc, "nullspace", &space)); 222 if (space) CHKERRQ(PetscObjectCompose((PetscObject) *is, "nullspace", space)); 223 CHKERRQ(PetscObjectQuery(disc, "nearnullspace", &space)); 224 if (space) CHKERRQ(PetscObjectCompose((PetscObject) *is, "nearnullspace", space)); 225 CHKERRQ(PetscObjectQuery(disc, "pmat", &pmat)); 226 if (pmat) CHKERRQ(PetscObjectCompose((PetscObject) *is, "pmat", pmat)); 227 } 228 /* Check if DSes record their DM fields */ 229 if (dm->probs[0].fields) { 230 PetscInt d, e; 231 232 for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) { 233 const PetscInt Nf = dm->probs[d].ds->Nf; 234 const PetscInt *fld; 235 PetscInt f, g; 236 237 CHKERRQ(ISGetIndices(dm->probs[d].fields, &fld)); 238 for (f = 0; f < Nf; ++f) { 239 for (g = 0; g < numFields; ++g) if (fld[f] == fields[g]) break; 240 if (g < numFields) break; 241 } 242 CHKERRQ(ISRestoreIndices(dm->probs[d].fields, &fld)); 243 if (f == Nf) continue; 244 CHKERRQ(PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds)); 245 CHKERRQ(PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds)); 246 /* Translate DM fields to DS fields */ 247 { 248 IS infields, dsfields; 249 const PetscInt *fld, *ofld; 250 PetscInt *fidx; 251 PetscInt onf, nf, f, g; 252 253 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields)); 254 CHKERRQ(ISIntersect(infields, dm->probs[d].fields, &dsfields)); 255 CHKERRQ(ISDestroy(&infields)); 256 CHKERRQ(ISGetLocalSize(dsfields, &nf)); 257 PetscCheck(nf,PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields"); 258 CHKERRQ(ISGetIndices(dsfields, &fld)); 259 CHKERRQ(ISGetLocalSize(dm->probs[d].fields, &onf)); 260 CHKERRQ(ISGetIndices(dm->probs[d].fields, &ofld)); 261 CHKERRQ(PetscMalloc1(nf, &fidx)); 262 for (f = 0, g = 0; f < onf && g < nf; ++f) { 263 if (ofld[f] == fld[g]) fidx[g++] = f; 264 } 265 CHKERRQ(ISRestoreIndices(dm->probs[d].fields, &ofld)); 266 CHKERRQ(ISRestoreIndices(dsfields, &fld)); 267 CHKERRQ(ISDestroy(&dsfields)); 268 CHKERRQ(PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds)); 269 CHKERRQ(PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds)); 270 CHKERRQ(PetscFree(fidx)); 271 } 272 ++e; 273 } 274 } else { 275 CHKERRQ(PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds)); 276 CHKERRQ(PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds)); 277 CHKERRQ(PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds)); 278 CHKERRQ(PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds)); 279 } 280 } 281 if (haveNull && is) { 282 MatNullSpace nullSpace; 283 284 CHKERRQ((*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace)); 285 CHKERRQ(PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace)); 286 CHKERRQ(MatNullSpaceDestroy(&nullSpace)); 287 } 288 if (dm->coarseMesh) { 289 CHKERRQ(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh)); 290 } 291 } 292 PetscFunctionReturn(0); 293 } 294 295 /*@C 296 DMCreateSectionSuperDM - Returns an arrays of ISes and DM+Section encapsulating a superproblem defined by the DM+Sections passed in. 297 298 Not collective 299 300 Input Parameters: 301 + dms - The DM objects 302 - len - The number of DMs 303 304 Output Parameters: 305 + is - The global indices for the subproblem, or NULL 306 - superdm - The DM for the superproblem, which must already have be cloned 307 308 Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes, 309 such as Plex and Forest. 310 311 Level: intermediate 312 313 .seealso DMCreateSuperDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView() 314 @*/ 315 PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm) 316 { 317 MPI_Comm comm; 318 PetscSection supersection, *sections, *sectionGlobals; 319 PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i; 320 PetscBool haveNull = PETSC_FALSE; 321 322 PetscFunctionBegin; 323 CHKERRQ(PetscObjectGetComm((PetscObject)dms[0], &comm)); 324 /* Pull out local and global sections */ 325 CHKERRQ(PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals)); 326 for (i = 0 ; i < len; ++i) { 327 CHKERRQ(DMGetLocalSection(dms[i], §ions[i])); 328 CHKERRQ(DMGetGlobalSection(dms[i], §ionGlobals[i])); 329 PetscCheck(sections[i],comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 330 PetscCheck(sectionGlobals[i],comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 331 CHKERRQ(PetscSectionGetNumFields(sections[i], &Nfs[i])); 332 Nf += Nfs[i]; 333 } 334 /* Create the supersection */ 335 CHKERRQ(PetscSectionCreateSupersection(sections, len, &supersection)); 336 CHKERRQ(DMSetLocalSection(*superdm, supersection)); 337 /* Create ISes */ 338 if (is) { 339 PetscSection supersectionGlobal; 340 PetscInt bs = -1, startf = 0; 341 342 CHKERRQ(PetscMalloc1(len, is)); 343 CHKERRQ(DMGetGlobalSection(*superdm, &supersectionGlobal)); 344 for (i = 0 ; i < len; startf += Nfs[i], ++i) { 345 PetscInt *subIndices; 346 PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy; 347 348 CHKERRQ(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd)); 349 CHKERRQ(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize)); 350 CHKERRQ(PetscMalloc1(subSize, &subIndices)); 351 for (p = pStart, subOff = 0; p < pEnd; ++p) { 352 PetscInt gdof, gcdof, gtdof, d; 353 354 CHKERRQ(PetscSectionGetDof(sectionGlobals[i], p, &gdof)); 355 CHKERRQ(PetscSectionGetConstraintDof(sections[i], p, &gcdof)); 356 gtdof = gdof - gcdof; 357 if (gdof > 0 && gtdof) { 358 if (bs < 0) {bs = gtdof;} 359 else if (bs != gtdof) {bs = 1;} 360 CHKERRQ(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy)); 361 CHKERRQ(DMGetGlobalFieldOffset_Private(*superdm, p, startf+Nfs[i]-1, &dummy, &end)); 362 PetscCheck(end-start == gtdof,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid number of global dofs %D != %D for dm %D on point %D", end-start, gtdof, i, p); 363 for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d; 364 } 365 } 366 CHKERRQ(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i])); 367 /* Must have same blocksize on all procs (some might have no points) */ 368 { 369 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 370 371 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 372 CHKERRQ(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax)); 373 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 374 else {bs = bsMinMax[0];} 375 CHKERRQ(ISSetBlockSize((*is)[i], bs)); 376 } 377 } 378 } 379 /* Preserve discretizations */ 380 if (len && dms[0]->probs) { 381 CHKERRQ(DMSetNumFields(*superdm, Nf)); 382 for (i = 0, supf = 0; i < len; ++i) { 383 for (f = 0; f < Nfs[i]; ++f, ++supf) { 384 PetscObject disc; 385 386 CHKERRQ(DMGetField(dms[i], f, NULL, &disc)); 387 CHKERRQ(DMSetField(*superdm, supf, NULL, disc)); 388 } 389 } 390 CHKERRQ(DMCreateDS(*superdm)); 391 } 392 /* Preserve nullspaces */ 393 for (i = 0, supf = 0; i < len; ++i) { 394 for (f = 0; f < Nfs[i]; ++f, ++supf) { 395 (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f]; 396 if ((*superdm)->nullspaceConstructors[supf]) { 397 haveNull = PETSC_TRUE; 398 nullf = supf; 399 oldf = f; 400 } 401 } 402 } 403 /* Attach nullspace to IS */ 404 if (haveNull && is) { 405 MatNullSpace nullSpace; 406 407 CHKERRQ((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace)); 408 CHKERRQ(PetscObjectCompose((PetscObject) (*is)[nullf], "nullspace", (PetscObject) nullSpace)); 409 CHKERRQ(MatNullSpaceDestroy(&nullSpace)); 410 } 411 CHKERRQ(PetscSectionDestroy(&supersection)); 412 CHKERRQ(PetscFree3(Nfs, sections, sectionGlobals)); 413 PetscFunctionReturn(0); 414 } 415