43#ifndef __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
44#define __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
48#ifdef PANZER_HAVE_EPETRA_STACK
55#include "Thyra_DefaultBlockedLinearOp.hpp"
56#include "Thyra_DefaultProductVectorSpace.hpp"
57#include "Thyra_SpmdVectorBase.hpp"
58#include "Thyra_TpetraLinearOp.hpp"
59#include "Thyra_TpetraThyraWrappers.hpp"
62#include "Tpetra_CrsMatrix.hpp"
63#include "Tpetra_MultiVector.hpp"
64#include "Tpetra_Vector.hpp"
74template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
77 const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
78 : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
80 for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
81 gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
90template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
93 const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
94 : gidProviders_(gidProviders), comm_(comm)
99template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
107template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
108Teuchos::RCP<LinearObjContainer>
112 std::vector<Teuchos::RCP<const MapType> > blockMaps;
113 std::size_t blockDim = gidProviders_.size();
114 for(std::size_t i=0;i<blockDim;i++)
115 blockMaps.push_back(getMap(i));
117 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
118 container->setMapsForBlocks(blockMaps);
123template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
124Teuchos::RCP<LinearObjContainer>
128 std::vector<Teuchos::RCP<const MapType> > blockMaps;
129 std::size_t blockDim = gidProviders_.size();
130 for(std::size_t i=0;i<blockDim;i++)
131 blockMaps.push_back(getGhostedMap(i));
133 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
134 container->setMapsForBlocks(blockMaps);
139template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
143 using Teuchos::is_null;
147 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
148 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
152 if ( !is_null(b_in.
get_x()) && !is_null(b_out.
get_x()) && ((mem & LOC::X)==LOC::X))
153 globalToGhostThyraVector(b_in.
get_x(),b_out.
get_x());
155 if ( !is_null(b_in.
get_dxdt()) && !is_null(b_out.
get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
158 if ( !is_null(b_in.
get_f()) && !is_null(b_out.
get_f()) && ((mem & LOC::F)==LOC::F))
159 globalToGhostThyraVector(b_in.
get_f(),b_out.
get_f());
162template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
166 using Teuchos::is_null;
170 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
171 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
175 if ( !is_null(b_in.
get_x()) && !is_null(b_out.
get_x()) && ((mem & LOC::X)==LOC::X))
176 ghostToGlobalThyraVector(b_in.
get_x(),b_out.
get_x());
178 if ( !is_null(b_in.
get_f()) && !is_null(b_out.
get_f()) && ((mem & LOC::F)==LOC::F))
179 ghostToGlobalThyraVector(b_in.
get_f(),b_out.
get_f());
181 if ( !is_null(b_in.
get_A()) && !is_null(b_out.
get_A()) && ((mem & LOC::Mat)==LOC::Mat))
182 ghostToGlobalThyraMatrix(*b_in.
get_A(),*b_out.
get_A());
185template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
190 bool zeroVectorRows,
bool adjustX)
const
193 using Teuchos::rcp_dynamic_cast;
195 using Thyra::PhysicallyBlockedLinearOpBase;
199 std::size_t blockDim = gidProviders_.size();
202 const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
203 const BTLOC & b_globalBCRows = Teuchos::dyn_cast<const BTLOC>(globalBCRows);
204 BTLOC & b_ghosted = Teuchos::dyn_cast<BTLOC>(ghostedObjs);
206 TEUCHOS_ASSERT(b_localBCRows.
get_f()!=Teuchos::null);
207 TEUCHOS_ASSERT(b_globalBCRows.
get_f()!=Teuchos::null);
210 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.
get_A());
211 RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.
get_f());
212 RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.
get_f(),
true);
213 RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.
get_f(),
true);
215 if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.
get_x());
218 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(
int) blockDim);
219 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(
int) blockDim);
220 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(
int) blockDim);
221 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(
int) blockDim);
222 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(
int) blockDim);
224 for(std::size_t i=0;i<blockDim;i++) {
226 RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
227 RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
230 RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
232 if(th_f==Teuchos::null)
235 t_f = rcp_dynamic_cast<ThyraVector>(th_f,
true)->getTpetraVector();
237 for(std::size_t j=0;j<blockDim;j++) {
238 RCP<const MapType> map_i = getGhostedMap(i);
239 RCP<const MapType> map_j = getGhostedMap(j);
242 RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
245 RCP<CrsMatrixType> t_A;
246 if(th_A==Teuchos::null)
249 RCP<OperatorType> t_A_op = rcp_dynamic_cast<ThyraLinearOp>(th_A,
true)->getTpetraOperator();
250 t_A = rcp_dynamic_cast<CrsMatrixType>(t_A_op,
true);
254 if(t_A!=Teuchos::null) {
258 adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
260 if(t_A!=Teuchos::null) {
269template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
273 const Teuchos::Ptr<VectorType> & f,
274 const Teuchos::Ptr<CrsMatrixType> & A,
275 bool zeroVectorRows)
const
277 if(f==Teuchos::null && A==Teuchos::null)
280 Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
282 Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
283 Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
285 TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
286 for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
287 if(global_bcs_array[i]==0.0)
290 if(local_bcs_array[i]==0.0 || zeroVectorRows) {
294 if(!Teuchos::is_null(f))
296 if(!Teuchos::is_null(A)) {
297 std::size_t numEntries = 0;
298 std::size_t sz = A->getNumEntriesInLocalRow(i);
299 typename CrsMatrixType::nonconst_local_inds_host_view_type indices(
"indices", sz);
300 typename CrsMatrixType::nonconst_values_host_view_type values(
"values", sz);
302 A->getLocalRowCopy(i,indices,values,numEntries);
304 for(std::size_t c=0;c<numEntries;c++)
307 A->replaceLocalValues(i,indices,values);
313 ScalarT scaleFactor = global_bcs_array[i];
316 if(!Teuchos::is_null(f))
317 f_array[i] /= scaleFactor;
318 if(!Teuchos::is_null(A)) {
319 std::size_t numEntries = 0;
320 std::size_t sz = A->getNumEntriesInLocalRow(i);
321 typename CrsMatrixType::nonconst_local_inds_host_view_type indices(
"indices", sz);
322 typename CrsMatrixType::nonconst_values_host_view_type values(
"values", sz);
324 A->getLocalRowCopy(i,indices,values,numEntries);
326 for(std::size_t c=0;c<numEntries;c++)
327 values(c) /= scaleFactor;
329 A->replaceLocalValues(i,indices,values);
335template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
340 TEUCHOS_ASSERT(
false);
348template<
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
349 typename GlobalOrdinalT,
typename NodeT>
350Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
353buildReadOnlyDomainContainer()
const
360 LocalOrdinalT, GlobalOrdinalT, NodeT>;
361 vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
362 for (
int i(0); i < getBlockColCount(); ++i)
364 auto tvroged = rcp(
new TVROGED);
365 tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
366 gedBlocks.push_back(tvroged);
368 auto ged = rcp(
new BVROGED);
369 ged->initialize(getGhostedThyraDomainSpace(), getThyraDomainSpace(),
374#ifdef PANZER_HAVE_EPETRA_STACK
380template<
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
381 typename GlobalOrdinalT,
typename NodeT>
382Teuchos::RCP<WriteVector_GlobalEvaluationData>
385buildWriteDomainContainer()
const
387 using std::logic_error;
390 auto ged = rcp(
new EVWGED);
391 TEUCHOS_TEST_FOR_EXCEPTION(
true, logic_error,
"NOT YET IMPLEMENTED")
396template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
403template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
407 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
408 initializeContainer(mem,bloc);
411template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
415 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
416 initializeGhostedContainer(mem,bloc);
422template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
430 if((mem & LOC::X) == LOC::X)
431 loc.
set_x(getThyraDomainVector());
433 if((mem & LOC::DxDt) == LOC::DxDt)
434 loc.
set_dxdt(getThyraDomainVector());
436 if((mem & LOC::F) == LOC::F)
437 loc.
set_f(getThyraRangeVector());
439 if((mem & LOC::Mat) == LOC::Mat)
440 loc.
set_A(getThyraMatrix());
443template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
451 if((mem & LOC::X) == LOC::X)
452 loc.
set_x(getGhostedThyraDomainVector());
454 if((mem & LOC::DxDt) == LOC::DxDt)
455 loc.
set_dxdt(getGhostedThyraDomainVector());
457 if((mem & LOC::F) == LOC::F) {
458 loc.
set_f(getGhostedThyraRangeVector());
462 if((mem & LOC::Mat) == LOC::Mat) {
463 loc.
set_A(getGhostedThyraMatrix());
468template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
472 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
475template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
479 for(std::size_t i=0;i<exPairs.size();i++)
480 excludedPairs_.insert(exPairs[i]);
483template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
484Teuchos::RCP<const GlobalIndexer>
488 return gidProviders_[i];
491template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
495 maps_.resize(blockCnt);
496 ghostedMaps_.resize(blockCnt);
497 importers_.resize(blockCnt);
498 exporters_.resize(blockCnt);
504template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
505Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
509 if(domainSpace_==Teuchos::null) {
511 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
512 for(std::size_t i=0;i<gidProviders_.size();i++)
513 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
515 domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
521template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
522Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
526 if(rangeSpace_==Teuchos::null) {
528 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
529 for(std::size_t i=0;i<gidProviders_.size();i++)
530 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
532 rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
538template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
539Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
543 if(domainSpace_==Teuchos::null) {
544 getThyraDomainSpace();
547 return domainSpace_->getBlock(blk);
550template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
551Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
555 if(rangeSpace_==Teuchos::null) {
556 getThyraRangeSpace();
559 return rangeSpace_->getBlock(blk);
562template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
563Teuchos::RCP<Thyra::VectorBase<ScalarT> >
567 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
568 Thyra::createMember<ScalarT>(*getThyraDomainSpace());
569 Thyra::assign(vec.ptr(),0.0);
571 Teuchos::RCP<Thyra::ProductVectorBase<ScalarT> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<ScalarT> >(vec);
572 for(std::size_t i=0;i<gidProviders_.size();i++) {
573 TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
574 Teuchos::as<int>(getMap(i)->getLocalNumElements()));
580template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
581Teuchos::RCP<Thyra::VectorBase<ScalarT> >
585 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
586 Thyra::createMember<ScalarT>(*getThyraRangeSpace());
587 Thyra::assign(vec.ptr(),0.0);
592template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
593Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
597 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
600 std::size_t blockDim = gidProviders_.size();
603 blockedOp->beginBlockFill(blockDim,blockDim);
606 for(std::size_t i=0;i<blockDim;i++) {
607 for(std::size_t j=0;j<blockDim;j++) {
608 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
610 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
611 blockedOp->setNonconstBlock(i,j,block);
617 blockedOp->endBlockFill();
622template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
623Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
627 if(ghostedDomainSpace_==Teuchos::null) {
629 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
630 for(std::size_t i=0;i<gidProviders_.size();i++)
631 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
633 ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
636 return ghostedDomainSpace_;
639template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
640Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
644 if(ghostedRangeSpace_==Teuchos::null) {
646 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
647 for(std::size_t i=0;i<gidProviders_.size();i++)
648 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
650 ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
653 return ghostedRangeSpace_;
656template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
657Teuchos::RCP<Thyra::VectorBase<ScalarT> >
661 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
662 Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
663 Thyra::assign(vec.ptr(),0.0);
668template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
669Teuchos::RCP<Thyra::VectorBase<ScalarT> >
673 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
674 Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
675 Thyra::assign(vec.ptr(),0.0);
680template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
681Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
685 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
688 std::size_t blockDim = gidProviders_.size();
691 blockedOp->beginBlockFill(blockDim,blockDim);
694 for(std::size_t i=0;i<blockDim;i++) {
695 for(std::size_t j=0;j<blockDim;j++) {
696 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
698 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
699 blockedOp->setNonconstBlock(i,j,block);
705 blockedOp->endBlockFill();
710template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
716 using Teuchos::rcp_dynamic_cast;
719 std::size_t blockDim = gidProviders_.size();
722 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,
true);
723 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
725 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(
int) blockDim);
726 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(
int) blockDim);
728 for(std::size_t i=0;i<blockDim;i++) {
730 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
731 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
734 ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
738template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
743 using Teuchos::rcp_dynamic_cast;
744 using Teuchos::dyn_cast;
746 using Thyra::PhysicallyBlockedLinearOpBase;
748 std::size_t blockDim = gidProviders_.size();
751 const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
752 PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
754 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(
int) blockDim);
755 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(
int) blockDim);
756 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(
int) blockDim);
757 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(
int) blockDim);
759 for(std::size_t i=0;i<blockDim;i++) {
760 for(std::size_t j=0;j<blockDim;j++) {
761 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
763 RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
764 RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
767 TEUCHOS_ASSERT(th_in!=Teuchos::null);
768 TEUCHOS_ASSERT(th_out!=Teuchos::null);
771 RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<const ThyraLinearOp>(th_in,
true)->getConstTpetraOperator();
772 RCP<OperatorType> tp_op_out = rcp_dynamic_cast<ThyraLinearOp>(th_out,
true)->getTpetraOperator();
774 RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<const CrsMatrixType>(tp_op_in,
true);
775 RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<CrsMatrixType>(tp_op_out,
true);
778 ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
784template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
790 using Teuchos::rcp_dynamic_cast;
793 std::size_t blockDim = gidProviders_.size();
796 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,
true);
797 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
799 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(
int) blockDim);
800 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(
int) blockDim);
802 for(std::size_t i=0;i<blockDim;i++) {
804 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
805 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
808 globalToGhostTpetraVector(i,*tp_in,*tp_out);
815template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
822 RCP<const ExportType> exporter = getGhostedExport(i);
824 out.doExport(in,*exporter,Tpetra::ADD);
827template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
833 RCP<const MapType> map_i = out.getRangeMap();
834 RCP<const MapType> map_j = out.getDomainMap();
837 RCP<const ExportType> exporter = getGhostedExport(blockRow);
840 out.setAllToScalar(0.0);
841 out.doExport(in,*exporter,Tpetra::ADD);
842 out.fillComplete(map_j,map_i);
845template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
852 RCP<const ImportType> importer = getGhostedImport(i);
854 out.doImport(in,*importer,Tpetra::INSERT);
858template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
859Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
863 if(maps_[i]==Teuchos::null)
864 maps_[i] = buildTpetraMap(i);
869template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
870Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
874 if(ghostedMaps_[i]==Teuchos::null)
875 ghostedMaps_[i] = buildTpetraGhostedMap(i);
877 return ghostedMaps_[i];
881template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
882Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
886 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
888 typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
889 Teuchos::RCP<const CrsGraphType> graph;
890 if(itr==graphs_.end()) {
891 graph = buildTpetraGraph(i,j);
892 graphs_[std::make_pair(i,j)] = graph;
897 TEUCHOS_ASSERT(graph!=Teuchos::null);
901template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
902Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
906 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
908 typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
909 Teuchos::RCP<const CrsGraphType> ghostedGraph;
910 if(itr==ghostedGraphs_.end()) {
911 ghostedGraph = buildTpetraGhostedGraph(i,j);
912 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
915 ghostedGraph = itr->second;
917 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
921template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
922Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
926 if(importers_[i]==Teuchos::null)
927 importers_[i] = Teuchos::rcp(
new ImportType(getMap(i),getGhostedMap(i)));
929 return importers_[i];
932template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
933Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
937 if(exporters_[i]==Teuchos::null)
938 exporters_[i] = Teuchos::rcp(
new ExportType(getGhostedMap(i),getMap(i)));
940 return exporters_[i];
943template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
944Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
948 std::vector<GlobalOrdinalT> indices;
951 getGlobalIndexer(i)->getOwnedIndices(indices);
953 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
957template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
958Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
962 std::vector<GlobalOrdinalT> indices;
965 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
967 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
971template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
972Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
981 RCP<const MapType> map_i = getMap(i);
982 RCP<const MapType> map_j = getMap(j);
984 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,0));
985 RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
988 RCP<const ExportType> exporter = getGhostedExport(i);
989 graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
990 graph->fillComplete(map_j,map_i);
995template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
996Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
1005 RCP<const MapType> map_i = getGhostedMap(i);
1006 RCP<const MapType> map_j = getGhostedMap(j);
1008 std::vector<std::string> elementBlockIds;
1010 Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
1012 rowProvider = getGlobalIndexer(i);
1013 colProvider = getGlobalIndexer(j);
1015 gidProviders_[0]->getElementBlockIds(elementBlockIds);
1019 std::vector<size_t> nEntriesPerRow(map_i->getLocalNumElements(), 0);
1020 std::vector<std::string>::const_iterator blockItr;
1021 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1022 std::string blockId = *blockItr;
1024 const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId);
1028 std::vector<GlobalOrdinalT> row_gids;
1029 std::vector<GlobalOrdinalT> col_gids;
1032 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1034 rowProvider->getElementGIDs(elements[elmt],row_gids);
1035 colProvider->getElementGIDs(elements[elmt],col_gids);
1036 for(std::size_t row=0;row<row_gids.size();row++) {
1037 LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
1038 nEntriesPerRow[lid] += col_gids.size();
1042 Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
1043 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,map_j, nEntriesPerRowView));
1048 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1049 std::string blockId = *blockItr;
1052 const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId);
1056 std::vector<GlobalOrdinalT> row_gids;
1057 std::vector<GlobalOrdinalT> col_gids;
1060 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1062 rowProvider->getElementGIDs(elements[elmt],row_gids);
1063 colProvider->getElementGIDs(elements[elmt],col_gids);
1064 for(std::size_t row=0;row<row_gids.size();row++)
1065 graph->insertGlobalIndices(row_gids[row],col_gids);
1071 graph->fillComplete(getMap(j),getMap(i));
1076template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1077Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1081 Teuchos::RCP<const MapType> map_i = getMap(i);
1082 Teuchos::RCP<const MapType> map_j = getMap(j);
1084 Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
1085 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
1086 mat->fillComplete(map_j,map_i);
1091template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1092Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1096 Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1097 Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1099 Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1100 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
1101 mat->fillComplete(getMap(j),getMap(i));
1106template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1107Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1111 Teuchos::RCP<const MapType> tMap = getMap(i);
1115template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1116Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1120 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1124template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1125Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1129 Teuchos::RCP<const MapType> tMap = getMap(i);
1133template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1134Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1138 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1142template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1147 return gidProviders_.size();
1150template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1155 return gidProviders_.size();
1158template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1162 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1163 if(tloc.
get_A()!=Teuchos::null)
1167template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1171 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1172 if(tloc.
get_A()!=Teuchos::null)
Teuchos::RCP< CrsMatrixType > get_A() const
Teuchos::RCP< VectorType > get_f() const
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
void set_f(const Teuchos::RCP< VectorType > &in)
void set_x(const Teuchos::RCP< VectorType > &in)
Teuchos::RCP< VectorType > get_x() const
Teuchos::RCP< VectorType > get_dxdt() const
void set_dxdt(const Teuchos::RCP< VectorType > &in)
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraRangeVector() const
Get a range vector.
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
std::vector< Teuchos::RCP< const GlobalIndexer > > gidProviders_
void initializeGhostedContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraDomainVector() const
Get a domain vector.
Teuchos::RCP< VectorType > getTpetraDomainVector(int i) const
virtual void endFill(LinearObjContainer &loc) const
void initializeContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
Teuchos::MpiComm< int > getComm() const
virtual Teuchos::RCP< const CrsGraphType > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
virtual Teuchos::RCP< const MapType > buildTpetraGhostedMap(int i) const
virtual Teuchos::RCP< const MapType > buildTpetraMap(int i) const
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGraph(int i, int j) const
Teuchos::RCP< VectorType > getGhostedTpetraRangeVector(int i) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void makeRoomForBlocks(std::size_t blockCnt)
Allocate the space in the std::vector objects so we can fill with appropriate Tpetra data.
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraRangeVector() const
Get a range vector.
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGhostedGraph(int i, int j) const
virtual Teuchos::RCP< const MapType > getGhostedMap(int i) const
get the ghosted map from the matrix
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
Teuchos::RCP< VectorType > getGhostedTpetraDomainVector(int i) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
virtual void beginFill(LinearObjContainer &loc) const
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
int getBlockColCount() const
how many block columns
void globalToGhostTpetraVector(int i, const VectorType &in, VectorType &out) const
BlockedTpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const BlockedDOFManager > &gidProvider)
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
void ghostToGlobalTpetraVector(int i, const VectorType &in, VectorType &out) const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< ScalarT > &in, Thyra::LinearOpBase< ScalarT > &out) const
void ghostToGlobalTpetraMatrix(int blockRow, const CrsMatrixType &in, CrsMatrixType &out) const
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< VectorType > getTpetraRangeVector(int i) const
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) 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
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a Thyra operator.
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraDomainVector() const
Get a domain vector.
virtual Teuchos::RCP< const CrsGraphType > getGraph(int i, int j) const
get the graph of the crs matrix
virtual ~BlockedTpetraLinearObjFactory()
Teuchos::RCP< const panzer::BlockedDOFManager > getGlobalIndexer() const
int getBlockRowCount() const
how many block rows
Teuchos::RCP< Thyra::BlockedLinearOpBase< ScalarT > > getGhostedThyraMatrix() const
Get a Thyra operator.
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors.
This class provides a boundary exchange communication mechanism for vectors.
void setRequiresDirichletAdjustment(bool b)
void buildGatherScatterEvaluators(const BuilderT &builder)