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 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
44 #define __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
45 
46 // Panzer
48 #include "Panzer_EpetraVector_Write_GlobalEvaluationData.hpp" // JMG: Remove this eventually.
50 #include "Panzer_GlobalIndexer.hpp"
51 
52 // Thyra
53 #include "Thyra_DefaultBlockedLinearOp.hpp"
54 #include "Thyra_DefaultProductVectorSpace.hpp"
55 #include "Thyra_SpmdVectorBase.hpp"
56 #include "Thyra_TpetraLinearOp.hpp"
57 #include "Thyra_TpetraThyraWrappers.hpp"
58 
59 // Tpetra
60 #include "Tpetra_CrsMatrix.hpp"
61 #include "Tpetra_MultiVector.hpp"
62 #include "Tpetra_Vector.hpp"
63 
64 namespace panzer {
65 
66 using Teuchos::RCP;
67 
68 // ************************************************************
69 // class BlockedTpetraLinearObjFactory
70 // ************************************************************
71 
72 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
75  const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
76  : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
77 {
78  for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
79  gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
80 
82 
83  // build and register the gather/scatter evaluators with
84  // the base class.
85  this->buildGatherScatterEvaluators(*this);
86 }
87 
88 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
91  const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
92  : gidProviders_(gidProviders), comm_(comm)
93 {
95 }
96 
97 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
100 { }
101 
102 // LinearObjectFactory functions
104 
105 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
109 {
110  std::vector<Teuchos::RCP<const MapType> > blockMaps;
111  std::size_t blockDim = gidProviders_.size();
112  for(std::size_t i=0;i<blockDim;i++)
113  blockMaps.push_back(getMap(i));
114 
115  Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
116  container->setMapsForBlocks(blockMaps);
117 
118  return container;
119 }
120 
121 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
125 {
126  std::vector<Teuchos::RCP<const MapType> > blockMaps;
127  std::size_t blockDim = gidProviders_.size();
128  for(std::size_t i=0;i<blockDim;i++)
129  blockMaps.push_back(getGhostedMap(i));
130 
131  Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
132  container->setMapsForBlocks(blockMaps);
133 
134  return container;
135 }
136 
137 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
140 {
141  using Teuchos::is_null;
142 
143  typedef LinearObjContainer LOC;
144 
145  const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
146  BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
147 
148  // Operations occur if the GLOBAL container has the correct targets!
149  // Users set the GLOBAL continer arguments
150  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
151  globalToGhostThyraVector(b_in.get_x(),b_out.get_x());
152 
153  if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
154  globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt());
155 
156  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
157  globalToGhostThyraVector(b_in.get_f(),b_out.get_f());
158 }
159 
160 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
163 {
164  using Teuchos::is_null;
165 
166  typedef LinearObjContainer LOC;
167 
168  const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
169  BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
170 
171  // Operations occur if the GLOBAL container has the correct targets!
172  // Users set the GLOBAL continer arguments
173  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
174  ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x());
175 
176  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
177  ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f());
178 
179  if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
180  ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
181 }
182 
183 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
186  const LinearObjContainer & globalBCRows,
187  LinearObjContainer & ghostedObjs,
188  bool zeroVectorRows, bool adjustX) const
189 {
190  using Teuchos::RCP;
191  using Teuchos::rcp_dynamic_cast;
192  using Thyra::LinearOpBase;
193  using Thyra::PhysicallyBlockedLinearOpBase;
194  using Thyra::VectorBase;
196 
197  std::size_t blockDim = gidProviders_.size();
198 
199  // first cast to block LOCs
200  const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
201  const BTLOC & b_globalBCRows = Teuchos::dyn_cast<const BTLOC>(globalBCRows);
202  BTLOC & b_ghosted = Teuchos::dyn_cast<BTLOC>(ghostedObjs);
203 
204  TEUCHOS_ASSERT(b_localBCRows.get_f()!=Teuchos::null);
205  TEUCHOS_ASSERT(b_globalBCRows.get_f()!=Teuchos::null);
206 
207  // cast each component as needed to their product form
208  RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.get_A());
209  RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_f());
210  RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.get_f(),true);
211  RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.get_f(),true);
212 
213  if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
214 
215  // sanity check!
216  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
217  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
218  if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
219  TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
220  TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
221 
222  for(std::size_t i=0;i<blockDim;i++) {
223  // grab epetra vector
224  RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),true)->getConstTpetraVector();
225  RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),true)->getConstTpetraVector();
226 
227  // pull out epetra values
228  RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
229  RCP<VectorType> t_f;
230  if(th_f==Teuchos::null)
231  t_f = Teuchos::null;
232  else
233  t_f = rcp_dynamic_cast<ThyraVector>(th_f,true)->getTpetraVector();
234 
235  for(std::size_t j=0;j<blockDim;j++) {
236  RCP<const MapType> map_i = getGhostedMap(i);
237  RCP<const MapType> map_j = getGhostedMap(j);
238 
239  // pull out epetra values
240  RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
241 
242  // don't do anyting if opertor is null
243  RCP<CrsMatrixType> t_A;
244  if(th_A==Teuchos::null)
245  t_A = Teuchos::null;
246  else {
247  RCP<OperatorType> t_A_op = rcp_dynamic_cast<ThyraLinearOp>(th_A,true)->getTpetraOperator();
248  t_A = rcp_dynamic_cast<CrsMatrixType>(t_A_op,true);
249  }
250 
251  // adjust Block operator
252  if(t_A!=Teuchos::null) {
253  t_A->resumeFill();
254  }
255 
256  adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
257 
258  if(t_A!=Teuchos::null) {
259  //t_A->fillComplete(map_j,map_i);
260  }
261 
262  t_f = Teuchos::null; // this is so we only adjust it once on the first pass
263  }
264  }
265 }
266 
267 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
270  const VectorType & global_bcs,
271  const Teuchos::Ptr<VectorType> & f,
272  const Teuchos::Ptr<CrsMatrixType> & A,
273  bool zeroVectorRows) const
274 {
275  if(f==Teuchos::null && A==Teuchos::null)
276  return;
277 
278  Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
279 
280  Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
281  Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
282 
283  TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
284  for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
285  if(global_bcs_array[i]==0.0)
286  continue;
287 
288  if(local_bcs_array[i]==0.0 || zeroVectorRows) {
289  // this boundary condition was NOT set by this processor
290 
291  // if they exist put 0.0 in each entry
292  if(!Teuchos::is_null(f))
293  f_array[i] = 0.0;
294  if(!Teuchos::is_null(A)) {
295  std::size_t numEntries = 0;
296  std::size_t sz = A->getNumEntriesInLocalRow(i);
297  Teuchos::Array<LocalOrdinalT> indices(sz);
298  Teuchos::Array<ScalarT> values(sz);
299 
300  A->getLocalRowCopy(i,indices,values,numEntries);
301 
302  for(std::size_t c=0;c<numEntries;c++)
303  values[c] = 0.0;
304 
305  A->replaceLocalValues(i,indices,values);
306  }
307  }
308  else {
309  // this boundary condition was set by this processor
310 
311  ScalarT scaleFactor = global_bcs_array[i];
312 
313  // if they exist scale linear objects by scale factor
314  if(!Teuchos::is_null(f))
315  f_array[i] /= scaleFactor;
316  if(!Teuchos::is_null(A)) {
317  std::size_t numEntries = 0;
318  std::size_t sz = A->getNumEntriesInLocalRow(i);
319  Teuchos::Array<LocalOrdinalT> indices(sz);
320  Teuchos::Array<ScalarT> values(sz);
321 
322  A->getLocalRowCopy(i,indices,values,numEntries);
323 
324  for(std::size_t c=0;c<numEntries;c++)
325  values[c] /= scaleFactor;
326 
327  A->replaceLocalValues(i,indices,values);
328  }
329  }
330  }
331 }
332 
333 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
335 applyDirichletBCs(const LinearObjContainer & /* counter */,
336  LinearObjContainer & /* result */) const
337 {
338  TEUCHOS_ASSERT(false); // not yet implemented
339 }
340 
342 //
343 // buildReadOnlyDomainContainer()
344 //
346 template<typename Traits, typename ScalarT, typename LocalOrdinalT,
347  typename GlobalOrdinalT, typename NodeT>
349 BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
350  NodeT>::
351 buildReadOnlyDomainContainer() const
352 {
353  using std::vector;
354  using Teuchos::RCP;
355  using Teuchos::rcp;
358  LocalOrdinalT, GlobalOrdinalT, NodeT>;
359  vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
360  for (int i(0); i < getBlockColCount(); ++i)
361  {
362  auto tvroged = rcp(new TVROGED);
363  tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
364  gedBlocks.push_back(tvroged);
365  }
366  auto ged = rcp(new BVROGED);
367  ged->initialize(getGhostedThyraDomainSpace(), getThyraDomainSpace(),
368  gedBlocks);
369  return ged;
370 } // end of buildReadOnlyDomainContainer()
371 
373 //
374 // buildWriteDomainContainer()
375 //
377 template<typename Traits, typename ScalarT, typename LocalOrdinalT,
378  typename GlobalOrdinalT, typename NodeT>
380 BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
381  NodeT>::
382 buildWriteDomainContainer() const
383 {
384  using std::logic_error;
385  using Teuchos::rcp;
387  auto ged = rcp(new EVWGED);
388  TEUCHOS_TEST_FOR_EXCEPTION(true, logic_error, "NOT YET IMPLEMENTED")
389  return ged;
390 } // end of buildWriteDomainContainer()
391 
392 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
394 getComm() const
395 {
396  return *comm_;
397 }
398 
399 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
402 {
403  BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
404  initializeContainer(mem,bloc);
405 }
406 
407 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
410 {
411  BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
412  initializeGhostedContainer(mem,bloc);
413 }
414 
415 // Generic methods
417 
418 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
420 initializeContainer(int mem,BTLOC & loc) const
421 {
422  typedef LinearObjContainer LOC;
423 
424  loc.clear();
425 
426  if((mem & LOC::X) == LOC::X)
427  loc.set_x(getThyraDomainVector());
428 
429  if((mem & LOC::DxDt) == LOC::DxDt)
430  loc.set_dxdt(getThyraDomainVector());
431 
432  if((mem & LOC::F) == LOC::F)
433  loc.set_f(getThyraRangeVector());
434 
435  if((mem & LOC::Mat) == LOC::Mat)
436  loc.set_A(getThyraMatrix());
437 }
438 
439 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
441 initializeGhostedContainer(int mem,BTLOC & loc) const
442 {
443  typedef LinearObjContainer LOC;
444 
445  loc.clear();
446 
447  if((mem & LOC::X) == LOC::X)
448  loc.set_x(getGhostedThyraDomainVector());
449 
450  if((mem & LOC::DxDt) == LOC::DxDt)
451  loc.set_dxdt(getGhostedThyraDomainVector());
452 
453  if((mem & LOC::F) == LOC::F) {
454  loc.set_f(getGhostedThyraRangeVector());
456  }
457 
458  if((mem & LOC::Mat) == LOC::Mat) {
459  loc.set_A(getGhostedThyraMatrix());
461  }
462 }
463 
464 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
466 addExcludedPair(int rowBlock,int colBlock)
467 {
468  excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
469 }
470 
471 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
473 addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs)
474 {
475  for(std::size_t i=0;i<exPairs.size();i++)
476  excludedPairs_.insert(exPairs[i]);
477 }
478 
479 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
482 getGlobalIndexer(int i) const
483 {
484  return gidProviders_[i];
485 }
486 
487 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
489 makeRoomForBlocks(std::size_t blockCnt)
490 {
491  maps_.resize(blockCnt);
492  ghostedMaps_.resize(blockCnt);
493  importers_.resize(blockCnt);
494  exporters_.resize(blockCnt);
495 }
496 
497 // Thyra methods
499 
500 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
504 {
505  if(domainSpace_==Teuchos::null) {
506  // loop over all vectors and build the vector space
507  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
508  for(std::size_t i=0;i<gidProviders_.size();i++)
509  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
510 
511  domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
512  }
513 
514  return domainSpace_;
515 }
516 
517 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
521 {
522  if(rangeSpace_==Teuchos::null) {
523  // loop over all vectors and build the vector space
524  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
525  for(std::size_t i=0;i<gidProviders_.size();i++)
526  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
527 
528  rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
529  }
530 
531  return rangeSpace_;
532 }
533 
534 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
537 getThyraDomainSpace(int blk) const
538 {
539  if(domainSpace_==Teuchos::null) {
540  getThyraDomainSpace();
541  }
542 
543  return domainSpace_->getBlock(blk);
544 }
545 
546 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
549 getThyraRangeSpace(int blk) const
550 {
551  if(rangeSpace_==Teuchos::null) {
552  getThyraRangeSpace();
553  }
554 
555  return rangeSpace_->getBlock(blk);
556 }
557 
558 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
562 {
564  Thyra::createMember<ScalarT>(*getThyraDomainSpace());
565  Thyra::assign(vec.ptr(),0.0);
566 
568  for(std::size_t i=0;i<gidProviders_.size();i++) {
569  TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
570  Teuchos::as<int>(getMap(i)->getNodeNumElements()));
571  }
572 
573  return vec;
574 }
575 
576 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
580 {
582  Thyra::createMember<ScalarT>(*getThyraRangeSpace());
583  Thyra::assign(vec.ptr(),0.0);
584 
585  return vec;
586 }
587 
588 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
592 {
593  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
594 
595  // get the block dimension
596  std::size_t blockDim = gidProviders_.size();
597 
598  // this operator will be square
599  blockedOp->beginBlockFill(blockDim,blockDim);
600 
601  // loop over each block
602  for(std::size_t i=0;i<blockDim;i++) {
603  for(std::size_t j=0;j<blockDim;j++) {
604  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
605  // build (i,j) block matrix and add it to blocked operator
606  Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
607  blockedOp->setNonconstBlock(i,j,block);
608  }
609  }
610  }
611 
612  // all done
613  blockedOp->endBlockFill();
614 
615  return blockedOp;
616 }
617 
618 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
622 {
623  if(ghostedDomainSpace_==Teuchos::null) {
624  // loop over all vectors and build the vector space
625  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
626  for(std::size_t i=0;i<gidProviders_.size();i++)
627  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
628 
629  ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
630  }
631 
632  return ghostedDomainSpace_;
633 }
634 
635 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
639 {
640  if(ghostedRangeSpace_==Teuchos::null) {
641  // loop over all vectors and build the vector space
642  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
643  for(std::size_t i=0;i<gidProviders_.size();i++)
644  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
645 
646  ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
647  }
648 
649  return ghostedRangeSpace_;
650 }
651 
652 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
656 {
658  Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
659  Thyra::assign(vec.ptr(),0.0);
660 
661  return vec;
662 }
663 
664 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
668 {
670  Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
671  Thyra::assign(vec.ptr(),0.0);
672 
673  return vec;
674 }
675 
676 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
680 {
681  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
682 
683  // get the block dimension
684  std::size_t blockDim = gidProviders_.size();
685 
686  // this operator will be square
687  blockedOp->beginBlockFill(blockDim,blockDim);
688 
689  // loop over each block
690  for(std::size_t i=0;i<blockDim;i++) {
691  for(std::size_t j=0;j<blockDim;j++) {
692  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
693  // build (i,j) block matrix and add it to blocked operator
694  Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
695  blockedOp->setNonconstBlock(i,j,block);
696  }
697  }
698  }
699 
700  // all done
701  blockedOp->endBlockFill();
702 
703  return blockedOp;
704 }
705 
706 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
709  const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
710 {
711  using Teuchos::RCP;
712  using Teuchos::rcp_dynamic_cast;
714 
715  std::size_t blockDim = gidProviders_.size();
716 
717  // get product vectors
718  RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
719  RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
720 
721  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
722  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
723 
724  for(std::size_t i=0;i<blockDim;i++) {
725  // first get each Tpetra vector
726  RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
727  RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
728 
729  // use Tpetra to do global communication
730  ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
731  }
732 }
733 
734 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
737 {
738  using Teuchos::RCP;
739  using Teuchos::rcp_dynamic_cast;
740  using Teuchos::dyn_cast;
741  using Thyra::LinearOpBase;
742  using Thyra::PhysicallyBlockedLinearOpBase;
743 
744  std::size_t blockDim = gidProviders_.size();
745 
746  // get product vectors
747  const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
748  PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
749 
750  TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
751  TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
752  TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
753  TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
754 
755  for(std::size_t i=0;i<blockDim;i++) {
756  for(std::size_t j=0;j<blockDim;j++) {
757  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
758  // extract the blocks
759  RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
760  RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
761 
762  // sanity check
763  TEUCHOS_ASSERT(th_in!=Teuchos::null);
764  TEUCHOS_ASSERT(th_out!=Teuchos::null);
765 
766  // get the epetra version of the blocks
767  RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<const ThyraLinearOp>(th_in,true)->getConstTpetraOperator();
768  RCP<OperatorType> tp_op_out = rcp_dynamic_cast<ThyraLinearOp>(th_out,true)->getTpetraOperator();
769 
770  RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<const CrsMatrixType>(tp_op_in,true);
771  RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<CrsMatrixType>(tp_op_out,true);
772 
773  // use Tpetra to do global communication
774  ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
775  }
776  }
777  }
778 }
779 
780 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
783  const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
784 {
785  using Teuchos::RCP;
786  using Teuchos::rcp_dynamic_cast;
788 
789  std::size_t blockDim = gidProviders_.size();
790 
791  // get product vectors
792  RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
793  RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
794 
795  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
796  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
797 
798  for(std::size_t i=0;i<blockDim;i++) {
799  // first get each Tpetra vector
800  RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
801  RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
802 
803  // use Tpetra to do global communication
804  globalToGhostTpetraVector(i,*tp_in,*tp_out);
805  }
806 }
807 
808 // Tpetra methods
810 
811 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
813 ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out) const
814 {
815  using Teuchos::RCP;
816 
817  // do the global distribution
818  RCP<const ExportType> exporter = getGhostedExport(i);
819  out.putScalar(0.0);
820  out.doExport(in,*exporter,Tpetra::ADD);
821 }
822 
823 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
825 ghostToGlobalTpetraMatrix(int blockRow,const CrsMatrixType & in,CrsMatrixType & out) const
826 {
827  using Teuchos::RCP;
828 
829  RCP<const MapType> map_i = out.getRangeMap();
830  RCP<const MapType> map_j = out.getDomainMap();
831 
832  // do the global distribution
833  RCP<const ExportType> exporter = getGhostedExport(blockRow);
834 
835  out.resumeFill();
836  out.setAllToScalar(0.0);
837  out.doExport(in,*exporter,Tpetra::ADD);
838  out.fillComplete(map_j,map_i);
839 }
840 
841 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
843 globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out) const
844 {
845  using Teuchos::RCP;
846 
847  // do the global distribution
848  RCP<const ImportType> importer = getGhostedImport(i);
849  out.putScalar(0.0);
850  out.doImport(in,*importer,Tpetra::INSERT);
851 }
852 
853 // get the map from the matrix
854 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
857 getMap(int i) const
858 {
859  if(maps_[i]==Teuchos::null)
860  maps_[i] = buildTpetraMap(i);
861 
862  return maps_[i];
863 }
864 
865 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
868 getGhostedMap(int i) const
869 {
870  if(ghostedMaps_[i]==Teuchos::null)
871  ghostedMaps_[i] = buildTpetraGhostedMap(i);
872 
873  return ghostedMaps_[i];
874 }
875 
876 // get the graph of the crs matrix
877 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
880 getGraph(int i,int j) const
881 {
882  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
883 
884  typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
885  Teuchos::RCP<const CrsGraphType> graph;
886  if(itr==graphs_.end()) {
887  graph = buildTpetraGraph(i,j);
888  graphs_[std::make_pair(i,j)] = graph;
889  }
890  else
891  graph = itr->second;
892 
893  TEUCHOS_ASSERT(graph!=Teuchos::null);
894  return graph;
895 }
896 
897 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
900 getGhostedGraph(int i,int j) const
901 {
902  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
903 
904  typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
905  Teuchos::RCP<const CrsGraphType> ghostedGraph;
906  if(itr==ghostedGraphs_.end()) {
907  ghostedGraph = buildTpetraGhostedGraph(i,j);
908  ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
909  }
910  else
911  ghostedGraph = itr->second;
912 
913  TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
914  return ghostedGraph;
915 }
916 
917 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
920 getGhostedImport(int i) const
921 {
922  if(importers_[i]==Teuchos::null)
923  importers_[i] = Teuchos::rcp(new ImportType(getMap(i),getGhostedMap(i)));
924 
925  return importers_[i];
926 }
927 
928 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
931 getGhostedExport(int i) const
932 {
933  if(exporters_[i]==Teuchos::null)
934  exporters_[i] = Teuchos::rcp(new ExportType(getGhostedMap(i),getMap(i)));
935 
936  return exporters_[i];
937 }
938 
939 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
942 buildTpetraMap(int i) const
943 {
944  std::vector<GlobalOrdinalT> indices;
945 
946  // get the global indices
947  getGlobalIndexer(i)->getOwnedIndices(indices);
948 
950 }
951 
952 // build the ghosted map
953 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
957 {
958  std::vector<GlobalOrdinalT> indices;
959 
960  // get the global indices
961  getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
962 
964 }
965 
966 // get the graph of the crs matrix
967 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
970 buildTpetraGraph(int i,int j) const
971 {
972  using Teuchos::RCP;
973  using Teuchos::rcp;
974 
975  // build the map and allocate the space for the graph and
976  // grab the ghosted graph
977  RCP<const MapType> map_i = getMap(i);
978  RCP<const MapType> map_j = getMap(j);
979 
980  RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,0));
981  RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
982 
983  // perform the communication to finish building graph
984  RCP<const ExportType> exporter = getGhostedExport(i);
985  graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
986  graph->fillComplete(map_j,map_i);
987 
988  return graph;
989 }
990 
991 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
994 buildTpetraGhostedGraph(int i,int j) const
995 {
996  using Teuchos::RCP;
997  using Teuchos::rcp;
998 
999  // build the map and allocate the space for the graph and
1000  // grab the ghosted graph
1001  RCP<const MapType> map_i = getGhostedMap(i);
1002  RCP<const MapType> map_j = getGhostedMap(j);
1003 
1004  std::vector<std::string> elementBlockIds;
1005 
1006  Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
1007 
1008  rowProvider = getGlobalIndexer(i);
1009  colProvider = getGlobalIndexer(j);
1010 
1011  gidProviders_[0]->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
1012  // same element blocks
1013 
1014  // Count number of entries in each row of graph; needed for graph constructor
1015  std::vector<size_t> nEntriesPerRow(map_i->getNodeNumElements(), 0);
1016  std::vector<std::string>::const_iterator blockItr;
1017  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1018  std::string blockId = *blockItr;
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  LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
1034  nEntriesPerRow[lid] += col_gids.size();
1035  }
1036  }
1037  }
1038  Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
1039  RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,map_j, nEntriesPerRowView,
1040  Tpetra::StaticProfile));
1041 
1042 
1043 
1044  // graph information about the mesh
1045  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1046  std::string blockId = *blockItr;
1047 
1048  // grab elements for this block
1049  const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the
1050  // same elements in each element block
1051 
1052  // get information about number of indicies
1053  std::vector<GlobalOrdinalT> row_gids;
1054  std::vector<GlobalOrdinalT> col_gids;
1055 
1056  // loop over the elemnts
1057  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1058 
1059  rowProvider->getElementGIDs(elements[elmt],row_gids);
1060  colProvider->getElementGIDs(elements[elmt],col_gids);
1061  for(std::size_t row=0;row<row_gids.size();row++)
1062  graph->insertGlobalIndices(row_gids[row],col_gids);
1063  }
1064  }
1065 
1066  // finish filling the graph: Make sure the colmap and row maps coincide to
1067  // minimize calls to LID lookups
1068  graph->fillComplete(getMap(j),getMap(i));
1069 
1070  return graph;
1071 }
1072 
1073 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1076 getTpetraMatrix(int i,int j) const
1077 {
1078  Teuchos::RCP<const MapType> map_i = getMap(i);
1079  Teuchos::RCP<const MapType> map_j = getMap(j);
1080 
1081  Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
1083  mat->fillComplete(map_j,map_i);
1084 
1085  return mat;
1086 }
1087 
1088 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1091 getGhostedTpetraMatrix(int i,int j) const
1092 {
1093  Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1094  Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1095 
1096  Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1098  mat->fillComplete(getMap(j),getMap(i));
1099 
1100  return mat;
1101 }
1102 
1103 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1107 {
1108  Teuchos::RCP<const MapType> tMap = getMap(i);
1109  return Teuchos::rcp(new VectorType(tMap));
1110 }
1111 
1112 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1116 {
1117  Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1118  return Teuchos::rcp(new VectorType(tMap));
1119 }
1120 
1121 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1125 {
1126  Teuchos::RCP<const MapType> tMap = getMap(i);
1127  return Teuchos::rcp(new VectorType(tMap));
1128 }
1129 
1130 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1134 {
1135  Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1136  return Teuchos::rcp(new VectorType(tMap));
1137 }
1138 
1139 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1140 int
1143 {
1144  return gidProviders_.size();
1145 }
1146 
1147 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1148 int
1151 {
1152  return gidProviders_.size();
1153 }
1154 
1155 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1158 {
1159  BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1160  if(tloc.get_A()!=Teuchos::null)
1161  tloc.beginFill();
1162 }
1163 
1164 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1167 {
1168  BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1169  if(tloc.get_A()!=Teuchos::null)
1170  tloc.endFill();
1171 }
1172 
1173 }
1174 
1175 #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)