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