46#ifndef MUELU_NULLSPACEFACTORY_DEF_HPP
47#define MUELU_NULLSPACEFACTORY_DEF_HPP
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_MultiVectorFactory.hpp>
51#include <Xpetra_BlockedMultiVector.hpp>
52#include <Xpetra_VectorFactory.hpp>
61 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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;
89 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 const ParameterList & pL = GetParameterList();
93 std::string nspName = pL.get<std::string>(
"Fine level nullspace");
99 Input(currentLevel,
"A");
103 pL.get<
bool>(
"nullspace: calculate rotations") ) {
104 calculateRotations_ =
true;
105 Input(currentLevel,
"Coordinates");
114 TEUCHOS_TEST_FOR_EXCEPTION(GetFactory(nspName)==Teuchos::null,
Exceptions::RuntimeError,
"MueLu::NullspaceFactory::DeclareInput(): You must declare an existing factory which produces the variable \"Nullspace\" in the NullspaceFactory (e.g. a TentativePFactory).");
115 currentLevel.
DeclareInput(
"Nullspace", GetFactory(nspName).get(),
this);
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
250 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
#define SET_VALID_ENTRY(name)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
void Build(Level ¤tLevel) const
Build an object with this factory.
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.
RCP< const ParameterList > GetValidParameterList() const
Define valid parameters for internal factory parameters.
void DeclareInput(Level ¤tLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Namespace for MueLu classes and methods.
@ Runtime1
Description of what is happening (more verbose)