Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_BlockedTpetraLinearObjFactory_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Panzer: A partial differential equation assembly
4 // engine for strongly coupled complex multiphysics systems
5 //
6 // Copyright 2011 NTESS and the Panzer contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
11 #ifndef __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
12 #define __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
13 
14 // Panzer
16 #ifdef PANZER_HAVE_EPETRA_STACK
17 #include "Panzer_EpetraVector_Write_GlobalEvaluationData.hpp" // JMG: Remove this eventually.
18 #endif
20 #include "Panzer_GlobalIndexer.hpp"
21 
22 // Thyra
23 #include "Thyra_DefaultBlockedLinearOp.hpp"
24 #include "Thyra_DefaultProductVectorSpace.hpp"
25 #include "Thyra_SpmdVectorBase.hpp"
26 #include "Thyra_TpetraLinearOp.hpp"
27 #include "Thyra_TpetraThyraWrappers.hpp"
28 
29 // Tpetra
30 #include "Tpetra_CrsMatrix.hpp"
31 #include "Tpetra_MultiVector.hpp"
32 #include "Tpetra_Vector.hpp"
33 
34 namespace panzer {
35 
36 using Teuchos::RCP;
37 
38 // ************************************************************
39 // class BlockedTpetraLinearObjFactory
40 // ************************************************************
41 
42 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
45  const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
46  : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
47 {
48  for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
49  gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
50 
52 
53  // build and register the gather/scatter evaluators with
54  // the base class.
55  this->buildGatherScatterEvaluators(*this);
56 }
57 
58 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
61  const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
62  : gidProviders_(gidProviders), comm_(comm)
63 {
65 }
66 
67 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
70 { }
71 
72 // LinearObjectFactory functions
74 
75 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
79 {
80  std::vector<Teuchos::RCP<const MapType> > blockMaps;
81  std::size_t blockDim = gidProviders_.size();
82  for(std::size_t i=0;i<blockDim;i++)
83  blockMaps.push_back(getMap(i));
84 
85  Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
86  container->setMapsForBlocks(blockMaps);
87 
88  return container;
89 }
90 
91 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
95 {
96  std::vector<Teuchos::RCP<const MapType> > blockMaps;
97  std::size_t blockDim = gidProviders_.size();
98  for(std::size_t i=0;i<blockDim;i++)
99  blockMaps.push_back(getGhostedMap(i));
100 
101  Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
102  container->setMapsForBlocks(blockMaps);
103 
104  return container;
105 }
106 
107 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
110 {
111  using Teuchos::is_null;
112 
113  typedef LinearObjContainer LOC;
114 
115  const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
116  BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
117 
118  // Operations occur if the GLOBAL container has the correct targets!
119  // Users set the GLOBAL continer arguments
120  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
121  globalToGhostThyraVector(b_in.get_x(),b_out.get_x());
122 
123  if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
124  globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt());
125 
126  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
127  globalToGhostThyraVector(b_in.get_f(),b_out.get_f());
128 }
129 
130 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
133 {
134  using Teuchos::is_null;
135 
136  typedef LinearObjContainer LOC;
137 
138  const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
139  BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
140 
141  // Operations occur if the GLOBAL container has the correct targets!
142  // Users set the GLOBAL continer arguments
143  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
144  ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x());
145 
146  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
147  ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f());
148 
149  if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
150  ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
151 }
152 
153 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
156  const LinearObjContainer & globalBCRows,
157  LinearObjContainer & ghostedObjs,
158  bool zeroVectorRows, bool adjustX) const
159 {
160  using Teuchos::RCP;
161  using Teuchos::rcp_dynamic_cast;
162  using Thyra::LinearOpBase;
163  using Thyra::PhysicallyBlockedLinearOpBase;
164  using Thyra::VectorBase;
166 
167  std::size_t blockDim = gidProviders_.size();
168 
169  // first cast to block LOCs
170  const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
171  const BTLOC & b_globalBCRows = Teuchos::dyn_cast<const BTLOC>(globalBCRows);
172  BTLOC & b_ghosted = Teuchos::dyn_cast<BTLOC>(ghostedObjs);
173 
174  TEUCHOS_ASSERT(b_localBCRows.get_f()!=Teuchos::null);
175  TEUCHOS_ASSERT(b_globalBCRows.get_f()!=Teuchos::null);
176 
177  // cast each component as needed to their product form
178  RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.get_A());
179  RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_f());
180  RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.get_f(),true);
181  RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.get_f(),true);
182 
183  if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
184 
185  // sanity check!
186  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
187  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
188  if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
189  TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
190  TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
191 
192  for(std::size_t i=0;i<blockDim;i++) {
193  // grab epetra vector
194  RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),true)->getConstTpetraVector();
195  RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),true)->getConstTpetraVector();
196 
197  // pull out epetra values
198  RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
199  RCP<VectorType> t_f;
200  if(th_f==Teuchos::null)
201  t_f = Teuchos::null;
202  else
203  t_f = rcp_dynamic_cast<ThyraVector>(th_f,true)->getTpetraVector();
204 
205  for(std::size_t j=0;j<blockDim;j++) {
206  RCP<const MapType> map_i = getGhostedMap(i);
207  RCP<const MapType> map_j = getGhostedMap(j);
208 
209  // pull out epetra values
210  RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
211 
212  // don't do anyting if opertor is null
213  RCP<CrsMatrixType> t_A;
214  if(th_A==Teuchos::null)
215  t_A = Teuchos::null;
216  else {
217  RCP<OperatorType> t_A_op = rcp_dynamic_cast<ThyraLinearOp>(th_A,true)->getTpetraOperator();
218  t_A = rcp_dynamic_cast<CrsMatrixType>(t_A_op,true);
219  }
220 
221  // adjust Block operator
222  if(t_A!=Teuchos::null) {
223  t_A->resumeFill();
224  }
225 
226  adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
227 
228  if(t_A!=Teuchos::null) {
229  //t_A->fillComplete(map_j,map_i);
230  }
231 
232  t_f = Teuchos::null; // this is so we only adjust it once on the first pass
233  }
234  }
235 }
236 
237 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
240  const VectorType & global_bcs,
241  const Teuchos::Ptr<VectorType> & f,
242  const Teuchos::Ptr<CrsMatrixType> & A,
243  bool zeroVectorRows) const
244 {
245  if(f==Teuchos::null && A==Teuchos::null)
246  return;
247 
248  Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
249 
250  Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
251  Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
252 
253  TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
254  for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
255  if(global_bcs_array[i]==0.0)
256  continue;
257 
258  if(local_bcs_array[i]==0.0 || zeroVectorRows) {
259  // this boundary condition was NOT set by this processor
260 
261  // if they exist put 0.0 in each entry
262  if(!Teuchos::is_null(f))
263  f_array[i] = 0.0;
264  if(!Teuchos::is_null(A)) {
265  std::size_t numEntries = 0;
266  std::size_t sz = A->getNumEntriesInLocalRow(i);
267  typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
268  typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
269 
270  A->getLocalRowCopy(i,indices,values,numEntries);
271 
272  for(std::size_t c=0;c<numEntries;c++)
273  values(c) = 0.0;
274 
275  A->replaceLocalValues(i,indices,values);
276  }
277  }
278  else {
279  // this boundary condition was set by this processor
280 
281  ScalarT scaleFactor = global_bcs_array[i];
282 
283  // if they exist scale linear objects by scale factor
284  if(!Teuchos::is_null(f))
285  f_array[i] /= scaleFactor;
286  if(!Teuchos::is_null(A)) {
287  std::size_t numEntries = 0;
288  std::size_t sz = A->getNumEntriesInLocalRow(i);
289  typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
290  typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
291 
292  A->getLocalRowCopy(i,indices,values,numEntries);
293 
294  for(std::size_t c=0;c<numEntries;c++)
295  values(c) /= scaleFactor;
296 
297  A->replaceLocalValues(i,indices,values);
298  }
299  }
300  }
301 }
302 
303 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
305 applyDirichletBCs(const LinearObjContainer & /* counter */,
306  LinearObjContainer & /* result */) const
307 {
308  TEUCHOS_ASSERT(false); // not yet implemented
309 }
310 
312 //
313 // buildReadOnlyDomainContainer()
314 //
316 template<typename Traits, typename ScalarT, typename LocalOrdinalT,
317  typename GlobalOrdinalT, typename NodeT>
319 BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
320  NodeT>::
321 buildReadOnlyDomainContainer() const
322 {
323  using std::vector;
324  using Teuchos::RCP;
325  using Teuchos::rcp;
328  LocalOrdinalT, GlobalOrdinalT, NodeT>;
329  vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
330  for (int i(0); i < getBlockColCount(); ++i)
331  {
332  auto tvroged = rcp(new TVROGED);
333  tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
334  gedBlocks.push_back(tvroged);
335  }
336  auto ged = rcp(new BVROGED);
337  ged->initialize(getGhostedThyraDomainSpace(), getThyraDomainSpace(),
338  gedBlocks);
339  return ged;
340 } // end of buildReadOnlyDomainContainer()
341 
342 #ifdef PANZER_HAVE_EPETRA_STACK
343 //
345 // buildWriteDomainContainer()
346 //
348 template<typename Traits, typename ScalarT, typename LocalOrdinalT,
349  typename GlobalOrdinalT, typename NodeT>
351 BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
352  NodeT>::
353 buildWriteDomainContainer() const
354 {
355  using std::logic_error;
356  using Teuchos::rcp;
358  auto ged = rcp(new EVWGED);
359  TEUCHOS_TEST_FOR_EXCEPTION(true, logic_error, "NOT YET IMPLEMENTED")
360  return ged;
361 } // end of buildWriteDomainContainer()
362 #endif // PANZER_HAVE_EPETRA_STACK
363 
364 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
366 getComm() const
367 {
368  return *comm_;
369 }
370 
371 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
374 {
375  BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
376  initializeContainer(mem,bloc);
377 }
378 
379 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
382 {
383  BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
384  initializeGhostedContainer(mem,bloc);
385 }
386 
387 // Generic methods
389 
390 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
392 initializeContainer(int mem,BTLOC & loc) const
393 {
394  typedef LinearObjContainer LOC;
395 
396  loc.clear();
397 
398  if((mem & LOC::X) == LOC::X)
399  loc.set_x(getThyraDomainVector());
400 
401  if((mem & LOC::DxDt) == LOC::DxDt)
402  loc.set_dxdt(getThyraDomainVector());
403 
404  if((mem & LOC::F) == LOC::F)
405  loc.set_f(getThyraRangeVector());
406 
407  if((mem & LOC::Mat) == LOC::Mat)
408  loc.set_A(getThyraMatrix());
409 }
410 
411 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
413 initializeGhostedContainer(int mem,BTLOC & loc) const
414 {
415  typedef LinearObjContainer LOC;
416 
417  loc.clear();
418 
419  if((mem & LOC::X) == LOC::X)
420  loc.set_x(getGhostedThyraDomainVector());
421 
422  if((mem & LOC::DxDt) == LOC::DxDt)
423  loc.set_dxdt(getGhostedThyraDomainVector());
424 
425  if((mem & LOC::F) == LOC::F) {
426  loc.set_f(getGhostedThyraRangeVector());
428  }
429 
430  if((mem & LOC::Mat) == LOC::Mat) {
431  loc.set_A(getGhostedThyraMatrix());
433  }
434 }
435 
436 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
438 addExcludedPair(int rowBlock,int colBlock)
439 {
440  excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
441 }
442 
443 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
445 addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs)
446 {
447  for(std::size_t i=0;i<exPairs.size();i++)
448  excludedPairs_.insert(exPairs[i]);
449 }
450 
451 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
454 getGlobalIndexer(int i) const
455 {
456  return gidProviders_[i];
457 }
458 
459 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
461 makeRoomForBlocks(std::size_t blockCnt)
462 {
463  maps_.resize(blockCnt);
464  ghostedMaps_.resize(blockCnt);
465  importers_.resize(blockCnt);
466  exporters_.resize(blockCnt);
467 }
468 
469 // Thyra methods
471 
472 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
476 {
477  if(domainSpace_==Teuchos::null) {
478  // loop over all vectors and build the vector space
479  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
480  for(std::size_t i=0;i<gidProviders_.size();i++)
481  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
482 
483  domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
484  }
485 
486  return domainSpace_;
487 }
488 
489 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
493 {
494  if(rangeSpace_==Teuchos::null) {
495  // loop over all vectors and build the vector space
496  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
497  for(std::size_t i=0;i<gidProviders_.size();i++)
498  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
499 
500  rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
501  }
502 
503  return rangeSpace_;
504 }
505 
506 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
509 getThyraDomainSpace(int blk) const
510 {
511  if(domainSpace_==Teuchos::null) {
512  getThyraDomainSpace();
513  }
514 
515  return domainSpace_->getBlock(blk);
516 }
517 
518 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
521 getThyraRangeSpace(int blk) const
522 {
523  if(rangeSpace_==Teuchos::null) {
524  getThyraRangeSpace();
525  }
526 
527  return rangeSpace_->getBlock(blk);
528 }
529 
530 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
534 {
536  Thyra::createMember<ScalarT>(*getThyraDomainSpace());
537  Thyra::assign(vec.ptr(),0.0);
538 
540  for(std::size_t i=0;i<gidProviders_.size();i++) {
541  TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
542  Teuchos::as<int>(getMap(i)->getLocalNumElements()));
543  }
544 
545  return vec;
546 }
547 
548 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
552 {
554  Thyra::createMember<ScalarT>(*getThyraRangeSpace());
555  Thyra::assign(vec.ptr(),0.0);
556 
557  return vec;
558 }
559 
560 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
564 {
565  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
566 
567  // get the block dimension
568  std::size_t blockDim = gidProviders_.size();
569 
570  // this operator will be square
571  blockedOp->beginBlockFill(blockDim,blockDim);
572 
573  // loop over each block
574  for(std::size_t i=0;i<blockDim;i++) {
575  for(std::size_t j=0;j<blockDim;j++) {
576  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
577  // build (i,j) block matrix and add it to blocked operator
578  Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
579  blockedOp->setNonconstBlock(i,j,block);
580  }
581  }
582  }
583 
584  // all done
585  blockedOp->endBlockFill();
586 
587  return blockedOp;
588 }
589 
590 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
594 {
595  if(ghostedDomainSpace_==Teuchos::null) {
596  // loop over all vectors and build the vector space
597  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
598  for(std::size_t i=0;i<gidProviders_.size();i++)
599  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
600 
601  ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
602  }
603 
604  return ghostedDomainSpace_;
605 }
606 
607 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
611 {
612  if(ghostedRangeSpace_==Teuchos::null) {
613  // loop over all vectors and build the vector space
614  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
615  for(std::size_t i=0;i<gidProviders_.size();i++)
616  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
617 
618  ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
619  }
620 
621  return ghostedRangeSpace_;
622 }
623 
624 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
628 {
630  Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
631  Thyra::assign(vec.ptr(),0.0);
632 
633  return vec;
634 }
635 
636 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
640 {
642  Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
643  Thyra::assign(vec.ptr(),0.0);
644 
645  return vec;
646 }
647 
648 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
652 {
653  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
654 
655  // get the block dimension
656  std::size_t blockDim = gidProviders_.size();
657 
658  // this operator will be square
659  blockedOp->beginBlockFill(blockDim,blockDim);
660 
661  // loop over each block
662  for(std::size_t i=0;i<blockDim;i++) {
663  for(std::size_t j=0;j<blockDim;j++) {
664  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
665  // build (i,j) block matrix and add it to blocked operator
666  Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
667  blockedOp->setNonconstBlock(i,j,block);
668  }
669  }
670  }
671 
672  // all done
673  blockedOp->endBlockFill();
674 
675  return blockedOp;
676 }
677 
678 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
681  const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
682 {
683  using Teuchos::RCP;
684  using Teuchos::rcp_dynamic_cast;
686 
687  std::size_t blockDim = gidProviders_.size();
688 
689  // get product vectors
690  RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
691  RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
692 
693  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
694  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
695 
696  for(std::size_t i=0;i<blockDim;i++) {
697  // first get each Tpetra vector
698  RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
699  RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
700 
701  // use Tpetra to do global communication
702  ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
703  }
704 }
705 
706 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
709 {
710  using Teuchos::RCP;
711  using Teuchos::rcp_dynamic_cast;
712  using Teuchos::dyn_cast;
713  using Thyra::LinearOpBase;
714  using Thyra::PhysicallyBlockedLinearOpBase;
715 
716  std::size_t blockDim = gidProviders_.size();
717 
718  // get product vectors
719  const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
720  PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
721 
722  TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
723  TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
724  TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
725  TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
726 
727  for(std::size_t i=0;i<blockDim;i++) {
728  for(std::size_t j=0;j<blockDim;j++) {
729  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
730  // extract the blocks
731  RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
732  RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
733 
734  // sanity check
735  TEUCHOS_ASSERT(th_in!=Teuchos::null);
736  TEUCHOS_ASSERT(th_out!=Teuchos::null);
737 
738  // get the epetra version of the blocks
739  RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<const ThyraLinearOp>(th_in,true)->getConstTpetraOperator();
740  RCP<OperatorType> tp_op_out = rcp_dynamic_cast<ThyraLinearOp>(th_out,true)->getTpetraOperator();
741 
742  RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<const CrsMatrixType>(tp_op_in,true);
743  RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<CrsMatrixType>(tp_op_out,true);
744 
745  // use Tpetra to do global communication
746  ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
747  }
748  }
749  }
750 }
751 
752 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
755  const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
756 {
757  using Teuchos::RCP;
758  using Teuchos::rcp_dynamic_cast;
760 
761  std::size_t blockDim = gidProviders_.size();
762 
763  // get product vectors
764  RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
765  RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
766 
767  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
768  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
769 
770  for(std::size_t i=0;i<blockDim;i++) {
771  // first get each Tpetra vector
772  RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
773  RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
774 
775  // use Tpetra to do global communication
776  globalToGhostTpetraVector(i,*tp_in,*tp_out);
777  }
778 }
779 
780 // Tpetra methods
782 
783 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
785 ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out) const
786 {
787  using Teuchos::RCP;
788 
789  // do the global distribution
790  RCP<const ExportType> exporter = getGhostedExport(i);
791  out.putScalar(0.0);
792  out.doExport(in,*exporter,Tpetra::ADD);
793 }
794 
795 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
797 ghostToGlobalTpetraMatrix(int blockRow,const CrsMatrixType & in,CrsMatrixType & out) const
798 {
799  using Teuchos::RCP;
800 
801  RCP<const MapType> map_i = out.getRangeMap();
802  RCP<const MapType> map_j = out.getDomainMap();
803 
804  // do the global distribution
805  RCP<const ExportType> exporter = getGhostedExport(blockRow);
806 
807  out.resumeFill();
808  out.setAllToScalar(0.0);
809  out.doExport(in,*exporter,Tpetra::ADD);
810  out.fillComplete(map_j,map_i);
811 }
812 
813 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
815 globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out) const
816 {
817  using Teuchos::RCP;
818 
819  // do the global distribution
820  RCP<const ImportType> importer = getGhostedImport(i);
821  out.putScalar(0.0);
822  out.doImport(in,*importer,Tpetra::INSERT);
823 }
824 
825 // get the map from the matrix
826 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
829 getMap(int i) const
830 {
831  if(maps_[i]==Teuchos::null)
832  maps_[i] = buildTpetraMap(i);
833 
834  return maps_[i];
835 }
836 
837 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
840 getGhostedMap(int i) const
841 {
842  if(ghostedMaps_[i]==Teuchos::null)
843  ghostedMaps_[i] = buildTpetraGhostedMap(i);
844 
845  return ghostedMaps_[i];
846 }
847 
848 // get the graph of the crs matrix
849 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
852 getGraph(int i,int j) const
853 {
854  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
855 
856  typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
857  Teuchos::RCP<const CrsGraphType> graph;
858  if(itr==graphs_.end()) {
859  graph = buildTpetraGraph(i,j);
860  graphs_[std::make_pair(i,j)] = graph;
861  }
862  else
863  graph = itr->second;
864 
865  TEUCHOS_ASSERT(graph!=Teuchos::null);
866  return graph;
867 }
868 
869 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
872 getGhostedGraph(int i,int j) const
873 {
874  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
875 
876  typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
877  Teuchos::RCP<const CrsGraphType> ghostedGraph;
878  if(itr==ghostedGraphs_.end()) {
879  ghostedGraph = buildTpetraGhostedGraph(i,j);
880  ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
881  }
882  else
883  ghostedGraph = itr->second;
884 
885  TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
886  return ghostedGraph;
887 }
888 
889 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
892 getGhostedImport(int i) const
893 {
894  if(importers_[i]==Teuchos::null)
895  importers_[i] = Teuchos::rcp(new ImportType(getMap(i),getGhostedMap(i)));
896 
897  return importers_[i];
898 }
899 
900 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
903 getGhostedExport(int i) const
904 {
905  if(exporters_[i]==Teuchos::null)
906  exporters_[i] = Teuchos::rcp(new ExportType(getGhostedMap(i),getMap(i)));
907 
908  return exporters_[i];
909 }
910 
911 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
914 buildTpetraMap(int i) const
915 {
916  std::vector<GlobalOrdinalT> indices;
917 
918  // get the global indices
919  getGlobalIndexer(i)->getOwnedIndices(indices);
920 
922 }
923 
924 // build the ghosted map
925 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
929 {
930  std::vector<GlobalOrdinalT> indices;
931 
932  // get the global indices
933  getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
934 
936 }
937 
938 // get the graph of the crs matrix
939 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
942 buildTpetraGraph(int i,int j) const
943 {
944  using Teuchos::RCP;
945  using Teuchos::rcp;
946 
947  // build the map and allocate the space for the graph and
948  // grab the ghosted graph
949  RCP<const MapType> map_i = getMap(i);
950  RCP<const MapType> map_j = getMap(j);
951 
952  RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,0));
953  RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
954 
955  // perform the communication to finish building graph
956  RCP<const ExportType> exporter = getGhostedExport(i);
957  graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
958  graph->fillComplete(map_j,map_i);
959 
960  return graph;
961 }
962 
963 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
966 buildTpetraGhostedGraph(int i,int j) const
967 {
968  using Teuchos::RCP;
969  using Teuchos::rcp;
970 
971  // build the map and allocate the space for the graph and
972  // grab the ghosted graph
973  RCP<const MapType> map_i = getGhostedMap(i);
974  RCP<const MapType> map_j = getGhostedMap(j);
975 
976  std::vector<std::string> elementBlockIds;
977 
978  Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
979 
980  rowProvider = getGlobalIndexer(i);
981  colProvider = getGlobalIndexer(j);
982 
983  gidProviders_[0]->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
984  // same element blocks
985 
986  // Count number of entries in each row of graph; needed for graph constructor
987  std::vector<size_t> nEntriesPerRow(map_i->getLocalNumElements(), 0);
988  std::vector<std::string>::const_iterator blockItr;
989  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
990  std::string blockId = *blockItr;
991  // grab elements for this block
992  const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the
993  // same elements in each element block
994 
995  // get information about number of indicies
996  std::vector<GlobalOrdinalT> row_gids;
997  std::vector<GlobalOrdinalT> col_gids;
998 
999  // loop over the elemnts
1000  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1001 
1002  rowProvider->getElementGIDs(elements[elmt],row_gids);
1003  colProvider->getElementGIDs(elements[elmt],col_gids);
1004  for(std::size_t row=0;row<row_gids.size();row++) {
1005  LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
1006  nEntriesPerRow[lid] += col_gids.size();
1007  }
1008  }
1009  }
1010  Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
1011  RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,map_j, nEntriesPerRowView));
1012 
1013 
1014 
1015  // graph information about the mesh
1016  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1017  std::string blockId = *blockItr;
1018 
1019  // grab elements for this block
1020  const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the
1021  // same elements in each element block
1022 
1023  // get information about number of indicies
1024  std::vector<GlobalOrdinalT> row_gids;
1025  std::vector<GlobalOrdinalT> col_gids;
1026 
1027  // loop over the elemnts
1028  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1029 
1030  rowProvider->getElementGIDs(elements[elmt],row_gids);
1031  colProvider->getElementGIDs(elements[elmt],col_gids);
1032  for(std::size_t row=0;row<row_gids.size();row++)
1033  graph->insertGlobalIndices(row_gids[row],col_gids);
1034  }
1035  }
1036 
1037  // finish filling the graph: Make sure the colmap and row maps coincide to
1038  // minimize calls to LID lookups
1039  graph->fillComplete(getMap(j),getMap(i));
1040 
1041  return graph;
1042 }
1043 
1044 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1047 getTpetraMatrix(int i,int j) const
1048 {
1049  Teuchos::RCP<const MapType> map_i = getMap(i);
1050  Teuchos::RCP<const MapType> map_j = getMap(j);
1051 
1052  Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
1054  mat->fillComplete(map_j,map_i);
1055 
1056  return mat;
1057 }
1058 
1059 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1062 getGhostedTpetraMatrix(int i,int j) const
1063 {
1064  Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1065  Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1066 
1067  Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1069  mat->fillComplete(getMap(j),getMap(i));
1070 
1071  return mat;
1072 }
1073 
1074 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1078 {
1079  Teuchos::RCP<const MapType> tMap = getMap(i);
1080  return Teuchos::rcp(new VectorType(tMap));
1081 }
1082 
1083 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1087 {
1088  Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1089  return Teuchos::rcp(new VectorType(tMap));
1090 }
1091 
1092 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1096 {
1097  Teuchos::RCP<const MapType> tMap = getMap(i);
1098  return Teuchos::rcp(new VectorType(tMap));
1099 }
1100 
1101 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1105 {
1106  Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1107  return Teuchos::rcp(new VectorType(tMap));
1108 }
1109 
1110 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1111 int
1114 {
1115  return gidProviders_.size();
1116 }
1117 
1118 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1119 int
1122 {
1123  return gidProviders_.size();
1124 }
1125 
1126 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1129 {
1130  BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1131  if(tloc.get_A()!=Teuchos::null)
1132  tloc.beginFill();
1133 }
1134 
1135 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1138 {
1139  BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1140  if(tloc.get_A()!=Teuchos::null)
1141  tloc.endFill();
1142 }
1143 
1144 }
1145 
1146 #endif // __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
void initializeContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
Thyra::TpetraVector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > ThyraVector
bool is_null(const boost::shared_ptr< T > &p)
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a &quot;normal&quot; object
Thyra::TpetraLinearOp< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > ThyraLinearOp
void set_dxdt(const Teuchos::RCP< VectorType > &in)
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
bool is_null(const std::shared_ptr< T > &p)
Teuchos::RCP< VectorType > getTpetraDomainVector(int i) const
void initializeGhostedContainer(int, LinearObjContainer &loc) const
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraDomainVector() const
Get a domain vector.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::vector< Teuchos::RCP< const GlobalIndexer > > gidProviders_
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraDomainVector() const
Get a domain vector.
T_To & dyn_cast(T_From &from)
Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > ImportType
Teuchos::RCP< VectorType > getGhostedTpetraRangeVector(int i) const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
Teuchos::RCP< const panzer::BlockedDOFManager > getGlobalIndexer() const
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraRangeVector() const
Get a range vector.
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
virtual Teuchos::RCP< const ImportType > getGhostedImport(int i) const
get importer for converting an overalapped object to a &quot;normal&quot; object
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range vector space (f)
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void ghostToGlobalTpetraMatrix(int blockRow, const CrsMatrixType &in, CrsMatrixType &out) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
virtual Teuchos::RCP< const MapType > buildTpetraMap(int i) const
void globalToGhostTpetraVector(int i, const VectorType &in, VectorType &out) const
virtual Teuchos::RCP< const MapType > buildTpetraGhostedMap(int i) const
BlockedTpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const BlockedDOFManager > &gidProvider)
Ptr< T > ptr() const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< ScalarT > &in, Thyra::LinearOpBase< ScalarT > &out) const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
Teuchos::RCP< Thyra::BlockedLinearOpBase< ScalarT > > getGhostedThyraMatrix() const
Get a Thyra operator.
void makeRoomForBlocks(std::size_t blockCnt)
Allocate the space in the std::vector objects so we can fill with appropriate Tpetra data...
virtual Teuchos::RCP< const MapType > getGhostedMap(int i) const
get the ghosted map from the matrix
Teuchos::RCP< VectorType > getGhostedTpetraDomainVector(int i) const
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
This class provides a boundary exchange communication mechanism for vectors.
Teuchos::RCP< VectorType > getTpetraRangeVector(int i) const
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
virtual void endFill(LinearObjContainer &loc) const
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
virtual Teuchos::RCP< const CrsGraphType > getGraph(int i, int j) const
get the graph of the crs matrix
virtual Teuchos::RCP< const CrsGraphType > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGraph(int i, int j) const
virtual void beginFill(LinearObjContainer &loc) const
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGhostedGraph(int i, int j) const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
virtual void getElementGIDs(panzer::LocalOrdinal localElmtId, std::vector< panzer::GlobalOrdinal > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
void ghostToGlobalTpetraVector(int i, const VectorType &in, VectorType &out) const
void buildGatherScatterEvaluators(const BuilderT &builder)
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraRangeVector() const
Get a range vector.
#define TEUCHOS_ASSERT(assertion_test)
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a Thyra operator.
void set_A(const Teuchos::RCP< CrsMatrixType > &in)