43 #ifndef __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__ 44 #define __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__ 53 #include "Thyra_DefaultBlockedLinearOp.hpp" 54 #include "Thyra_DefaultProductVectorSpace.hpp" 55 #include "Thyra_SpmdVectorBase.hpp" 56 #include "Thyra_TpetraLinearOp.hpp" 57 #include "Thyra_TpetraThyraWrappers.hpp" 60 #include "Tpetra_CrsMatrix.hpp" 61 #include "Tpetra_MultiVector.hpp" 62 #include "Tpetra_Vector.hpp" 72 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
75 const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
76 : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
78 for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
79 gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
88 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
91 const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
92 : gidProviders_(gidProviders), comm_(comm)
97 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
105 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
106 Teuchos::RCP<LinearObjContainer>
110 std::vector<Teuchos::RCP<const MapType> > blockMaps;
111 std::size_t blockDim = gidProviders_.size();
112 for(std::size_t i=0;i<blockDim;i++)
113 blockMaps.push_back(getMap(i));
115 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
116 container->setMapsForBlocks(blockMaps);
121 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
122 Teuchos::RCP<LinearObjContainer>
126 std::vector<Teuchos::RCP<const MapType> > blockMaps;
127 std::size_t blockDim = gidProviders_.size();
128 for(std::size_t i=0;i<blockDim;i++)
129 blockMaps.push_back(getGhostedMap(i));
131 Teuchos::RCP<BTLOC> container = Teuchos::rcp(
new BTLOC);
132 container->setMapsForBlocks(blockMaps);
137 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
141 using Teuchos::is_null;
145 const BTLOC & b_in = Teuchos::dyn_cast<
const BTLOC>(in);
146 BTLOC & b_out = Teuchos::dyn_cast<
BTLOC>(out);
150 if ( !is_null(b_in.
get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
151 globalToGhostThyraVector(b_in.
get_x(),b_out.get_x());
153 if ( !is_null(b_in.
get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
154 globalToGhostThyraVector(b_in.
get_dxdt(),b_out.get_dxdt());
156 if ( !is_null(b_in.
get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
157 globalToGhostThyraVector(b_in.
get_f(),b_out.get_f());
160 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
164 using Teuchos::is_null;
168 const BTLOC & b_in = Teuchos::dyn_cast<
const BTLOC>(in);
169 BTLOC & b_out = Teuchos::dyn_cast<
BTLOC>(out);
173 if ( !is_null(b_in.
get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
174 ghostToGlobalThyraVector(b_in.
get_x(),b_out.get_x());
176 if ( !is_null(b_in.
get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
177 ghostToGlobalThyraVector(b_in.
get_f(),b_out.get_f());
179 if ( !is_null(b_in.
get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
180 ghostToGlobalThyraMatrix(*b_in.
get_A(),*b_out.get_A());
183 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
188 bool zeroVectorRows,
bool adjustX)
const 191 using Teuchos::rcp_dynamic_cast;
193 using Thyra::PhysicallyBlockedLinearOpBase;
197 std::size_t blockDim = gidProviders_.size();
200 const BTLOC & b_localBCRows = Teuchos::dyn_cast<
const BTLOC>(localBCRows);
201 const BTLOC & b_globalBCRows = Teuchos::dyn_cast<
const BTLOC>(globalBCRows);
202 BTLOC & b_ghosted = Teuchos::dyn_cast<
BTLOC>(ghostedObjs);
204 TEUCHOS_ASSERT(b_localBCRows.
get_f()!=Teuchos::null);
205 TEUCHOS_ASSERT(b_globalBCRows.get_f()!=Teuchos::null);
208 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.get_A());
209 RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_f());
210 RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.
get_f(),
true);
211 RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.get_f(),
true);
213 if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
216 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
217 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
218 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
219 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
220 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
222 for(std::size_t i=0;i<blockDim;i++) {
224 RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<
const ThyraVector>(local_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
225 RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<
const ThyraVector>(global_bcs->getVectorBlock(i),
true)->getConstTpetraVector();
228 RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
230 if(th_f==Teuchos::null)
233 t_f = rcp_dynamic_cast<
ThyraVector>(th_f,
true)->getTpetraVector();
235 for(std::size_t j=0;j<blockDim;j++) {
236 RCP<const MapType> map_i = getGhostedMap(i);
237 RCP<const MapType> map_j = getGhostedMap(j);
240 RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
243 RCP<CrsMatrixType> t_A;
244 if(th_A==Teuchos::null)
247 RCP<OperatorType> t_A_op = rcp_dynamic_cast<
ThyraLinearOp>(th_A,
true)->getTpetraOperator();
252 if(t_A!=Teuchos::null) {
256 adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
258 if(t_A!=Teuchos::null) {
267 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
271 const Teuchos::Ptr<VectorType> & f,
272 const Teuchos::Ptr<CrsMatrixType> & A,
273 bool zeroVectorRows)
const 275 if(f==Teuchos::null && A==Teuchos::null)
278 Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
280 Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
281 Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
283 TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
284 for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
285 if(global_bcs_array[i]==0.0)
288 if(local_bcs_array[i]==0.0 || zeroVectorRows) {
292 if(!Teuchos::is_null(f))
294 if(!Teuchos::is_null(A)) {
295 std::size_t numEntries = 0;
296 std::size_t sz = A->getNumEntriesInLocalRow(i);
297 typename CrsMatrixType::nonconst_local_inds_host_view_type indices(
"indices", sz);
298 typename CrsMatrixType::nonconst_values_host_view_type values(
"values", sz);
300 A->getLocalRowCopy(i,indices,values,numEntries);
302 for(std::size_t c=0;c<numEntries;c++)
305 A->replaceLocalValues(i,indices,values);
311 ScalarT scaleFactor = global_bcs_array[i];
314 if(!Teuchos::is_null(f))
315 f_array[i] /= scaleFactor;
316 if(!Teuchos::is_null(A)) {
317 std::size_t numEntries = 0;
318 std::size_t sz = A->getNumEntriesInLocalRow(i);
319 typename CrsMatrixType::nonconst_local_inds_host_view_type indices(
"indices", sz);
320 typename CrsMatrixType::nonconst_values_host_view_type values(
"values", sz);
322 A->getLocalRowCopy(i,indices,values,numEntries);
324 for(std::size_t c=0;c<numEntries;c++)
325 values(c) /= scaleFactor;
327 A->replaceLocalValues(i,indices,values);
333 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
338 TEUCHOS_ASSERT(
false);
346 template<
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
347 typename GlobalOrdinalT,
typename NodeT>
348 Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
351 buildReadOnlyDomainContainer()
const 358 LocalOrdinalT, GlobalOrdinalT, NodeT>;
359 vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
360 for (
int i(0); i < getBlockColCount(); ++i)
362 auto tvroged = rcp(
new TVROGED);
363 tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
364 gedBlocks.push_back(tvroged);
366 auto ged = rcp(
new BVROGED);
367 ged->initialize(getGhostedThyraDomainSpace(), getThyraDomainSpace(),
377 template<
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
378 typename GlobalOrdinalT,
typename NodeT>
379 Teuchos::RCP<WriteVector_GlobalEvaluationData>
382 buildWriteDomainContainer()
const 384 using std::logic_error;
387 auto ged = rcp(
new EVWGED);
388 TEUCHOS_TEST_FOR_EXCEPTION(
true, logic_error,
"NOT YET IMPLEMENTED")
392 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
399 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
404 initializeContainer(mem,bloc);
407 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
412 initializeGhostedContainer(mem,bloc);
418 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
426 if((mem & LOC::X) == LOC::X)
427 loc.
set_x(getThyraDomainVector());
429 if((mem & LOC::DxDt) == LOC::DxDt)
430 loc.
set_dxdt(getThyraDomainVector());
432 if((mem & LOC::F) == LOC::F)
433 loc.
set_f(getThyraRangeVector());
435 if((mem & LOC::Mat) == LOC::Mat)
436 loc.
set_A(getThyraMatrix());
439 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
447 if((mem & LOC::X) == LOC::X)
448 loc.
set_x(getGhostedThyraDomainVector());
450 if((mem & LOC::DxDt) == LOC::DxDt)
451 loc.
set_dxdt(getGhostedThyraDomainVector());
453 if((mem & LOC::F) == LOC::F) {
454 loc.
set_f(getGhostedThyraRangeVector());
458 if((mem & LOC::Mat) == LOC::Mat) {
459 loc.
set_A(getGhostedThyraMatrix());
464 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
468 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
471 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
475 for(std::size_t i=0;i<exPairs.size();i++)
476 excludedPairs_.insert(exPairs[i]);
479 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
480 Teuchos::RCP<const GlobalIndexer>
484 return gidProviders_[i];
487 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
491 maps_.resize(blockCnt);
492 ghostedMaps_.resize(blockCnt);
493 importers_.resize(blockCnt);
494 exporters_.resize(blockCnt);
500 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
501 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
505 if(domainSpace_==Teuchos::null) {
507 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
508 for(std::size_t i=0;i<gidProviders_.size();i++)
509 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
511 domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
517 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
518 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
522 if(rangeSpace_==Teuchos::null) {
524 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
525 for(std::size_t i=0;i<gidProviders_.size();i++)
526 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
528 rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
534 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
535 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
539 if(domainSpace_==Teuchos::null) {
540 getThyraDomainSpace();
543 return domainSpace_->getBlock(blk);
546 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
547 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
551 if(rangeSpace_==Teuchos::null) {
552 getThyraRangeSpace();
555 return rangeSpace_->getBlock(blk);
558 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
559 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
563 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
564 Thyra::createMember<ScalarT>(*getThyraDomainSpace());
565 Thyra::assign(vec.ptr(),0.0);
568 for(std::size_t i=0;i<gidProviders_.size();i++) {
569 TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
570 Teuchos::as<int>(getMap(i)->getNodeNumElements()));
576 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
577 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
581 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
582 Thyra::createMember<ScalarT>(*getThyraRangeSpace());
583 Thyra::assign(vec.ptr(),0.0);
588 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
589 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
593 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
596 std::size_t blockDim = gidProviders_.size();
599 blockedOp->beginBlockFill(blockDim,blockDim);
602 for(std::size_t i=0;i<blockDim;i++) {
603 for(std::size_t j=0;j<blockDim;j++) {
604 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
606 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
607 blockedOp->setNonconstBlock(i,j,block);
613 blockedOp->endBlockFill();
618 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
619 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
623 if(ghostedDomainSpace_==Teuchos::null) {
625 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
626 for(std::size_t i=0;i<gidProviders_.size();i++)
627 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
629 ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
632 return ghostedDomainSpace_;
635 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
636 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
640 if(ghostedRangeSpace_==Teuchos::null) {
642 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
643 for(std::size_t i=0;i<gidProviders_.size();i++)
644 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
646 ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
649 return ghostedRangeSpace_;
652 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
653 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
657 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
658 Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
659 Thyra::assign(vec.ptr(),0.0);
664 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
665 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
669 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
670 Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
671 Thyra::assign(vec.ptr(),0.0);
676 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
677 Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
681 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
684 std::size_t blockDim = gidProviders_.size();
687 blockedOp->beginBlockFill(blockDim,blockDim);
690 for(std::size_t i=0;i<blockDim;i++) {
691 for(std::size_t j=0;j<blockDim;j++) {
692 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
694 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
695 blockedOp->setNonconstBlock(i,j,block);
701 blockedOp->endBlockFill();
706 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
712 using Teuchos::rcp_dynamic_cast;
715 std::size_t blockDim = gidProviders_.size();
718 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<
const ProductVectorBase<ScalarT> >(in,
true);
719 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
721 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
722 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
724 for(std::size_t i=0;i<blockDim;i++) {
726 RCP<const VectorType> tp_in = rcp_dynamic_cast<
const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
727 RCP<VectorType> tp_out = rcp_dynamic_cast<
ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
730 ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
734 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
739 using Teuchos::rcp_dynamic_cast;
740 using Teuchos::dyn_cast;
742 using Thyra::PhysicallyBlockedLinearOpBase;
744 std::size_t blockDim = gidProviders_.size();
747 const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<
const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
748 PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
750 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
751 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
752 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
753 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
755 for(std::size_t i=0;i<blockDim;i++) {
756 for(std::size_t j=0;j<blockDim;j++) {
757 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
759 RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
760 RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
763 TEUCHOS_ASSERT(th_in!=Teuchos::null);
764 TEUCHOS_ASSERT(th_out!=Teuchos::null);
767 RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<
const ThyraLinearOp>(th_in,
true)->getConstTpetraOperator();
768 RCP<OperatorType> tp_op_out = rcp_dynamic_cast<
ThyraLinearOp>(th_out,
true)->getTpetraOperator();
770 RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<
const CrsMatrixType>(tp_op_in,
true);
771 RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<
CrsMatrixType>(tp_op_out,
true);
774 ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
780 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
786 using Teuchos::rcp_dynamic_cast;
789 std::size_t blockDim = gidProviders_.size();
792 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<
const ProductVectorBase<ScalarT> >(in,
true);
793 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,
true);
795 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
796 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
798 for(std::size_t i=0;i<blockDim;i++) {
800 RCP<const VectorType> tp_in = rcp_dynamic_cast<
const ThyraVector>(prod_in->getVectorBlock(i),
true)->getConstTpetraVector();
801 RCP<VectorType> tp_out = rcp_dynamic_cast<
ThyraVector>(prod_out->getNonconstVectorBlock(i),
true)->getTpetraVector();
804 globalToGhostTpetraVector(i,*tp_in,*tp_out);
811 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
818 RCP<const ExportType> exporter = getGhostedExport(i);
820 out.doExport(in,*exporter,Tpetra::ADD);
823 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
829 RCP<const MapType> map_i = out.getRangeMap();
830 RCP<const MapType> map_j = out.getDomainMap();
833 RCP<const ExportType> exporter = getGhostedExport(blockRow);
836 out.setAllToScalar(0.0);
837 out.doExport(in,*exporter,Tpetra::ADD);
838 out.fillComplete(map_j,map_i);
841 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
848 RCP<const ImportType> importer = getGhostedImport(i);
850 out.doImport(in,*importer,Tpetra::INSERT);
854 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
855 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
859 if(maps_[i]==Teuchos::null)
860 maps_[i] = buildTpetraMap(i);
865 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
866 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
870 if(ghostedMaps_[i]==Teuchos::null)
871 ghostedMaps_[i] = buildTpetraGhostedMap(i);
873 return ghostedMaps_[i];
877 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
878 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
882 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
884 typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
885 Teuchos::RCP<const CrsGraphType> graph;
886 if(itr==graphs_.end()) {
887 graph = buildTpetraGraph(i,j);
888 graphs_[std::make_pair(i,j)] = graph;
893 TEUCHOS_ASSERT(graph!=Teuchos::null);
897 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
898 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
902 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,
panzer::pair_hash> GraphMap;
904 typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
905 Teuchos::RCP<const CrsGraphType> ghostedGraph;
906 if(itr==ghostedGraphs_.end()) {
907 ghostedGraph = buildTpetraGhostedGraph(i,j);
908 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
911 ghostedGraph = itr->second;
913 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
917 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
918 Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
922 if(importers_[i]==Teuchos::null)
923 importers_[i] = Teuchos::rcp(
new ImportType(getMap(i),getGhostedMap(i)));
925 return importers_[i];
928 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
929 Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
933 if(exporters_[i]==Teuchos::null)
934 exporters_[i] = Teuchos::rcp(
new ExportType(getGhostedMap(i),getMap(i)));
936 return exporters_[i];
939 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
940 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
944 std::vector<GlobalOrdinalT> indices;
947 getGlobalIndexer(i)->getOwnedIndices(indices);
949 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
953 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
954 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
958 std::vector<GlobalOrdinalT> indices;
961 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
963 return Teuchos::rcp(
new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
967 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
968 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
977 RCP<const MapType> map_i = getMap(i);
978 RCP<const MapType> map_j = getMap(j);
980 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,0));
981 RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
984 RCP<const ExportType> exporter = getGhostedExport(i);
985 graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
986 graph->fillComplete(map_j,map_i);
991 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
992 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
1001 RCP<const MapType> map_i = getGhostedMap(i);
1002 RCP<const MapType> map_j = getGhostedMap(j);
1004 std::vector<std::string> elementBlockIds;
1006 Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
1008 rowProvider = getGlobalIndexer(i);
1009 colProvider = getGlobalIndexer(j);
1011 gidProviders_[0]->getElementBlockIds(elementBlockIds);
1015 std::vector<size_t> nEntriesPerRow(map_i->getNodeNumElements(), 0);
1016 std::vector<std::string>::const_iterator blockItr;
1017 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1018 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 LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
1034 nEntriesPerRow[lid] += col_gids.size();
1038 Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
1039 RCP<CrsGraphType> graph = rcp(
new CrsGraphType(map_i,map_j, nEntriesPerRowView,
1040 Tpetra::StaticProfile));
1045 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1046 std::string blockId = *blockItr;
1049 const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId);
1053 std::vector<GlobalOrdinalT> row_gids;
1054 std::vector<GlobalOrdinalT> col_gids;
1057 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1059 rowProvider->getElementGIDs(elements[elmt],row_gids);
1060 colProvider->getElementGIDs(elements[elmt],col_gids);
1061 for(std::size_t row=0;row<row_gids.size();row++)
1062 graph->insertGlobalIndices(row_gids[row],col_gids);
1068 graph->fillComplete(getMap(j),getMap(i));
1073 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1074 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1078 Teuchos::RCP<const MapType> map_i = getMap(i);
1079 Teuchos::RCP<const MapType> map_j = getMap(j);
1081 Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
1082 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
1083 mat->fillComplete(map_j,map_i);
1088 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1089 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1093 Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1094 Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1096 Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1097 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(
new CrsMatrixType(tGraph));
1098 mat->fillComplete(getMap(j),getMap(i));
1103 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1104 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1108 Teuchos::RCP<const MapType> tMap = getMap(i);
1112 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1113 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1117 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1121 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1122 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1126 Teuchos::RCP<const MapType> tMap = getMap(i);
1130 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1131 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1135 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1139 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1144 return gidProviders_.size();
1147 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1152 return gidProviders_.size();
1155 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1159 BTLOC & tloc = Teuchos::dyn_cast<
BTLOC>(loc);
1160 if(tloc.
get_A()!=Teuchos::null)
1164 template <
typename Traits,
typename ScalarT,
typename LocalOrdinalT,
typename GlobalOrdinalT,
typename NodeT>
1168 BTLOC & tloc = Teuchos::dyn_cast<
BTLOC>(loc);
1169 if(tloc.
get_A()!=Teuchos::null)
1175 #endif // __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
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< VectorType > getTpetraDomainVector(int i) 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
Teuchos::RCP< const panzer::BlockedDOFManager > getGlobalIndexer() const
std::vector< Teuchos::RCP< const GlobalIndexer > > 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
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
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
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.
BlockedTpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const BlockedDOFManager > &gidProvider)
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)
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...
Teuchos::RCP< VectorType > getGhostedTpetraRangeVector(int i) const
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
This class provides a boundary exchange communication mechanism for vectors.
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraDomainVector() const
Get a domain vector.
virtual void endFill(LinearObjContainer &loc) const
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
Teuchos::RCP< VectorType > getTpetraRangeVector(int i) 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< VectorType > getGhostedTpetraDomainVector(int i) const
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