11#ifndef __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
12#define __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
16#ifdef PANZER_HAVE_EPETRA_STACK
23#include "Thyra_DefaultBlockedLinearOp.hpp"
24#include "Thyra_DefaultProductVectorSpace.hpp"
25#include "Thyra_SpmdVectorBase.hpp"
26#include "Thyra_TpetraLinearOp.hpp"
27#include "Thyra_TpetraThyraWrappers.hpp"
30#include "Tpetra_CrsMatrix.hpp"
31#include "Tpetra_MultiVector.hpp"
32#include "Tpetra_Vector.hpp"
42template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
45 const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
46 : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
48 for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
49 gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
58template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
61 const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
62 : gidProviders_(gidProviders), comm_(comm)
67template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
75template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
76Teuchos::RCP<LinearObjContainer>
80 std::vector<Teuchos::RCP<const MapType> > blockMaps;
81 std::size_t blockDim = gidProviders_.size();
82 for(std::size_t i=0;i<blockDim;i++)
83 blockMaps.push_back(getMap(i));
85 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
86 container->setMapsForBlocks(blockMaps);
91template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
92Teuchos::RCP<LinearObjContainer>
96 std::vector<Teuchos::RCP<const MapType> > blockMaps;
97 std::size_t blockDim = gidProviders_.size();
98 for(std::size_t i=0;i<blockDim;i++)
99 blockMaps.push_back(getGhostedMap(i));
101 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
102 container->setMapsForBlocks(blockMaps);
107template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
111 using Teuchos::is_null;
115 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
116 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
120 if ( !is_null(b_in.
get_x()) && !is_null(b_out.
get_x()) && ((mem & LOC::X)==LOC::X))
121 globalToGhostThyraVector(b_in.
get_x(),b_out.
get_x());
123 if ( !is_null(b_in.
get_dxdt()) && !is_null(b_out.
get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
126 if ( !is_null(b_in.
get_f()) && !is_null(b_out.
get_f()) && ((mem & LOC::F)==LOC::F))
127 globalToGhostThyraVector(b_in.
get_f(),b_out.
get_f());
130template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
134 using Teuchos::is_null;
138 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
139 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
143 if ( !is_null(b_in.
get_x()) && !is_null(b_out.
get_x()) && ((mem & LOC::X)==LOC::X))
144 ghostToGlobalThyraVector(b_in.
get_x(),b_out.
get_x());
146 if ( !is_null(b_in.
get_f()) && !is_null(b_out.
get_f()) && ((mem & LOC::F)==LOC::F))
147 ghostToGlobalThyraVector(b_in.
get_f(),b_out.
get_f());
149 if ( !is_null(b_in.
get_A()) && !is_null(b_out.
get_A()) && ((mem & LOC::Mat)==LOC::Mat))
150 ghostToGlobalThyraMatrix(*b_in.
get_A(),*b_out.
get_A());
153template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
158 bool zeroVectorRows,
bool adjustX)
const
161 using Teuchos::rcp_dynamic_cast;
163 using Thyra::PhysicallyBlockedLinearOpBase;
167 std::size_t blockDim = gidProviders_.size();
170 const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
171 const BTLOC & b_globalBCRows = Teuchos::dyn_cast<const BTLOC>(globalBCRows);
172 BTLOC & b_ghosted = Teuchos::dyn_cast<BTLOC>(ghostedObjs);
174 TEUCHOS_ASSERT(b_localBCRows.
get_f()!=Teuchos::null);
175 TEUCHOS_ASSERT(b_globalBCRows.
get_f()!=Teuchos::null);
178 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.
get_A());
179 RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.
get_f());
180 RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.
get_f(),
true);
181 RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.
get_f(),
true);
183 if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.
get_x());
186 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(
int) blockDim);
187 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(
int) blockDim);
188 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(
int) blockDim);
189 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(
int) blockDim);
190 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(
int) blockDim);
192 for(std::size_t i=0;i<blockDim;i++) {
194 RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
195 RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
198 RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
200 if(th_f==Teuchos::null)
203 t_f = rcp_dynamic_cast<ThyraVector>(th_f,
true)->getTpetraVector();
205 for(std::size_t j=0;j<blockDim;j++) {
206 RCP<const MapType> map_i = getGhostedMap(i);
207 RCP<const MapType> map_j = getGhostedMap(j);
210 RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
213 RCP<CrsMatrixType> t_A;
214 if(th_A==Teuchos::null)
217 RCP<OperatorType> t_A_op = rcp_dynamic_cast<ThyraLinearOp>(th_A,
true)->getTpetraOperator();
218 t_A = rcp_dynamic_cast<CrsMatrixType>(t_A_op,
true);
222 if(t_A!=Teuchos::null) {
226 adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
228 if(t_A!=Teuchos::null) {
237template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
241 const Teuchos::Ptr<VectorType> & f,
242 const Teuchos::Ptr<CrsMatrixType> & A,
243 bool zeroVectorRows)
const
245 if(f==Teuchos::null && A==Teuchos::null)
248 Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
250 Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
251 Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
253 TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
254 for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
255 if(global_bcs_array[i]==0.0)
258 if(local_bcs_array[i]==0.0 || zeroVectorRows) {
262 if(!Teuchos::is_null(f))
264 if(!Teuchos::is_null(A)) {
265 std::size_t numEntries = 0;
266 std::size_t sz = A->getNumEntriesInLocalRow(i);
267 typename CrsMatrixType::nonconst_local_inds_host_view_type indices(
"indices", sz);
268 typename CrsMatrixType::nonconst_values_host_view_type values(
"values", sz);
270 A->getLocalRowCopy(i,indices,values,numEntries);
272 for(std::size_t c=0;c<numEntries;c++)
275 A->replaceLocalValues(i,indices,values);
281 ScalarT scaleFactor = global_bcs_array[i];
284 if(!Teuchos::is_null(f))
285 f_array[i] /= scaleFactor;
286 if(!Teuchos::is_null(A)) {
287 std::size_t numEntries = 0;
288 std::size_t sz = A->getNumEntriesInLocalRow(i);
289 typename CrsMatrixType::nonconst_local_inds_host_view_type indices(
"indices", sz);
290 typename CrsMatrixType::nonconst_values_host_view_type values(
"values", sz);
292 A->getLocalRowCopy(i,indices,values,numEntries);
294 for(std::size_t c=0;c<numEntries;c++)
295 values(c) /= scaleFactor;
297 A->replaceLocalValues(i,indices,values);
303template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
308 TEUCHOS_ASSERT(
false);
316template<
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
317 typename GlobalOrdinalT,
typename NodeT>
318Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
321buildReadOnlyDomainContainer()
const
328 LocalOrdinalT, GlobalOrdinalT, NodeT>;
329 vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
330 for (
int i(0); i < getBlockColCount(); ++i)
332 auto tvroged = rcp(
new TVROGED);
333 tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
334 gedBlocks.push_back(tvroged);
336 auto ged = rcp(
new BVROGED);
337 ged->initialize(getGhostedThyraDomainSpace(), getThyraDomainSpace(),
342#ifdef PANZER_HAVE_EPETRA_STACK
348template<
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
349 typename GlobalOrdinalT,
typename NodeT>
350Teuchos::RCP<WriteVector_GlobalEvaluationData>
353buildWriteDomainContainer()
const
355 using std::logic_error;
358 auto ged = rcp(
new EVWGED);
359 TEUCHOS_TEST_FOR_EXCEPTION(
true, logic_error,
"NOT YET IMPLEMENTED")
364template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
371template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
375 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
376 initializeContainer(mem,bloc);
379template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
383 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
384 initializeGhostedContainer(mem,bloc);
390template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
398 if((mem & LOC::X) == LOC::X)
399 loc.
set_x(getThyraDomainVector());
401 if((mem & LOC::DxDt) == LOC::DxDt)
402 loc.
set_dxdt(getThyraDomainVector());
404 if((mem & LOC::F) == LOC::F)
405 loc.
set_f(getThyraRangeVector());
407 if((mem & LOC::Mat) == LOC::Mat)
408 loc.
set_A(getThyraMatrix());
411template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
419 if((mem & LOC::X) == LOC::X)
420 loc.
set_x(getGhostedThyraDomainVector());
422 if((mem & LOC::DxDt) == LOC::DxDt)
423 loc.
set_dxdt(getGhostedThyraDomainVector());
425 if((mem & LOC::F) == LOC::F) {
426 loc.
set_f(getGhostedThyraRangeVector());
430 if((mem & LOC::Mat) == LOC::Mat) {
431 loc.
set_A(getGhostedThyraMatrix());
436template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
440 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
443template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
447 for(std::size_t i=0;i<exPairs.size();i++)
448 excludedPairs_.insert(exPairs[i]);
451template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
452Teuchos::RCP<const GlobalIndexer>
456 return gidProviders_[i];
459template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
463 maps_.resize(blockCnt);
464 ghostedMaps_.resize(blockCnt);
465 importers_.resize(blockCnt);
466 exporters_.resize(blockCnt);
472template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
473Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
477 if(domainSpace_==Teuchos::null) {
479 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
480 for(std::size_t i=0;i<gidProviders_.size();i++)
481 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
483 domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
489template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
490Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
494 if(rangeSpace_==Teuchos::null) {
496 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
497 for(std::size_t i=0;i<gidProviders_.size();i++)
498 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
500 rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
506template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
507Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
511 if(domainSpace_==Teuchos::null) {
512 getThyraDomainSpace();
515 return domainSpace_->getBlock(blk);
518template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
519Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
523 if(rangeSpace_==Teuchos::null) {
524 getThyraRangeSpace();
527 return rangeSpace_->getBlock(blk);
530template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
531Teuchos::RCP<Thyra::VectorBase<ScalarT> >
535 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
536 Thyra::createMember<ScalarT>(*getThyraDomainSpace());
537 Thyra::assign(vec.ptr(),0.0);
539 Teuchos::RCP<Thyra::ProductVectorBase<ScalarT> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<ScalarT> >(vec);
540 for(std::size_t i=0;i<gidProviders_.size();i++) {
541 TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
542 Teuchos::as<int>(getMap(i)->getLocalNumElements()));
548template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
549Teuchos::RCP<Thyra::VectorBase<ScalarT> >
553 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
554 Thyra::createMember<ScalarT>(*getThyraRangeSpace());
555 Thyra::assign(vec.ptr(),0.0);
560template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
561Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
565 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
568 std::size_t blockDim = gidProviders_.size();
571 blockedOp->beginBlockFill(blockDim,blockDim);
574 for(std::size_t i=0;i<blockDim;i++) {
575 for(std::size_t j=0;j<blockDim;j++) {
576 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
578 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
579 blockedOp->setNonconstBlock(i,j,block);
585 blockedOp->endBlockFill();
590template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
591Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
595 if(ghostedDomainSpace_==Teuchos::null) {
597 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
598 for(std::size_t i=0;i<gidProviders_.size();i++)
599 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
601 ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
604 return ghostedDomainSpace_;
607template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
608Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
612 if(ghostedRangeSpace_==Teuchos::null) {
614 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
615 for(std::size_t i=0;i<gidProviders_.size();i++)
616 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
618 ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
621 return ghostedRangeSpace_;
624template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
625Teuchos::RCP<Thyra::VectorBase<ScalarT> >
629 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
630 Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
631 Thyra::assign(vec.ptr(),0.0);
636template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
637Teuchos::RCP<Thyra::VectorBase<ScalarT> >
641 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
642 Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
643 Thyra::assign(vec.ptr(),0.0);
648template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
649Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
653 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
656 std::size_t blockDim = gidProviders_.size();
659 blockedOp->beginBlockFill(blockDim,blockDim);
662 for(std::size_t i=0;i<blockDim;i++) {
663 for(std::size_t j=0;j<blockDim;j++) {
664 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
666 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
667 blockedOp->setNonconstBlock(i,j,block);
673 blockedOp->endBlockFill();
678template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
684 using Teuchos::rcp_dynamic_cast;
687 std::size_t blockDim = gidProviders_.size();
690 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,
true);
691 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
693 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(
int) blockDim);
694 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(
int) blockDim);
696 for(std::size_t i=0;i<blockDim;i++) {
698 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
699 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
702 ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
706template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
711 using Teuchos::rcp_dynamic_cast;
712 using Teuchos::dyn_cast;
714 using Thyra::PhysicallyBlockedLinearOpBase;
716 std::size_t blockDim = gidProviders_.size();
719 const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
720 PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
722 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(
int) blockDim);
723 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(
int) blockDim);
724 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(
int) blockDim);
725 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(
int) blockDim);
727 for(std::size_t i=0;i<blockDim;i++) {
728 for(std::size_t j=0;j<blockDim;j++) {
729 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
731 RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
732 RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
735 TEUCHOS_ASSERT(th_in!=Teuchos::null);
736 TEUCHOS_ASSERT(th_out!=Teuchos::null);
739 RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<const ThyraLinearOp>(th_in,
true)->getConstTpetraOperator();
740 RCP<OperatorType> tp_op_out = rcp_dynamic_cast<ThyraLinearOp>(th_out,
true)->getTpetraOperator();
742 RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<const CrsMatrixType>(tp_op_in,
true);
743 RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<CrsMatrixType>(tp_op_out,
true);
746 ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
752template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
758 using Teuchos::rcp_dynamic_cast;
761 std::size_t blockDim = gidProviders_.size();
764 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,
true);
765 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
767 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(
int) blockDim);
768 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(
int) blockDim);
770 for(std::size_t i=0;i<blockDim;i++) {
772 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
773 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
776 globalToGhostTpetraVector(i,*tp_in,*tp_out);
783template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
790 RCP<const ExportType> exporter = getGhostedExport(i);
792 out.doExport(in,*exporter,Tpetra::ADD);
795template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
801 RCP<const MapType> map_i = out.getRangeMap();
802 RCP<const MapType> map_j = out.getDomainMap();
805 RCP<const ExportType> exporter = getGhostedExport(blockRow);
808 out.setAllToScalar(0.0);
809 out.doExport(in,*exporter,Tpetra::ADD);
810 out.fillComplete(map_j,map_i);
813template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
820 RCP<const ImportType> importer = getGhostedImport(i);
822 out.doImport(in,*importer,Tpetra::INSERT);
826template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
827Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
831 if(maps_[i]==Teuchos::null)
832 maps_[i] = buildTpetraMap(i);
837template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
838Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
842 if(ghostedMaps_[i]==Teuchos::null)
843 ghostedMaps_[i] = buildTpetraGhostedMap(i);
845 return ghostedMaps_[i];
849template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
850Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
854 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
856 typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
857 Teuchos::RCP<const CrsGraphType> graph;
858 if(itr==graphs_.end()) {
859 graph = buildTpetraGraph(i,j);
860 graphs_[std::make_pair(i,j)] = graph;
865 TEUCHOS_ASSERT(graph!=Teuchos::null);
869template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
870Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
874 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
876 typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
877 Teuchos::RCP<const CrsGraphType> ghostedGraph;
878 if(itr==ghostedGraphs_.end()) {
879 ghostedGraph = buildTpetraGhostedGraph(i,j);
880 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
883 ghostedGraph = itr->second;
885 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
889template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
890Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
894 if(importers_[i]==Teuchos::null)
895 importers_[i] = Teuchos::rcp(
new ImportType(getMap(i),getGhostedMap(i)));
897 return importers_[i];
900template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
901Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
905 if(exporters_[i]==Teuchos::null)
906 exporters_[i] = Teuchos::rcp(
new ExportType(getGhostedMap(i),getMap(i)));
908 return exporters_[i];
911template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
912Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
916 std::vector<GlobalOrdinalT> indices;
919 getGlobalIndexer(i)->getOwnedIndices(indices);
921 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
925template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
926Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
930 std::vector<GlobalOrdinalT> indices;
933 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
935 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
939template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
940Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
949 RCP<const MapType> map_i = getMap(i);
950 RCP<const MapType> map_j = getMap(j);
952 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,0));
953 RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
956 RCP<const ExportType> exporter = getGhostedExport(i);
957 graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
958 graph->fillComplete(map_j,map_i);
963template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
964Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
973 RCP<const MapType> map_i = getGhostedMap(i);
974 RCP<const MapType> map_j = getGhostedMap(j);
976 std::vector<std::string> elementBlockIds;
978 Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
980 rowProvider = getGlobalIndexer(i);
981 colProvider = getGlobalIndexer(j);
983 gidProviders_[0]->getElementBlockIds(elementBlockIds);
987 std::vector<size_t> nEntriesPerRow(map_i->getLocalNumElements(), 0);
988 std::vector<std::string>::const_iterator blockItr;
989 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
990 std::string blockId = *blockItr;
992 const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId);
996 std::vector<GlobalOrdinalT> row_gids;
997 std::vector<GlobalOrdinalT> col_gids;
1000 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1002 rowProvider->getElementGIDs(elements[elmt],row_gids);
1003 colProvider->getElementGIDs(elements[elmt],col_gids);
1004 for(std::size_t row=0;row<row_gids.size();row++) {
1005 LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
1006 nEntriesPerRow[lid] += col_gids.size();
1010 Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
1011 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,map_j, nEntriesPerRowView));
1016 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1017 std::string blockId = *blockItr;
1020 const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId);
1024 std::vector<GlobalOrdinalT> row_gids;
1025 std::vector<GlobalOrdinalT> col_gids;
1028 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1030 rowProvider->getElementGIDs(elements[elmt],row_gids);
1031 colProvider->getElementGIDs(elements[elmt],col_gids);
1032 for(std::size_t row=0;row<row_gids.size();row++)
1033 graph->insertGlobalIndices(row_gids[row],col_gids);
1039 graph->fillComplete(getMap(j),getMap(i));
1044template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1045Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1049 Teuchos::RCP<const MapType> map_i = getMap(i);
1050 Teuchos::RCP<const MapType> map_j = getMap(j);
1052 Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
1053 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
1054 mat->fillComplete(map_j,map_i);
1059template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1060Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1064 Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1065 Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1067 Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1068 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
1069 mat->fillComplete(getMap(j),getMap(i));
1074template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1075Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1079 Teuchos::RCP<const MapType> tMap = getMap(i);
1083template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1084Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1088 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1092template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1093Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1097 Teuchos::RCP<const MapType> tMap = getMap(i);
1101template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1102Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1106 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1110template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1115 return gidProviders_.size();
1118template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1123 return gidProviders_.size();
1126template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1130 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1131 if(tloc.
get_A()!=Teuchos::null)
1135template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1139 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1140 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)