43 #ifndef PANZER_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP 44 #define PANZER_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP 52 #include "Thyra_SpmdVectorBase.hpp" 54 #include "Epetra_MultiVector.h" 55 #include "Epetra_Vector.h" 56 #include "Epetra_CrsMatrix.h" 57 #include "Epetra_MpiComm.h" 70 template <
typename Traits,
typename LocalOrdinalT>
72 const Teuchos::RCP<
const UniqueGlobalIndexer<LocalOrdinalT,int> > & gidProvider,
73 bool useDiscreteAdjoint)
74 : comm_(
Teuchos::null), gidProvider_(gidProvider), rawMpiComm_(
comm->getRawMpiComm()), useDiscreteAdjoint_(useDiscreteAdjoint)
76 comm_ = Teuchos::rcp(
new Epetra_MpiComm(*rawMpiComm_));
77 hasColProvider_ = colGidProvider_!=Teuchos::null;
81 this->buildGatherScatterEvaluators(*
this);
84 template <
typename Traits,
typename LocalOrdinalT>
86 const Teuchos::RCP<
const UniqueGlobalIndexer<LocalOrdinalT,int> > & gidProvider,
87 const Teuchos::RCP<
const UniqueGlobalIndexer<LocalOrdinalT,int> > & colGidProvider,
88 bool useDiscreteAdjoint)
89 : comm_(
Teuchos::null), gidProvider_(gidProvider), colGidProvider_(colGidProvider), rawMpiComm_(
comm->getRawMpiComm()), useDiscreteAdjoint_(useDiscreteAdjoint)
91 comm_ = Teuchos::rcp(
new Epetra_MpiComm(*rawMpiComm_));
92 hasColProvider_ = colGidProvider_!=Teuchos::null;
96 this->buildGatherScatterEvaluators(*
this);
99 template <
typename Traits,
typename LocalOrdinalT>
101 : comm_(
Teuchos::null), gidProvider_(gidProvider), useDiscreteAdjoint_(false)
103 hasColProvider_ = colGidProvider_!=Teuchos::null;
107 this->buildGatherScatterEvaluators(*
this);
110 template <
typename Traits,
typename LocalOrdinalT>
111 EpetraLinearObjFactory<Traits,LocalOrdinalT>::~EpetraLinearObjFactory()
117 template <
typename Traits,
typename LocalOrdinalT>
119 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
120 readVector(
const std::string & identifier,LinearObjContainer & loc,
int id)
const 123 using Teuchos::rcp_dynamic_cast;
124 using Teuchos::dyn_cast;
126 EpetraLinearObjContainer & eloc = dyn_cast<EpetraLinearObjContainer>(loc);
129 RCP<Thyra::VectorBase<double> > vec;
131 case LinearObjContainer::X:
132 vec = eloc.get_x_th();
134 case LinearObjContainer::DxDt:
135 vec = eloc.get_dxdt_th();
137 case LinearObjContainer::F:
138 vec = eloc.get_f_th();
141 TEUCHOS_ASSERT(
false);
146 RCP<Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(),vec);
149 std::stringstream ss;
150 ss << identifier <<
".mm";
153 Epetra_Vector * ptr_ex = 0;
160 template <
typename Traits,
typename LocalOrdinalT>
162 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
163 writeVector(
const std::string & identifier,
const LinearObjContainer & loc,
int id)
const 166 using Teuchos::rcp_dynamic_cast;
167 using Teuchos::dyn_cast;
169 const EpetraLinearObjContainer & eloc = dyn_cast<
const EpetraLinearObjContainer>(loc);
172 RCP<const Thyra::VectorBase<double> > vec;
174 case LinearObjContainer::X:
175 vec = eloc.get_x_th();
177 case LinearObjContainer::DxDt:
178 vec = eloc.get_dxdt_th();
180 case LinearObjContainer::F:
181 vec = eloc.get_f_th();
184 TEUCHOS_ASSERT(
false);
189 RCP<const Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(),vec);
192 std::stringstream ss;
193 ss << identifier <<
".mm";
199 template <
typename Traits,
typename LocalOrdinalT>
200 Teuchos::RCP<LinearObjContainer> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildLinearObjContainer()
const 202 Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(
new EpetraLinearObjContainer(getColMap(),getMap()));
207 template <
typename Traits,
typename LocalOrdinalT>
208 Teuchos::RCP<LinearObjContainer> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedLinearObjContainer()
const 210 Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(
new EpetraLinearObjContainer(getGhostedColMap(),getGhostedMap()));
215 template <
typename Traits,
typename LocalOrdinalT>
216 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::globalToGhostContainer(
const LinearObjContainer & in,
217 LinearObjContainer & out,
int mem)
const 219 using Teuchos::is_null;
221 typedef LinearObjContainer LOC;
222 const EpetraLinearObjContainer & e_in = Teuchos::dyn_cast<
const EpetraLinearObjContainer>(in);
223 EpetraLinearObjContainer & e_out = Teuchos::dyn_cast<EpetraLinearObjContainer>(out);
227 if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
228 globalToGhostEpetraVector(*e_in.get_x(),*e_out.get_x(),
true);
230 if ( !is_null(e_in.get_dxdt()) && !is_null(e_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
231 globalToGhostEpetraVector(*e_in.get_dxdt(),*e_out.get_dxdt(),
true);
233 if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
234 globalToGhostEpetraVector(*e_in.get_f(),*e_out.get_f(),
false);
237 template <
typename Traits,
typename LocalOrdinalT>
238 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::ghostToGlobalContainer(
const LinearObjContainer & in,
239 LinearObjContainer & out,
int mem)
const 241 using Teuchos::is_null;
243 typedef LinearObjContainer LOC;
244 const EpetraLinearObjContainer & e_in = Teuchos::dyn_cast<
const EpetraLinearObjContainer>(in);
245 EpetraLinearObjContainer & e_out = Teuchos::dyn_cast<EpetraLinearObjContainer>(out);
249 if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
250 ghostToGlobalEpetraVector(*e_in.get_x(),*e_out.get_x(),
true);
252 if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
253 ghostToGlobalEpetraVector(*e_in.get_f(),*e_out.get_f(),
false);
255 if ( !is_null(e_in.get_A()) && !is_null(e_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
256 ghostToGlobalEpetraMatrix(*e_in.get_A(),*e_out.get_A());
259 template <
typename Traits,
typename LocalOrdinalT>
260 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::ghostToGlobalEpetraVector(
const Epetra_Vector & in,Epetra_Vector & out,
bool col)
const 265 RCP<Epetra_Export> exporter = col ? getGhostedColExport() : getGhostedExport();
266 TEUCHOS_ASSERT(out.PutScalar(0.0)==0);
268 int errCode = out.Export(in,*exporter,Add);
269 TEUCHOS_TEST_FOR_EXCEPTION(errCode!=0,std::runtime_error,
270 "Epetra_Vector::Export returned an error code of " << errCode <<
"!");
273 template <
typename Traits,
typename LocalOrdinalT>
274 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::ghostToGlobalEpetraMatrix(
const Epetra_CrsMatrix & in,Epetra_CrsMatrix & out)
const 279 RCP<Epetra_Export> exporter = getGhostedExport();
281 out.Export(in,*exporter,Add);
284 template <
typename Traits,
typename LocalOrdinalT>
285 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::globalToGhostEpetraVector(
const Epetra_Vector & in,Epetra_Vector & out,
bool col)
const 290 RCP<Epetra_Import> importer = col ? getGhostedColImport() : getGhostedImport();
292 out.Import(in,*importer,Insert);
299 template <
typename Traits,
typename LocalOrdinalT>
300 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::
301 adjustForDirichletConditions(
const LinearObjContainer & localBCRows,
302 const LinearObjContainer & globalBCRows,
303 LinearObjContainer & ghostedObjs,
304 bool zeroVectorRows,
bool adjustX)
const 307 const EpetraLinearObjContainer & e_localBCRows = Teuchos::dyn_cast<
const EpetraLinearObjContainer>(localBCRows);
308 const EpetraLinearObjContainer & e_globalBCRows = Teuchos::dyn_cast<
const EpetraLinearObjContainer>(globalBCRows);
309 EpetraLinearObjContainer & e_ghosted = Teuchos::dyn_cast<EpetraLinearObjContainer>(ghostedObjs);
311 TEUCHOS_ASSERT(!Teuchos::is_null(e_localBCRows.get_f()));
312 TEUCHOS_ASSERT(!Teuchos::is_null(e_globalBCRows.get_f()));
315 Teuchos::RCP<Epetra_CrsMatrix> A = e_ghosted.get_A();
316 Teuchos::RCP<Epetra_Vector> f = e_ghosted.get_f();
317 if(adjustX) f = e_ghosted.get_x();
319 const Epetra_Vector & local_bcs = *(e_localBCRows.get_f());
320 const Epetra_Vector & global_bcs = *(e_globalBCRows.get_f());
322 TEUCHOS_ASSERT(local_bcs.MyLength()==global_bcs.MyLength());
323 for(
int i=0;i<local_bcs.MyLength();i++) {
324 if(global_bcs[i]==0.0)
331 if(local_bcs[i]==0.0 || zeroVectorRows) {
336 if(!Teuchos::is_null(f))
338 if(!Teuchos::is_null(A)) {
339 A->ExtractMyRowView(i,numEntries,
values,indices);
340 for(
int c=0;c<numEntries;c++)
347 double scaleFactor = global_bcs[i];
350 if(!Teuchos::is_null(f))
351 (*f)[i] /= scaleFactor;
352 if(!Teuchos::is_null(A)) {
353 A->ExtractMyRowView(i,numEntries,
values,indices);
354 for(
int c=0;c<numEntries;c++)
361 template <
typename Traits,
typename LocalOrdinalT>
362 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::
363 applyDirichletBCs(
const LinearObjContainer & counter,
364 LinearObjContainer & result)
const 367 using Teuchos::rcp_dynamic_cast;
368 using Teuchos::dyn_cast;
370 const ThyraObjContainer<double> & th_counter = dyn_cast<
const ThyraObjContainer<double> >(counter);
371 ThyraObjContainer<double> & th_result = dyn_cast<ThyraObjContainer<double> >(result);
373 RCP<const Thyra::VectorBase<double> > count = th_counter.get_f_th();
374 RCP<const Thyra::VectorBase<double> > f_in = th_counter.get_f_th();
375 RCP<Thyra::VectorBase<double> > f_out = th_result.get_f_th();
377 Teuchos::ArrayRCP<const double> count_array,f_in_array;
378 Teuchos::ArrayRCP<double> f_out_array;
380 rcp_dynamic_cast<
const Thyra::SpmdVectorBase<double> >(count,
true)->getLocalData(Teuchos::ptrFromRef(count_array));
381 rcp_dynamic_cast<
const Thyra::SpmdVectorBase<double> >(f_in,
true)->getLocalData(Teuchos::ptrFromRef(f_in_array));
382 rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(f_out,
true)->getNonconstLocalData(Teuchos::ptrFromRef(f_out_array));
384 TEUCHOS_ASSERT(count_array.size()==f_in_array.size());
385 TEUCHOS_ASSERT(count_array.size()==f_out_array.size());
387 for(Teuchos::ArrayRCP<double>::size_type i=0;i<count_array.size();i++) {
388 if(count_array[i]!=0.0)
389 f_out_array[i] = f_in_array[i];
393 template <
typename Traits,
typename LocalOrdinalT>
394 Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
395 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
396 buildDomainContainer()
const 398 Teuchos::RCP<EpetraVector_ReadOnly_GlobalEvaluationData> vec_ged
399 = Teuchos::rcp(
new EpetraVector_ReadOnly_GlobalEvaluationData);
400 vec_ged->initialize(getGhostedImport(),getGhostedColMap(),getColMap());
405 template <
typename Traits,
typename LocalOrdinalT>
406 Teuchos::MpiComm<int> EpetraLinearObjFactory<Traits,LocalOrdinalT>::
409 return Teuchos::MpiComm<int>(Teuchos::opaqueWrapper(dynamic_cast<const Epetra_MpiComm &>(*getEpetraComm()).Comm()));
412 template <
typename Traits,
typename LocalOrdinalT>
413 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
414 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
415 getThyraDomainSpace()
const 417 if(domainSpace_ == Teuchos::null) {
421 domainSpace_ = getThyraRangeSpace();
423 domainSpace_ = Thyra::create_VectorSpace(getColMap());
429 template <
typename Traits,
typename LocalOrdinalT>
430 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
431 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
432 getThyraRangeSpace()
const 434 if(rangeSpace_ == Teuchos::null)
435 rangeSpace_ = Thyra::create_VectorSpace(getMap());
440 template <
typename Traits,
typename LocalOrdinalT>
441 Teuchos::RCP<Thyra::LinearOpBase<double> >
442 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
443 getThyraMatrix()
const 445 return Thyra::nonconstEpetraLinearOp(getEpetraMatrix());
451 template <
typename Traits,
typename LocalOrdinalT>
452 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeContainer(
int mem,LinearObjContainer & loc)
const 454 EpetraLinearObjContainer & eloc = Teuchos::dyn_cast<EpetraLinearObjContainer>(loc);
455 initializeContainer(mem,eloc);
458 template <
typename Traits,
typename LocalOrdinalT>
459 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeContainer(
int mem,EpetraLinearObjContainer & loc)
const 461 typedef EpetraLinearObjContainer ELOC;
465 if((mem & ELOC::X) == ELOC::X)
466 loc.set_x(getEpetraColVector());
468 if((mem & ELOC::DxDt) == ELOC::DxDt)
469 loc.set_dxdt(getEpetraColVector());
471 if((mem & ELOC::F) == ELOC::F)
472 loc.set_f(getEpetraVector());
474 if((mem & ELOC::Mat) == ELOC::Mat)
475 loc.set_A(getEpetraMatrix());
478 template <
typename Traits,
typename LocalOrdinalT>
479 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeGhostedContainer(
int mem,LinearObjContainer & loc)
const 481 EpetraLinearObjContainer & eloc = Teuchos::dyn_cast<EpetraLinearObjContainer>(loc);
482 initializeGhostedContainer(mem,eloc);
485 template <
typename Traits,
typename LocalOrdinalT>
486 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeGhostedContainer(
int mem,EpetraLinearObjContainer & loc)
const 488 typedef EpetraLinearObjContainer ELOC;
492 if((mem & ELOC::X) == ELOC::X)
493 loc.set_x(getGhostedEpetraColVector());
495 if((mem & ELOC::DxDt) == ELOC::DxDt)
496 loc.set_dxdt(getGhostedEpetraColVector());
498 if((mem & ELOC::F) == ELOC::F) {
499 loc.set_f(getGhostedEpetraVector());
500 loc.setRequiresDirichletAdjustment(
true);
503 if((mem & ELOC::Mat) == ELOC::Mat) {
504 loc.set_A(getGhostedEpetraMatrix());
505 loc.setRequiresDirichletAdjustment(
true);
513 template <
typename Traits,
typename LocalOrdinalT>
514 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getMap()
const 516 if(map_==Teuchos::null) map_ = buildMap();
522 template <
typename Traits,
typename LocalOrdinalT>
523 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getColMap()
const 525 if(cMap_==Teuchos::null) cMap_ = buildColMap();
530 template <
typename Traits,
typename LocalOrdinalT>
531 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedMap()
const 533 if(ghostedMap_==Teuchos::null) ghostedMap_ = buildGhostedMap();
538 template <
typename Traits,
typename LocalOrdinalT>
539 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedColMap()
const 541 if(cGhostedMap_==Teuchos::null) cGhostedMap_ = buildGhostedColMap();
547 template <
typename Traits,
typename LocalOrdinalT>
548 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGraph()
const 550 if(graph_==Teuchos::null) graph_ = buildGraph();
555 template <
typename Traits,
typename LocalOrdinalT>
556 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedGraph()
const 558 if(ghostedGraph_==Teuchos::null) ghostedGraph_ = buildGhostedGraph(
true);
560 return ghostedGraph_;
563 template <
typename Traits,
typename LocalOrdinalT>
564 const Teuchos::RCP<Epetra_Import> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedImport()
const 566 if(importer_==Teuchos::null)
567 importer_ = Teuchos::rcp(
new Epetra_Import(*getGhostedMap(),*getMap()));
572 template <
typename Traits,
typename LocalOrdinalT>
573 const Teuchos::RCP<Epetra_Import> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedColImport()
const 576 colImporter_ = getGhostedImport();
578 if(colImporter_==Teuchos::null)
579 colImporter_ = Teuchos::rcp(
new Epetra_Import(*getGhostedColMap(),*getColMap()));
584 template <
typename Traits,
typename LocalOrdinalT>
585 const Teuchos::RCP<Epetra_Export> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedExport()
const 587 if(exporter_==Teuchos::null)
588 exporter_ = Teuchos::rcp(
new Epetra_Export(*getGhostedMap(),*getMap()));
593 template <
typename Traits,
typename LocalOrdinalT>
594 const Teuchos::RCP<Epetra_Export> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedColExport()
const 597 colExporter_ = getGhostedExport();
599 if(colExporter_==Teuchos::null)
600 colExporter_ = Teuchos::rcp(
new Epetra_Export(*getGhostedColMap(),*getColMap()));
608 template <
typename Traits,
typename LocalOrdinalT>
609 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildMap()
const 611 Teuchos::RCP<Epetra_Map> map;
612 std::vector<int> indices;
615 gidProvider_->getOwnedIndices(indices);
617 map = Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
622 template <
typename Traits,
typename LocalOrdinalT>
623 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildColMap()
const 628 std::vector<int> indices;
631 colGidProvider_->getOwnedIndices(indices);
633 return Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
637 template <
typename Traits,
typename LocalOrdinalT>
638 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedMap()
const 640 std::vector<int> indices;
643 gidProvider_->getOwnedAndGhostedIndices(indices);
645 return Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
649 template <
typename Traits,
typename LocalOrdinalT>
650 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedColMap()
const 653 return buildGhostedMap();
655 std::vector<int> indices;
658 colGidProvider_->getOwnedAndGhostedIndices(indices);
660 return Teuchos::rcp(
new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
664 template <
typename Traits,
typename LocalOrdinalT>
665 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGraph()
const 672 RCP<Epetra_Map> rMap = getMap();
673 RCP<Epetra_Map> cMap = getColMap();
674 RCP<Epetra_CrsGraph> graph = rcp(
new Epetra_CrsGraph(Copy,*rMap,0));
676 RCP<Epetra_CrsGraph> oGraph = buildFilteredGhostedGraph();
681 RCP<Epetra_Export> exporter = getGhostedExport();
682 graph->Export( *oGraph, *exporter, Insert );
683 graph->FillComplete(*cMap,*rMap);
684 graph->OptimizeStorage();
688 template <
typename Traits,
typename LocalOrdinalT>
689 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedGraph(
bool optimizeStorage)
const 692 Teuchos::RCP<Epetra_Map> rMap = getGhostedMap();
693 Teuchos::RCP<Epetra_Map> cMap = getGhostedColMap();
694 Teuchos::RCP<Epetra_CrsGraph> graph = Teuchos::rcp(
new Epetra_CrsGraph(Copy,*rMap,*cMap,0));
696 std::vector<std::string> elementBlockIds;
697 gidProvider_->getElementBlockIds(elementBlockIds);
699 const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,int> >
700 colGidProvider = hasColProvider_ ? colGidProvider_ : gidProvider_;
701 const Teuchos::RCP<const ConnManagerBase<LocalOrdinalT> > conn_mgr = colGidProvider->getConnManagerBase();
702 const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors();
705 std::vector<std::string>::const_iterator blockItr;
706 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
707 std::string blockId = *blockItr;
710 const std::vector<LocalOrdinalT> & elements = gidProvider_->getElementBlock(blockId);
713 std::vector<int> gids;
714 std::vector<int> col_gids;
717 for(std::size_t i=0;i<elements.size();i++) {
718 gidProvider_->getElementGIDs(elements[i],gids);
720 colGidProvider->getElementGIDs(elements[i],col_gids);
722 const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[i]);
723 for (
typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
724 eit != aes.end(); ++eit) {
725 std::vector<int> other_col_gids;
726 colGidProvider->getElementGIDs(*eit, other_col_gids);
727 col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
731 for(std::size_t j=0;j<gids.size();j++)
732 graph->InsertGlobalIndices(gids[j],col_gids.size(),&col_gids[0]);
737 graph->FillComplete(*cMap,*rMap);
738 if(optimizeStorage) {
739 graph->OptimizeStorage();
745 template <
typename Traits,
typename LocalOrdinalT>
746 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildFilteredGhostedGraph()
const 749 using Teuchos::rcp_dynamic_cast;
752 RCP<const Filtered_UniqueGlobalIndexer<LocalOrdinalT,int> > filtered_ugi
753 = rcp_dynamic_cast<
const Filtered_UniqueGlobalIndexer<LocalOrdinalT,int> >(getDomainGlobalIndexer());
756 if(filtered_ugi==Teuchos::null)
757 return getGhostedGraph();
760 std::vector<int> ghostedActive;
761 filtered_ugi->getOwnedAndGhostedNotFilteredIndicator(ghostedActive);
765 Teuchos::RCP<Epetra_CrsGraph> filteredGraph = buildGhostedGraph(
false);
769 for(
int i=0;i<filteredGraph->NumMyRows();i++) {
770 std::vector<int> removedIndices;
773 TEUCHOS_ASSERT(filteredGraph->ExtractMyRowView(i,numIndices,indices)==0);
775 for(
int j=0;j<numIndices;j++) {
776 if(ghostedActive[indices[j]]==0)
777 removedIndices.push_back(indices[j]);
780 TEUCHOS_ASSERT(filteredGraph->RemoveMyIndices(i,Teuchos::as<int>(removedIndices.size()),&removedIndices[0])==0);
784 Teuchos::RCP<Epetra_Map> rMap = getGhostedMap();
785 Teuchos::RCP<Epetra_Map> cMap = getGhostedColMap();
787 TEUCHOS_ASSERT(filteredGraph->FillComplete(*cMap,*rMap)==0);
788 TEUCHOS_ASSERT(filteredGraph->OptimizeStorage()==0);
790 return filteredGraph;
793 template <
typename Traits,
typename LocalOrdinalT>
794 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedEpetraColVector()
const 796 Teuchos::RCP<const Epetra_Map> eMap = getGhostedColMap();
797 return Teuchos::rcp(
new Epetra_Vector(*eMap));
800 template <
typename Traits,
typename LocalOrdinalT>
801 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedEpetraVector()
const 803 Teuchos::RCP<const Epetra_Map> eMap = getGhostedMap();
804 return Teuchos::rcp(
new Epetra_Vector(*eMap));
807 template <
typename Traits,
typename LocalOrdinalT>
808 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraColVector()
const 810 Teuchos::RCP<const Epetra_Map> eMap = getColMap();
811 return Teuchos::rcp(
new Epetra_Vector(*eMap));
814 template <
typename Traits,
typename LocalOrdinalT>
815 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraVector()
const 817 Teuchos::RCP<const Epetra_Map> eMap = getMap();
818 return Teuchos::rcp(
new Epetra_Vector(*eMap));
821 template <
typename Traits,
typename LocalOrdinalT>
822 Teuchos::RCP<Epetra_CrsMatrix> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraMatrix()
const 824 Teuchos::RCP<Epetra_CrsGraph> eGraph = getGraph();
825 return Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *eGraph));
828 template <
typename Traits,
typename LocalOrdinalT>
829 Teuchos::RCP<Epetra_CrsMatrix> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedEpetraMatrix()
const 831 Teuchos::RCP<Epetra_CrsGraph> eGraph = getGhostedGraph();
832 return Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *eGraph));
835 template <
typename Traits,
typename LocalOrdinalT>
836 const Teuchos::RCP<const Epetra_Comm> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraComm()
const
BlockedEpetraLinearObjFactory< Traits, LocalOrdinalT > EpetraLinearObjFactory
std::vector< ScalarT > values
int VectorToMatrixMarketFile(const char *filename, const Epetra_Vector &A, const char *matrixName=0, const char *matrixDescription=0, bool writeHeader=true)
Teuchos::RCP< const Teuchos::Comm< int > > comm
int MatrixMarketFileToVector(const char *filename, const Epetra_BlockMap &map, Epetra_Vector *&A)