63 RCP<ParameterList> validParamList = rcp(
new ParameterList());
65#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
68 validParamList->set< std::string >(
"Fine level nullspace",
"Nullspace",
"Variable name which is used to store null space multi vector on the finest level (default=\"Nullspace\"). For block matrices also \"Nullspace1\" to \"Nullspace9\" are accepted to describe the null space vectors for the (i,i) block (i=1..9).");
70 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the fine level matrix (only needed if default null space is generated)");
71 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the fine level null space");
72 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory of the coordinates");
76 validParamList->set< RCP<const FactoryBase> >(
"Nullspace1", Teuchos::null,
"Generating factory of the fine level null space associated with the (1,1) block in your n x n block matrix.");
77 validParamList->set< RCP<const FactoryBase> >(
"Nullspace2", Teuchos::null,
"Generating factory of the fine level null space associated with the (2,2) block in your n x n block matrix.");
78 validParamList->set< RCP<const FactoryBase> >(
"Nullspace3", Teuchos::null,
"Generating factory of the fine level null space associated with the (3,3) block in your n x n block matrix.");
79 validParamList->set< RCP<const FactoryBase> >(
"Nullspace4", Teuchos::null,
"Generating factory of the fine level null space associated with the (4,4) block in your n x n block matrix.");
80 validParamList->set< RCP<const FactoryBase> >(
"Nullspace5", Teuchos::null,
"Generating factory of the fine level null space associated with the (5,5) block in your n x n block matrix.");
81 validParamList->set< RCP<const FactoryBase> >(
"Nullspace6", Teuchos::null,
"Generating factory of the fine level null space associated with the (6,6) block in your n x n block matrix.");
82 validParamList->set< RCP<const FactoryBase> >(
"Nullspace7", Teuchos::null,
"Generating factory of the fine level null space associated with the (7,7) block in your n x n block matrix.");
83 validParamList->set< RCP<const FactoryBase> >(
"Nullspace8", Teuchos::null,
"Generating factory of the fine level null space associated with the (8,8) block in your n x n block matrix.");
84 validParamList->set< RCP<const FactoryBase> >(
"Nullspace9", Teuchos::null,
"Generating factory of the fine level null space associated with the (9,9) block in your n x n block matrix.");
86 return validParamList;
123 RCP<MultiVector> nullspace;
126 const ParameterList & pL = GetParameterList();
127 std::string nspName = pL.get<std::string>(
"Fine level nullspace");
131 typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType coordinate_type;
132 typedef Xpetra::MultiVector<coordinate_type,LO,GO,NO> RealValuedMultiVector;
133 RCP<RealValuedMultiVector> Coords;
134 ArrayRCP<const coordinate_type> xvals, yvals, zvals;
137 cx = 0.; cy = 0.; cz = 0.;
138 if(calculateRotations_) {
139 Coords = Get< RCP<RealValuedMultiVector> >(currentLevel,
"Coordinates");
141 cx = Coords->getVector(0)->meanValue();
142 if (Coords->getNumVectors() > 1)
143 cy = Coords->getVector(1)->meanValue();
144 if (Coords->getNumVectors() > 2)
145 cz = Coords->getVector(2)->meanValue();
147 xvals = Coords->getData(0);
148 if (Coords->getNumVectors() > 1)
149 yvals = Coords->getData(1);
150 if (Coords->getNumVectors() > 2)
151 zvals = Coords->getData(2);
160 GetOStream(
Runtime1) <<
"Use user-given nullspace " << nspName <<
": nullspace dimension=" << nullspace->getNumVectors() <<
" nullspace length=" << nullspace->getGlobalLength() << std::endl;
163 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
167 if(A->IsView(
"stridedMaps")==
true) {
168 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
169 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap()) == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
170 numPDEs = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap())->getFixedBlockSize();
171 oldView = A->SwitchToView(oldView);
174 LO nullspaceDim = numPDEs;
176 if(calculateRotations_) {
177 if (Coords->getNumVectors() > 1) nullspaceDim++;
178 if (Coords->getNumVectors() > 2) nullspaceDim += 2;
179 GetOStream(
Runtime1) <<
"Generating nullspace with rotations: dimension = " << nullspaceDim << std::endl;
181 else GetOStream(
Runtime1) <<
"Generating canonical nullspace: dimension = " << numPDEs << std::endl;
183 nullspace = MultiVectorFactory::Build(A->getDomainMap(), nullspaceDim);
185 RCP<BlockedMultiVector> bnsp = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(nullspace);
186 if(bnsp.is_null() ==
true) {
187 for (
int i=0; i<numPDEs; ++i) {
188 ArrayRCP<Scalar> nsValues = nullspace->getDataNonConst(i);
189 int numBlocks = nsValues.size() / numPDEs;
190 for (
int j=0; j< numBlocks; ++j) {
191 nsValues[j*numPDEs + i] = 1.0;
194 if (( (
int) nullspaceDim > numPDEs ) && ((
int) numPDEs > 1)) {
196 ArrayRCP<Scalar> nsValues = nullspace->getDataNonConst(numPDEs);
197 int numBlocks = nsValues.size() / numPDEs;
198 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks != (
int) xvals.size(),
Exceptions::RuntimeError,
"MueLu::NullspaceFactory::Build(): number of coordinates does not match ndofs/numPDEs.");
199 for (
int j=0; j< numBlocks; ++j) {
200 nsValues[j*numPDEs + 0] = -(yvals[j]-cy);
201 nsValues[j*numPDEs + 1] = (xvals[j]-cx);
204 if (( (
int) nullspaceDim == numPDEs+3 ) && ((
int) numPDEs > 2)) {
206 ArrayRCP<Scalar> nsValues = nullspace->getDataNonConst(numPDEs+1);
207 int numBlocks = nsValues.size() / numPDEs;
208 for (
int j=0; j< numBlocks; ++j) {
209 nsValues[j*numPDEs + 1] = -(zvals[j]-cz);
210 nsValues[j*numPDEs + 2] = (yvals[j]-cy);
213 nsValues = nullspace->getDataNonConst(numPDEs+2);
214 numBlocks = nsValues.size() / numPDEs;
215 for (
int j=0; j< numBlocks; ++j) {
216 nsValues[j*numPDEs + 0] = -(zvals[j]-cz);
217 nsValues[j*numPDEs + 2] = (xvals[j]-cx);
234 fillNullspaceVector(bnsp,numPDEs, xvals, yvals, zvals, nullspaceDim, cx, cy, cz);
241 nullspace = currentLevel.
Get< RCP<MultiVector> >(
"Nullspace", GetFactory(nspName).get());
246 Set(currentLevel,
"Nullspace", nullspace);
251 void NullspaceFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::fillNullspaceVector(
const RCP<BlockedMultiVector>& nsp,
LocalOrdinal numPDEs, ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xvals, ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yvals, ArrayRCP<
const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zvals,
LocalOrdinal nullspaceDim,
Scalar cx,
Scalar cy,
Scalar cz)
const {
252 RCP< const BlockedMap> bmap = nsp->getBlockedMap();
254 for(
size_t r = 0; r < bmap->getNumMaps(); r++) {
255 Teuchos::RCP<MultiVector> part = nsp->getMultiVector(r);
256 Teuchos::RCP<BlockedMultiVector> bpart = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(part);
257 if(bpart.is_null() ==
true) {
258 for (
int i=0; i<numPDEs; ++i) {
259 ArrayRCP<Scalar> nsValues = part->getDataNonConst(i);
260 int numBlocks = nsValues.size() / numPDEs;
261 for (
int j=0; j< numBlocks; ++j) {
262 nsValues[j*numPDEs + i] = 1.0;
265 if ( (
int) nullspaceDim > numPDEs ) {
267 ArrayRCP<Scalar> nsValues = part->getDataNonConst(numPDEs);
268 int numBlocks = nsValues.size() / numPDEs;
269 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks != (
int) xvals.size(),
Exceptions::RuntimeError,
"MueLu::NullspaceFactory::fillNullspaceVector(): number of coordinates does not match ndofs/numPDEs.");
270 for (
int j=0; j< numBlocks; ++j) {
271 nsValues[j*numPDEs + 0] = -(yvals[j]-cy);
272 nsValues[j*numPDEs + 1] = (xvals[j]-cx);
275 if ( (
int) nullspaceDim == numPDEs+3 ) {
277 ArrayRCP<Scalar> nsValues = part->getDataNonConst(numPDEs+1);
278 int numBlocks = nsValues.size() / numPDEs;
279 for (
int j=0; j< numBlocks; ++j) {
280 nsValues[j*numPDEs + 1] = -(zvals[j]-cz);
281 nsValues[j*numPDEs + 2] = (yvals[j]-cy);
284 nsValues = part->getDataNonConst(numPDEs+2);
285 numBlocks = nsValues.size() / numPDEs;
286 for (
int j=0; j< numBlocks; ++j) {
287 nsValues[j*numPDEs + 0] = -(zvals[j]-cz);
288 nsValues[j*numPDEs + 2] = (xvals[j]-cx);
306 fillNullspaceVector(bpart,numPDEs,xvals, yvals, zvals, nullspaceDim, cx, cy, cz);
void fillNullspaceVector(const RCP< BlockedMultiVector > &nsp, LocalOrdinal numPDEs, ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > xvals, ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > yvals, ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > zvals, LocalOrdinal nullspaceDim, Scalar cx, Scalar cy, Scalar cz) const
Helper function to recursively fill BlockedMultiVector with default null space vectors.