Panzer Version of the Day
Loading...
Searching...
No Matches
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
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
34namespace panzer {
35
36using Teuchos::RCP;
37
38// ************************************************************
39// class BlockedTpetraLinearObjFactory
40// ************************************************************
41
42template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
44BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
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.
56}
57
58template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
60BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
61 const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
62 : gidProviders_(gidProviders), comm_(comm)
63{
65}
66
67template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
71
72// LinearObjectFactory functions
74
75template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
76Teuchos::RCP<LinearObjContainer>
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
91template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
92Teuchos::RCP<LinearObjContainer>
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
107template <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
130template <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
153template <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;
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
237template <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
303template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
305applyDirichletBCs(const LinearObjContainer & /* counter */,
306 LinearObjContainer & /* result */) const
307{
308 TEUCHOS_ASSERT(false); // not yet implemented
309}
310
312//
313// buildReadOnlyDomainContainer()
314//
316template<typename Traits, typename ScalarT, typename LocalOrdinalT,
317 typename GlobalOrdinalT, typename NodeT>
318Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
319BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
320 NodeT>::
321buildReadOnlyDomainContainer() 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
344//
345// buildWriteDomainContainer()
346//
348template<typename Traits, typename ScalarT, typename LocalOrdinalT,
349 typename GlobalOrdinalT, typename NodeT>
350Teuchos::RCP<WriteVector_GlobalEvaluationData>
351BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
352 NodeT>::
353buildWriteDomainContainer() 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
364template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
370
371template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
373initializeContainer(int mem,LinearObjContainer & loc) const
374{
375 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
376 initializeContainer(mem,bloc);
377}
378
379template <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
390template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
392initializeContainer(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
411template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
413initializeGhostedContainer(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
436template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
438addExcludedPair(int rowBlock,int colBlock)
439{
440 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
441}
442
443template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
445addExcludedPairs(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
451template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
452Teuchos::RCP<const GlobalIndexer>
458
459template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
461makeRoomForBlocks(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
472template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
473Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
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
489template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
490Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
492getThyraRangeSpace() const
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
506template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
507Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
509getThyraDomainSpace(int blk) const
510{
511 if(domainSpace_==Teuchos::null) {
512 getThyraDomainSpace();
513 }
514
515 return domainSpace_->getBlock(blk);
516}
517
518template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
519Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
521getThyraRangeSpace(int blk) const
522{
523 if(rangeSpace_==Teuchos::null) {
524 getThyraRangeSpace();
525 }
526
527 return rangeSpace_->getBlock(blk);
528}
529
530template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
531Teuchos::RCP<Thyra::VectorBase<ScalarT> >
534{
535 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
536 Thyra::createMember<ScalarT>(*getThyraDomainSpace());
537 Thyra::assign(vec.ptr(),0.0);
538
539 Teuchos::RCP<Thyra::ProductVectorBase<ScalarT> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<ScalarT> >(vec);
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
548template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
549Teuchos::RCP<Thyra::VectorBase<ScalarT> >
552{
553 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
554 Thyra::createMember<ScalarT>(*getThyraRangeSpace());
555 Thyra::assign(vec.ptr(),0.0);
556
557 return vec;
558}
559
560template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
561Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
563getThyraMatrix() const
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
590template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
591Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
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
607template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
608Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
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
624template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
625Teuchos::RCP<Thyra::VectorBase<ScalarT> >
628{
629 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
630 Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
631 Thyra::assign(vec.ptr(),0.0);
632
633 return vec;
634}
635
636template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
637Teuchos::RCP<Thyra::VectorBase<ScalarT> >
640{
641 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
642 Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
643 Thyra::assign(vec.ptr(),0.0);
644
645 return vec;
646}
647
648template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
649Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
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
678template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
680ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
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
706template <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;
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
752template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
754globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
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
783template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
785ghostToGlobalTpetraVector(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
795template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
797ghostToGlobalTpetraMatrix(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
813template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
815globalToGhostTpetraVector(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
826template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
827Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
829getMap(int i) const
830{
831 if(maps_[i]==Teuchos::null)
832 maps_[i] = buildTpetraMap(i);
833
834 return maps_[i];
835}
836
837template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
838Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
840getGhostedMap(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
849template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
850Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
852getGraph(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
869template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
870Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
872getGhostedGraph(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
889template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
890Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
892getGhostedImport(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
900template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
901Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
903getGhostedExport(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
911template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
912Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
914buildTpetraMap(int i) const
915{
916 std::vector<GlobalOrdinalT> indices;
917
918 // get the global indices
919 getGlobalIndexer(i)->getOwnedIndices(indices);
920
921 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
922}
923
924// build the ghosted map
925template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
926Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
928buildTpetraGhostedMap(int i) const
929{
930 std::vector<GlobalOrdinalT> indices;
931
932 // get the global indices
933 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
934
935 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
936}
937
938// get the graph of the crs matrix
939template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
940Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
942buildTpetraGraph(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
963template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
964Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
966buildTpetraGhostedGraph(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
1044template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1045Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1047getTpetraMatrix(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);
1053 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
1054 mat->fillComplete(map_j,map_i);
1055
1056 return mat;
1057}
1058
1059template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1060Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1062getGhostedTpetraMatrix(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);
1068 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
1069 mat->fillComplete(getMap(j),getMap(i));
1070
1071 return mat;
1072}
1073
1074template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1075Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1077getTpetraDomainVector(int i) const
1078{
1079 Teuchos::RCP<const MapType> tMap = getMap(i);
1080 return Teuchos::rcp(new VectorType(tMap));
1081}
1082
1083template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1084Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1087{
1088 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1089 return Teuchos::rcp(new VectorType(tMap));
1090}
1091
1092template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1093Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1095getTpetraRangeVector(int i) const
1096{
1097 Teuchos::RCP<const MapType> tMap = getMap(i);
1098 return Teuchos::rcp(new VectorType(tMap));
1099}
1100
1101template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1102Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1104getGhostedTpetraRangeVector(int i) const
1105{
1106 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1107 return Teuchos::rcp(new VectorType(tMap));
1108}
1109
1110template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1111int
1117
1118template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1119int
1125
1126template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1128beginFill(LinearObjContainer & loc) const
1129{
1130 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1131 if(tloc.get_A()!=Teuchos::null)
1132 tloc.beginFill();
1133}
1134
1135template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1137endFill(LinearObjContainer & loc) const
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 set_A(const Teuchos::RCP< CrsMatrixType > &in)
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraRangeVector() const
Get a range vector.
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
std::vector< Teuchos::RCP< const GlobalIndexer > > gidProviders_
void initializeGhostedContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getThyraDomainVector() const
Get a domain vector.
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
virtual Teuchos::RCP< const CrsGraphType > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
virtual Teuchos::RCP< const MapType > buildTpetraGhostedMap(int i) const
virtual Teuchos::RCP< const MapType > buildTpetraMap(int i) const
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGraph(int i, int j) const
Teuchos::RCP< VectorType > getGhostedTpetraRangeVector(int i) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void makeRoomForBlocks(std::size_t blockCnt)
Allocate the space in the std::vector objects so we can fill with appropriate Tpetra data.
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraRangeVector() const
Get a range vector.
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
virtual Teuchos::RCP< const CrsGraphType > buildTpetraGhostedGraph(int i, int j) const
virtual Teuchos::RCP< const MapType > getGhostedMap(int i) const
get the ghosted map from the matrix
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
Teuchos::RCP< VectorType > getGhostedTpetraDomainVector(int i) const
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< ScalarT > > &in, const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &out) const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
void globalToGhostTpetraVector(int i, const VectorType &in, VectorType &out) const
BlockedTpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const BlockedDOFManager > &gidProvider)
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
void ghostToGlobalTpetraVector(int i, const VectorType &in, VectorType &out) const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< ScalarT > &in, Thyra::LinearOpBase< ScalarT > &out) const
void ghostToGlobalTpetraMatrix(int blockRow, const CrsMatrixType &in, CrsMatrixType &out) const
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) 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
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a Thyra operator.
Teuchos::RCP< Thyra::VectorBase< ScalarT > > getGhostedThyraDomainVector() const
Get a domain vector.
virtual Teuchos::RCP< const CrsGraphType > getGraph(int i, int j) const
get the graph of the crs matrix
Teuchos::RCP< const panzer::BlockedDOFManager > getGlobalIndexer() const
Teuchos::RCP< Thyra::BlockedLinearOpBase< ScalarT > > getGhostedThyraMatrix() const
Get a Thyra operator.
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors.
This class provides a boundary exchange communication mechanism for vectors.
void buildGatherScatterEvaluators(const BuilderT &builder)