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" 64 #include "Tpetra_MultiVector_decl.hpp" 65 #include "MatrixMarket_Tpetra.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 =
88 Tpetra::Operator<>::scalar_type,
90 typename Tpetra::Operator<Scalar>::local_ordinal_type,
92 typename Tpetra::Operator<Scalar, LocalOrdinal>::global_ordinal_type,
94 typename Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
95 class RefMaxwell :
public Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
97 #undef MUELU_REFMAXWELL_SHORT 102 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>
TMap;
103 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
TCRS;
104 typedef Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
TROW;
105 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
TMV;
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;
150 bool ComputePrec =
true)
152 initialize(D0_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
177 initialize(D0_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
196 bool ComputePrec =
true)
198 initialize(D0_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
221 initialize(D0_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
260 void apply(
const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
261 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
269 template <
class NewNode>
280 std::vector<LocalOrdinal>& dirichletRows) {
281 dirichletRows.resize(0);
282 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
285 A->getLocalRowView(i,indices,values);
287 for (
int j=0; j<indices.
size(); j++) {
297 if (nnz == 1 || nnz == 2) {
298 dirichletRows.push_back(i);
304 std::vector<LocalOrdinal>& dirichletRows,
305 std::vector<LocalOrdinal>& dirichletCols) {
309 = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
312 myColsToZero->putScalar((Scalar)0.0);
313 globalColsToZero->putScalar((Scalar)0.0);
314 for(
size_t i=0; i<dirichletRows.size(); i++) {
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);
324 dirichletCols.resize(colMap->getNodeNumElements());
325 for(
size_t i=0; i<colMap->getNodeNumElements(); i++) {
334 std::vector<LocalOrdinal>& dirichletRows) {
335 for(
size_t i=0; i<dirichletRows.size(); i++) {
338 A->getLocalRowView(dirichletRows[i],indices,values);
339 std::vector<Scalar> vec;
340 vec.resize(indices.
size());
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++) {
353 A->getLocalRowView(i,indices,values);
354 std::vector<Scalar> vec;
355 vec.resize(indices.
size());
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);
371 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
374 A->getLocalRowView(i,indices,values);
376 for (
int j=0; j<indices.
size(); j++) {
381 Scalar one = (Scalar)1.0;
382 Scalar zero = (Scalar)0.0;
383 GlobalOrdinal row = rowMap->getGlobalElement(i);
385 DiagMatrix->insertGlobalValues(row,
390 DiagMatrix->insertGlobalValues(row,
395 DiagMatrix->fillComplete();
399 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*DiagMatrix,
false,(Scalar)1.0,*A,
false,(Scalar)1.0,NewMatrix,*out);
400 NewMatrix->fillComplete();
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
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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_
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.
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.