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 = DMGetSection(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, const 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 = DMGetSection(dm, §ion);CHKERRQ(ierr); 95 ierr = DMGetGlobalSection(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, bsLocal[2], bsMinMax[2]; 102 103 for (f = 0, bs = 0; f < numFields; ++f) { 104 PetscInt Nc; 105 106 ierr = PetscSectionGetFieldComponents(section, fields[f], &Nc);CHKERRQ(ierr); 107 bs += Nc; 108 } 109 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 110 for (p = pStart; p < pEnd; ++p) { 111 PetscInt gdof, pSubSize = 0; 112 113 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 114 if (gdof > 0) { 115 for (f = 0; f < numFields; ++f) { 116 PetscInt fdof, fcdof; 117 118 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 119 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 120 pSubSize += fdof-fcdof; 121 } 122 subSize += pSubSize; 123 if (pSubSize && bs != pSubSize) { 124 /* Layout does not admit a pointwise block size */ 125 bs = 1; 126 } 127 } 128 } 129 /* Must have same blocksize on all procs (some might have no points) */ 130 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 131 ierr = PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);CHKERRQ(ierr); 132 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 133 else {bs = bsMinMax[0];} 134 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 135 for (p = pStart; p < pEnd; ++p) { 136 PetscInt gdof, goff; 137 138 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 139 if (gdof > 0) { 140 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); 141 for (f = 0; f < numFields; ++f) { 142 PetscInt fdof, fcdof, fc, f2, poff = 0; 143 144 /* Can get rid of this loop by storing field information in the global section */ 145 for (f2 = 0; f2 < fields[f]; ++f2) { 146 ierr = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr); 147 ierr = PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);CHKERRQ(ierr); 148 poff += fdof-fcdof; 149 } 150 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 151 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 152 for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) { 153 subIndices[subOff] = goff+poff+fc; 154 } 155 } 156 } 157 } 158 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 159 if (bs > 1) { 160 /* We need to check that the block size does not come from non-contiguous fields */ 161 PetscInt i, j, set = 1; 162 for (i = 0; i < subSize; i += bs) { 163 for (j = 0; j < bs; ++j) { 164 if (subIndices[i+j] != subIndices[i]+j) {set = 0; break;} 165 } 166 } 167 if (set) {ierr = ISSetBlockSize(*is, bs);CHKERRQ(ierr);} 168 } 169 } 170 if (subdm) { 171 PetscSection subsection; 172 PetscBool haveNull = PETSC_FALSE; 173 PetscInt f, nf = 0; 174 175 ierr = PetscSectionCreateSubsection(section, numFields, fields, &subsection);CHKERRQ(ierr); 176 ierr = DMSetSection(*subdm, subsection);CHKERRQ(ierr); 177 ierr = PetscSectionDestroy(&subsection);CHKERRQ(ierr); 178 for (f = 0; f < numFields; ++f) { 179 (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]]; 180 if ((*subdm)->nullspaceConstructors[f]) { 181 haveNull = PETSC_TRUE; 182 nf = f; 183 } 184 } 185 if (haveNull && is) { 186 MatNullSpace nullSpace; 187 188 ierr = (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);CHKERRQ(ierr); 189 ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 190 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 191 } 192 if (dm->probs) { 193 ierr = DMSetNumFields(*subdm, numFields);CHKERRQ(ierr); 194 for (f = 0; f < numFields; ++f) { 195 PetscObject disc; 196 197 ierr = DMGetField(dm, fields[f], NULL, &disc);CHKERRQ(ierr); 198 ierr = DMSetField(*subdm, f, NULL, disc);CHKERRQ(ierr); 199 } 200 ierr = DMCreateDS(*subdm);CHKERRQ(ierr); 201 if (numFields == 1 && is) { 202 PetscObject disc, space, pmat; 203 204 ierr = DMGetField(*subdm, 0, NULL, &disc);CHKERRQ(ierr); 205 ierr = PetscObjectQuery(disc, "nullspace", &space);CHKERRQ(ierr); 206 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nullspace", space);CHKERRQ(ierr);} 207 ierr = PetscObjectQuery(disc, "nearnullspace", &space);CHKERRQ(ierr); 208 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nearnullspace", space);CHKERRQ(ierr);} 209 ierr = PetscObjectQuery(disc, "pmat", &pmat);CHKERRQ(ierr); 210 if (pmat) {ierr = PetscObjectCompose((PetscObject) *is, "pmat", pmat);CHKERRQ(ierr);} 211 } 212 ierr = PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds);CHKERRQ(ierr); 213 ierr = PetscDSCopyBoundary(dm->probs[0].ds, (*subdm)->probs[0].ds);CHKERRQ(ierr); 214 ierr = PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);CHKERRQ(ierr); 215 } 216 if (dm->coarseMesh) { 217 ierr = DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh);CHKERRQ(ierr); 218 } 219 } 220 PetscFunctionReturn(0); 221 } 222 223 /* This assumes that the DM has been cloned prior to the call */ 224 PetscErrorCode DMCreateSuperDM_Section_Private(DM dms[], PetscInt len, IS **is, DM *superdm) 225 { 226 MPI_Comm comm; 227 PetscSection supersection, *sections, *sectionGlobals; 228 PetscInt *Nfs, Nf = 0, f, supf, nullf = -1, i; 229 PetscBool haveNull = PETSC_FALSE; 230 PetscErrorCode ierr; 231 232 PetscFunctionBegin; 233 ierr = PetscObjectGetComm((PetscObject)dms[0], &comm);CHKERRQ(ierr); 234 /* Pull out local and global sections */ 235 ierr = PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals);CHKERRQ(ierr); 236 for (i = 0 ; i < len; ++i) { 237 ierr = DMGetSection(dms[i], §ions[i]);CHKERRQ(ierr); 238 ierr = DMGetGlobalSection(dms[i], §ionGlobals[i]);CHKERRQ(ierr); 239 if (!sections[i]) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 240 if (!sectionGlobals[i]) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 241 ierr = PetscSectionGetNumFields(sections[i], &Nfs[i]);CHKERRQ(ierr); 242 Nf += Nfs[i]; 243 } 244 /* Create the supersection */ 245 ierr = PetscSectionCreateSupersection(sections, len, &supersection);CHKERRQ(ierr); 246 ierr = DMSetSection(*superdm, supersection);CHKERRQ(ierr); 247 /* Create ISes */ 248 if (is) { 249 PetscSection supersectionGlobal; 250 PetscInt bs = -1, startf = 0; 251 252 ierr = PetscMalloc1(len, is);CHKERRQ(ierr); 253 ierr = DMGetDefaultGlobalSection(*superdm, &supersectionGlobal);CHKERRQ(ierr); 254 for (i = 0 ; i < len; startf += Nfs[i], ++i) { 255 PetscInt *subIndices; 256 PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy; 257 258 ierr = PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);CHKERRQ(ierr); 259 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);CHKERRQ(ierr); 260 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 261 for (p = pStart, subOff = 0; p < pEnd; ++p) { 262 PetscInt gdof, gcdof, gtdof, d; 263 264 ierr = PetscSectionGetDof(sectionGlobals[i], p, &gdof);CHKERRQ(ierr); 265 ierr = PetscSectionGetConstraintDof(sections[i], p, &gcdof);CHKERRQ(ierr); 266 gtdof = gdof - gcdof; 267 if (gdof > 0 && gtdof) { 268 if (bs < 0) {bs = gtdof;} 269 else if (bs != gtdof) {bs = 1;} 270 ierr = DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy);CHKERRQ(ierr); 271 ierr = DMGetGlobalFieldOffset_Private(*superdm, p, startf+Nfs[i]-1, &dummy, &end);CHKERRQ(ierr); 272 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); 273 for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d; 274 } 275 } 276 ierr = ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);CHKERRQ(ierr); 277 /* Must have same blocksize on all procs (some might have no points) */ 278 { 279 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 280 281 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 282 ierr = PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax);CHKERRQ(ierr); 283 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 284 else {bs = bsMinMax[0];} 285 ierr = ISSetBlockSize((*is)[i], bs);CHKERRQ(ierr); 286 } 287 } 288 } 289 /* Preserve discretizations */ 290 if (len && dms[0]->probs) { 291 ierr = DMSetNumFields(*superdm, Nf);CHKERRQ(ierr); 292 for (i = 0, supf = 0; i < len; ++i) { 293 for (f = 0; f < Nfs[i]; ++f, ++supf) { 294 PetscObject disc; 295 296 ierr = DMGetField(dms[i], f, NULL, &disc);CHKERRQ(ierr); 297 ierr = DMSetField(*superdm, supf, NULL, disc);CHKERRQ(ierr); 298 } 299 } 300 ierr = DMCreateDS(*superdm);CHKERRQ(ierr); 301 } 302 /* Preserve nullspaces */ 303 for (i = 0, supf = 0; i < len; ++i) { 304 for (f = 0; f < Nfs[i]; ++f, ++supf) { 305 (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f]; 306 if ((*superdm)->nullspaceConstructors[supf]) { 307 haveNull = PETSC_TRUE; 308 nullf = supf; 309 } 310 } 311 } 312 /* Attach nullspace to IS */ 313 if (haveNull && is) { 314 MatNullSpace nullSpace; 315 316 ierr = (*(*superdm)->nullspaceConstructors[nullf])(*superdm, nullf, &nullSpace);CHKERRQ(ierr); 317 ierr = PetscObjectCompose((PetscObject) (*is)[nullf], "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 318 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 319 } 320 ierr = PetscSectionDestroy(&supersection);CHKERRQ(ierr); 321 ierr = PetscFree3(Nfs, sections, sectionGlobals);CHKERRQ(ierr); 322 PetscFunctionReturn(0); 323 } 324