Panzer  Version of the Day
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_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
44 #define PANZER_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
45 
47 
48 #include "Tpetra_MultiVector.hpp"
49 #include "Tpetra_Vector.hpp"
50 #include "Tpetra_CrsMatrix.hpp"
51 
52 #include "Thyra_TpetraThyraWrappers.hpp"
53 #include "Thyra_DefaultProductVectorSpace.hpp"
54 #include "Thyra_DefaultBlockedLinearOp.hpp"
55 #include "Thyra_TpetraLinearOp.hpp"
56 #include "Thyra_SpmdVectorBase.hpp"
57 
60 
61 using Teuchos::RCP;
62 
63 namespace panzer {
64 
65 // ************************************************************
66 // class BlockedTpetraLinearObjFactory
67 // ************************************************************
68 
69 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
71 BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
72  const Teuchos::RCP<const BlockedDOFManager<LocalOrdinalT,GlobalOrdinalT> > & gidProvider)
73  : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
74 {
75  for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
76  gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
77 
79 
80  // build and register the gather/scatter evaluators with
81  // the base class.
82  this->buildGatherScatterEvaluators(*this);
83 }
84 
85 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
88 { }
89 
90 // LinearObjectFactory functions
92 
93 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
94 Teuchos::RCP<LinearObjContainer>
97 {
98  std::vector<Teuchos::RCP<const MapType> > blockMaps;
99  std::size_t blockDim = gidProviders_.size();
100  for(std::size_t i=0;i<blockDim;i++)
101  blockMaps.push_back(getMap(i));
102 
103  Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
104  container->setMapsForBlocks(blockMaps);
105 
106  return container;
107 }
108 
109 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
110 Teuchos::RCP<LinearObjContainer>
113 {
114  std::vector<Teuchos::RCP<const MapType> > blockMaps;
115  std::size_t blockDim = gidProviders_.size();
116  for(std::size_t i=0;i<blockDim;i++)
117  blockMaps.push_back(getGhostedMap(i));
118 
119  Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
120  container->setMapsForBlocks(blockMaps);
121 
122  return container;
123 }
124 
125 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
128 {
129  using Teuchos::is_null;
130 
131  typedef LinearObjContainer LOC;
132 
133  const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
134  BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
135 
136  // Operations occur if the GLOBAL container has the correct targets!
137  // Users set the GLOBAL continer arguments
138  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
139  globalToGhostThyraVector(b_in.get_x(),b_out.get_x());
140 
141  if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
142  globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt());
143 
144  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
145  globalToGhostThyraVector(b_in.get_f(),b_out.get_f());
146 }
147 
148 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
151 {
152  using Teuchos::is_null;
153 
154  typedef LinearObjContainer LOC;
155 
156  const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
157  BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
158 
159  // Operations occur if the GLOBAL container has the correct targets!
160  // Users set the GLOBAL continer arguments
161  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
162  ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x());
163 
164  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
165  ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f());
166 
167  if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
168  ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
169 }
170 
171 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
174  const LinearObjContainer & globalBCRows,
175  LinearObjContainer & ghostedObjs,
176  bool zeroVectorRows, bool adjustX) const
177 {
178  using Teuchos::RCP;
179  using Teuchos::rcp_dynamic_cast;
180  using Thyra::LinearOpBase;
181  using Thyra::PhysicallyBlockedLinearOpBase;
182  using Thyra::VectorBase;
184 
185  std::size_t blockDim = gidProviders_.size();
186 
187  // first cast to block LOCs
188  const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
189  const BTLOC & b_globalBCRows = Teuchos::dyn_cast<const BTLOC>(globalBCRows);
190  BTLOC & b_ghosted = Teuchos::dyn_cast<BTLOC>(ghostedObjs);
191 
192  TEUCHOS_ASSERT(b_localBCRows.get_f()!=Teuchos::null);
193  TEUCHOS_ASSERT(b_globalBCRows.get_f()!=Teuchos::null);
194 
195  // cast each component as needed to their product form
196  RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.get_A());
197  RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_f());
198  RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.get_f(),true);
199  RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.get_f(),true);
200 
201  if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
202 
203  // sanity check!
204  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
205  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
206  if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
207  TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
208  TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
209 
210  for(std::size_t i=0;i<blockDim;i++) {
211  // grab epetra vector
212  RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),true)->getConstTpetraVector();
213  RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),true)->getConstTpetraVector();
214 
215  // pull out epetra values
216  RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
217  RCP<VectorType> t_f;
218  if(th_f==Teuchos::null)
219  t_f = Teuchos::null;
220  else
221  t_f = rcp_dynamic_cast<ThyraVector>(th_f,true)->getTpetraVector();
222 
223  for(std::size_t j=0;j<blockDim;j++) {
224  RCP<const MapType> map_i = getGhostedMap(i);
225  RCP<const MapType> map_j = getGhostedMap(j);
226 
227  // pull out epetra values
228  RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
229 
230  // don't do anyting if opertor is null
231  RCP<CrsMatrixType> t_A;
232  if(th_A==Teuchos::null)
233  t_A = Teuchos::null;
234  else {
235  RCP<OperatorType> t_A_op = rcp_dynamic_cast<ThyraLinearOp>(th_A,true)->getTpetraOperator();
236  t_A = rcp_dynamic_cast<CrsMatrixType>(t_A_op,true);
237  }
238 
239  // adjust Block operator
240  if(t_A!=Teuchos::null) {
241  t_A->resumeFill();
242  }
243 
244  adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
245 
246  if(t_A!=Teuchos::null) {
247  //t_A->fillComplete(map_j,map_i);
248  }
249 
250  t_f = Teuchos::null; // this is so we only adjust it once on the first pass
251  }
252  }
253 }
254 
255 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
258  const VectorType & global_bcs,
259  const Teuchos::Ptr<VectorType> & f,
260  const Teuchos::Ptr<CrsMatrixType> & A,
261  bool zeroVectorRows) const
262 {
263  if(f==Teuchos::null && A==Teuchos::null)
264  return;
265 
266  Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
267 
268  Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
269  Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
270 
271  TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
272  for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
273  if(global_bcs_array[i]==0.0)
274  continue;
275 
276  if(local_bcs_array[i]==0.0 || zeroVectorRows) {
277  // this boundary condition was NOT set by this processor
278 
279  // if they exist put 0.0 in each entry
280  if(!Teuchos::is_null(f))
281  f_array[i] = 0.0;
282  if(!Teuchos::is_null(A)) {
283  std::size_t numEntries = 0;
284  std::size_t sz = A->getNumEntriesInLocalRow(i);
285  Teuchos::Array<LocalOrdinalT> indices(sz);
286  Teuchos::Array<ScalarT> values(sz);
287 
288  A->getLocalRowCopy(i,indices,values,numEntries);
289 
290  for(std::size_t c=0;c<numEntries;c++)
291  values[c] = 0.0;
292 
293  A->replaceLocalValues(i,indices,values);
294  }
295  }
296  else {
297  // this boundary condition was set by this processor
298 
299  ScalarT scaleFactor = global_bcs_array[i];
300 
301  // if they exist scale linear objects by scale factor
302  if(!Teuchos::is_null(f))
303  f_array[i] /= scaleFactor;
304  if(!Teuchos::is_null(A)) {
305  std::size_t numEntries = 0;
306  std::size_t sz = A->getNumEntriesInLocalRow(i);
307  Teuchos::Array<LocalOrdinalT> indices(sz);
308  Teuchos::Array<ScalarT> values(sz);
309 
310  A->getLocalRowCopy(i,indices,values,numEntries);
311 
312  for(std::size_t c=0;c<numEntries;c++)
313  values[c] /= scaleFactor;
314 
315  A->replaceLocalValues(i,indices,values);
316  }
317  }
318  }
319 }
320 
321 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
324  LinearObjContainer & result) const
325 {
326  TEUCHOS_ASSERT(false); // not yet implemented
327 }
328 
329 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
330 Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
333 {
334  using Teuchos::RCP;
335  using Teuchos::rcp;
336 
337  std::vector<RCP<ReadOnlyVector_GlobalEvaluationData> > gedBlocks;
338  for(int i=0;i<getBlockColCount();i++) {
339  RCP<TpetraVector_ReadOnly_GlobalEvaluationData<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > vec_ged
341  vec_ged->initialize(getGhostedImport(i),getGhostedMap(i),getMap(i));
342 
343  gedBlocks.push_back(vec_ged);
344  }
345 
346  RCP<BlockedVector_ReadOnly_GlobalEvaluationData> ged = rcp(new BlockedVector_ReadOnly_GlobalEvaluationData);
347  ged->initialize(getGhostedThyraDomainSpace(),getThyraDomainSpace(),gedBlocks);
348 
349  return ged;
350 }
351 
352 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
354 getComm() const
355 {
356  return *comm_;
357 }
358 
359 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
362 {
363  BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
364  initializeContainer(mem,bloc);
365 }
366 
367 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
370 {
371  BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
372  initializeGhostedContainer(mem,bloc);
373 }
374 
375 // Generic methods
377 
378 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
380 initializeContainer(int mem,BTLOC & loc) const
381 {
382  typedef LinearObjContainer LOC;
383 
384  loc.clear();
385 
386  if((mem & LOC::X) == LOC::X)
387  loc.set_x(getThyraDomainVector());
388 
389  if((mem & LOC::DxDt) == LOC::DxDt)
390  loc.set_dxdt(getThyraDomainVector());
391 
392  if((mem & LOC::F) == LOC::F)
393  loc.set_f(getThyraRangeVector());
394 
395  if((mem & LOC::Mat) == LOC::Mat)
396  loc.set_A(getThyraMatrix());
397 }
398 
399 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
401 initializeGhostedContainer(int mem,BTLOC & loc) const
402 {
403  typedef LinearObjContainer LOC;
404 
405  loc.clear();
406 
407  if((mem & LOC::X) == LOC::X)
408  loc.set_x(getGhostedThyraDomainVector());
409 
410  if((mem & LOC::DxDt) == LOC::DxDt)
411  loc.set_dxdt(getGhostedThyraDomainVector());
412 
413  if((mem & LOC::F) == LOC::F) {
414  loc.set_f(getGhostedThyraRangeVector());
416  }
417 
418  if((mem & LOC::Mat) == LOC::Mat) {
419  loc.set_A(getGhostedThyraMatrix());
421  }
422 }
423 
424 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
426 addExcludedPair(int rowBlock,int colBlock)
427 {
428  excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
429 }
430 
431 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
433 addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs)
434 {
435  for(std::size_t i=0;i<exPairs.size();i++)
436  excludedPairs_.insert(exPairs[i]);
437 }
438 
439 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
440 Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> >
442 getGlobalIndexer(int i) const
443 {
444  return gidProviders_[i];
445 }
446 
447 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
449 makeRoomForBlocks(std::size_t blockCnt)
450 {
451  maps_.resize(blockCnt);
452  ghostedMaps_.resize(blockCnt);
453  importers_.resize(blockCnt);
454  exporters_.resize(blockCnt);
455 }
456 
457 // Thyra methods
459 
460 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
461 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
464 {
465  if(domainSpace_==Teuchos::null) {
466  // loop over all vectors and build the vector space
467  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
468  for(std::size_t i=0;i<gidProviders_.size();i++)
469  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
470 
471  domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
472  }
473 
474  return domainSpace_;
475 }
476 
477 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
478 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
481 {
482  if(rangeSpace_==Teuchos::null) {
483  // loop over all vectors and build the vector space
484  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
485  for(std::size_t i=0;i<gidProviders_.size();i++)
486  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
487 
488  rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
489  }
490 
491  return rangeSpace_;
492 }
493 
494 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
495 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
498 {
499  Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
500  Thyra::createMember<ScalarT>(*getThyraDomainSpace());
501  Thyra::assign(vec.ptr(),0.0);
502 
503  Teuchos::RCP<Thyra::ProductVectorBase<ScalarT> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<ScalarT> >(vec);
504  for(std::size_t i=0;i<gidProviders_.size();i++) {
505  TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
506  Teuchos::as<int>(getMap(i)->getNodeNumElements()));
507  }
508 
509  return vec;
510 }
511 
512 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
513 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
516 {
517  Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
518  Thyra::createMember<ScalarT>(*getThyraRangeSpace());
519  Thyra::assign(vec.ptr(),0.0);
520 
521  return vec;
522 }
523 
524 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
525 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
528 {
529  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
530 
531  // get the block dimension
532  std::size_t blockDim = gidProviders_.size();
533 
534  // this operator will be square
535  blockedOp->beginBlockFill(blockDim,blockDim);
536 
537  // loop over each block
538  for(std::size_t i=0;i<blockDim;i++) {
539  for(std::size_t j=0;j<blockDim;j++) {
540  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
541  // build (i,j) block matrix and add it to blocked operator
542  Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
543  blockedOp->setNonconstBlock(i,j,block);
544  }
545  }
546  }
547 
548  // all done
549  blockedOp->endBlockFill();
550 
551  return blockedOp;
552 }
553 
554 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
555 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
558 {
559  if(ghostedDomainSpace_==Teuchos::null) {
560  // loop over all vectors and build the vector space
561  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
562  for(std::size_t i=0;i<gidProviders_.size();i++)
563  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
564 
565  ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
566  }
567 
568  return ghostedDomainSpace_;
569 }
570 
571 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
572 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
575 {
576  if(ghostedRangeSpace_==Teuchos::null) {
577  // loop over all vectors and build the vector space
578  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
579  for(std::size_t i=0;i<gidProviders_.size();i++)
580  vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
581 
582  ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
583  }
584 
585  return ghostedRangeSpace_;
586 }
587 
588 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
589 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
592 {
593  Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
594  Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
595  Thyra::assign(vec.ptr(),0.0);
596 
597  return vec;
598 }
599 
600 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
601 Teuchos::RCP<Thyra::VectorBase<ScalarT> >
604 {
605  Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
606  Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
607  Thyra::assign(vec.ptr(),0.0);
608 
609  return vec;
610 }
611 
612 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
613 Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
616 {
617  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
618 
619  // get the block dimension
620  std::size_t blockDim = gidProviders_.size();
621 
622  // this operator will be square
623  blockedOp->beginBlockFill(blockDim,blockDim);
624 
625  // loop over each block
626  for(std::size_t i=0;i<blockDim;i++) {
627  for(std::size_t j=0;j<blockDim;j++) {
628  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
629  // build (i,j) block matrix and add it to blocked operator
630  Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
631  blockedOp->setNonconstBlock(i,j,block);
632  }
633  }
634  }
635 
636  // all done
637  blockedOp->endBlockFill();
638 
639  return blockedOp;
640 }
641 
642 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
644 ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
645  const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
646 {
647  using Teuchos::RCP;
648  using Teuchos::rcp_dynamic_cast;
650 
651  std::size_t blockDim = gidProviders_.size();
652 
653  // get product vectors
654  RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
655  RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
656 
657  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
658  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
659 
660  for(std::size_t i=0;i<blockDim;i++) {
661  // first get each Tpetra vector
662  RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
663  RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
664 
665  // use Tpetra to do global communication
666  ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
667  }
668 }
669 
670 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
673 {
674  using Teuchos::RCP;
675  using Teuchos::rcp_dynamic_cast;
676  using Teuchos::dyn_cast;
677  using Thyra::LinearOpBase;
678  using Thyra::PhysicallyBlockedLinearOpBase;
679 
680  std::size_t blockDim = gidProviders_.size();
681 
682  // get product vectors
683  const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
684  PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
685 
686  TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
687  TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
688  TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
689  TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
690 
691  for(std::size_t i=0;i<blockDim;i++) {
692  for(std::size_t j=0;j<blockDim;j++) {
693  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
694  // extract the blocks
695  RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
696  RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
697 
698  // sanity check
699  TEUCHOS_ASSERT(th_in!=Teuchos::null);
700  TEUCHOS_ASSERT(th_out!=Teuchos::null);
701 
702  // get the epetra version of the blocks
703  RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<const ThyraLinearOp>(th_in,true)->getConstTpetraOperator();
704  RCP<OperatorType> tp_op_out = rcp_dynamic_cast<ThyraLinearOp>(th_out,true)->getTpetraOperator();
705 
706  RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<const CrsMatrixType>(tp_op_in,true);
707  RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<CrsMatrixType>(tp_op_out,true);
708 
709  // use Tpetra to do global communication
710  ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
711  }
712  }
713  }
714 }
715 
716 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
718 globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
719  const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
720 {
721  using Teuchos::RCP;
722  using Teuchos::rcp_dynamic_cast;
724 
725  std::size_t blockDim = gidProviders_.size();
726 
727  // get product vectors
728  RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
729  RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
730 
731  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
732  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
733 
734  for(std::size_t i=0;i<blockDim;i++) {
735  // first get each Tpetra vector
736  RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
737  RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
738 
739  // use Tpetra to do global communication
740  globalToGhostTpetraVector(i,*tp_in,*tp_out);
741  }
742 }
743 
744 // Tpetra methods
746 
747 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
749 ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out) const
750 {
751  using Teuchos::RCP;
752 
753  // do the global distribution
754  RCP<const ExportType> exporter = getGhostedExport(i);
755  out.putScalar(0.0);
756  out.doExport(in,*exporter,Tpetra::ADD);
757 }
758 
759 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
761 ghostToGlobalTpetraMatrix(int blockRow,const CrsMatrixType & in,CrsMatrixType & out) const
762 {
763  using Teuchos::RCP;
764 
765  RCP<const MapType> map_i = out.getRangeMap();
766  RCP<const MapType> map_j = out.getDomainMap();
767 
768  // do the global distribution
769  RCP<const ExportType> exporter = getGhostedExport(blockRow);
770 
771  out.resumeFill();
772  out.setAllToScalar(0.0);
773  out.doExport(in,*exporter,Tpetra::ADD);
774  out.fillComplete(map_j,map_i);
775 }
776 
777 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
779 globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out) const
780 {
781  using Teuchos::RCP;
782 
783  // do the global distribution
784  RCP<const ImportType> importer = getGhostedImport(i);
785  out.putScalar(0.0);
786  out.doImport(in,*importer,Tpetra::INSERT);
787 }
788 
789 // get the map from the matrix
790 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
791 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
793 getMap(int i) const
794 {
795  if(maps_[i]==Teuchos::null)
796  maps_[i] = buildTpetraMap(i);
797 
798  return maps_[i];
799 }
800 
801 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
802 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
804 getGhostedMap(int i) const
805 {
806  if(ghostedMaps_[i]==Teuchos::null)
807  ghostedMaps_[i] = buildTpetraGhostedMap(i);
808 
809  return ghostedMaps_[i];
810 }
811 
812 // get the graph of the crs matrix
813 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
814 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
816 getGraph(int i,int j) const
817 {
818  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
819 
820  typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
821  Teuchos::RCP<const CrsGraphType> graph;
822  if(itr==graphs_.end()) {
823  graph = buildTpetraGraph(i,j);
824  graphs_[std::make_pair(i,j)] = graph;
825  }
826  else
827  graph = itr->second;
828 
829  TEUCHOS_ASSERT(graph!=Teuchos::null);
830  return graph;
831 }
832 
833 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
834 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
836 getGhostedGraph(int i,int j) const
837 {
838  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
839 
840  typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
841  Teuchos::RCP<const CrsGraphType> ghostedGraph;
842  if(itr==ghostedGraphs_.end()) {
843  ghostedGraph = buildTpetraGhostedGraph(i,j);
844  ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
845  }
846  else
847  ghostedGraph = itr->second;
848 
849  TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
850  return ghostedGraph;
851 }
852 
853 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
854 Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
856 getGhostedImport(int i) const
857 {
858  if(importers_[i]==Teuchos::null)
859  importers_[i] = Teuchos::rcp(new ImportType(getMap(i),getGhostedMap(i)));
860 
861  return importers_[i];
862 }
863 
864 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
865 Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
867 getGhostedExport(int i) const
868 {
869  if(exporters_[i]==Teuchos::null)
870  exporters_[i] = Teuchos::rcp(new ExportType(getGhostedMap(i),getMap(i)));
871 
872  return exporters_[i];
873 }
874 
875 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
876 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
878 buildTpetraMap(int i) const
879 {
880  std::vector<GlobalOrdinalT> indices;
881 
882  // get the global indices
883  getGlobalIndexer(i)->getOwnedIndices(indices);
884 
885  return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
886 }
887 
888 // build the ghosted map
889 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
890 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
893 {
894  std::vector<GlobalOrdinalT> indices;
895 
896  // get the global indices
897  getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
898 
899  return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
900 }
901 
902 // get the graph of the crs matrix
903 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
904 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
906 buildTpetraGraph(int i,int j) const
907 {
908  using Teuchos::RCP;
909  using Teuchos::rcp;
910 
911  // build the map and allocate the space for the graph and
912  // grab the ghosted graph
913  RCP<const MapType> map_i = getMap(i);
914  RCP<const MapType> map_j = getMap(j);
915 
916  RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,0));
917  RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
918 
919  // perform the communication to finish building graph
920  RCP<const ExportType> exporter = getGhostedExport(i);
921  graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
922  graph->fillComplete(map_j,map_i);
923 
924  return graph;
925 }
926 
927 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
928 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
930 buildTpetraGhostedGraph(int i,int j) const
931 {
932  using Teuchos::RCP;
933  using Teuchos::rcp;
934 
935  // build the map and allocate the space for the graph and
936  // grab the ghosted graph
937  RCP<const MapType> map_i = getGhostedMap(i);
938  RCP<const MapType> map_j = getGhostedMap(j);
939 
940  RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,map_j,0));
941 
942  std::vector<std::string> elementBlockIds;
943 
944  Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> > rowProvider, colProvider;
945 
946  rowProvider = getGlobalIndexer(i);
947  colProvider = getGlobalIndexer(j);
948 
949  blockProvider_->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
950  // same element blocks
951 
952  // graph information about the mesh
953  std::vector<std::string>::const_iterator blockItr;
954  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
955  std::string blockId = *blockItr;
956 
957  // grab elements for this block
958  const std::vector<LocalOrdinalT> & elements = blockProvider_->getElementBlock(blockId); // each sub provider "should" have the
959  // same elements in each element block
960 
961  // get information about number of indicies
962  std::vector<GlobalOrdinalT> row_gids;
963  std::vector<GlobalOrdinalT> col_gids;
964 
965  // loop over the elemnts
966  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
967 
968  rowProvider->getElementGIDs(elements[elmt],row_gids);
969  colProvider->getElementGIDs(elements[elmt],col_gids);
970  for(std::size_t row=0;row<row_gids.size();row++)
971  graph->insertGlobalIndices(row_gids[row],col_gids);
972  }
973  }
974 
975  // finish filling the graph: Make sure the colmap and row maps coincide to
976  // minimize calls to LID lookups
977  graph->fillComplete(map_j,map_i);
978 
979  return graph;
980 }
981 
982 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
983 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
985 getTpetraMatrix(int i,int j) const
986 {
987  Teuchos::RCP<const MapType> map_i = getMap(i);
988  Teuchos::RCP<const MapType> map_j = getMap(j);
989 
990  Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
991  Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
992  mat->fillComplete(map_j,map_i);
993 
994  return mat;
995 }
996 
997 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
998 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1000 getGhostedTpetraMatrix(int i,int j) const
1001 {
1002  Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1003  Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1004 
1005  Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1006  Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
1007  mat->fillComplete(map_j,map_i);
1008 
1009  return mat;
1010 }
1011 
1012 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1013 int
1016 {
1017  return gidProviders_.size();
1018 }
1019 
1020 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1021 int
1024 {
1025  return gidProviders_.size();
1026 }
1027 
1028 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1031 {
1032  BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1033  if(tloc.get_A()!=Teuchos::null)
1034  tloc.beginFill();
1035 }
1036 
1037 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1040 {
1041  BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1042  if(tloc.get_A()!=Teuchos::null)
1043  tloc.endFill();
1044 }
1045 
1046 }
1047 
1048 #endif
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
Thyra::TpetraVector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > ThyraVector
void ghostToGlobalTpetraMatrix(int blockRow, const CrsMatrixType &in, CrsMatrixType &out) const
Teuchos::RCP< const panzer::BlockedDOFManager< int, GlobalOrdinalT > > getGlobalIndexer() const
virtual Teuchos::RCP< const MapType > buildTpetraMap(int i) const
Thyra::TpetraLinearOp< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > ThyraLinearOp
void set_dxdt(const Teuchos::RCP< VectorType > &in)
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
std::vector< Teuchos::RCP< const UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > gidProviders_
virtual Teuchos::RCP< const MapType > buildTpetraGhostedMap(int i) const
virtual Teuchos::RCP< const ImportType > getGhostedImport(int i) const
get importer for converting an overalapped object to a "normal" object
Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > ImportType
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGraph(int i, int j) const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a Thyra operator.
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
PHX::MDField< ScalarT > vector
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
Teuchos::RCP< Thyra::BlockedLinearOpBase< ScalarT > > getGhostedThyraMatrix() const
Get a Thyra operator.
void initializeGhostedContainer(int, LinearObjContainer &loc) const
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
virtual Teuchos::RCP< const CrsGraphType > getGraph(int i, int j) const
get the graph of the crs matrix
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraRangeVector() const
Get a range vector.
void ghostToGlobalTpetraVector(int i, const VectorType &in, VectorType &out) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
std::vector< ScalarT > values
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraRangeVector() const
Get a range vector.
virtual Teuchos::RCP< const CrsGraphType > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
void makeRoomForBlocks(std::size_t blockCnt)
Allocate the space in the std::vector objects so we can fill with appropriate Tpetra data...
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraDomainVector() const
Get a domain vector.
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildDomainContainer() const
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraDomainVector() const
Get a domain vector.
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range vector space (f)
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
virtual Teuchos::RCP< const MapType > getGhostedMap(int i) const
get the ghosted map from the matrix
void globalToGhostTpetraVector(int i, const VectorType &in, VectorType &out) const
Teuchos::RCP< const Teuchos::Comm< int > > comm
BlockedTpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const BlockedDOFManager< LocalOrdinalT, GlobalOrdinalT > > &gidProvider)
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void buildGatherScatterEvaluators(const BuilderT &builder)
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
void initializeContainer(int, LinearObjContainer &loc) const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< ScalarT > &in, Thyra::LinearOpBase< ScalarT > &out) const
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGhostedGraph(int i, int j) const
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const