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) :
49 Teko::Epetra::BlockedEpetraOperator(vars, content, label),
50 pressureMassMatrix_(Teuchos::null), gamma_(gamma)
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()));
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);
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>")