7 #include "Teko_BlockedEpetraOperator.hpp" 8 #include "Teko_BlockedMappingStrategy.hpp" 9 #include "Teko_ReorderedMappingStrategy.hpp" 11 #include "Teuchos_VerboseObject.hpp" 13 #include "Thyra_LinearOpBase.hpp" 14 #include "Thyra_EpetraLinearOp.hpp" 15 #include "Thyra_EpetraThyraWrappers.hpp" 16 #include "Thyra_DefaultProductMultiVector.hpp" 17 #include "Thyra_DefaultProductVectorSpace.hpp" 18 #include "Thyra_DefaultBlockedLinearOp.hpp" 19 #include "Thyra_get_Epetra_Operator.hpp" 21 #include "EpetraExt_MultiVectorOut.h" 22 #include "EpetraExt_RowMatrixOut.h" 26 #include "Teko_ALOperator.hpp" 35 const Teuchos::RCP<Epetra_Operator> & content,
36 LinearOp pressureMassMatrix,
double gamma,
37 const std::string & label) :
38 Teko::Epetra::BlockedEpetraOperator(vars, content, label),
39 pressureMassMatrix_(pressureMassMatrix), gamma_(gamma)
47 const Teuchos::RCP<Epetra_Operator> & content,
48 double gamma,
const std::string & label) :
60 if(pressureMassMatrix != Teuchos::null)
69 TEUCHOS_ASSERT(gamma > 0.0);
73 const Teuchos::RCP<const Epetra_Operator>
76 const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
77 = Teuchos::rcp_dynamic_cast< Thyra::BlockedLinearOpBase<double> >
78 (blockedOperator_,
true);
80 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i, j));
86 dim_ = vars.size() - 1;
87 TEUCHOS_ASSERT(
dim_ == 2 ||
dim_ == 3);
93 TEUCHOS_ASSERT(blockedMapping_ != Teuchos::null);
96 const Teuchos::RCP<const Epetra_CrsMatrix> crsContent
97 = Teuchos::rcp_dynamic_cast<
const Epetra_CrsMatrix>(fullContent_);
100 if(blockedOperator_ == Teuchos::null)
103 = blockedMapping_->buildBlockedThyraOp(crsContent, label_);
107 const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
108 = Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >
109 (blockedOperator_,
true);
110 blockedMapping_->rebuildBlockedThyraOp(crsContent, blkOp);
114 const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
115 = Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >
116 (blockedOperator_,
true);
118 Teuchos::RCP<const Thyra::LinearOpBase<double> > blockedOpBlocks[4][4];
119 for(
int i = 0; i <=
dim_; i++)
121 for(
int j = 0; j <=
dim_; j++)
123 blockedOpBlocks[i][j] = blkOp->getBlock(i, j);
135 std::cout <<
"Pressure mass matrix is null. Use identity." << std::endl;
137 = Thyra::identity<double>(blockedOpBlocks[
dim_][0]->range());
139 = Thyra::identity<double>(blockedOpBlocks[
dim_][0]->range());
143 Teuchos::RCP<Thyra::DefaultBlockedLinearOp<double> >
alOperator_ 144 = Thyra::defaultBlockedLinearOp<double>();
145 alOperator_->beginBlockFill(
dim_ + 1,
dim_ + 1);
148 for(
int i = 0; i <
dim_; i++)
150 for(
int j = 0; j <
dim_; j++)
153 alOperator_->setBlock(i, j,
154 Thyra::add(blockedOpBlocks[i][j],
156 Thyra::multiply(blockedOpBlocks[i][dim_],
162 for(
int j = 0; j <=
dim_; j++)
164 alOperator_->setBlock(dim_, j, blockedOpBlocks[dim_][j]);
168 for(
int i = 0; i <
dim_; i++)
170 alOperator_->setBlock(i, dim_,
171 Thyra::add(blockedOpBlocks[i][dim_],
173 Thyra::multiply(blockedOpBlocks[i][dim_],
177 alOperator_->endBlockFill();
180 SetOperator(alOperator_,
false);
183 Teuchos::RCP<Thyra::DefaultBlockedLinearOp<double> > alOpRhs_
184 = Thyra::defaultBlockedLinearOp<double>();
185 alOpRhs_->beginBlockFill(dim_ + 1, dim_ + 1);
188 for(
int i = 0; i <
dim_; i++)
191 alOpRhs_->setBlock(i, i,
192 Thyra::identity<double>(blockedOpBlocks[0][0]->range()));
194 alOpRhs_->setBlock(dim_, dim_,
195 Thyra::identity<double>(blockedOpBlocks[dim_][dim_]->range()));
198 for(
int i = 0; i <
dim_; i++)
200 alOpRhs_->setBlock(i, dim_,
205 alOpRhs_->endBlockFill();
210 if(reorderManager_ != Teuchos::null)
217 Teuchos::RCP<const Teko::Epetra::MappingStrategy> mapping
219 Teuchos::RCP<Thyra::MultiVectorBase<double> > bThyra
220 = Thyra::createMembers(thyraOp_->range(), b.NumVectors());
221 Teuchos::RCP<Thyra::MultiVectorBase<double> > bThyraAugmented
222 = Thyra::createMembers(thyraOp_->range(), b.NumVectors());
225 mapping->copyEpetraIntoThyra(b, bThyra.ptr());
227 alOperatorRhs_->apply(Thyra::NOTRANS, *bThyra, bThyraAugmented.ptr(), 1.0, 0.0);
229 mapping->copyThyraIntoEpetra(bThyraAugmented, bAugmented);
Teuchos::RCP< Thyra::LinearOpBase< double > > alOperatorRhs_
void augmentRHS(const Epetra_MultiVector &b, Epetra_MultiVector &bAugmented)
LinearOp invPressureMassMatrix_
void checkDim(const std::vector< std::vector< int > > &vars)
void Reorder(const BlockReorderManager &brm)
void setPressureMassMatrix(LinearOp pressureMassMatrix)
virtual void SetContent(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< const Epetra_Operator > &content)
const RCP< const MappingStrategy > getMapStrategy() const
Get the mapping strategy for this wrapper (translate between Thyra and Epetra)
LinearOp pressureMassMatrix_
BlockedEpetraOperator(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< const Epetra_Operator > &content, const std::string &label="<ANYM>")
const Teuchos::RCP< const Epetra_Operator > GetBlock(int i, int j) const
ALOperator(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< Epetra_Operator > &content, LinearOp pressureMassMatrix, double gamma=0.05, const std::string &label="<ANYM>")
void setGamma(double gamma)
Teuchos::RCP< Thyra::LinearOpBase< double > > alOperator_