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 = DMGetDefaultGlobalSection(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 = DMGetDefaultSection(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 /* This assumes that the DM has been cloned prior to the call */ 85 PetscErrorCode DMCreateSubDM_Section_Private(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 86 { 87 PetscSection section, sectionGlobal; 88 PetscInt *subIndices; 89 PetscInt subSize = 0, subOff = 0, nF, f, pStart, pEnd, p; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 if (!numFields) PetscFunctionReturn(0); 94 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 95 ierr = DMGetDefaultGlobalSection(dm, §ionGlobal);CHKERRQ(ierr); 96 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 97 if (!sectionGlobal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 98 ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); 99 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); 100 if (is) { 101 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 102 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 103 for (p = pStart; p < pEnd; ++p) { 104 PetscInt gdof, pSubSize = 0; 105 106 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 107 if (gdof > 0) { 108 for (f = 0; f < numFields; ++f) { 109 PetscInt fdof, fcdof; 110 111 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 112 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 113 pSubSize += fdof-fcdof; 114 } 115 subSize += pSubSize; 116 if (pSubSize) { 117 if (bs < 0) { 118 bs = pSubSize; 119 } else if (bs != pSubSize) { 120 /* Layout does not admit a pointwise block size */ 121 bs = 1; 122 } 123 } 124 } 125 } 126 /* Must have same blocksize on all procs (some might have no points) */ 127 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 128 ierr = PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);CHKERRQ(ierr); 129 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 130 else {bs = bsMinMax[0];} 131 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 132 for (p = pStart; p < pEnd; ++p) { 133 PetscInt gdof, goff; 134 135 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 136 if (gdof > 0) { 137 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); 138 for (f = 0; f < numFields; ++f) { 139 PetscInt fdof, fcdof, fc, f2, poff = 0; 140 141 /* Can get rid of this loop by storing field information in the global section */ 142 for (f2 = 0; f2 < fields[f]; ++f2) { 143 ierr = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr); 144 ierr = PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);CHKERRQ(ierr); 145 poff += fdof-fcdof; 146 } 147 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 148 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 149 for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) { 150 subIndices[subOff] = goff+poff+fc; 151 } 152 } 153 } 154 } 155 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 156 ierr = ISSetBlockSize(*is, bs);CHKERRQ(ierr); 157 } 158 if (subdm) { 159 PetscSection subsection; 160 PetscBool haveNull = PETSC_FALSE; 161 PetscInt f, nf = 0; 162 163 ierr = PetscSectionCreateSubsection(section, numFields, fields, &subsection);CHKERRQ(ierr); 164 ierr = DMSetDefaultSection(*subdm, subsection);CHKERRQ(ierr); 165 ierr = PetscSectionDestroy(&subsection);CHKERRQ(ierr); 166 for (f = 0; f < numFields; ++f) { 167 (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]]; 168 if ((*subdm)->nullspaceConstructors[f]) { 169 haveNull = PETSC_TRUE; 170 nf = f; 171 } 172 } 173 if (haveNull && is) { 174 MatNullSpace nullSpace; 175 176 ierr = (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);CHKERRQ(ierr); 177 ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 178 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 179 } 180 if (dm->prob) { 181 PetscInt Nf; 182 183 ierr = PetscDSGetNumFields(dm->prob, &Nf);CHKERRQ(ierr); 184 if (nF != Nf) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "The number of DM fields %d does not match the number of Section fields %d", Nf, nF); 185 ierr = DMSetNumFields(*subdm, numFields);CHKERRQ(ierr); 186 for (f = 0; f < numFields; ++f) { 187 PetscObject disc; 188 189 ierr = DMGetField(dm, fields[f], &disc);CHKERRQ(ierr); 190 ierr = DMSetField(*subdm, f, disc);CHKERRQ(ierr); 191 } 192 if (numFields == 1 && is) { 193 PetscObject disc, space, pmat; 194 195 ierr = DMGetField(*subdm, 0, &disc);CHKERRQ(ierr); 196 ierr = PetscObjectQuery(disc, "nullspace", &space);CHKERRQ(ierr); 197 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nullspace", space);CHKERRQ(ierr);} 198 ierr = PetscObjectQuery(disc, "nearnullspace", &space);CHKERRQ(ierr); 199 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nearnullspace", space);CHKERRQ(ierr);} 200 ierr = PetscObjectQuery(disc, "pmat", &pmat);CHKERRQ(ierr); 201 if (pmat) {ierr = PetscObjectCompose((PetscObject) *is, "pmat", pmat);CHKERRQ(ierr);} 202 } 203 } 204 } 205 #if 0 206 /* We need a way to filter the original SF for given fields: 207 - Keeping the original section around it too much I think 208 - We could keep the distributed section, and subset it 209 */ 210 if (dm->sfNatural) { 211 PetscSF sfNatural; 212 213 ierr = PetscSectionCreateSubsection(dm->originalSection, numFields, fields, &(*subdm)->originalSection);CHKERRQ(ierr); 214 ierr = DMPlexCreateGlobalToNaturalPetscSF(*subdm, &sfNatural);CHKERRQ(ierr); 215 ierr = DMPlexSetGlobalToNaturalPetscSF(*subdm, sfNatural);CHKERRQ(ierr); 216 } 217 #endif 218 PetscFunctionReturn(0); 219 } 220 221 /* This assumes that the DM has been cloned prior to the call */ 222 PetscErrorCode DMCreateSuperDM_Section_Private(DM dms[], PetscInt len, IS **is, DM *superdm) 223 { 224 PetscSection *sections, *sectionGlobals; 225 PetscInt *Nfs, Nf = 0, *subIndices, i; 226 PetscErrorCode ierr; 227 228 PetscFunctionBegin; 229 ierr = PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals);CHKERRQ(ierr); 230 for (i = 0 ; i < len; ++i) { 231 ierr = DMGetDefaultSection(dms[i], §ions[i]);CHKERRQ(ierr); 232 ierr = DMGetDefaultGlobalSection(dms[i], §ionGlobals[i]);CHKERRQ(ierr); 233 if (!sections[i]) SETERRQ(PetscObjectComm((PetscObject)dms[0]), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 234 if (!sectionGlobals[i]) SETERRQ(PetscObjectComm((PetscObject)dms[0]), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 235 ierr = PetscSectionGetNumFields(sections[i], &Nfs[i]);CHKERRQ(ierr); 236 Nf += Nfs[i]; 237 } 238 if (is) { 239 PetscInt *offs, *globalOffs, iOff = 0; 240 241 ierr = PetscMalloc1(len, is);CHKERRQ(ierr); 242 ierr = PetscCalloc2(len+1, &offs, len+1, &globalOffs);CHKERRQ(ierr); 243 for (i = 0 ; i < len; ++i) { 244 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &offs[i]);CHKERRQ(ierr); 245 offs[len] += offs[i]; 246 } 247 ierr = MPI_Scan(offs, globalOffs, len+1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dms[0]));CHKERRQ(ierr); 248 for (i = 0 ; i <= len; ++i) globalOffs[i] -= offs[i]; 249 for (i = 0 ; i < len; ++i, iOff += offs[i-1]) { 250 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 251 PetscInt subSize = 0, subOff = 0, gtoff = globalOffs[len] - globalOffs[i], pStart, pEnd, p; 252 253 ierr = PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);CHKERRQ(ierr); 254 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);CHKERRQ(ierr); 255 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 256 for (p = pStart; p < pEnd; ++p) { 257 PetscInt gdof, gcdof, gtdof, goff, d; 258 259 ierr = PetscSectionGetDof(sectionGlobals[i], p, &gdof);CHKERRQ(ierr); 260 ierr = PetscSectionGetConstraintDof(sections[i], p, &gcdof);CHKERRQ(ierr); 261 gtdof = gdof-gcdof; 262 if (gdof > 0 && gtdof) { 263 if (bs < 0) {bs = gtdof;} 264 else if (bs != gtdof) {bs = 1;} 265 ierr = PetscSectionGetOffset(sectionGlobals[i], p, &goff);CHKERRQ(ierr); 266 for (d = 0; d < gtdof; ++d, ++subOff) { 267 subIndices[subOff] = goff+gtoff+d+iOff; 268 } 269 } 270 } 271 ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dms[0]), subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);CHKERRQ(ierr); 272 /* Must have same blocksize on all procs (some might have no points) */ 273 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 274 ierr = PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dms[0]), bsLocal, bsMinMax);CHKERRQ(ierr); 275 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 276 else {bs = bsMinMax[0];} 277 ierr = ISSetBlockSize((*is)[i], bs);CHKERRQ(ierr); 278 } 279 ierr = PetscFree2(offs, globalOffs);CHKERRQ(ierr); 280 } 281 if (superdm) { 282 PetscSection supersection; 283 PetscBool haveNull = PETSC_FALSE; 284 PetscInt field, f, nf = 0; 285 286 ierr = PetscSectionCreateSupersection(sections, len, &supersection);CHKERRQ(ierr); 287 ierr = DMSetDefaultSection(*superdm, supersection);CHKERRQ(ierr); 288 ierr = PetscSectionDestroy(&supersection);CHKERRQ(ierr); 289 for (i = 0, field = 0; i < len; ++i) { 290 for (f = 0; f < Nfs[i]; ++f, ++field) { 291 (*superdm)->nullspaceConstructors[field] = dms[i]->nullspaceConstructors[f]; 292 if ((*superdm)->nullspaceConstructors[field]) { 293 haveNull = PETSC_TRUE; 294 nf = field; 295 } 296 } 297 } 298 if (haveNull && is) { 299 MatNullSpace nullSpace; 300 301 ierr = (*(*superdm)->nullspaceConstructors[nf])(*superdm, nf, &nullSpace);CHKERRQ(ierr); 302 ierr = PetscObjectCompose((PetscObject) (*is)[nf], "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 303 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 304 } 305 if (len && dms[0]->prob) { 306 ierr = DMSetNumFields(*superdm, Nf);CHKERRQ(ierr); 307 for (i = 0, field = 0; i < len; ++i) { 308 for (f = 0; f < Nfs[i]; ++f, ++field) { 309 PetscObject disc; 310 311 ierr = DMGetField(dms[i], f, &disc);CHKERRQ(ierr); 312 ierr = DMSetField(*superdm, field, disc);CHKERRQ(ierr); 313 } 314 } 315 } 316 } 317 #if 0 318 /* We need a way to filter the original SF for given fields: 319 - Keeping the original section around is too much I think 320 - We could keep the distributed section, and subset it 321 */ 322 if (len && dms[0]->sfNatural) { 323 PetscSF sfNatural; 324 325 ierr = PetscSectionCreateSupersection(originalSections, len, &(*superdm)->originalSection);CHKERRQ(ierr); 326 ierr = DMPlexCreateGlobalToNaturalPetscSF(*superdm, &sfNatural);CHKERRQ(ierr); 327 ierr = DMPlexSetGlobalToNaturalPetscSF(*superdm, sfNatural);CHKERRQ(ierr); 328 } 329 #endif 330 ierr = PetscFree3(Nfs, sections, sectionGlobals);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333