43 #ifndef PANZER_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP 44 #define PANZER_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP 50 #include "Epetra_MultiVector.h" 51 #include "Epetra_Vector.h" 52 #include "Epetra_CrsMatrix.h" 53 #include "Epetra_MpiComm.h" 58 #include "Thyra_EpetraThyraWrappers.hpp" 59 #include "Thyra_DefaultProductVectorSpace.hpp" 60 #include "Thyra_DefaultProductVector.hpp" 61 #include "Thyra_DefaultBlockedLinearOp.hpp" 62 #include "Thyra_EpetraLinearOp.hpp" 63 #include "Thyra_SpmdVectorBase.hpp" 64 #include "Thyra_get_Epetra_Operator.hpp" 65 #include "Thyra_VectorStdOps.hpp" 78 template <
typename Traits,
typename LocalOrdinalT>
81 const Teuchos::RCP<const UniqueGlobalIndexerBase> & gidProvider,
82 bool useDiscreteAdjoint)
83 : useColGidProviders_(false), eComm_(
Teuchos::null)
84 , rawMpiComm_(
comm->getRawMpiComm())
85 , useDiscreteAdjoint_(useDiscreteAdjoint)
101 template <
typename Traits,
typename LocalOrdinalT>
104 const Teuchos::RCP<const UniqueGlobalIndexerBase> & gidProvider,
105 const Teuchos::RCP<const UniqueGlobalIndexerBase> & colGidProvider,
106 bool useDiscreteAdjoint)
108 , rawMpiComm_(
comm->getRawMpiComm())
109 , useDiscreteAdjoint_(useDiscreteAdjoint)
127 template <
typename Traits,
typename LocalOrdinalT>
134 template <
typename Traits,
typename LocalOrdinalT>
140 using Teuchos::rcp_dynamic_cast;
141 using Teuchos::dyn_cast;
147 RCP<Thyra::VectorBase<double> > vec;
159 TEUCHOS_ASSERT(
false);
163 int blockRows = this->getBlockRowCount();
164 RCP<ProductVectorBase<double> > b_vec = Thyra::nonconstProductVectorBase(vec);
167 for(
int i=0;i<blockRows;i++) {
168 RCP<Thyra::VectorBase<double> > x = b_vec->getNonconstVectorBlock(i);
169 RCP<Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(i),x);
172 std::stringstream ss;
173 ss << identifier <<
"-" << i <<
".mm";
176 Epetra_Vector * ptr_ex = 0;
184 template <
typename Traits,
typename LocalOrdinalT>
190 using Teuchos::rcp_dynamic_cast;
191 using Teuchos::dyn_cast;
197 RCP<const Thyra::VectorBase<double> > vec;
209 TEUCHOS_ASSERT(
false);
213 int blockRows = this->getBlockRowCount();
214 RCP<const ProductVectorBase<double> > b_vec = Thyra::productVectorBase(vec);
217 for(
int i=0;i<blockRows;i++) {
218 RCP<const Thyra::VectorBase<double> > x = b_vec->getVectorBlock(i);
219 RCP<const Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(i),x);
222 std::stringstream ss;
223 ss << identifier <<
"-" << i <<
".mm";
230 template <
typename Traits,
typename LocalOrdinalT>
234 if(!rowDOFManagerContainer_->containsBlockedDOFManager() && !colDOFManagerContainer_->containsBlockedDOFManager()) {
235 Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(
new EpetraLinearObjContainer(getColMap(0),getMap(0)));
240 std::vector<Teuchos::RCP<const Epetra_Map> > blockMaps;
241 std::size_t blockDim = getBlockRowCount();
242 for(std::size_t i=0;i<blockDim;i++)
243 blockMaps.push_back(getMap(i));
246 container->setMapsForBlocks(blockMaps);
251 template <
typename Traits,
typename LocalOrdinalT>
255 if(!rowDOFManagerContainer_->containsBlockedDOFManager() && !colDOFManagerContainer_->containsBlockedDOFManager()) {
256 Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(
new EpetraLinearObjContainer(getGhostedColMap(0),getGhostedMap(0)));
261 std::vector<Teuchos::RCP<const Epetra_Map> > blockMaps;
262 std::size_t blockDim = getBlockRowCount();
263 for(std::size_t i=0;i<blockDim;i++)
264 blockMaps.push_back(getGhostedMap(i));
267 container->setMapsForBlocks(blockMaps);
272 template <
typename Traits,
typename LocalOrdinalT>
276 using Teuchos::is_null;
281 if( !rowDOFManagerContainer_->containsBlockedDOFManager()
282 && !colDOFManagerContainer_->containsBlockedDOFManager()) {
288 if ( !is_null(e_in.
get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
289 globalToGhostEpetraVector(0,*e_in.
get_x(),*e_out.get_x(),
true);
291 if ( !is_null(e_in.
get_dxdt()) && !is_null(e_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
292 globalToGhostEpetraVector(0,*e_in.
get_dxdt(),*e_out.get_dxdt(),
true);
294 if ( !is_null(e_in.
get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
295 globalToGhostEpetraVector(0,*e_in.
get_f(),*e_out.get_f(),
false);
298 const BLOC & b_in = Teuchos::dyn_cast<
const BLOC>(in);
299 BLOC & b_out = Teuchos::dyn_cast<BLOC>(out);
303 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
304 globalToGhostThyraVector(b_in.get_x(),b_out.get_x(),
true);
306 if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
307 globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt(),
true);
309 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
310 globalToGhostThyraVector(b_in.get_f(),b_out.get_f(),
false);
314 template <
typename Traits,
typename LocalOrdinalT>
318 using Teuchos::is_null;
323 if( !rowDOFManagerContainer_->containsBlockedDOFManager()
324 && !colDOFManagerContainer_->containsBlockedDOFManager()) {
330 if ( !is_null(e_in.
get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
331 ghostToGlobalEpetraVector(0,*e_in.
get_x(),*e_out.get_x(),
true);
333 if ( !is_null(e_in.
get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
334 ghostToGlobalEpetraVector(0,*e_in.
get_f(),*e_out.get_f(),
false);
336 if ( !is_null(e_in.
get_A()) && !is_null(e_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
337 ghostToGlobalEpetraMatrix(0,*e_in.
get_A(),*e_out.get_A());
340 const BLOC & b_in = Teuchos::dyn_cast<
const BLOC>(in);
341 BLOC & b_out = Teuchos::dyn_cast<BLOC>(out);
345 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
346 ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x(),
true);
348 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
349 ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f(),
false);
351 if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
352 ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
356 template <
typename Traits,
typename LocalOrdinalT>
361 bool zeroVectorRows,
bool adjustX)
const 366 using Teuchos::rcp_dynamic_cast;
368 using Thyra::PhysicallyBlockedLinearOpBase;
371 using Thyra::get_Epetra_Vector;
372 using Thyra::get_Epetra_Operator;
374 int rBlockDim = getBlockRowCount();
375 int cBlockDim = getBlockColCount();
378 const TOC & b_localBCRows = Teuchos::dyn_cast<
const TOC>(localBCRows);
379 const TOC & b_globalBCRows = Teuchos::dyn_cast<
const TOC>(globalBCRows);
380 TOC & b_ghosted = Teuchos::dyn_cast<TOC>(ghostedObjs);
382 TEUCHOS_ASSERT(b_localBCRows.get_f_th()!=Teuchos::null);
383 TEUCHOS_ASSERT(b_globalBCRows.get_f_th()!=Teuchos::null);
386 RCP<PhysicallyBlockedLinearOpBase<double> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<double> >(b_ghosted.get_A_th());
387 if(A==Teuchos::null && b_ghosted.get_A_th()!=Teuchos::null) {
389 A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<double> >(Thyra::nonconstBlock1x1(b_ghosted.get_A_th()));
392 RCP<ProductVectorBase<double> > f = b_ghosted.get_f_th()==Teuchos::null
394 : Thyra::castOrCreateNonconstProductVectorBase(b_ghosted.get_f_th());
395 RCP<ProductVectorBase<double> > local_bcs = b_localBCRows.get_f_th()==Teuchos::null
397 : Thyra::castOrCreateNonconstProductVectorBase(b_localBCRows.get_f_th());
398 RCP<ProductVectorBase<double> > global_bcs = b_globalBCRows.get_f_th()==Teuchos::null
400 : Thyra::castOrCreateNonconstProductVectorBase(b_globalBCRows.get_f_th());
402 if(adjustX) f = Thyra::castOrCreateNonconstProductVectorBase(b_ghosted.get_x_th());
405 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==rBlockDim);
406 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==cBlockDim);
407 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==rBlockDim);
408 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==rBlockDim);
409 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==rBlockDim);
411 for(
int i=0;i<rBlockDim;i++) {
413 RCP<const Epetra_Vector> e_local_bcs = get_Epetra_Vector(*getGhostedMap(i),local_bcs->getVectorBlock(i));
414 RCP<const Epetra_Vector> e_global_bcs = get_Epetra_Vector(*getGhostedMap(i),global_bcs->getVectorBlock(i));
417 RCP<VectorBase<double> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
418 RCP<Epetra_Vector> e_f;
419 if(th_f==Teuchos::null)
422 e_f = get_Epetra_Vector(*getGhostedMap(i),th_f);
424 for(
int j=0;j<cBlockDim;j++) {
427 RCP<LinearOpBase<double> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
430 RCP<Epetra_CrsMatrix> e_A;
431 if(th_A==Teuchos::null)
434 e_A = rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*th_A),
true);
437 adjustForDirichletConditions(*e_local_bcs,*e_global_bcs,e_f.ptr(),e_A.ptr(),zeroVectorRows);
444 template <
typename Traits,
typename LocalOrdinalT>
447 const Epetra_Vector & global_bcs,
448 const Teuchos::Ptr<Epetra_Vector> & f,
449 const Teuchos::Ptr<Epetra_CrsMatrix> & A,
450 bool zeroVectorRows)
const 452 if(f==Teuchos::null && A==Teuchos::null)
455 TEUCHOS_ASSERT(local_bcs.MyLength()==global_bcs.MyLength());
456 for(
int i=0;i<local_bcs.MyLength();i++) {
457 if(global_bcs[i]==0.0)
464 if(local_bcs[i]==0.0 || zeroVectorRows) {
468 if(!Teuchos::is_null(f))
470 if(!Teuchos::is_null(A)) {
471 A->ExtractMyRowView(i,numEntries,
values,indices);
472 for(
int c=0;c<numEntries;c++)
479 double scaleFactor = global_bcs[i];
482 if(!Teuchos::is_null(f))
483 (*f)[i] /= scaleFactor;
484 if(!Teuchos::is_null(A)) {
485 A->ExtractMyRowView(i,numEntries,
values,indices);
486 for(
int c=0;c<numEntries;c++)
493 template <
typename Traits,
typename LocalOrdinalT>
499 using Teuchos::rcp_dynamic_cast;
500 using Teuchos::dyn_cast;
507 RCP<const PVector> count = Thyra::castOrCreateProductVectorBase(th_counter.
get_f_th().getConst());
508 RCP<const PVector> f_in = Thyra::castOrCreateProductVectorBase(th_counter.
get_f_th().getConst());
509 RCP<PVector> f_out = Thyra::castOrCreateNonconstProductVectorBase(th_result.get_f_th());
511 int rBlockDim = getBlockRowCount();
512 for(
int i=0;i<rBlockDim;i++) {
514 Teuchos::ArrayRCP<const double> count_array,f_in_array;
515 Teuchos::ArrayRCP<double> f_out_array;
517 rcp_dynamic_cast<
const Thyra::SpmdVectorBase<double> >(count->getVectorBlock(i),
true)->getLocalData(Teuchos::ptrFromRef(count_array));
518 rcp_dynamic_cast<
const Thyra::SpmdVectorBase<double> >(f_in->getVectorBlock(i),
true)->getLocalData(Teuchos::ptrFromRef(f_in_array));
519 rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(f_out->getNonconstVectorBlock(i),
true)->getNonconstLocalData(Teuchos::ptrFromRef(f_out_array));
521 TEUCHOS_ASSERT(count_array.size()==f_in_array.size());
522 TEUCHOS_ASSERT(count_array.size()==f_out_array.size());
524 for(Teuchos::ArrayRCP<double>::size_type i=0;i<count_array.size();i++) {
525 if(count_array[i]!=0.0)
526 f_out_array[i] = f_in_array[i];
531 template <
typename Traits,
typename LocalOrdinalT>
532 Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
540 if(!colDOFManagerContainer_->containsBlockedDOFManager()) {
542 ged->initialize(getGhostedColImport(0),getGhostedColMap(0),getColMap(0));
547 std::vector<RCP<ReadOnlyVector_GlobalEvaluationData> > gedBlocks;
548 for(
int i=0;i<getBlockColCount();i++) {
550 vec_ged->initialize(getGhostedColImport(i),getGhostedColMap(i),getColMap(i));
552 gedBlocks.push_back(vec_ged);
556 ged->initialize(getGhostedThyraDomainSpace(),getThyraDomainSpace(),gedBlocks);
562 template <
typename Traits,
typename LocalOrdinalT>
569 template <
typename Traits,
typename LocalOrdinalT>
575 TOC & toc = Teuchos::dyn_cast<TOC>(loc);
576 initializeContainer_internal(mem,toc);
579 template <
typename Traits,
typename LocalOrdinalT>
586 TOC & toc = Teuchos::dyn_cast<TOC>(loc);
587 initializeGhostedContainer_internal(mem,toc);
589 if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
592 BLOC & bloc = Teuchos::dyn_cast<BLOC>(loc);
594 if((mem & LOC::F) == LOC::F)
597 if((mem & LOC::Mat) == LOC::Mat)
598 bloc.setRequiresDirichletAdjustment(
true);
603 if((mem & LOC::F) == LOC::F)
606 if((mem & LOC::Mat) == LOC::Mat)
614 template <
typename Traits,
typename LocalOrdinalT>
622 if((mem & LOC::X) == LOC::X)
623 loc.
set_x_th(getThyraDomainVector());
625 if((mem & LOC::DxDt) == LOC::DxDt)
628 if((mem & LOC::F) == LOC::F)
629 loc.
set_f_th(getThyraRangeVector());
631 if((mem & LOC::Mat) == LOC::Mat)
635 template <
typename Traits,
typename LocalOrdinalT>
643 if((mem & LOC::X) == LOC::X)
644 loc.
set_x_th(getGhostedThyraDomainVector());
646 if((mem & LOC::DxDt) == LOC::DxDt)
649 if((mem & LOC::F) == LOC::F)
650 loc.
set_f_th(getGhostedThyraRangeVector());
652 if((mem & LOC::Mat) == LOC::Mat)
653 loc.
set_A_th(getGhostedThyraMatrix());
656 template <
typename Traits,
typename LocalOrdinalT>
660 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
663 template <
typename Traits,
typename LocalOrdinalT>
667 for(std::size_t i=0;i<exPairs.size();i++)
668 excludedPairs_.insert(exPairs[i]);
671 template <
typename Traits,
typename LocalOrdinalT>
675 return rowDOFManagerContainer_->getFieldDOFManagers()[i];
678 template <
typename Traits,
typename LocalOrdinalT>
682 return colDOFManagerContainer_->getFieldDOFManagers()[i];
685 template <
typename Traits,
typename LocalOrdinalT>
689 maps_.resize(blockCnt);
690 ghostedMaps_.resize(blockCnt);
691 importers_.resize(blockCnt);
692 exporters_.resize(blockCnt);
695 colMaps_.resize(colBlockCnt);
696 colGhostedMaps_.resize(colBlockCnt);
697 colImporters_.resize(colBlockCnt);
698 colExporters_.resize(colBlockCnt);
705 template <
typename Traits,
typename LocalOrdinalT>
709 if(domainSpace_==Teuchos::null) {
710 if(colDOFManagerContainer_->containsBlockedDOFManager()) {
712 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
713 for(
int i=0;i<getBlockColCount();i++)
714 vsArray.push_back(Thyra::create_VectorSpace(getColMap(i)));
716 domainSpace_ = Thyra::productVectorSpace<double>(vsArray);
721 domainSpace_ = Thyra::create_VectorSpace(getColMap(0));
728 template <
typename Traits,
typename LocalOrdinalT>
732 if(rangeSpace_==Teuchos::null) {
733 if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
735 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
736 for(
int i=0;i<getBlockRowCount();i++)
737 vsArray.push_back(Thyra::create_VectorSpace(getMap(i)));
739 rangeSpace_ = Thyra::productVectorSpace<double>(vsArray);
744 rangeSpace_ = Thyra::create_VectorSpace(getMap(0));
751 template <
typename Traits,
typename LocalOrdinalT>
755 Teuchos::RCP<Thyra::VectorBase<double> > vec =
756 Thyra::createMember<double>(*getThyraDomainSpace());
757 Thyra::assign(vec.ptr(),0.0);
762 template <
typename Traits,
typename LocalOrdinalT>
766 Teuchos::RCP<Thyra::VectorBase<double> > vec =
767 Thyra::createMember<double>(*getThyraRangeSpace());
768 Thyra::assign(vec.ptr(),0.0);
773 template <
typename Traits,
typename LocalOrdinalT>
778 if(!rowDOFManagerContainer_->containsBlockedDOFManager() &&
779 !colDOFManagerContainer_->containsBlockedDOFManager()) {
780 return Thyra::nonconstEpetraLinearOp(getEpetraMatrix(0,0));
783 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockedOp = Thyra::defaultBlockedLinearOp<double>();
786 std::size_t rBlockDim = getBlockRowCount();
787 std::size_t cBlockDim = getBlockColCount();
790 blockedOp->beginBlockFill(rBlockDim,cBlockDim);
793 for(std::size_t i=0;i<rBlockDim;i++) {
794 for(std::size_t j=0;j<cBlockDim;j++) {
795 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
797 Teuchos::RCP<Thyra::LinearOpBase<double> > block = Thyra::nonconstEpetraLinearOp(getEpetraMatrix(i,j));
798 blockedOp->setNonconstBlock(i,j,block);
804 blockedOp->endBlockFill();
809 template <
typename Traits,
typename LocalOrdinalT>
813 if(ghostedDomainSpace_==Teuchos::null) {
814 if(colDOFManagerContainer_->containsBlockedDOFManager()) {
816 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
817 for(
int i=0;i<getBlockColCount();i++)
818 vsArray.push_back(Thyra::create_VectorSpace(getGhostedColMap(i)));
820 ghostedDomainSpace_ = Thyra::productVectorSpace<double>(vsArray);
825 ghostedDomainSpace_ = Thyra::create_VectorSpace(getGhostedColMap(0));
829 return ghostedDomainSpace_;
832 template <
typename Traits,
typename LocalOrdinalT>
836 if(ghostedRangeSpace_==Teuchos::null) {
837 if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
839 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
840 for(
int i=0;i<getBlockRowCount();i++)
841 vsArray.push_back(Thyra::create_VectorSpace(getGhostedMap(i)));
843 ghostedRangeSpace_ = Thyra::productVectorSpace<double>(vsArray);
848 ghostedRangeSpace_ = Thyra::create_VectorSpace(getGhostedMap(0));
852 return ghostedRangeSpace_;
855 template <
typename Traits,
typename LocalOrdinalT>
859 Teuchos::RCP<Thyra::VectorBase<double> > vec =
860 Thyra::createMember<double>(*getGhostedThyraDomainSpace());
861 Thyra::assign(vec.ptr(),0.0);
866 template <
typename Traits,
typename LocalOrdinalT>
870 Teuchos::RCP<Thyra::VectorBase<double> > vec =
871 Thyra::createMember<double>(*getGhostedThyraRangeSpace());
872 Thyra::assign(vec.ptr(),0.0);
877 template <
typename Traits,
typename LocalOrdinalT>
882 if(!rowDOFManagerContainer_->containsBlockedDOFManager() &&
883 !colDOFManagerContainer_->containsBlockedDOFManager()) {
884 return Thyra::nonconstEpetraLinearOp(getGhostedEpetraMatrix(0,0));
887 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockedOp = Thyra::defaultBlockedLinearOp<double>();
890 std::size_t rBlockDim = getBlockRowCount();
891 std::size_t cBlockDim = getBlockColCount();
894 blockedOp->beginBlockFill(rBlockDim,cBlockDim);
897 for(std::size_t i=0;i<rBlockDim;i++) {
898 for(std::size_t j=0;j<cBlockDim;j++) {
899 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
901 Teuchos::RCP<Thyra::LinearOpBase<double> > block = Thyra::nonconstEpetraLinearOp(getGhostedEpetraMatrix(i,j));
902 blockedOp->setNonconstBlock(i,j,block);
908 blockedOp->endBlockFill();
913 template <
typename Traits,
typename LocalOrdinalT>
919 using Teuchos::rcp_dynamic_cast;
921 using Thyra::get_Epetra_Vector;
923 std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
926 RCP<const ProductVectorBase<double> > prod_in = Thyra::castOrCreateProductVectorBase(in);
927 RCP<ProductVectorBase<double> > prod_out = Thyra::castOrCreateNonconstProductVectorBase(out);
929 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
930 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
932 for(std::size_t i=0;i<blockDim;i++) {
934 RCP<const Epetra_Vector> ep_in;
935 RCP<Epetra_Vector> ep_out;
938 ep_in = get_Epetra_Vector(*getGhostedMap(i),prod_in->getVectorBlock(i));
939 ep_out = get_Epetra_Vector(*getMap(i),prod_out->getNonconstVectorBlock(i));
941 ep_in = get_Epetra_Vector(*getGhostedColMap(i),prod_in->getVectorBlock(i));
942 ep_out = get_Epetra_Vector(*getColMap(i),prod_out->getNonconstVectorBlock(i));
946 ghostToGlobalEpetraVector(i,*ep_in,*ep_out,col);
950 template <
typename Traits,
typename LocalOrdinalT>
955 using Teuchos::rcp_dynamic_cast;
956 using Teuchos::dyn_cast;
958 using Thyra::PhysicallyBlockedLinearOpBase;
959 using Thyra::get_Epetra_Operator;
962 std::size_t rBlockDim = getBlockRowCount();
963 std::size_t cBlockDim = getBlockColCount();
966 const PhysicallyBlockedLinearOpBase<double> & prod_in = dyn_cast<
const PhysicallyBlockedLinearOpBase<double> >(in);
967 PhysicallyBlockedLinearOpBase<double> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<double> >(out);
969 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) rBlockDim);
970 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) cBlockDim);
971 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) rBlockDim);
972 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) cBlockDim);
974 for(std::size_t i=0;i<rBlockDim;i++) {
975 for(std::size_t j=0;j<cBlockDim;j++) {
976 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
978 RCP<const LinearOpBase<double> > th_in = prod_in.getBlock(i,j);
979 RCP<LinearOpBase<double> > th_out = prod_out.getNonconstBlock(i,j);
982 TEUCHOS_ASSERT(th_in!=Teuchos::null);
983 TEUCHOS_ASSERT(th_out!=Teuchos::null);
986 RCP<const Epetra_CrsMatrix> ep_in = rcp_dynamic_cast<
const Epetra_CrsMatrix>(get_Epetra_Operator(*th_in),
true);
987 RCP<Epetra_CrsMatrix> ep_out = rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*th_out),
true);
990 ghostToGlobalEpetraMatrix(i,*ep_in,*ep_out);
996 template <
typename Traits,
typename LocalOrdinalT>
1002 using Teuchos::rcp_dynamic_cast;
1004 using Thyra::get_Epetra_Vector;
1006 std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
1009 RCP<const ProductVectorBase<double> > prod_in = Thyra::castOrCreateProductVectorBase(in);
1010 RCP<ProductVectorBase<double> > prod_out = Thyra::castOrCreateNonconstProductVectorBase(out);
1012 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
1013 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
1015 for(std::size_t i=0;i<blockDim;i++) {
1017 RCP<const Epetra_Vector> ep_in;
1018 RCP<Epetra_Vector> ep_out;
1021 ep_in = get_Epetra_Vector(*getMap(i),prod_in->getVectorBlock(i));
1022 ep_out = get_Epetra_Vector(*getGhostedMap(i),prod_out->getNonconstVectorBlock(i));
1025 ep_in = get_Epetra_Vector(*getColMap(i),prod_in->getVectorBlock(i));
1026 ep_out = get_Epetra_Vector(*getGhostedColMap(i),prod_out->getNonconstVectorBlock(i));
1030 globalToGhostEpetraVector(i,*ep_in,*ep_out,col);
1037 template <
typename Traits,
typename LocalOrdinalT>
1044 RCP<Epetra_Export> exporter = col ? getGhostedColExport(i) : getGhostedExport(i);
1046 int err = out.Export(in,*exporter,Add);
1047 TEUCHOS_ASSERT_EQUALITY(err,0);
1050 template <
typename Traits,
typename LocalOrdinalT>
1057 RCP<Epetra_Export> exporter = getGhostedExport(blockRow);
1059 int err = out.Export(in,*exporter,Add);
1060 TEUCHOS_ASSERT_EQUALITY(err,0);
1063 template <
typename Traits,
typename LocalOrdinalT>
1070 RCP<Epetra_Import> importer = col ? getGhostedColImport(i) : getGhostedImport(i);
1072 int err = out.Import(in,*importer,Insert);
1073 TEUCHOS_ASSERT_EQUALITY(err,0);
1077 template <
typename Traits,
typename LocalOrdinalT>
1081 if(maps_[i]==Teuchos::null)
1082 maps_[i] = buildMap(i);
1088 template <
typename Traits,
typename LocalOrdinalT>
1092 if(not useColGidProviders_)
1095 if(colMaps_[i]==Teuchos::null)
1096 colMaps_[i] = buildColMap(i);
1101 template <
typename Traits,
typename LocalOrdinalT>
1105 if(ghostedMaps_[i]==Teuchos::null)
1106 ghostedMaps_[i] = buildGhostedMap(i);
1108 return ghostedMaps_[i];
1111 template <
typename Traits,
typename LocalOrdinalT>
1115 if(not useColGidProviders_)
1116 return getGhostedMap(i);
1118 if(colGhostedMaps_[i]==Teuchos::null)
1119 colGhostedMaps_[i] = buildColGhostedMap(i);
1121 return colGhostedMaps_[i];
1125 template <
typename Traits,
typename LocalOrdinalT>
1129 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,
panzer::pair_hash> GraphMap;
1131 GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
1132 Teuchos::RCP<Epetra_CrsGraph> graph;
1133 if(itr==graphs_.end()) {
1134 graph = buildGraph(i,j);
1135 graphs_[std::make_pair(i,j)] = graph;
1138 graph = itr->second;
1140 TEUCHOS_ASSERT(graph!=Teuchos::null);
1144 template <
typename Traits,
typename LocalOrdinalT>
1148 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,
panzer::pair_hash> GraphMap;
1150 GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
1151 Teuchos::RCP<Epetra_CrsGraph> ghostedGraph;
1152 if(itr==ghostedGraphs_.end()) {
1153 ghostedGraph = buildGhostedGraph(i,j,
true);
1154 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
1157 ghostedGraph = itr->second;
1159 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
1160 return ghostedGraph;
1163 template <
typename Traits,
typename LocalOrdinalT>
1167 if(importers_[i]==Teuchos::null)
1168 importers_[i] = Teuchos::rcp(
new Epetra_Import(*getGhostedMap(i),*getMap(i)));
1170 return importers_[i];
1173 template <
typename Traits,
typename LocalOrdinalT>
1177 if(not useColGidProviders_)
1178 return getGhostedImport(i);
1180 if(colImporters_[i]==Teuchos::null)
1181 colImporters_[i] = Teuchos::rcp(
new Epetra_Import(*getGhostedColMap(i),*getColMap(i)));
1183 return colImporters_[i];
1186 template <
typename Traits,
typename LocalOrdinalT>
1190 if(exporters_[i]==Teuchos::null)
1191 exporters_[i] = Teuchos::rcp(
new Epetra_Export(*getGhostedMap(i),*getMap(i)));
1193 return exporters_[i];
1196 template <
typename Traits,
typename LocalOrdinalT>
1200 if(not useColGidProviders_)
1201 return getGhostedExport(i);
1203 if(colExporters_[i]==Teuchos::null)
1204 colExporters_[i] = Teuchos::rcp(
new Epetra_Export(*getGhostedColMap(i),*getColMap(i)));
1206 return colExporters_[i];
1209 template <
typename Traits,
typename LocalOrdinalT>
1213 std::vector<int> indices;
1216 getGlobalIndexer(i)->getOwnedIndices(indices);
1218 return Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1221 template <
typename Traits,
typename LocalOrdinalT>
1225 if(not useColGidProviders_)
1228 std::vector<int> indices;
1231 getColGlobalIndexer(i)->getOwnedIndices(indices);
1233 return Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1237 template <
typename Traits,
typename LocalOrdinalT>
1241 std::vector<int> indices;
1244 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
1246 return Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1250 template <
typename Traits,
typename LocalOrdinalT>
1254 if(not useColGidProviders_)
1255 return buildGhostedMap(i);
1257 std::vector<int> indices;
1260 getColGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
1262 return Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1266 template <
typename Traits,
typename LocalOrdinalT>
1275 RCP<Epetra_Map> map_i = getMap(i);
1276 RCP<Epetra_Map> map_j = getColMap(j);
1278 TEUCHOS_ASSERT(map_i!=Teuchos::null);
1279 TEUCHOS_ASSERT(map_j!=Teuchos::null);
1281 RCP<Epetra_CrsGraph> graph = rcp(
new Epetra_CrsGraph(Copy,*map_i,0));
1282 RCP<Epetra_CrsGraph> oGraph = buildFilteredGhostedGraph(i,j);
1287 RCP<Epetra_Export> exporter = getGhostedExport(i);
1288 int err = graph->Export( *oGraph, *exporter, Insert );
1289 TEUCHOS_ASSERT_EQUALITY(err,0);
1290 graph->FillComplete(*map_j,*map_i);
1291 graph->OptimizeStorage();
1296 template <
typename Traits,
typename LocalOrdinalT>
1301 Teuchos::RCP<Epetra_Map> rowMap = getGhostedMap(i);
1302 Teuchos::RCP<Epetra_Map> colMap = getGhostedColMap(j);
1303 Teuchos::RCP<Epetra_CrsGraph> graph = Teuchos::rcp(
new Epetra_CrsGraph(Copy,*rowMap,*colMap,0));
1305 std::vector<std::string> elementBlockIds;
1307 Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,int> > rowProvider, colProvider;
1309 rowProvider = getGlobalIndexer(i);
1310 colProvider = getColGlobalIndexer(j);
1312 rowProvider->getElementBlockIds(elementBlockIds);
1315 const Teuchos::RCP<const ConnManagerBase<LocalOrdinalT> > conn_mgr = colProvider->getConnManagerBase();
1316 const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors();
1319 std::vector<std::string>::const_iterator blockItr;
1320 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1321 std::string blockId = *blockItr;
1324 const std::vector<LocalOrdinalT> & elements = rowProvider->getElementBlock(blockId);
1328 std::vector<int> row_gids;
1329 std::vector<int> col_gids;
1332 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1334 rowProvider->getElementGIDs(elements[elmt],row_gids);
1335 colProvider->getElementGIDs(elements[elmt],col_gids);
1338 const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[elmt]);
1339 for (
typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
1340 eit != aes.end(); ++eit) {
1341 std::vector<int> other_col_gids;
1342 colProvider->getElementGIDs(*eit, other_col_gids);
1343 col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
1347 for(std::size_t row=0;row<row_gids.size();row++)
1348 graph->InsertGlobalIndices(row_gids[row],col_gids.size(),&col_gids[0]);
1354 graph->FillComplete(*colMap,*rowMap);
1356 graph->OptimizeStorage();
1361 template <
typename Traits,
typename LocalOrdinalT>
1367 using Teuchos::rcp_dynamic_cast;
1370 RCP<const Filtered_UniqueGlobalIndexer<LocalOrdinalT,int> > filtered_ugi
1374 if(filtered_ugi==Teuchos::null)
1375 return buildGhostedGraph(i,j,
true);
1378 std::vector<int> ghostedActive;
1379 filtered_ugi->getOwnedAndGhostedNotFilteredIndicator(ghostedActive);
1382 Teuchos::RCP<Epetra_CrsGraph> filteredGraph = buildGhostedGraph(i,j,
false);
1386 for(
int i=0;i<filteredGraph->NumMyRows();i++) {
1387 std::vector<int> removedIndices;
1390 TEUCHOS_ASSERT(filteredGraph->ExtractMyRowView(i,numIndices,indices)==0);
1392 for(
int j=0;j<numIndices;j++) {
1393 if(ghostedActive[indices[j]]==0)
1394 removedIndices.push_back(indices[j]);
1397 TEUCHOS_ASSERT(filteredGraph->RemoveMyIndices(i,Teuchos::as<int>(removedIndices.size()),&removedIndices[0])==0);
1401 Teuchos::RCP<Epetra_Map> rowMap = getGhostedMap(i);
1402 Teuchos::RCP<Epetra_Map> colMap = getGhostedColMap(j);
1404 TEUCHOS_ASSERT(filteredGraph->FillComplete(*colMap,*rowMap)==0);
1405 TEUCHOS_ASSERT(filteredGraph->OptimizeStorage()==0);
1407 return filteredGraph;
1410 template <
typename Traits,
typename LocalOrdinalT>
1414 Teuchos::RCP<Epetra_CrsGraph> eGraph = getGraph(i,j);
1415 Teuchos::RCP<Epetra_CrsMatrix> mat = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *eGraph));
1416 TEUCHOS_ASSERT(mat->Filled());
1420 template <
typename Traits,
typename LocalOrdinalT>
1424 Teuchos::RCP<Epetra_CrsGraph> eGraph = getGhostedGraph(i,j);
1425 Teuchos::RCP<Epetra_CrsMatrix> mat = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *eGraph));
1426 TEUCHOS_ASSERT(mat->Filled());
1430 template <
typename Traits,
typename LocalOrdinalT>
1437 template <
typename Traits,
typename LocalOrdinalT>
1441 return rowDOFManagerContainer_->getFieldBlocks();
1444 template <
typename Traits,
typename LocalOrdinalT>
1448 return colDOFManagerContainer_->getFieldBlocks();
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraRangeVector() const
Get a range vector.
void ghostToGlobalEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
virtual const Teuchos::RCP< Epetra_Map > buildMap(int i) const
Teuchos::RCP< Epetra_CrsMatrix > getGhostedEpetraMatrix(int i, int j) const
virtual const Teuchos::RCP< Epetra_Export > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< Thyra::LinearOpBase< double > > getThyraMatrix() const
Get a Thyra operator.
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
virtual void readVector(const std::string &identifier, LinearObjContainer &loc, int id) const
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
Teuchos::RCP< const Epetra_Comm > eComm_
Teuchos::RCP< Thyra::VectorBase< double > > getThyraDomainVector() const
Get a domain vector.
virtual const Teuchos::RCP< Epetra_Map > buildColMap(int i) const
void globalToGhostEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
Teuchos::RCP< const UniqueGlobalIndexer< LocalOrdinalT, int > > getColGlobalIndexer(int i) const
int getBlockColCount() const
how many block columns
virtual void set_x_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
virtual const Teuchos::RCP< Epetra_CrsGraph > buildGhostedGraph(int i, int j, bool optimizeStorage) const
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildDomainContainer() const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
virtual Teuchos::RCP< Thyra::VectorBase< ScalarT > > get_f_th() const =0
virtual void set_dxdt_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
virtual const Teuchos::RCP< Epetra_Map > getMap(int i) const
get the map from the matrix
Teuchos::RCP< Thyra::VectorBase< double > > getThyraRangeVector() const
Get a range vector.
void ghostToGlobalEpetraMatrix(int blockRow, const Epetra_CrsMatrix &in, Epetra_CrsMatrix &out) const
virtual const Teuchos::RCP< Epetra_Import > getGhostedImport(int i) const
get importer for converting an overalapped object to a "normal" object
PHX::MDField< ScalarT > vector
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< double > &in, Thyra::LinearOpBase< double > &out) const
virtual const Teuchos::RCP< Epetra_CrsGraph > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
virtual const Teuchos::RCP< Epetra_CrsGraph > buildFilteredGhostedGraph(int i, int j) const
virtual const Teuchos::RCP< Epetra_Map > getColMap(int i) const
get the map from the matrix
virtual ~BlockedEpetraLinearObjFactory()
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
Teuchos::RCP< const DOFManagerContainer > colDOFManagerContainer_
Teuchos::RCP< const Teuchos::OpaqueWrapper< MPI_Comm > > rawMpiComm_
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
const Teuchos::RCP< Epetra_Vector > get_x() const
virtual void set_f_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
void setRequiresDirichletAdjustment(bool b)
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraDomainVector() const
Get a domain vector.
std::vector< ScalarT > values
virtual const Teuchos::RCP< const Epetra_Comm > getEpetraComm() const
get exporter for converting an overalapped object to a "normal" object
void makeRoomForBlocks(std::size_t blockCnt, std::size_t colBlockCnt=0)
Allocate the space in the std::vector objects so we can fill with appropriate Epetra data...
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
Teuchos::RCP< VectorType > get_dxdt() const
void initializeGhostedContainer(int, LinearObjContainer &loc) const
virtual Teuchos::MpiComm< int > getComm() const
virtual const Teuchos::RCP< Epetra_Map > getGhostedMap(int i) const
get the ghosted map from the matrix
Teuchos::RCP< Epetra_CrsMatrix > getEpetraMatrix(int i, int j) const
Teuchos::RCP< VectorType > get_f() const
virtual const Teuchos::RCP< Epetra_CrsGraph > buildGraph(int i, int j) const
void initializeContainer(int, LinearObjContainer &loc) const
virtual const Teuchos::RCP< Epetra_Map > buildColGhostedMap(int i) const
int VectorToMatrixMarketFile(const char *filename, const Epetra_Vector &A, const char *matrixName=0, const char *matrixDescription=0, bool writeHeader=true)
virtual const Teuchos::RCP< Epetra_Map > buildGhostedMap(int i) const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void initializeGhostedContainer_internal(int mem, ThyraObjContainer< double > &loc) const
virtual void set_A_th(const Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > &in)=0
Teuchos::RCP< VectorType > get_x() const
Teuchos::RCP< const Teuchos::Comm< int > > comm
virtual const Teuchos::RCP< Epetra_Import > getGhostedColImport(int i) const
get importer for converting an overalapped object to a "normal" object
virtual const Teuchos::RCP< Epetra_CrsGraph > getGraph(int i, int j) const
get the graph of the crs matrix
virtual void writeVector(const std::string &identifier, const LinearObjContainer &loc, int id) const
int MatrixMarketFileToVector(const char *filename, const Epetra_BlockMap &map, Epetra_Vector *&A)
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
Teuchos::RCP< Thyra::LinearOpBase< double > > getGhostedThyraMatrix() const
Get a Thyra operator.
virtual const Teuchos::RCP< Epetra_Map > getGhostedColMap(int i) const
get the ghosted map from the matrix
Teuchos::RCP< const panzer::BlockedDOFManager< int, int > > getGlobalIndexer() const
void buildGatherScatterEvaluators(const BuilderT &builder)
void initializeContainer_internal(int mem, ThyraObjContainer< double > &loc) const
int getBlockRowCount() const
how many block rows
const Teuchos::RCP< Epetra_Vector > get_f() const
Teuchos::RCP< Teuchos::MpiComm< int > > tComm_
Teuchos::RCP< const DOFManagerContainer > rowDOFManagerContainer_
virtual const Teuchos::RCP< Epetra_Export > getGhostedColExport(int j) const
get exporter for converting an overalapped object to a "normal" object
const Teuchos::RCP< Epetra_Vector > get_dxdt() const
BlockedEpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const UniqueGlobalIndexerBase > &gidProvider, bool useDiscreteAdjoint=false)
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)