Panzer Version of the Day
Loading...
Searching...
No Matches
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#ifdef PANZER_HAVE_EPETRA_STACK
49#include "Panzer_EpetraVector_Write_GlobalEvaluationData.hpp" // JMG: Remove this eventually.
50#endif
53
54// Thyra
55#include "Thyra_DefaultBlockedLinearOp.hpp"
56#include "Thyra_DefaultProductVectorSpace.hpp"
57#include "Thyra_SpmdVectorBase.hpp"
58#include "Thyra_TpetraLinearOp.hpp"
59#include "Thyra_TpetraThyraWrappers.hpp"
60
61// Tpetra
62#include "Tpetra_CrsMatrix.hpp"
63#include "Tpetra_MultiVector.hpp"
64#include "Tpetra_Vector.hpp"
65
66namespace panzer {
67
68using Teuchos::RCP;
69
70// ************************************************************
71// class BlockedTpetraLinearObjFactory
72// ************************************************************
73
74template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
76BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
77 const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
78 : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
79{
80 for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
81 gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
82
84
85 // build and register the gather/scatter evaluators with
86 // the base class.
88}
89
90template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
92BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
93 const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
94 : gidProviders_(gidProviders), comm_(comm)
95{
97}
98
99template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
102{ }
103
104// LinearObjectFactory functions
106
107template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
108Teuchos::RCP<LinearObjContainer>
111{
112 std::vector<Teuchos::RCP<const MapType> > blockMaps;
113 std::size_t blockDim = gidProviders_.size();
114 for(std::size_t i=0;i<blockDim;i++)
115 blockMaps.push_back(getMap(i));
116
117 Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
118 container->setMapsForBlocks(blockMaps);
119
120 return container;
121}
122
123template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
124Teuchos::RCP<LinearObjContainer>
127{
128 std::vector<Teuchos::RCP<const MapType> > blockMaps;
129 std::size_t blockDim = gidProviders_.size();
130 for(std::size_t i=0;i<blockDim;i++)
131 blockMaps.push_back(getGhostedMap(i));
132
133 Teuchos::RCP<BTLOC> container = Teuchos::rcp(new BTLOC);
134 container->setMapsForBlocks(blockMaps);
135
136 return container;
137}
138
139template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
142{
143 using Teuchos::is_null;
144
145 typedef LinearObjContainer LOC;
146
147 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
148 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
149
150 // Operations occur if the GLOBAL container has the correct targets!
151 // Users set the GLOBAL continer arguments
152 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
153 globalToGhostThyraVector(b_in.get_x(),b_out.get_x());
154
155 if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
156 globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt());
157
158 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
159 globalToGhostThyraVector(b_in.get_f(),b_out.get_f());
160}
161
162template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
165{
166 using Teuchos::is_null;
167
168 typedef LinearObjContainer LOC;
169
170 const BTLOC & b_in = Teuchos::dyn_cast<const BTLOC>(in);
171 BTLOC & b_out = Teuchos::dyn_cast<BTLOC>(out);
172
173 // Operations occur if the GLOBAL container has the correct targets!
174 // Users set the GLOBAL continer arguments
175 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
176 ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x());
177
178 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
179 ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f());
180
181 if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
182 ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
183}
184
185template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
188 const LinearObjContainer & globalBCRows,
189 LinearObjContainer & ghostedObjs,
190 bool zeroVectorRows, bool adjustX) const
191{
192 using Teuchos::RCP;
193 using Teuchos::rcp_dynamic_cast;
195 using Thyra::PhysicallyBlockedLinearOpBase;
196 using Thyra::VectorBase;
198
199 std::size_t blockDim = gidProviders_.size();
200
201 // first cast to block LOCs
202 const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
203 const BTLOC & b_globalBCRows = Teuchos::dyn_cast<const BTLOC>(globalBCRows);
204 BTLOC & b_ghosted = Teuchos::dyn_cast<BTLOC>(ghostedObjs);
205
206 TEUCHOS_ASSERT(b_localBCRows.get_f()!=Teuchos::null);
207 TEUCHOS_ASSERT(b_globalBCRows.get_f()!=Teuchos::null);
208
209 // cast each component as needed to their product form
210 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(b_ghosted.get_A());
211 RCP<ProductVectorBase<ScalarT> > f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_f());
212 RCP<ProductVectorBase<ScalarT> > local_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_localBCRows.get_f(),true);
213 RCP<ProductVectorBase<ScalarT> > global_bcs = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_globalBCRows.get_f(),true);
214
215 if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
216
217 // sanity check!
218 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
219 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
220 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
221 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
222 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
223
224 for(std::size_t i=0;i<blockDim;i++) {
225 // grab epetra vector
226 RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),true)->getConstTpetraVector();
227 RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),true)->getConstTpetraVector();
228
229 // pull out epetra values
230 RCP<VectorBase<ScalarT> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
231 RCP<VectorType> t_f;
232 if(th_f==Teuchos::null)
233 t_f = Teuchos::null;
234 else
235 t_f = rcp_dynamic_cast<ThyraVector>(th_f,true)->getTpetraVector();
236
237 for(std::size_t j=0;j<blockDim;j++) {
238 RCP<const MapType> map_i = getGhostedMap(i);
239 RCP<const MapType> map_j = getGhostedMap(j);
240
241 // pull out epetra values
242 RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
243
244 // don't do anyting if opertor is null
245 RCP<CrsMatrixType> t_A;
246 if(th_A==Teuchos::null)
247 t_A = Teuchos::null;
248 else {
249 RCP<OperatorType> t_A_op = rcp_dynamic_cast<ThyraLinearOp>(th_A,true)->getTpetraOperator();
250 t_A = rcp_dynamic_cast<CrsMatrixType>(t_A_op,true);
251 }
252
253 // adjust Block operator
254 if(t_A!=Teuchos::null) {
255 t_A->resumeFill();
256 }
257
258 adjustForDirichletConditions(*t_local_bcs,*t_global_bcs,t_f.ptr(),t_A.ptr(),zeroVectorRows);
259
260 if(t_A!=Teuchos::null) {
261 //t_A->fillComplete(map_j,map_i);
262 }
263
264 t_f = Teuchos::null; // this is so we only adjust it once on the first pass
265 }
266 }
267}
268
269template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
272 const VectorType & global_bcs,
273 const Teuchos::Ptr<VectorType> & f,
274 const Teuchos::Ptr<CrsMatrixType> & A,
275 bool zeroVectorRows) const
276{
277 if(f==Teuchos::null && A==Teuchos::null)
278 return;
279
280 Teuchos::ArrayRCP<ScalarT> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
281
282 Teuchos::ArrayRCP<const ScalarT> local_bcs_array = local_bcs.get1dView();
283 Teuchos::ArrayRCP<const ScalarT> global_bcs_array = global_bcs.get1dView();
284
285 TEUCHOS_ASSERT(local_bcs.getLocalLength()==global_bcs.getLocalLength());
286 for(std::size_t i=0;i<local_bcs.getLocalLength();i++) {
287 if(global_bcs_array[i]==0.0)
288 continue;
289
290 if(local_bcs_array[i]==0.0 || zeroVectorRows) {
291 // this boundary condition was NOT set by this processor
292
293 // if they exist put 0.0 in each entry
294 if(!Teuchos::is_null(f))
295 f_array[i] = 0.0;
296 if(!Teuchos::is_null(A)) {
297 std::size_t numEntries = 0;
298 std::size_t sz = A->getNumEntriesInLocalRow(i);
299 typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
300 typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
301
302 A->getLocalRowCopy(i,indices,values,numEntries);
303
304 for(std::size_t c=0;c<numEntries;c++)
305 values(c) = 0.0;
306
307 A->replaceLocalValues(i,indices,values);
308 }
309 }
310 else {
311 // this boundary condition was set by this processor
312
313 ScalarT scaleFactor = global_bcs_array[i];
314
315 // if they exist scale linear objects by scale factor
316 if(!Teuchos::is_null(f))
317 f_array[i] /= scaleFactor;
318 if(!Teuchos::is_null(A)) {
319 std::size_t numEntries = 0;
320 std::size_t sz = A->getNumEntriesInLocalRow(i);
321 typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
322 typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
323
324 A->getLocalRowCopy(i,indices,values,numEntries);
325
326 for(std::size_t c=0;c<numEntries;c++)
327 values(c) /= scaleFactor;
328
329 A->replaceLocalValues(i,indices,values);
330 }
331 }
332 }
333}
334
335template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
337applyDirichletBCs(const LinearObjContainer & /* counter */,
338 LinearObjContainer & /* result */) const
339{
340 TEUCHOS_ASSERT(false); // not yet implemented
341}
342
344//
345// buildReadOnlyDomainContainer()
346//
348template<typename Traits, typename ScalarT, typename LocalOrdinalT,
349 typename GlobalOrdinalT, typename NodeT>
350Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
351BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
352 NodeT>::
353buildReadOnlyDomainContainer() const
354{
355 using std::vector;
356 using Teuchos::RCP;
357 using Teuchos::rcp;
360 LocalOrdinalT, GlobalOrdinalT, NodeT>;
361 vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
362 for (int i(0); i < getBlockColCount(); ++i)
363 {
364 auto tvroged = rcp(new TVROGED);
365 tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
366 gedBlocks.push_back(tvroged);
367 }
368 auto ged = rcp(new BVROGED);
369 ged->initialize(getGhostedThyraDomainSpace(), getThyraDomainSpace(),
370 gedBlocks);
371 return ged;
372} // end of buildReadOnlyDomainContainer()
373
374#ifdef PANZER_HAVE_EPETRA_STACK
376//
377// buildWriteDomainContainer()
378//
380template<typename Traits, typename ScalarT, typename LocalOrdinalT,
381 typename GlobalOrdinalT, typename NodeT>
382Teuchos::RCP<WriteVector_GlobalEvaluationData>
383BlockedTpetraLinearObjFactory<Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT,
384 NodeT>::
385buildWriteDomainContainer() const
386{
387 using std::logic_error;
388 using Teuchos::rcp;
390 auto ged = rcp(new EVWGED);
391 TEUCHOS_TEST_FOR_EXCEPTION(true, logic_error, "NOT YET IMPLEMENTED")
392 return ged;
393} // end of buildWriteDomainContainer()
394#endif // PANZER_HAVE_EPETRA_STACK
395
396template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
398getComm() const
399{
400 return *comm_;
401}
402
403template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
405initializeContainer(int mem,LinearObjContainer & loc) const
406{
407 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
408 initializeContainer(mem,bloc);
409}
410
411template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
414{
415 BTLOC & bloc = Teuchos::dyn_cast<BTLOC>(loc);
416 initializeGhostedContainer(mem,bloc);
417}
418
419// Generic methods
421
422template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
424initializeContainer(int mem,BTLOC & loc) const
425{
426 typedef LinearObjContainer LOC;
427
428 loc.clear();
429
430 if((mem & LOC::X) == LOC::X)
431 loc.set_x(getThyraDomainVector());
432
433 if((mem & LOC::DxDt) == LOC::DxDt)
434 loc.set_dxdt(getThyraDomainVector());
435
436 if((mem & LOC::F) == LOC::F)
437 loc.set_f(getThyraRangeVector());
438
439 if((mem & LOC::Mat) == LOC::Mat)
440 loc.set_A(getThyraMatrix());
441}
442
443template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
445initializeGhostedContainer(int mem,BTLOC & loc) const
446{
447 typedef LinearObjContainer LOC;
448
449 loc.clear();
450
451 if((mem & LOC::X) == LOC::X)
452 loc.set_x(getGhostedThyraDomainVector());
453
454 if((mem & LOC::DxDt) == LOC::DxDt)
455 loc.set_dxdt(getGhostedThyraDomainVector());
456
457 if((mem & LOC::F) == LOC::F) {
458 loc.set_f(getGhostedThyraRangeVector());
460 }
461
462 if((mem & LOC::Mat) == LOC::Mat) {
463 loc.set_A(getGhostedThyraMatrix());
465 }
466}
467
468template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
470addExcludedPair(int rowBlock,int colBlock)
471{
472 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
473}
474
475template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
477addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs)
478{
479 for(std::size_t i=0;i<exPairs.size();i++)
480 excludedPairs_.insert(exPairs[i]);
481}
482
483template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
484Teuchos::RCP<const GlobalIndexer>
486getGlobalIndexer(int i) const
487{
488 return gidProviders_[i];
489}
490
491template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
493makeRoomForBlocks(std::size_t blockCnt)
494{
495 maps_.resize(blockCnt);
496 ghostedMaps_.resize(blockCnt);
497 importers_.resize(blockCnt);
498 exporters_.resize(blockCnt);
499}
500
501// Thyra methods
503
504template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
505Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
508{
509 if(domainSpace_==Teuchos::null) {
510 // loop over all vectors and build the vector space
511 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
512 for(std::size_t i=0;i<gidProviders_.size();i++)
513 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
514
515 domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
516 }
517
518 return domainSpace_;
519}
520
521template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
522Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
524getThyraRangeSpace() const
525{
526 if(rangeSpace_==Teuchos::null) {
527 // loop over all vectors and build the vector space
528 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
529 for(std::size_t i=0;i<gidProviders_.size();i++)
530 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
531
532 rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
533 }
534
535 return rangeSpace_;
536}
537
538template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
539Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
541getThyraDomainSpace(int blk) const
542{
543 if(domainSpace_==Teuchos::null) {
544 getThyraDomainSpace();
545 }
546
547 return domainSpace_->getBlock(blk);
548}
549
550template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
551Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
553getThyraRangeSpace(int blk) const
554{
555 if(rangeSpace_==Teuchos::null) {
556 getThyraRangeSpace();
557 }
558
559 return rangeSpace_->getBlock(blk);
560}
561
562template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
563Teuchos::RCP<Thyra::VectorBase<ScalarT> >
566{
567 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
568 Thyra::createMember<ScalarT>(*getThyraDomainSpace());
569 Thyra::assign(vec.ptr(),0.0);
570
571 Teuchos::RCP<Thyra::ProductVectorBase<ScalarT> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<ScalarT> >(vec);
572 for(std::size_t i=0;i<gidProviders_.size();i++) {
573 TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
574 Teuchos::as<int>(getMap(i)->getLocalNumElements()));
575 }
576
577 return vec;
578}
579
580template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
581Teuchos::RCP<Thyra::VectorBase<ScalarT> >
584{
585 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
586 Thyra::createMember<ScalarT>(*getThyraRangeSpace());
587 Thyra::assign(vec.ptr(),0.0);
588
589 return vec;
590}
591
592template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
593Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
595getThyraMatrix() const
596{
597 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
598
599 // get the block dimension
600 std::size_t blockDim = gidProviders_.size();
601
602 // this operator will be square
603 blockedOp->beginBlockFill(blockDim,blockDim);
604
605 // loop over each block
606 for(std::size_t i=0;i<blockDim;i++) {
607 for(std::size_t j=0;j<blockDim;j++) {
608 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
609 // build (i,j) block matrix and add it to blocked operator
610 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
611 blockedOp->setNonconstBlock(i,j,block);
612 }
613 }
614 }
615
616 // all done
617 blockedOp->endBlockFill();
618
619 return blockedOp;
620}
621
622template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
623Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
626{
627 if(ghostedDomainSpace_==Teuchos::null) {
628 // loop over all vectors and build the vector space
629 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
630 for(std::size_t i=0;i<gidProviders_.size();i++)
631 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
632
633 ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
634 }
635
636 return ghostedDomainSpace_;
637}
638
639template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
640Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
643{
644 if(ghostedRangeSpace_==Teuchos::null) {
645 // loop over all vectors and build the vector space
646 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
647 for(std::size_t i=0;i<gidProviders_.size();i++)
648 vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
649
650 ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
651 }
652
653 return ghostedRangeSpace_;
654}
655
656template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
657Teuchos::RCP<Thyra::VectorBase<ScalarT> >
660{
661 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
662 Thyra::createMember<ScalarT>(*getGhostedThyraDomainSpace());
663 Thyra::assign(vec.ptr(),0.0);
664
665 return vec;
666}
667
668template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
669Teuchos::RCP<Thyra::VectorBase<ScalarT> >
672{
673 Teuchos::RCP<Thyra::VectorBase<ScalarT> > vec =
674 Thyra::createMember<ScalarT>(*getGhostedThyraRangeSpace());
675 Thyra::assign(vec.ptr(),0.0);
676
677 return vec;
678}
679
680template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
681Teuchos::RCP<Thyra::BlockedLinearOpBase<ScalarT> >
684{
685 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
686
687 // get the block dimension
688 std::size_t blockDim = gidProviders_.size();
689
690 // this operator will be square
691 blockedOp->beginBlockFill(blockDim,blockDim);
692
693 // loop over each block
694 for(std::size_t i=0;i<blockDim;i++) {
695 for(std::size_t j=0;j<blockDim;j++) {
696 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
697 // build (i,j) block matrix and add it to blocked operator
698 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
699 blockedOp->setNonconstBlock(i,j,block);
700 }
701 }
702 }
703
704 // all done
705 blockedOp->endBlockFill();
706
707 return blockedOp;
708}
709
710template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
712ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
713 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
714{
715 using Teuchos::RCP;
716 using Teuchos::rcp_dynamic_cast;
718
719 std::size_t blockDim = gidProviders_.size();
720
721 // get product vectors
722 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
723 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
724
725 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
726 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
727
728 for(std::size_t i=0;i<blockDim;i++) {
729 // first get each Tpetra vector
730 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
731 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
732
733 // use Tpetra to do global communication
734 ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
735 }
736}
737
738template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
741{
742 using Teuchos::RCP;
743 using Teuchos::rcp_dynamic_cast;
744 using Teuchos::dyn_cast;
746 using Thyra::PhysicallyBlockedLinearOpBase;
747
748 std::size_t blockDim = gidProviders_.size();
749
750 // get product vectors
751 const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
752 PhysicallyBlockedLinearOpBase<ScalarT> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
753
754 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
755 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
756 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
757 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
758
759 for(std::size_t i=0;i<blockDim;i++) {
760 for(std::size_t j=0;j<blockDim;j++) {
761 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
762 // extract the blocks
763 RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
764 RCP<LinearOpBase<ScalarT> > th_out = prod_out.getNonconstBlock(i,j);
765
766 // sanity check
767 TEUCHOS_ASSERT(th_in!=Teuchos::null);
768 TEUCHOS_ASSERT(th_out!=Teuchos::null);
769
770 // get the epetra version of the blocks
771 RCP<const OperatorType> tp_op_in = rcp_dynamic_cast<const ThyraLinearOp>(th_in,true)->getConstTpetraOperator();
772 RCP<OperatorType> tp_op_out = rcp_dynamic_cast<ThyraLinearOp>(th_out,true)->getTpetraOperator();
773
774 RCP<const CrsMatrixType> tp_in = rcp_dynamic_cast<const CrsMatrixType>(tp_op_in,true);
775 RCP<CrsMatrixType> tp_out = rcp_dynamic_cast<CrsMatrixType>(tp_op_out,true);
776
777 // use Tpetra to do global communication
778 ghostToGlobalTpetraMatrix(i,*tp_in,*tp_out);
779 }
780 }
781 }
782}
783
784template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
786globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
787 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
788{
789 using Teuchos::RCP;
790 using Teuchos::rcp_dynamic_cast;
792
793 std::size_t blockDim = gidProviders_.size();
794
795 // get product vectors
796 RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
797 RCP<ProductVectorBase<ScalarT> > prod_out = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(out,true);
798
799 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
800 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
801
802 for(std::size_t i=0;i<blockDim;i++) {
803 // first get each Tpetra vector
804 RCP<const VectorType> tp_in = rcp_dynamic_cast<const ThyraVector>(prod_in->getVectorBlock(i),true)->getConstTpetraVector();
805 RCP<VectorType> tp_out = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
806
807 // use Tpetra to do global communication
808 globalToGhostTpetraVector(i,*tp_in,*tp_out);
809 }
810}
811
812// Tpetra methods
814
815template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
817ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out) const
818{
819 using Teuchos::RCP;
820
821 // do the global distribution
822 RCP<const ExportType> exporter = getGhostedExport(i);
823 out.putScalar(0.0);
824 out.doExport(in,*exporter,Tpetra::ADD);
825}
826
827template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
829ghostToGlobalTpetraMatrix(int blockRow,const CrsMatrixType & in,CrsMatrixType & out) const
830{
831 using Teuchos::RCP;
832
833 RCP<const MapType> map_i = out.getRangeMap();
834 RCP<const MapType> map_j = out.getDomainMap();
835
836 // do the global distribution
837 RCP<const ExportType> exporter = getGhostedExport(blockRow);
838
839 out.resumeFill();
840 out.setAllToScalar(0.0);
841 out.doExport(in,*exporter,Tpetra::ADD);
842 out.fillComplete(map_j,map_i);
843}
844
845template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
847globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out) const
848{
849 using Teuchos::RCP;
850
851 // do the global distribution
852 RCP<const ImportType> importer = getGhostedImport(i);
853 out.putScalar(0.0);
854 out.doImport(in,*importer,Tpetra::INSERT);
855}
856
857// get the map from the matrix
858template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
859Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
861getMap(int i) const
862{
863 if(maps_[i]==Teuchos::null)
864 maps_[i] = buildTpetraMap(i);
865
866 return maps_[i];
867}
868
869template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
870Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
872getGhostedMap(int i) const
873{
874 if(ghostedMaps_[i]==Teuchos::null)
875 ghostedMaps_[i] = buildTpetraGhostedMap(i);
876
877 return ghostedMaps_[i];
878}
879
880// get the graph of the crs matrix
881template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
882Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
884getGraph(int i,int j) const
885{
886 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
887
888 typename GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
889 Teuchos::RCP<const CrsGraphType> graph;
890 if(itr==graphs_.end()) {
891 graph = buildTpetraGraph(i,j);
892 graphs_[std::make_pair(i,j)] = graph;
893 }
894 else
895 graph = itr->second;
896
897 TEUCHOS_ASSERT(graph!=Teuchos::null);
898 return graph;
899}
900
901template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
902Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
904getGhostedGraph(int i,int j) const
905{
906 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> GraphMap;
907
908 typename GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
909 Teuchos::RCP<const CrsGraphType> ghostedGraph;
910 if(itr==ghostedGraphs_.end()) {
911 ghostedGraph = buildTpetraGhostedGraph(i,j);
912 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
913 }
914 else
915 ghostedGraph = itr->second;
916
917 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
918 return ghostedGraph;
919}
920
921template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
922Teuchos::RCP<const Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
924getGhostedImport(int i) const
925{
926 if(importers_[i]==Teuchos::null)
927 importers_[i] = Teuchos::rcp(new ImportType(getMap(i),getGhostedMap(i)));
928
929 return importers_[i];
930}
931
932template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
933Teuchos::RCP<const Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
935getGhostedExport(int i) const
936{
937 if(exporters_[i]==Teuchos::null)
938 exporters_[i] = Teuchos::rcp(new ExportType(getGhostedMap(i),getMap(i)));
939
940 return exporters_[i];
941}
942
943template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
944Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
946buildTpetraMap(int i) const
947{
948 std::vector<GlobalOrdinalT> indices;
949
950 // get the global indices
951 getGlobalIndexer(i)->getOwnedIndices(indices);
952
953 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
954}
955
956// build the ghosted map
957template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
958Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
960buildTpetraGhostedMap(int i) const
961{
962 std::vector<GlobalOrdinalT> indices;
963
964 // get the global indices
965 getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
966
967 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
968}
969
970// get the graph of the crs matrix
971template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
972Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
974buildTpetraGraph(int i,int j) const
975{
976 using Teuchos::RCP;
977 using Teuchos::rcp;
978
979 // build the map and allocate the space for the graph and
980 // grab the ghosted graph
981 RCP<const MapType> map_i = getMap(i);
982 RCP<const MapType> map_j = getMap(j);
983
984 RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,0));
985 RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
986
987 // perform the communication to finish building graph
988 RCP<const ExportType> exporter = getGhostedExport(i);
989 graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
990 graph->fillComplete(map_j,map_i);
991
992 return graph;
993}
994
995template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
996Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
998buildTpetraGhostedGraph(int i,int j) const
999{
1000 using Teuchos::RCP;
1001 using Teuchos::rcp;
1002
1003 // build the map and allocate the space for the graph and
1004 // grab the ghosted graph
1005 RCP<const MapType> map_i = getGhostedMap(i);
1006 RCP<const MapType> map_j = getGhostedMap(j);
1007
1008 std::vector<std::string> elementBlockIds;
1009
1010 Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
1011
1012 rowProvider = getGlobalIndexer(i);
1013 colProvider = getGlobalIndexer(j);
1014
1015 gidProviders_[0]->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
1016 // same element blocks
1017
1018 // Count number of entries in each row of graph; needed for graph constructor
1019 std::vector<size_t> nEntriesPerRow(map_i->getLocalNumElements(), 0);
1020 std::vector<std::string>::const_iterator blockItr;
1021 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1022 std::string blockId = *blockItr;
1023 // grab elements for this block
1024 const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the
1025 // same elements in each element block
1026
1027 // get information about number of indicies
1028 std::vector<GlobalOrdinalT> row_gids;
1029 std::vector<GlobalOrdinalT> col_gids;
1030
1031 // loop over the elemnts
1032 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1033
1034 rowProvider->getElementGIDs(elements[elmt],row_gids);
1035 colProvider->getElementGIDs(elements[elmt],col_gids);
1036 for(std::size_t row=0;row<row_gids.size();row++) {
1037 LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
1038 nEntriesPerRow[lid] += col_gids.size();
1039 }
1040 }
1041 }
1042 Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
1043 RCP<CrsGraphType> graph = rcp(new CrsGraphType(map_i,map_j, nEntriesPerRowView));
1044
1045
1046
1047 // graph information about the mesh
1048 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1049 std::string blockId = *blockItr;
1050
1051 // grab elements for this block
1052 const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the
1053 // same elements in each element block
1054
1055 // get information about number of indicies
1056 std::vector<GlobalOrdinalT> row_gids;
1057 std::vector<GlobalOrdinalT> col_gids;
1058
1059 // loop over the elemnts
1060 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1061
1062 rowProvider->getElementGIDs(elements[elmt],row_gids);
1063 colProvider->getElementGIDs(elements[elmt],col_gids);
1064 for(std::size_t row=0;row<row_gids.size();row++)
1065 graph->insertGlobalIndices(row_gids[row],col_gids);
1066 }
1067 }
1068
1069 // finish filling the graph: Make sure the colmap and row maps coincide to
1070 // minimize calls to LID lookups
1071 graph->fillComplete(getMap(j),getMap(i));
1072
1073 return graph;
1074}
1075
1076template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1077Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1079getTpetraMatrix(int i,int j) const
1080{
1081 Teuchos::RCP<const MapType> map_i = getMap(i);
1082 Teuchos::RCP<const MapType> map_j = getMap(j);
1083
1084 Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
1085 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
1086 mat->fillComplete(map_j,map_i);
1087
1088 return mat;
1089}
1090
1091template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1092Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1094getGhostedTpetraMatrix(int i,int j) const
1095{
1096 Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
1097 Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
1098
1099 Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
1100 Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
1101 mat->fillComplete(getMap(j),getMap(i));
1102
1103 return mat;
1104}
1105
1106template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1107Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1109getTpetraDomainVector(int i) const
1110{
1111 Teuchos::RCP<const MapType> tMap = getMap(i);
1112 return Teuchos::rcp(new VectorType(tMap));
1113}
1114
1115template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1116Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1119{
1120 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1121 return Teuchos::rcp(new VectorType(tMap));
1122}
1123
1124template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1125Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1127getTpetraRangeVector(int i) const
1128{
1129 Teuchos::RCP<const MapType> tMap = getMap(i);
1130 return Teuchos::rcp(new VectorType(tMap));
1131}
1132
1133template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1134Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
1136getGhostedTpetraRangeVector(int i) const
1137{
1138 Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
1139 return Teuchos::rcp(new VectorType(tMap));
1140}
1141
1142template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1143int
1145getBlockRowCount() const
1146{
1147 return gidProviders_.size();
1148}
1149
1150template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1151int
1153getBlockColCount() const
1154{
1155 return gidProviders_.size();
1156}
1157
1158template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1160beginFill(LinearObjContainer & loc) const
1161{
1162 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1163 if(tloc.get_A()!=Teuchos::null)
1164 tloc.beginFill();
1165}
1166
1167template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
1169endFill(LinearObjContainer & loc) const
1170{
1171 BTLOC & tloc = Teuchos::dyn_cast<BTLOC>(loc);
1172 if(tloc.get_A()!=Teuchos::null)
1173 tloc.endFill();
1174}
1175
1176}
1177
1178#endif // __Panzer_BlockedTpetraLinearObjFactory_impl_hpp__
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
void set_dxdt(const Teuchos::RCP< VectorType > &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.
void initializeContainer(int, LinearObjContainer &loc) const
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)