43 #ifndef PANZER_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP 44 #define PANZER_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP 48 #include "Tpetra_MultiVector.hpp" 49 #include "Tpetra_Vector.hpp" 50 #include "Tpetra_CrsMatrix.hpp" 52 #include "Thyra_TpetraThyraWrappers.hpp" 53 #include "Thyra_DefaultProductVectorSpace.hpp" 54 #include "Thyra_DefaultBlockedLinearOp.hpp" 55 #include "Thyra_TpetraLinearOp.hpp" 56 #include "Thyra_SpmdVectorBase.hpp" 69 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
73 : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(
comm)
75 for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
76 gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
85 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
93 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
94 Teuchos::RCP<LinearObjContainer>
98 std::vector<Teuchos::RCP<const MapType> > blockMaps;
99 std::size_t blockDim = gidProviders_.size();
100 for(std::size_t i=0;i<blockDim;i++)
101 blockMaps.push_back(getMap(i));
103 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
104 container->setMapsForBlocks(blockMaps);
109 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
110 Teuchos::RCP<LinearObjContainer>
114 std::vector<Teuchos::RCP<const MapType> > blockMaps;
115 std::size_t blockDim = gidProviders_.size();
116 for(std::size_t i=0;i<blockDim;i++)
117 blockMaps.push_back(getGhostedMap(i));
119 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
120 container->setMapsForBlocks(blockMaps);
125 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
129 using Teuchos::is_null;
133 const BTLOC & b_in = Teuchos::dyn_cast<
const BTLOC>(in);
134 BTLOC & b_out = Teuchos::dyn_cast<
BTLOC>(out);
138 if ( !is_null(b_in.
get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
139 globalToGhostThyraVector(b_in.
get_x(),b_out.get_x());
141 if ( !is_null(b_in.
get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
142 globalToGhostThyraVector(b_in.
get_dxdt(),b_out.get_dxdt());
144 if ( !is_null(b_in.
get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
145 globalToGhostThyraVector(b_in.
get_f(),b_out.get_f());
148 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
152 using Teuchos::is_null;
156 const BTLOC & b_in = Teuchos::dyn_cast<
const BTLOC>(in);
157 BTLOC & b_out = Teuchos::dyn_cast<
BTLOC>(out);
161 if ( !is_null(b_in.
get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
162 ghostToGlobalThyraVector(b_in.
get_x(),b_out.get_x());
164 if ( !is_null(b_in.
get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
165 ghostToGlobalThyraVector(b_in.
get_f(),b_out.get_f());
167 if ( !is_null(b_in.
get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
168 ghostToGlobalThyraMatrix(*b_in.
get_A(),*b_out.get_A());
171 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
176 bool zeroVectorRows,
bool adjustX)
const 179 using Teuchos::rcp_dynamic_cast;
181 using Thyra::PhysicallyBlockedLinearOpBase;
185 std::size_t blockDim = gidProviders_.size();
188 const BTLOC & b_localBCRows = Teuchos::dyn_cast<
const BTLOC>(localBCRows);
189 const BTLOC & b_globalBCRows = Teuchos::dyn_cast<
const BTLOC>(globalBCRows);
190 BTLOC & b_ghosted = Teuchos::dyn_cast<
BTLOC>(ghostedObjs);
192 TEUCHOS_ASSERT(b_localBCRows.
get_f()!=Teuchos::null);
193 TEUCHOS_ASSERT(b_globalBCRows.get_f()!=Teuchos::null);
196 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.get_A());
197 RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_f());
198 RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.
get_f(),
true);
199 RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.get_f(),
true);
201 if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
204 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
205 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
206 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
207 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
208 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
210 for(std::size_t i=0;i<blockDim;i++) {
212 RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<
const ThyraVector>(local_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
213 RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<
const ThyraVector>(global_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
216 RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
218 if(th_f==Teuchos::null)
221 t_f = rcp_dynamic_cast<
ThyraVector>(th_f,
true)->getTpetraVector();
223 for(std::size_t j=0;j<blockDim;j++) {
224 RCP<const MapType> map_i = getGhostedMap(i);
225 RCP<const MapType> map_j = getGhostedMap(j);
228 RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
231 RCP<CrsMatrixType> t_A;
232 if(th_A==Teuchos::null)
235 RCP<OperatorType> t_A_op = rcp_dynamic_cast<
ThyraLinearOp>(th_A,
true)->getTpetraOperator();
240 if(t_A!=Teuchos::null) {
244 adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
246 if(t_A!=Teuchos::null) {
255 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
259 const Teuchos::Ptr<VectorType> & f,
260 const Teuchos::Ptr<CrsMatrixType> & A,
261 bool zeroVectorRows)
const 263 if(f==Teuchos::null && A==Teuchos::null)
266 Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
268 Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
269 Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
271 TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
272 for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
273 if(global_bcs_array[i]==0.0)
276 if(local_bcs_array[i]==0.0 || zeroVectorRows) {
280 if(!Teuchos::is_null(f))
282 if(!Teuchos::is_null(A)) {
283 std::size_t numEntries = 0;
284 std::size_t sz = A->getNumEntriesInLocalRow(i);
285 Teuchos::Array<LocalOrdinalT> indices(sz);
286 Teuchos::Array<ScalarT>
values(sz);
288 A->getLocalRowCopy(i,indices,
values,numEntries);
290 for(std::size_t c=0;c<numEntries;c++)
293 A->replaceLocalValues(i,indices,
values);
299 ScalarT scaleFactor = global_bcs_array[i];
302 if(!Teuchos::is_null(f))
303 f_array[i] /= scaleFactor;
304 if(!Teuchos::is_null(A)) {
305 std::size_t numEntries = 0;
306 std::size_t sz = A->getNumEntriesInLocalRow(i);
307 Teuchos::Array<LocalOrdinalT> indices(sz);
308 Teuchos::Array<ScalarT>
values(sz);
310 A->getLocalRowCopy(i,indices,
values,numEntries);
312 for(std::size_t c=0;c<numEntries;c++)
315 A->replaceLocalValues(i,indices,
values);
321 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
326 TEUCHOS_ASSERT(
false);
329 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
330 Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
337 std::vector<RCP<ReadOnlyVector_GlobalEvaluationData> > gedBlocks;
338 for(
int i=0;i<getBlockColCount();i++) {
339 RCP<TpetraVector_ReadOnly_GlobalEvaluationData<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > vec_ged
341 vec_ged->initialize(getGhostedImport(i),getGhostedMap(i),getMap(i));
343 gedBlocks.push_back(vec_ged);
347 ged->initialize(getGhostedThyraDomainSpace(),getThyraDomainSpace(),gedBlocks);
352 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
359 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
364 initializeContainer(mem,bloc);
367 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
372 initializeGhostedContainer(mem,bloc);
378 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
386 if((mem & LOC::X) == LOC::X)
387 loc.
set_x(getThyraDomainVector());
389 if((mem & LOC::DxDt) == LOC::DxDt)
390 loc.
set_dxdt(getThyraDomainVector());
392 if((mem & LOC::F) == LOC::F)
393 loc.
set_f(getThyraRangeVector());
395 if((mem & LOC::Mat) == LOC::Mat)
396 loc.
set_A(getThyraMatrix());
399 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
407 if((mem & LOC::X) == LOC::X)
408 loc.
set_x(getGhostedThyraDomainVector());
410 if((mem & LOC::DxDt) == LOC::DxDt)
411 loc.
set_dxdt(getGhostedThyraDomainVector());
413 if((mem & LOC::F) == LOC::F) {
414 loc.
set_f(getGhostedThyraRangeVector());
418 if((mem & LOC::Mat) == LOC::Mat) {
419 loc.
set_A(getGhostedThyraMatrix());
424 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
428 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
431 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
435 for(std::size_t i=0;i<exPairs.size();i++)
436 excludedPairs_.insert(exPairs[i]);
439 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
440 Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> >
444 return gidProviders_[i];
447 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
451 maps_.resize(blockCnt);
452 ghostedMaps_.resize(blockCnt);
453 importers_.resize(blockCnt);
454 exporters_.resize(blockCnt);
460 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
461 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
465 if(domainSpace_==Teuchos::null) {
467 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
468 for(std::size_t i=0;i<gidProviders_.size();i++)
469 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
471 domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
477 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
478 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
482 if(rangeSpace_==Teuchos::null) {
484 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
485 for(std::size_t i=0;i<gidProviders_.size();i++)
486 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
488 rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
494 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
495 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
499 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
500 Thyra::createMember<ScalarT>(*getThyraDomainSpace());
501 Thyra::assign(vec.ptr(),0.0);
504 for(std::size_t i=0;i<gidProviders_.size();i++) {
505 TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
506 Teuchos::as<int>(getMap(i)->getNodeNumElements()));
512 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
513 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
517 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
518 Thyra::createMember<ScalarT>(*getThyraRangeSpace());
519 Thyra::assign(vec.ptr(),0.0);
524 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
525 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
529 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
532 std::size_t blockDim = gidProviders_.size();
535 blockedOp->beginBlockFill(blockDim,blockDim);
538 for(std::size_t i=0;i<blockDim;i++) {
539 for(std::size_t j=0;j<blockDim;j++) {
540 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
542 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
543 blockedOp->setNonconstBlock(i,j,block);
549 blockedOp->endBlockFill();
554 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
555 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
559 if(ghostedDomainSpace_==Teuchos::null) {
561 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
562 for(std::size_t i=0;i<gidProviders_.size();i++)
563 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
565 ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
568 return ghostedDomainSpace_;
571 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
572 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
576 if(ghostedRangeSpace_==Teuchos::null) {
578 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
579 for(std::size_t i=0;i<gidProviders_.size();i++)
580 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
582 ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
585 return ghostedRangeSpace_;
588 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
589 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
593 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
594 Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
595 Thyra::assign(vec.ptr(),0.0);
600 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
601 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
605 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
606 Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
607 Thyra::assign(vec.ptr(),0.0);
612 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
613 Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
617 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
620 std::size_t blockDim = gidProviders_.size();
623 blockedOp->beginBlockFill(blockDim,blockDim);
626 for(std::size_t i=0;i<blockDim;i++) {
627 for(std::size_t j=0;j<blockDim;j++) {
628 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
630 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
631 blockedOp->setNonconstBlock(i,j,block);
637 blockedOp->endBlockFill();
642 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
648 using Teuchos::rcp_dynamic_cast;
651 std::size_t blockDim = gidProviders_.size();
654 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<
const ProductVectorBase<ScalarT> >(in,
true);
655 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
657 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
658 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
660 for(std::size_t i=0;i<blockDim;i++) {
662 RCP<const VectorType> tp_in = rcp_dynamic_cast<
const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
663 RCP<VectorType> tp_out = rcp_dynamic_cast<
ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
666 ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
670 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
675 using Teuchos::rcp_dynamic_cast;
676 using Teuchos::dyn_cast;
678 using Thyra::PhysicallyBlockedLinearOpBase;
680 std::size_t blockDim = gidProviders_.size();
683 const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<
const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
684 PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
686 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
687 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
688 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
689 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
691 for(std::size_t i=0;i<blockDim;i++) {
692 for(std::size_t j=0;j<blockDim;j++) {
693 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
695 RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
696 RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
699 TEUCHOS_ASSERT(th_in!=Teuchos::null);
700 TEUCHOS_ASSERT(th_out!=Teuchos::null);
703 RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<
const ThyraLinearOp>(th_in,
true)->getConstTpetraOperator();
704 RCP<OperatorType> tp_op_out = rcp_dynamic_cast<
ThyraLinearOp>(th_out,
true)->getTpetraOperator();
706 RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<
const CrsMatrixType>(tp_op_in,
true);
707 RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<
CrsMatrixType>(tp_op_out,
true);
710 ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
716 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
722 using Teuchos::rcp_dynamic_cast;
725 std::size_t blockDim = gidProviders_.size();
728 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<
const ProductVectorBase<ScalarT> >(in,
true);
729 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
731 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
732 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
734 for(std::size_t i=0;i<blockDim;i++) {
736 RCP<const VectorType> tp_in = rcp_dynamic_cast<
const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
737 RCP<VectorType> tp_out = rcp_dynamic_cast<
ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
740 globalToGhostTpetraVector(i,*tp_in,*tp_out);
747 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
754 RCP<const ExportType> exporter = getGhostedExport(i);
756 out.doExport(in,*exporter,Tpetra::ADD);
759 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
765 RCP<const MapType> map_i = out.getRangeMap();
766 RCP<const MapType> map_j = out.getDomainMap();
769 RCP<const ExportType> exporter = getGhostedExport(blockRow);
772 out.setAllToScalar(0.0);
773 out.doExport(in,*exporter,Tpetra::ADD);
774 out.fillComplete(map_j,map_i);
777 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
784 RCP<const ImportType> importer = getGhostedImport(i);
786 out.doImport(in,*importer,Tpetra::INSERT);
790 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
791 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
795 if(maps_[i]==Teuchos::null)
796 maps_[i] = buildTpetraMap(i);
801 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
802 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
806 if(ghostedMaps_[i]==Teuchos::null)
807 ghostedMaps_[i] = buildTpetraGhostedMap(i);
809 return ghostedMaps_[i];
813 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
814 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
818 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
820 typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
821 Teuchos::RCP<const CrsGraphType> graph;
822 if(itr==graphs_.end()) {
823 graph = buildTpetraGraph(i,j);
824 graphs_[std::make_pair(i,j)] = graph;
829 TEUCHOS_ASSERT(graph!=Teuchos::null);
833 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
834 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
838 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
840 typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
841 Teuchos::RCP<const CrsGraphType> ghostedGraph;
842 if(itr==ghostedGraphs_.end()) {
843 ghostedGraph = buildTpetraGhostedGraph(i,j);
844 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
847 ghostedGraph = itr->second;
849 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
853 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
854 Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
858 if(importers_[i]==Teuchos::null)
859 importers_[i] = Teuchos::rcp(
new ImportType(getMap(i),getGhostedMap(i)));
861 return importers_[i];
864 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
865 Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
869 if(exporters_[i]==Teuchos::null)
870 exporters_[i] = Teuchos::rcp(
new ExportType(getGhostedMap(i),getMap(i)));
872 return exporters_[i];
875 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
876 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
880 std::vector<GlobalOrdinalT> indices;
883 getGlobalIndexer(i)->getOwnedIndices(indices);
885 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
889 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
890 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
894 std::vector<GlobalOrdinalT> indices;
897 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
899 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
903 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
904 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
913 RCP<const MapType> map_i = getMap(i);
914 RCP<const MapType> map_j = getMap(j);
916 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,0));
917 RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
920 RCP<const ExportType> exporter = getGhostedExport(i);
921 graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
922 graph->fillComplete(map_j,map_i);
927 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
928 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
937 RCP<const MapType> map_i = getGhostedMap(i);
938 RCP<const MapType> map_j = getGhostedMap(j);
940 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,map_j,0));
942 std::vector<std::string> elementBlockIds;
944 Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> > rowProvider, colProvider;
946 rowProvider = getGlobalIndexer(i);
947 colProvider = getGlobalIndexer(j);
949 blockProvider_->getElementBlockIds(elementBlockIds);
953 std::vector<std::string>::const_iterator blockItr;
954 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
955 std::string blockId = *blockItr;
958 const std::vector<LocalOrdinalT> & elements = blockProvider_->getElementBlock(blockId);
962 std::vector<GlobalOrdinalT> row_gids;
963 std::vector<GlobalOrdinalT> col_gids;
966 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
968 rowProvider->getElementGIDs(elements[elmt],row_gids);
969 colProvider->getElementGIDs(elements[elmt],col_gids);
970 for(std::size_t row=0;row<row_gids.size();row++)
971 graph->insertGlobalIndices(row_gids[row],col_gids);
977 graph->fillComplete(map_j,map_i);
982 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
983 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
987 Teuchos::RCP<const MapType> map_i = getMap(i);
988 Teuchos::RCP<const MapType> map_j = getMap(j);
990 Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
991 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
992 mat->fillComplete(map_j,map_i);
997 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
998 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1002 Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1003 Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1005 Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1006 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
1007 mat->fillComplete(map_j,map_i);
1012 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1017 return gidProviders_.size();
1020 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1025 return gidProviders_.size();
1028 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1032 BTLOC & tloc = Teuchos::dyn_cast<
BTLOC>(loc);
1033 if(tloc.
get_A()!=Teuchos::null)
1037 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1041 BTLOC & tloc = Teuchos::dyn_cast<
BTLOC>(loc);
1042 if(tloc.
get_A()!=Teuchos::null)
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
Thyra::TpetraVector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > ThyraVector
void ghostToGlobalTpetraMatrix(int blockRow, const CrsMatrixType &in, CrsMatrixType &out) const
Teuchos::RCP< const panzer::BlockedDOFManager< int, GlobalOrdinalT > > getGlobalIndexer() const
virtual ~BlockedTpetraLinearObjFactory()
virtual Teuchos::RCP< const MapType > buildTpetraMap(int i) const
Thyra::TpetraLinearOp< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > ThyraLinearOp
void set_dxdt(const Teuchos::RCP< VectorType > &in)
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
std::vector< Teuchos::RCP< const UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > gidProviders_
virtual Teuchos::RCP< const MapType > buildTpetraGhostedMap(int i) const
virtual Teuchos::RCP< const ImportType > getGhostedImport(int i) const
get importer for converting an overalapped object to a "normal" object
Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > ImportType
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGraph(int i, int j) const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a Thyra operator.
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
Teuchos::MpiComm< int > getComm() const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
Teuchos::RCP< VectorType > get_f() const
PHX::MDField< ScalarT > vector
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
Teuchos::RCP< Thyra::BlockedLinearOpBase< ScalarT > > getGhostedThyraMatrix() const
Get a Thyra operator.
Teuchos::RCP< VectorType > get_dxdt() const
int getBlockColCount() const
how many block columns
void initializeGhostedContainer(int, LinearObjContainer &loc) const
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
virtual Teuchos::RCP< const CrsGraphType > getGraph(int i, int j) const
get the graph of the crs matrix
Teuchos::RCP< CrsMatrixType > get_A() const
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraRangeVector() const
Get a range vector.
void ghostToGlobalTpetraVector(int i, const VectorType &in, VectorType &out) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void setRequiresDirichletAdjustment(bool b)
std::vector< ScalarT > values
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraRangeVector() const
Get a range vector.
virtual Teuchos::RCP< const CrsGraphType > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
void makeRoomForBlocks(std::size_t blockCnt)
Allocate the space in the std::vector objects so we can fill with appropriate Tpetra data...
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraDomainVector() const
Get a domain vector.
virtual void endFill(LinearObjContainer &loc) const
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildDomainContainer() const
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraDomainVector() const
Get a domain vector.
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range vector space (f)
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
virtual Teuchos::RCP< const MapType > getGhostedMap(int i) const
get the ghosted map from the matrix
Teuchos::RCP< VectorType > get_x() const
void globalToGhostTpetraVector(int i, const VectorType &in, VectorType &out) const
int getBlockRowCount() const
how many block rows
void set_x(const Teuchos::RCP< VectorType > &in)
Teuchos::RCP< const Teuchos::Comm< int > > comm
BlockedTpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const BlockedDOFManager< LocalOrdinalT, GlobalOrdinalT > > &gidProvider)
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void buildGatherScatterEvaluators(const BuilderT &builder)
virtual void beginFill(LinearObjContainer &loc) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
void initializeContainer(int, LinearObjContainer &loc) const
void set_f(const Teuchos::RCP< VectorType > &in)
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< ScalarT > &in, Thyra::LinearOpBase< ScalarT > &out) const
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGhostedGraph(int i, int j) const
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const