46 #ifndef MUELU_REFMAXWELL_DECL_HPP 47 #define MUELU_REFMAXWELL_DECL_HPP 52 #include "MueLu_Utilities.hpp" 53 #include "MueLu_SaPFactory.hpp" 54 #include "MueLu_TentativePFactory.hpp" 55 #include "MueLu_SmootherFactory.hpp" 56 #include "MueLu_CoalesceDropFactory.hpp" 57 #include "MueLu_UncoupledAggregationFactory.hpp" 58 #include "MueLu_TrilinosSmoother.hpp" 60 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2) 62 #include "Tpetra_Operator.hpp" 63 #include "Tpetra_CrsMatrix.hpp" 66 #include "Xpetra_Matrix.hpp" 67 #include "Xpetra_MatrixFactory.hpp" 68 #include "Xpetra_CrsMatrixWrap.hpp" 69 #include "Xpetra_BlockedCrsMatrix.hpp" 70 #include "Xpetra_TpetraMultiVector.hpp" 71 #include "Xpetra_ExportFactory.hpp" 87 template <
class Scalar =
95 class RefMaxwell :
public Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
97 #undef MUELU_REFMAXWELL_SHORT 106 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>
XMap;
107 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
XMV;
108 typedef Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
XTMV;
109 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
XCRS;
110 typedef Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
XTCRS;
111 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
XMat;
112 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>
XCrsWrap;
124 RefMaxwell(Teuchos::RCP<Hierarchy> H11, Teuchos::RCP<Hierarchy> H22) :
144 const Teuchos::RCP<TCRS> & D0_Matrix,
145 const Teuchos::RCP<TCRS> & M0inv_Matrix,
146 const Teuchos::RCP<TCRS> & M1_Matrix,
147 const Teuchos::RCP<TMV> & Nullspace,
148 const Teuchos::RCP<TMV> & Coords,
149 Teuchos::ParameterList& List,
150 bool ComputePrec =
true)
152 initialize(D0_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
171 const Teuchos::RCP<TCRS> & M0inv_Matrix,
172 const Teuchos::RCP<TCRS> & M1_Matrix,
173 const Teuchos::RCP<TMV> & Nullspace,
174 const Teuchos::RCP<TMV> & Coords,
175 Teuchos::ParameterList& List)
177 initialize(D0_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
191 const Teuchos::RCP<TCRS> & D0_Matrix,
192 const Teuchos::RCP<TCRS> & M1_Matrix,
193 const Teuchos::RCP<TMV> & Nullspace,
194 const Teuchos::RCP<TMV> & Coords,
195 Teuchos::ParameterList& List,
196 bool ComputePrec =
true)
198 initialize(D0_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
216 const Teuchos::RCP<TCRS> & M1_Matrix,
217 const Teuchos::RCP<TMV> & Nullspace,
218 const Teuchos::RCP<TMV> & Coords,
219 Teuchos::ParameterList& List)
221 initialize(D0_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
228 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
getDomainMap()
const;
231 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
getRangeMap()
const;
246 void resetMatrix(Teuchos::RCP<TCRS> SM_Matrix_new);
262 Teuchos::ETransp mode = Teuchos::NO_TRANS,
263 Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
264 Scalar beta = Teuchos::ScalarTraits<Scalar>::one())
const;
269 template <
class NewNode>
270 Teuchos::RCP< RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, NewNode> >
271 clone (
const RCP<NewNode>& new_node)
const {
280 std::vector<LocalOrdinal>& dirichletRows) {
281 dirichletRows.resize(0);
282 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
283 Teuchos::ArrayView<const LocalOrdinal> indices;
284 Teuchos::ArrayView<const Scalar> values;
285 A->getLocalRowView(i,indices,values);
287 for (
int j=0; j<indices.size(); j++) {
293 if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > 1.0e-16) {
297 if (nnz == 1 || nnz == 2) {
298 dirichletRows.push_back(i);
304 std::vector<LocalOrdinal>& dirichletRows,
305 std::vector<LocalOrdinal>& dirichletCols) {
306 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A->getDomainMap();
307 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A->getColMap();
308 Teuchos::RCP< Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter
309 = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
310 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> ::Build(colMap,1);
311 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> ::Build(domMap,1);
312 myColsToZero->putScalar((Scalar)0.0);
313 globalColsToZero->putScalar((Scalar)0.0);
314 for(
size_t i=0; i<dirichletRows.size(); i++) {
315 Teuchos::ArrayView<const LocalOrdinal> indices;
316 Teuchos::ArrayView<const Scalar> values;
317 A->getLocalRowView(dirichletRows[i],indices,values);
318 for(
int j=0; j<indices.size(); j++)
319 myColsToZero->replaceLocalValue(indices[j],0,(Scalar)1.0);
321 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
322 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
323 Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
324 dirichletCols.resize(colMap->getNodeNumElements());
325 for(
size_t i=0; i<colMap->getNodeNumElements(); i++) {
326 if(Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>0.0)
334 std::vector<LocalOrdinal>& dirichletRows) {
335 for(
size_t i=0; i<dirichletRows.size(); i++) {
336 Teuchos::ArrayView<const LocalOrdinal> indices;
337 Teuchos::ArrayView<const Scalar> values;
338 A->getLocalRowView(dirichletRows[i],indices,values);
339 std::vector<Scalar> vec;
340 vec.resize(indices.size());
341 Teuchos::ArrayView<Scalar> zerovalues(vec);
342 for(
int j=0; j<indices.size(); j++)
343 zerovalues[j]=(Scalar)1.0e-32;
344 A->replaceLocalValues(dirichletRows[i],indices,zerovalues);
349 std::vector<LocalOrdinal>& dirichletCols) {
350 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
351 Teuchos::ArrayView<const LocalOrdinal> indices;
352 Teuchos::ArrayView<const Scalar> values;
353 A->getLocalRowView(i,indices,values);
354 std::vector<Scalar> vec;
355 vec.resize(indices.size());
356 Teuchos::ArrayView<Scalar> zerovalues(vec);
357 for(
int j=0; j<indices.size(); j++) {
358 if(dirichletCols[indices[j]]==1)
359 zerovalues[j]=(Scalar)1.0e-32;
361 zerovalues[j]=values[j];
363 A->replaceLocalValues(i,indices,zerovalues);
367 void Remove_Zeroed_Rows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
double tol=1.0e-14) {
368 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A->getRowMap();
369 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > DiagMatrix = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,1);
370 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > NewMatrix = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,1);
371 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
372 Teuchos::ArrayView<const LocalOrdinal> indices;
373 Teuchos::ArrayView<const Scalar> values;
374 A->getLocalRowView(i,indices,values);
376 for (
int j=0; j<indices.size(); j++) {
377 if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > tol) {
381 Scalar one = (Scalar)1.0;
382 Scalar zero = (Scalar)0.0;
383 GlobalOrdinal row = rowMap->getGlobalElement(i);
385 DiagMatrix->insertGlobalValues(row,
386 Teuchos::ArrayView<GlobalOrdinal>(&row,1),
387 Teuchos::ArrayView<Scalar>(&one,1));
390 DiagMatrix->insertGlobalValues(row,
391 Teuchos::ArrayView<GlobalOrdinal>(&row,1),
392 Teuchos::ArrayView<Scalar>(&zero,1));
395 DiagMatrix->fillComplete();
398 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
399 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*DiagMatrix,
false,(Scalar)1.0,*A,
false,(Scalar)1.0,NewMatrix,*out);
400 NewMatrix->fillComplete();
414 void initialize(
const Teuchos::RCP<TCRS> & D0_Matrix,
415 const Teuchos::RCP<TCRS> & M0inv_Matrix,
416 const Teuchos::RCP<TCRS> & M1_Matrix,
417 const Teuchos::RCP<TMV> & Nullspace,
418 const Teuchos::RCP<TMV> & Coords,
419 Teuchos::ParameterList& List);
442 #endif //ifdef HAVE_MUELU_TPETRA 444 #define MUELU_REFMAXWELL_SHORT 445 #endif // MUELU_REFMAXWELL_DECL_HPP RefMaxwell(Teuchos::RCP< Hierarchy > H11, Teuchos::RCP< Hierarchy > H22)
Constructor with Hierarchies.
Teuchos::RCP< XMat > Ms_Matrix_
RefMaxwell(const Teuchos::RCP< TCRS > &D0_Matrix, const Teuchos::RCP< TCRS > &M0inv_Matrix, const Teuchos::RCP< TCRS > &M1_Matrix, const Teuchos::RCP< TMV > &Nullspace, const Teuchos::RCP< TMV > &Coords, Teuchos::ParameterList &List)
Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > XCRS
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > XMV
Teuchos::RCP< XMat > M0inv_Matrix_
Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > TROW
void applyInverse121(const XTMV &RHS, XTMV &X) const
apply 1-2-1 algorithm for 2x2 solve
std::vector< LocalOrdinal > BCcols_
Teuchos::ParameterList smootherList_
Teuchos::RCP< XMat > TMT_Matrix_
RefMaxwell(const Teuchos::RCP< TCRS > &SM_Matrix, const Teuchos::RCP< TCRS > &D0_Matrix, const Teuchos::RCP< TCRS > &M0inv_Matrix, const Teuchos::RCP< TCRS > &M1_Matrix, const Teuchos::RCP< TMV > &Nullspace, const Teuchos::RCP< TMV > &Coords, Teuchos::ParameterList &List, bool ComputePrec=true)
std::vector< LocalOrdinal > BCrows_
Vectors for BCs.
void Remove_Zeroed_Rows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, double tol=1.0e-14)
Teuchos::RCP< Hierarchy > HierarchySmoother_
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > XMat
RefMaxwell(const Teuchos::RCP< TCRS > &D0_Matrix, const Teuchos::RCP< TCRS > &M1_Matrix, const Teuchos::RCP< TMV > &Nullspace, const Teuchos::RCP< TMV > &Coords, Teuchos::ParameterList &List)
Namespace for MueLu classes and methods.
Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > TMap
Teuchos::RCP< XMat > A22_
void setParameters(Teuchos::ParameterList &list)
Set parameters.
void buildProlongator()
Setup the prolongator for the (1,1)-block.
Teuchos::RCP< XMat > D0_Matrix_
Xpetra::TpetraMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > XTMV
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Preconditioner (wrapped as a Tpetra::Operator) for Maxwell's equations in curl-curl form...
void applyInverse212(const XTMV &RHS, XTMV &X) const
apply 2-1-2 algorithm for 2x2 solve
Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > XMap
void findDirichletCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > A, std::vector< LocalOrdinal > &dirichletRows, std::vector< LocalOrdinal > &dirichletCols)
Teuchos::ParameterList precList11_
Xpetra::CrsMatrixWrap< Scalar, LocalOrdinal, GlobalOrdinal, Node > XCrsWrap
void compute()
Setup the preconditioner.
Teuchos::RCP< XMat > SM_Matrix_
Various matrices.
Teuchos::RCP< XMat > A11_
Teuchos::RCP< Hierarchy > Hierarchy22_
Teuchos::RCP< XMat > TMT_Agg_Matrix_
GlobalOrdinal global_ordinal_type
Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > TCRS
void applyInverseAdditive(const XTMV &RHS, XTMV &X) const
apply additive algorithm for 2x2 solve
void Apply_BCsToMatrixCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletCols)
Teuchos::RCP< Hierarchy > Hierarchy11_
Two hierarchies: one for the (1,1)-block, another for the (2,2)-block.
bool disable_addon_
Some options.
Teuchos::RCP< Level > TopLevel_
Top Level.
virtual ~RefMaxwell()
Destructor.
Teuchos::ParameterList parameterList_
Parameter lists.
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
void findDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > A, std::vector< LocalOrdinal > &dirichletRows)
void Apply_BCsToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows)
Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > TMV
Teuchos::RCP< RefMaxwell< Scalar, LocalOrdinal, GlobalOrdinal, NewNode > > clone(const RCP< NewNode > &new_node) const
Teuchos::RCP< XMat > P11_
void initialize(const Teuchos::RCP< TCRS > &D0_Matrix, const Teuchos::RCP< TCRS > &M0inv_Matrix, const Teuchos::RCP< TCRS > &M1_Matrix, const Teuchos::RCP< TMV > &Nullspace, const Teuchos::RCP< TMV > &Coords, Teuchos::ParameterList &List)
Teuchos::RCP< XMV > Coords_
RefMaxwell(const Teuchos::RCP< TCRS > &SM_Matrix, const Teuchos::RCP< TCRS > &D0_Matrix, const Teuchos::RCP< TCRS > &M1_Matrix, const Teuchos::RCP< TMV > &Nullspace, const Teuchos::RCP< TMV > &Coords, Teuchos::ParameterList &List, bool ComputePrec=true)
Xpetra::TpetraCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > XTCRS
void apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::one()) const
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
LocalOrdinal local_ordinal_type
Teuchos::RCP< XMat > M1_Matrix_
void resetMatrix(Teuchos::RCP< TCRS > SM_Matrix_new)
Reset system matrix.
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Teuchos::ParameterList precList22_
Teuchos::RCP< XMV > Nullspace_
Nullspace.