1 #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMCreateGlobalVector_Section_Private" 5 PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm,Vec *vec) 6 { 7 PetscSection gSection; 8 PetscInt localSize, bs, blockSize = -1, pStart, pEnd, p; 9 PetscErrorCode ierr; 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 if ((blockSize < 0) && (dof > 0) && (dof-cdof > 0)) blockSize = dof-cdof; 20 if ((dof > 0) && (dof-cdof != blockSize)) { 21 blockSize = 1; 22 break; 23 } 24 } 25 if (blockSize < 0) blockSize = 1; 26 ierr = MPI_Allreduce(&blockSize, &bs, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 27 ierr = PetscSectionGetConstrainedStorageSize(gSection, &localSize);CHKERRQ(ierr); 28 if (localSize%blockSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %d and local storage size %d", blockSize, localSize); 29 ierr = VecCreate(PetscObjectComm((PetscObject)dm), vec);CHKERRQ(ierr); 30 ierr = VecSetSizes(*vec, localSize, PETSC_DETERMINE);CHKERRQ(ierr); 31 ierr = VecSetBlockSize(*vec, bs);CHKERRQ(ierr); 32 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr); 33 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr); 34 /* ierr = VecSetLocalToGlobalMapping(*vec, dm->ltogmap);CHKERRQ(ierr); */ 35 /* ierr = VecSetLocalToGlobalMappingBlock(*vec, dm->ltogmapb);CHKERRQ(ierr); */ 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "DMCreateLocalVector_Section_Private" 41 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec) 42 { 43 PetscSection section; 44 PetscInt localSize, blockSize = -1, pStart, pEnd, p; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 49 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 50 for (p = pStart; p < pEnd; ++p) { 51 PetscInt dof; 52 53 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 54 if ((blockSize < 0) && (dof > 0)) blockSize = dof; 55 if ((dof > 0) && (dof != blockSize)) { 56 blockSize = 1; 57 break; 58 } 59 } 60 ierr = PetscSectionGetStorageSize(section, &localSize);CHKERRQ(ierr); 61 ierr = VecCreate(PETSC_COMM_SELF, vec);CHKERRQ(ierr); 62 ierr = VecSetSizes(*vec, localSize, localSize);CHKERRQ(ierr); 63 ierr = VecSetBlockSize(*vec, blockSize);CHKERRQ(ierr); 64 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr); 65 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "DMCreateSubDM_Section_Private" 71 /* This assumes that the DM has been cloned prior to the call */ 72 PetscErrorCode DMCreateSubDM_Section_Private(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 73 { 74 PetscSection section, sectionGlobal; 75 PetscInt *subIndices; 76 PetscInt subSize = 0, subOff = 0, nF, f, pStart, pEnd, p; 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 if (!numFields) PetscFunctionReturn(0); 81 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 82 ierr = DMGetDefaultGlobalSection(dm, §ionGlobal);CHKERRQ(ierr); 83 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 84 if (!sectionGlobal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 85 ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); 86 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); 87 if (is) { 88 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 89 for (p = pStart; p < pEnd; ++p) { 90 PetscInt gdof; 91 92 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 93 if (gdof > 0) { 94 for (f = 0; f < numFields; ++f) { 95 PetscInt fdof, fcdof; 96 97 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 98 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 99 subSize += fdof-fcdof; 100 } 101 } 102 } 103 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 104 for (p = pStart; p < pEnd; ++p) { 105 PetscInt gdof, goff; 106 107 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 108 if (gdof > 0) { 109 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); 110 for (f = 0; f < numFields; ++f) { 111 PetscInt fdof, fcdof, fc, f2, poff = 0; 112 113 /* Can get rid of this loop by storing field information in the global section */ 114 for (f2 = 0; f2 < fields[f]; ++f2) { 115 ierr = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr); 116 ierr = PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);CHKERRQ(ierr); 117 poff += fdof-fcdof; 118 } 119 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 120 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 121 for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) { 122 subIndices[subOff] = goff+poff+fc; 123 } 124 } 125 } 126 } 127 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 128 } 129 if (subdm) { 130 PetscSection subsection; 131 PetscBool haveNull = PETSC_FALSE; 132 PetscInt f, nf = 0; 133 134 ierr = PetscSectionCreateSubsection(section, numFields, fields, &subsection);CHKERRQ(ierr); 135 ierr = DMSetDefaultSection(*subdm, subsection);CHKERRQ(ierr); 136 ierr = PetscSectionDestroy(&subsection);CHKERRQ(ierr); 137 for (f = 0; f < numFields; ++f) { 138 (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]]; 139 if ((*subdm)->nullspaceConstructors[f]) { 140 haveNull = PETSC_TRUE; 141 nf = f; 142 } 143 } 144 if (haveNull) { 145 MatNullSpace nullSpace; 146 147 ierr = (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);CHKERRQ(ierr); 148 ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 149 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 150 } 151 if (dm->fields) { 152 if (nF != dm->numFields) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "The number of DM fields %d does not match the number of Section fields %d", dm->numFields, nF); 153 ierr = DMSetNumFields(*subdm, numFields);CHKERRQ(ierr); 154 for (f = 0; f < numFields; ++f) { 155 ierr = PetscObjectListDuplicate(dm->fields[fields[f]]->olist, &(*subdm)->fields[f]->olist);CHKERRQ(ierr); 156 } 157 if (numFields == 1) { 158 MatNullSpace space; 159 Mat pmat; 160 161 ierr = PetscObjectQuery((*subdm)->fields[0], "nullspace", (PetscObject*) &space);CHKERRQ(ierr); 162 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) space);CHKERRQ(ierr);} 163 ierr = PetscObjectQuery((*subdm)->fields[0], "nearnullspace", (PetscObject*) &space);CHKERRQ(ierr); 164 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nearnullspace", (PetscObject) space);CHKERRQ(ierr);} 165 ierr = PetscObjectQuery((*subdm)->fields[0], "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 166 if (pmat) {ierr = PetscObjectCompose((PetscObject) *is, "pmat", (PetscObject) pmat);CHKERRQ(ierr);} 167 } 168 } 169 } 170 PetscFunctionReturn(0); 171 } 172