Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_BlockedEpetraLinearObjFactory_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_BlockedEpetraLinearObjFactory_impl_hpp__
44#define __Panzer_BlockedEpetraLinearObjFactory_impl_hpp__
45
46
47// Epetra
48#include "Epetra_CrsMatrix.h"
49#include "Epetra_MpiComm.h"
50#include "Epetra_MultiVector.h"
51#include "Epetra_Vector.h"
52
53// EpetraExt
54#include "EpetraExt_VectorIn.h"
55#include "EpetraExt_VectorOut.h"
56
57// Panzer
63#include "Panzer_HashUtils.hpp"
65
66// Thyra
67#include "Thyra_DefaultBlockedLinearOp.hpp"
68#include "Thyra_DefaultProductVector.hpp"
69#include "Thyra_DefaultProductVectorSpace.hpp"
70#include "Thyra_EpetraLinearOp.hpp"
71#include "Thyra_EpetraThyraWrappers.hpp"
72#include "Thyra_get_Epetra_Operator.hpp"
73#include "Thyra_SpmdVectorBase.hpp"
74#include "Thyra_VectorStdOps.hpp"
75
76using Teuchos::RCP;
77
78namespace panzer {
79
80// ************************************************************
81// class BlockedEpetraLinearObjFactory
82// ************************************************************
83
84template <typename Traits,typename LocalOrdinalT>
86BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
87 const Teuchos::RCP<const GlobalIndexer> & gidProvider,
88 bool useDiscreteAdjoint)
89 : useColGidProviders_(false), eComm_(Teuchos::null)
90 , rawMpiComm_(comm->getRawMpiComm())
91 , useDiscreteAdjoint_(useDiscreteAdjoint)
92{
93 rowDOFManagerContainer_ = Teuchos::rcp(new DOFManagerContainer(gidProvider));
95
96 eComm_ = Teuchos::rcp(new Epetra_MpiComm(*rawMpiComm_));
97
99
100 // build and register the gather/scatter evaluators with
101 // the base class.
102 this->buildGatherScatterEvaluators(*this);
103
104 tComm_ = Teuchos::rcp(new Teuchos::MpiComm<int>(rawMpiComm_));
105}
106
107template <typename Traits,typename LocalOrdinalT>
109BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
110 const Teuchos::RCP<const GlobalIndexer> & gidProvider,
111 const Teuchos::RCP<const GlobalIndexer> & colGidProvider,
112 bool useDiscreteAdjoint)
113 : eComm_(Teuchos::null)
114 , rawMpiComm_(comm->getRawMpiComm())
115 , useDiscreteAdjoint_(useDiscreteAdjoint)
116{
117 rowDOFManagerContainer_ = Teuchos::rcp(new DOFManagerContainer(gidProvider));
118 colDOFManagerContainer_ = Teuchos::rcp(new DOFManagerContainer(colGidProvider));
119
120 eComm_ = Teuchos::rcp(new Epetra_MpiComm(*rawMpiComm_));
121
122 useColGidProviders_ = true;
123
124 makeRoomForBlocks(rowDOFManagerContainer_->getFieldBlocks(),colDOFManagerContainer_->getFieldBlocks());
125
126 // build and register the gather/scatter evaluators with
127 // the base class.
128 this->buildGatherScatterEvaluators(*this);
129
130 tComm_ = Teuchos::rcp(new Teuchos::MpiComm<int>(rawMpiComm_));
131}
132
133template <typename Traits,typename LocalOrdinalT>
136
137// LinearObjectFactory functions
139
140template <typename Traits,typename LocalOrdinalT>
141void
143readVector(const std::string & identifier,LinearObjContainer & loc,int id) const
144{
145 using Teuchos::RCP;
146 using Teuchos::rcp_dynamic_cast;
147 using Teuchos::dyn_cast;
149
150 BlockedEpetraLinearObjContainer & eloc = dyn_cast<BlockedEpetraLinearObjContainer>(loc);
151
152 // extract the vector from linear object container
153 RCP<Thyra::VectorBase<double> > vec;
154 switch(id) {
156 vec = eloc.get_x();
157 break;
159 vec = eloc.get_dxdt();
160 break;
162 vec = eloc.get_f();
163 break;
164 default:
165 TEUCHOS_ASSERT(false);
166 break;
167 };
168
169 int blockRows = this->getBlockRowCount();
170 RCP<ProductVectorBase<double> > b_vec = Thyra::nonconstProductVectorBase(vec);
171
172 // convert to Epetra then write out each vector to file
173 for(int i=0;i<blockRows;i++) {
174 RCP<Thyra::VectorBase<double> > x = b_vec->getNonconstVectorBlock(i);
175 RCP<Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(i),x);
176
177 // build the file name from the identifier
178 std::stringstream ss;
179 ss << identifier << "-" << i << ".mm";
180
181 // read in vector (wow the MM to Vector is a poorly designed interface!)
182 Epetra_Vector * ptr_ex = 0;
183 TEUCHOS_ASSERT(0==EpetraExt::MatrixMarketFileToVector(ss.str().c_str(),*getMap(i),ptr_ex));
184
185 *ex = *ptr_ex;
186 delete ptr_ex;
187 }
188}
189
190template <typename Traits,typename LocalOrdinalT>
191void
193writeVector(const std::string & identifier,const LinearObjContainer & loc,int id) const
194{
195 using Teuchos::RCP;
196 using Teuchos::rcp_dynamic_cast;
197 using Teuchos::dyn_cast;
199
200 const BlockedEpetraLinearObjContainer & eloc = dyn_cast<const BlockedEpetraLinearObjContainer>(loc);
201
202 // extract the vector from linear object container
203 RCP<const Thyra::VectorBase<double> > vec;
204 switch(id) {
206 vec = eloc.get_x();
207 break;
209 vec = eloc.get_dxdt();
210 break;
212 vec = eloc.get_f();
213 break;
214 default:
215 TEUCHOS_ASSERT(false);
216 break;
217 };
218
219 int blockRows = this->getBlockRowCount();
220 RCP<const ProductVectorBase<double> > b_vec = Thyra::productVectorBase(vec);
221
222 // convert to Epetra then write out each vector to file
223 for(int i=0;i<blockRows;i++) {
224 RCP<const Thyra::VectorBase<double> > x = b_vec->getVectorBlock(i);
225 RCP<const Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(i),x);
226
227 // build the file name from the identifier
228 std::stringstream ss;
229 ss << identifier << "-" << i << ".mm";
230
231 // write out vector
232 TEUCHOS_ASSERT(0==EpetraExt::VectorToMatrixMarketFile(ss.str().c_str(),*ex));
233 }
234}
235
236template <typename Traits,typename LocalOrdinalT>
238{
239 // if a "single field" DOFManager is used
240 if(!rowDOFManagerContainer_->containsBlockedDOFManager() && !colDOFManagerContainer_->containsBlockedDOFManager()) {
241 Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(new EpetraLinearObjContainer(getColMap(0),getMap(0)));
242
243 return container;
244 }
245
246 std::vector<Teuchos::RCP<const Epetra_Map> > blockMaps;
247 std::size_t blockDim = getBlockRowCount();
248 for(std::size_t i=0;i<blockDim;i++)
249 blockMaps.push_back(getMap(i));
250
251 Teuchos::RCP<BlockedEpetraLinearObjContainer > container = Teuchos::rcp(new BlockedEpetraLinearObjContainer);
252 container->setMapsForBlocks(blockMaps);
253
254 return container;
255}
256
257template <typename Traits,typename LocalOrdinalT>
259{
260 // if a "single field" DOFManager is used
261 if(!rowDOFManagerContainer_->containsBlockedDOFManager() && !colDOFManagerContainer_->containsBlockedDOFManager()) {
262 Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(new EpetraLinearObjContainer(getGhostedColMap(0),getGhostedMap(0)));
263
264 return container;
265 }
266
267 std::vector<Teuchos::RCP<const Epetra_Map> > blockMaps;
268 std::size_t blockDim = getBlockRowCount();
269 for(std::size_t i=0;i<blockDim;i++)
270 blockMaps.push_back(getGhostedMap(i));
271
272 Teuchos::RCP<BlockedEpetraLinearObjContainer > container = Teuchos::rcp(new BlockedEpetraLinearObjContainer);
273 container->setMapsForBlocks(blockMaps);
274
275 return container;
276}
277
278template <typename Traits,typename LocalOrdinalT>
280 LinearObjContainer & out,int mem) const
281{
282 using Teuchos::is_null;
283
284 typedef LinearObjContainer LOC;
286
287 if( !rowDOFManagerContainer_->containsBlockedDOFManager()
288 && !colDOFManagerContainer_->containsBlockedDOFManager()) {
289 const EpetraLinearObjContainer & e_in = Teuchos::dyn_cast<const EpetraLinearObjContainer>(in);
290 EpetraLinearObjContainer & e_out = Teuchos::dyn_cast<EpetraLinearObjContainer>(out);
291
292 // Operations occur if the GLOBAL container has the correct targets!
293 // Users set the GLOBAL continer arguments
294 if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
295 globalToGhostEpetraVector(0,*e_in.get_x(),*e_out.get_x(),true);
296
297 if ( !is_null(e_in.get_dxdt()) && !is_null(e_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
298 globalToGhostEpetraVector(0,*e_in.get_dxdt(),*e_out.get_dxdt(),true);
299
300 if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
301 globalToGhostEpetraVector(0,*e_in.get_f(),*e_out.get_f(),false);
302 }
303 else {
304 const BLOC & b_in = Teuchos::dyn_cast<const BLOC>(in);
305 BLOC & b_out = Teuchos::dyn_cast<BLOC>(out);
306
307 // Operations occur if the GLOBAL container has the correct targets!
308 // Users set the GLOBAL continer arguments
309 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
310 globalToGhostThyraVector(b_in.get_x(),b_out.get_x(),true);
311
312 if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
313 globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt(),true);
314
315 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
316 globalToGhostThyraVector(b_in.get_f(),b_out.get_f(),false);
317 }
318}
319
320template <typename Traits,typename LocalOrdinalT>
322 LinearObjContainer & out,int mem) const
323{
324 using Teuchos::is_null;
325
326 typedef LinearObjContainer LOC;
328
329 if( !rowDOFManagerContainer_->containsBlockedDOFManager()
330 && !colDOFManagerContainer_->containsBlockedDOFManager()) {
331 const EpetraLinearObjContainer & e_in = Teuchos::dyn_cast<const EpetraLinearObjContainer>(in);
332 EpetraLinearObjContainer & e_out = Teuchos::dyn_cast<EpetraLinearObjContainer>(out);
333
334 // Operations occur if the GLOBAL container has the correct targets!
335 // Users set the GLOBAL continer arguments
336 if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
337 ghostToGlobalEpetraVector(0,*e_in.get_x(),*e_out.get_x(),true);
338
339 if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
340 ghostToGlobalEpetraVector(0,*e_in.get_f(),*e_out.get_f(),false);
341
342 if ( !is_null(e_in.get_A()) && !is_null(e_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
343 ghostToGlobalEpetraMatrix(0,*e_in.get_A(),*e_out.get_A());
344 }
345 else {
346 const BLOC & b_in = Teuchos::dyn_cast<const BLOC>(in);
347 BLOC & b_out = Teuchos::dyn_cast<BLOC>(out);
348
349 // Operations occur if the GLOBAL container has the correct targets!
350 // Users set the GLOBAL continer arguments
351 if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
352 ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x(),true);
353
354 if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
355 ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f(),false);
356
357 if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
358 ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
359 }
360}
361
362template <typename Traits,typename LocalOrdinalT>
365 const LinearObjContainer & globalBCRows,
366 LinearObjContainer & ghostedObjs,
367 bool zeroVectorRows, bool adjustX) const
368{
369 typedef ThyraObjContainer<double> TOC;
370
371 using Teuchos::RCP;
372 using Teuchos::rcp_dynamic_cast;
374 using Thyra::PhysicallyBlockedLinearOpBase;
375 using Thyra::VectorBase;
377 using Thyra::get_Epetra_Vector;
378 using Thyra::get_Epetra_Operator;
379
380 int rBlockDim = getBlockRowCount();
381 int cBlockDim = getBlockColCount();
382
383 // first cast to block LOCs
384 const TOC & b_localBCRows = Teuchos::dyn_cast<const TOC>(localBCRows);
385 const TOC & b_globalBCRows = Teuchos::dyn_cast<const TOC>(globalBCRows);
386 TOC & b_ghosted = Teuchos::dyn_cast<TOC>(ghostedObjs);
387
388 TEUCHOS_ASSERT(b_localBCRows.get_f_th()!=Teuchos::null);
389 TEUCHOS_ASSERT(b_globalBCRows.get_f_th()!=Teuchos::null);
390
391 // cast each component as needed to their product form
392 RCP<PhysicallyBlockedLinearOpBase<double> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<double> >(b_ghosted.get_A_th());
393 if(A==Teuchos::null && b_ghosted.get_A_th()!=Teuchos::null) {
394 // assume it isn't physically blocked, for convenience physically block it
395 A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<double> >(Thyra::nonconstBlock1x1(b_ghosted.get_A_th()));
396 }
397
398 RCP<ProductVectorBase<double> > f = b_ghosted.get_f_th()==Teuchos::null
399 ? Teuchos::null
400 : Thyra::castOrCreateNonconstProductVectorBase(b_ghosted.get_f_th());
401 RCP<ProductVectorBase<double> > local_bcs = b_localBCRows.get_f_th()==Teuchos::null
402 ? Teuchos::null
403 : Thyra::castOrCreateNonconstProductVectorBase(b_localBCRows.get_f_th());
404 RCP<ProductVectorBase<double> > global_bcs = b_globalBCRows.get_f_th()==Teuchos::null
405 ? Teuchos::null
406 : Thyra::castOrCreateNonconstProductVectorBase(b_globalBCRows.get_f_th());
407
408 if(adjustX) f = Thyra::castOrCreateNonconstProductVectorBase(b_ghosted.get_x_th());
409
410 // sanity check!
411 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==rBlockDim);
412 if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==cBlockDim);
413 if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==rBlockDim);
414 TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==rBlockDim);
415 TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==rBlockDim);
416
417 for(int i=0;i<rBlockDim;i++) {
418 // grab epetra vector
419 RCP<const Epetra_Vector> e_local_bcs = get_Epetra_Vector(*getGhostedMap(i),local_bcs->getVectorBlock(i));
420 RCP<const Epetra_Vector> e_global_bcs = get_Epetra_Vector(*getGhostedMap(i),global_bcs->getVectorBlock(i));
421
422 // pull out epetra values
423 RCP<VectorBase<double> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
424 RCP<Epetra_Vector> e_f;
425 if(th_f==Teuchos::null)
426 e_f = Teuchos::null;
427 else
428 e_f = get_Epetra_Vector(*getGhostedMap(i),th_f);
429
430 for(int j=0;j<cBlockDim;j++) {
431
432 // pull out epetra values
433 RCP<LinearOpBase<double> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
434
435 // don't do anyting if opertor is null
436 RCP<Epetra_CrsMatrix> e_A;
437 if(th_A==Teuchos::null)
438 e_A = Teuchos::null;
439 else
440 e_A = rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*th_A),true);
441
442 // adjust Block operator
443 adjustForDirichletConditions(*e_local_bcs,*e_global_bcs,e_f.ptr(),e_A.ptr(),zeroVectorRows);
444
445 e_f = Teuchos::null; // this is so we only adjust it once on the first pass
446 }
447 }
448}
449
450template <typename Traits,typename LocalOrdinalT>
453 const Epetra_Vector & global_bcs,
454 const Teuchos::Ptr<Epetra_Vector> & f,
455 const Teuchos::Ptr<Epetra_CrsMatrix> & A,
456 bool zeroVectorRows) const
457{
458 if(f==Teuchos::null && A==Teuchos::null)
459 return;
460
461 TEUCHOS_ASSERT(local_bcs.MyLength()==global_bcs.MyLength());
462 for(int i=0;i<local_bcs.MyLength();i++) {
463 if(global_bcs[i]==0.0)
464 continue;
465
466 int numEntries = 0;
467 double * values = 0;
468 int * indices = 0;
469
470 if(local_bcs[i]==0.0 || zeroVectorRows) {
471 // this boundary condition was NOT set by this processor
472
473 // if they exist put 0.0 in each entry
474 if(!Teuchos::is_null(f))
475 (*f)[i] = 0.0;
476 if(!Teuchos::is_null(A)) {
477 A->ExtractMyRowView(i,numEntries,values,indices);
478 for(int c=0;c<numEntries;c++)
479 values[c] = 0.0;
480 }
481 }
482 else {
483 // this boundary condition was set by this processor
484
485 double scaleFactor = global_bcs[i];
486
487 // if they exist scale linear objects by scale factor
488 if(!Teuchos::is_null(f))
489 (*f)[i] /= scaleFactor;
490 if(!Teuchos::is_null(A)) {
491 A->ExtractMyRowView(i,numEntries,values,indices);
492 for(int c=0;c<numEntries;c++)
493 values[c] /= scaleFactor;
494 }
495 }
496 }
497}
498
499template <typename Traits,typename LocalOrdinalT>
503{
504 using Teuchos::RCP;
505 using Teuchos::rcp_dynamic_cast;
506 using Teuchos::dyn_cast;
507
508 typedef Thyra::ProductVectorBase<double> PVector;
509
510 const ThyraObjContainer<double> & th_counter = dyn_cast<const ThyraObjContainer<double> >(counter);
511 ThyraObjContainer<double> & th_result = dyn_cast<ThyraObjContainer<double> >(result);
512
513 RCP<const PVector> count = Thyra::castOrCreateProductVectorBase(th_counter.get_f_th().getConst());
514 RCP<const PVector> f_in = Thyra::castOrCreateProductVectorBase(th_counter.get_f_th().getConst());
515 RCP<PVector> f_out = Thyra::castOrCreateNonconstProductVectorBase(th_result.get_f_th());
516
517 int rBlockDim = getBlockRowCount();
518 for(int i=0;i<rBlockDim;i++) {
519
520 Teuchos::ArrayRCP<const double> count_array,f_in_array;
521 Teuchos::ArrayRCP<double> f_out_array;
522
523 rcp_dynamic_cast<const Thyra::SpmdVectorBase<double> >(count->getVectorBlock(i),true)->getLocalData(Teuchos::ptrFromRef(count_array));
524 rcp_dynamic_cast<const Thyra::SpmdVectorBase<double> >(f_in->getVectorBlock(i),true)->getLocalData(Teuchos::ptrFromRef(f_in_array));
525 rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(f_out->getNonconstVectorBlock(i),true)->getNonconstLocalData(Teuchos::ptrFromRef(f_out_array));
526
527 TEUCHOS_ASSERT(count_array.size()==f_in_array.size());
528 TEUCHOS_ASSERT(count_array.size()==f_out_array.size());
529
530 for(Teuchos::ArrayRCP<double>::size_type j=0;j<count_array.size();++j) {
531 if(count_array[j]!=0.0)
532 f_out_array[j] = f_in_array[j];
533 }
534 }
535}
536
538//
539// buildReadOnlyDomainContainer()
540//
542template<typename Traits, typename LocalOrdinalT>
543Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
546{
547 using std::vector;
548 using Teuchos::RCP;
549 using Teuchos::rcp;
553
554 // If a "single field" DOFManager is used, return a single
555 // EpetraVector_ReadOnly_GlobalEvaluationData.
556 if (not colDOFManagerContainer_->containsBlockedDOFManager())
557 {
558 auto ged = rcp(new EVROGED);
559 ged->initialize(getGhostedColImport2(0), getGhostedColMap2(0),
560 getColMap(0));
561 return ged;
562 } // end if a "single field" DOFManager is used
563
564 // Otherwise, return a BlockedVector_ReadOnly_GlobalEvaluationData.
565 vector<RCP<ROVGED>> gedBlocks;
566 for (int i(0); i < getBlockColCount(); ++i)
567 {
568 auto vecGed = rcp(new EVROGED);
569 vecGed->initialize(getGhostedColImport2(i), getGhostedColMap2(i),
570 getColMap(i));
571 gedBlocks.push_back(vecGed);
572 } // end loop over the blocks
573 auto ged = rcp(new BVROGED);
574 ged->initialize(getGhostedThyraDomainSpace2(), getThyraDomainSpace(),
575 gedBlocks);
576 return ged;
577} // end of buildReadOnlyDomainContainer()
578
580//
581// buildWriteDomainContainer()
582//
584template<typename Traits, typename LocalOrdinalT>
585Teuchos::RCP<WriteVector_GlobalEvaluationData>
588{
589 using std::vector;
590 using Teuchos::RCP;
591 using Teuchos::rcp;
595
596 // If a "single field" DOFManager is used, return a single
597 // EpetraVector_Write_GlobalEvaluationData.
598 if (not colDOFManagerContainer_->containsBlockedDOFManager())
599 {
600 auto ged = rcp(new EVWGED);
601 ged->initialize(getGhostedColExport2(0), getGhostedColMap2(0),
602 getColMap(0));
603 return ged;
604 } // end if a "single field" DOFManager is used
605
606 // Otherwise, return a BlockedVector_Write_GlobalEvaluationData.
607 vector<RCP<WVGED>> gedBlocks;
608 for (int i(0); i < getBlockColCount(); ++i)
609 {
610 auto vecGed = rcp(new EVWGED);
611 vecGed->initialize(getGhostedColExport2(i), getGhostedColMap2(i),
612 getColMap(i));
613 gedBlocks.push_back(vecGed);
614 } // end loop over the blocks
615 auto ged = rcp(new BVWGED);
616 ged->initialize(getGhostedThyraDomainSpace2(), getThyraDomainSpace(),
617 gedBlocks);
618 return ged;
619} // end of buildWriteDomainContainer()
620
621template <typename Traits,typename LocalOrdinalT>
627
628template <typename Traits,typename LocalOrdinalT>
630initializeContainer(int mem,LinearObjContainer & loc) const
631{
632 typedef ThyraObjContainer<double> TOC;
633
634 TOC & toc = Teuchos::dyn_cast<TOC>(loc);
635 initializeContainer_internal(mem,toc);
636}
637
638template <typename Traits,typename LocalOrdinalT>
641{
642 typedef LinearObjContainer LOC;
643 typedef ThyraObjContainer<double> TOC;
644
645 TOC & toc = Teuchos::dyn_cast<TOC>(loc);
646 initializeGhostedContainer_internal(mem,toc);
647
648 if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
650
651 BLOC & bloc = Teuchos::dyn_cast<BLOC>(loc);
652
653 if((mem & LOC::F) == LOC::F)
654 bloc.setRequiresDirichletAdjustment(true);
655
656 if((mem & LOC::Mat) == LOC::Mat)
657 bloc.setRequiresDirichletAdjustment(true);
658 }
659 else {
660 EpetraLinearObjContainer & eloc = Teuchos::dyn_cast<EpetraLinearObjContainer>(loc);
661
662 if((mem & LOC::F) == LOC::F)
664
665 if((mem & LOC::Mat) == LOC::Mat)
667 }
668}
669
670// Generic methods
672
673template <typename Traits,typename LocalOrdinalT>
676{
677 typedef LinearObjContainer LOC;
678
679 loc.clear();
680
681 if((mem & LOC::X) == LOC::X)
682 loc.set_x_th(getThyraDomainVector());
683
684 if((mem & LOC::DxDt) == LOC::DxDt)
685 loc.set_dxdt_th(getThyraDomainVector());
686
687 if((mem & LOC::F) == LOC::F)
688 loc.set_f_th(getThyraRangeVector());
689
690 if((mem & LOC::Mat) == LOC::Mat)
691 loc.set_A_th(getThyraMatrix());
692}
693
694template <typename Traits,typename LocalOrdinalT>
697{
698 typedef LinearObjContainer LOC;
699
700 loc.clear();
701
702 if((mem & LOC::X) == LOC::X)
703 loc.set_x_th(getGhostedThyraDomainVector());
704
705 if((mem & LOC::DxDt) == LOC::DxDt)
706 loc.set_dxdt_th(getGhostedThyraDomainVector());
707
708 if((mem & LOC::F) == LOC::F)
709 loc.set_f_th(getGhostedThyraRangeVector());
710
711 if((mem & LOC::Mat) == LOC::Mat)
712 loc.set_A_th(getGhostedThyraMatrix());
713}
714
715template <typename Traits,typename LocalOrdinalT>
717addExcludedPair(int rowBlock,int colBlock)
718{
719 excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
720}
721
722template <typename Traits,typename LocalOrdinalT>
724addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs)
725{
726 for(std::size_t i=0;i<exPairs.size();i++)
727 excludedPairs_.insert(exPairs[i]);
728}
729
730template <typename Traits,typename LocalOrdinalT>
732getGlobalIndexer(int i) const
733{
734 return rowDOFManagerContainer_->getFieldDOFManagers()[i];
735}
736
737template <typename Traits,typename LocalOrdinalT>
739getColGlobalIndexer(int i) const
740{
741 return colDOFManagerContainer_->getFieldDOFManagers()[i];
742}
743
745//
746// makeRoomForBlocks()
747//
749template<typename Traits, typename LocalOrdinalT>
750void
753 std::size_t blockCnt,
754 std::size_t colBlockCnt)
755{
756 maps_.resize(blockCnt);
757 ghostedMaps_.resize(blockCnt);
758 ghostedMaps2_.resize(blockCnt);
759 importers_.resize(blockCnt);
760 importers2_.resize(blockCnt);
761 exporters_.resize(blockCnt);
762 if (colBlockCnt > 0)
763 {
764 colMaps_.resize(colBlockCnt);
765 colGhostedMaps_.resize(colBlockCnt);
766 colGhostedMaps2_.resize(colBlockCnt);
767 colImporters_.resize(colBlockCnt);
768 colImporters2_.resize(colBlockCnt);
769 colExporters_.resize(colBlockCnt);
770 } // end if (colBlockCnt > 0)
771} // end of makeRoomForBlocks()
772
773// Thyra methods
775
776template <typename Traits,typename LocalOrdinalT>
777Teuchos::RCP<const Thyra::VectorSpaceBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
779{
780 if(domainSpace_==Teuchos::null) {
781 if(colDOFManagerContainer_->containsBlockedDOFManager()) {
782 // loop over all vectors and build the vector space
783 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
784 for(int i=0;i<getBlockColCount();i++)
785 vsArray.push_back(Thyra::create_VectorSpace(getColMap(i)));
786
787 domainSpace_ = Thyra::productVectorSpace<double>(vsArray);
788 }
789 else {
790 // the domain space is not blocked (just an SPMD vector), build it from
791 // the zeroth column
792 domainSpace_ = Thyra::create_VectorSpace(getColMap(0));
793 }
794 }
795
796 return domainSpace_;
797}
798
799template <typename Traits,typename LocalOrdinalT>
800Teuchos::RCP<const Thyra::VectorSpaceBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
801getThyraRangeSpace() const
802{
803 if(rangeSpace_==Teuchos::null) {
804 if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
805 // loop over all vectors and build the vector space
806 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
807 for(int i=0;i<getBlockRowCount();i++)
808 vsArray.push_back(Thyra::create_VectorSpace(getMap(i)));
809
810 rangeSpace_ = Thyra::productVectorSpace<double>(vsArray);
811 }
812 else {
813 // the range space is not blocked (just an SPMD vector), build it from
814 // the zeroth row
815 rangeSpace_ = Thyra::create_VectorSpace(getMap(0));
816 }
817 }
818
819 return rangeSpace_;
820}
821
822template <typename Traits,typename LocalOrdinalT>
823Teuchos::RCP<Thyra::VectorBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
825{
826 Teuchos::RCP<Thyra::VectorBase<double> > vec =
827 Thyra::createMember<double>(*getThyraDomainSpace());
828 Thyra::assign(vec.ptr(),0.0);
829
830 return vec;
831}
832
833template <typename Traits,typename LocalOrdinalT>
834Teuchos::RCP<Thyra::VectorBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
836{
837 Teuchos::RCP<Thyra::VectorBase<double> > vec =
838 Thyra::createMember<double>(*getThyraRangeSpace());
839 Thyra::assign(vec.ptr(),0.0);
840
841 return vec;
842}
843
844template <typename Traits,typename LocalOrdinalT>
845Teuchos::RCP<Thyra::LinearOpBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
846getThyraMatrix() const
847{
848 // return a flat epetra matrix
849 if(!rowDOFManagerContainer_->containsBlockedDOFManager() &&
850 !colDOFManagerContainer_->containsBlockedDOFManager()) {
851 return Thyra::nonconstEpetraLinearOp(getEpetraMatrix(0,0));
852 }
853
854 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockedOp = Thyra::defaultBlockedLinearOp<double>();
855
856 // get the block dimension
857 std::size_t rBlockDim = getBlockRowCount();
858 std::size_t cBlockDim = getBlockColCount();
859
860 // this operator will be square
861 blockedOp->beginBlockFill(rBlockDim,cBlockDim);
862
863 // loop over each block
864 for(std::size_t i=0;i<rBlockDim;i++) {
865 for(std::size_t j=0;j<cBlockDim;j++) {
866 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
867 // build (i,j) block matrix and add it to blocked operator
868 Teuchos::RCP<Thyra::LinearOpBase<double> > block = Thyra::nonconstEpetraLinearOp(getEpetraMatrix(i,j));
869 blockedOp->setNonconstBlock(i,j,block);
870 }
871 }
872 }
873
874 // all done
875 blockedOp->endBlockFill();
876
877 return blockedOp;
878}
879
881//
882// getGhostedThyraDomainSpace()
883//
885template<typename Traits, typename LocalOrdinalT>
886Teuchos::RCP<const Thyra::VectorSpaceBase<double>>
889{
890 using std::vector;
891 using Teuchos::RCP;
892 using Thyra::create_VectorSpace;
893 using Thyra::productVectorSpace;
895 if (ghostedDomainSpace_.is_null())
896 {
897 if (colDOFManagerContainer_->containsBlockedDOFManager())
898 {
899 // Loop over all vectors and build the vector space.
900 vector<RCP<const VectorSpaceBase<double>>> vsArray;
901 for (int i(0); i < getBlockColCount(); ++i)
902 vsArray.push_back(create_VectorSpace(getGhostedColMap(i)));
903 ghostedDomainSpace_ = productVectorSpace<double>(vsArray);
904 }
905 else // if (not colDOFManagerContainer_->containsBlockedDOFManager())
906 {
907 // The domain space is not blocked (that is, we're just dealing with a
908 // SPMD vector), so build it from the zeroth column.
909 ghostedDomainSpace_ = create_VectorSpace(getGhostedColMap(0));
910 } // end if (colDOFManagerContainer_->containsBlockedDOFManager()) or not
911 } // end if (ghostedDomainSpace_.is_null())
912 return ghostedDomainSpace_;
913} // end of getGhostedThyraDomainSpace()
914
916//
917// getGhostedThyraDomainSpace2()
918//
920template<typename Traits, typename LocalOrdinalT>
921Teuchos::RCP<const Thyra::VectorSpaceBase<double>>
924{
925 using std::vector;
926 using Teuchos::RCP;
927 using Thyra::create_VectorSpace;
928 using Thyra::productVectorSpace;
930 if (ghostedDomainSpace_.is_null())
931 {
932 if (colDOFManagerContainer_->containsBlockedDOFManager())
933 {
934 // Loop over all vectors and build the vector space.
935 vector<RCP<const VectorSpaceBase<double>>> vsArray;
936 for (int i(0); i < getBlockColCount(); ++i)
937 vsArray.push_back(create_VectorSpace(getGhostedColMap2(i)));
938 ghostedDomainSpace_ = productVectorSpace<double>(vsArray);
939 }
940 else // if (not colDOFManagerContainer_->containsBlockedDOFManager())
941 {
942 // The domain space is not blocked (that is, we're just dealing with a
943 // SPMD vector), so build it from the zeroth column.
944 ghostedDomainSpace_ = create_VectorSpace(getGhostedColMap2(0));
945 } // end if (colDOFManagerContainer_->containsBlockedDOFManager()) or not
946 } // end if (ghostedDomainSpace_.is_null())
947 return ghostedDomainSpace_;
948} // end of getGhostedThyraDomainSpace2()
949
950template <typename Traits,typename LocalOrdinalT>
951Teuchos::RCP<const Thyra::VectorSpaceBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
953{
954 if(ghostedRangeSpace_==Teuchos::null) {
955 if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
956 // loop over all vectors and build the vector space
957 std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
958 for(int i=0;i<getBlockRowCount();i++)
959 vsArray.push_back(Thyra::create_VectorSpace(getGhostedMap(i)));
960
961 ghostedRangeSpace_ = Thyra::productVectorSpace<double>(vsArray);
962 }
963 else {
964 // the range space is not blocked (just an SPMD vector), build it from
965 // the zeroth row
966 ghostedRangeSpace_ = Thyra::create_VectorSpace(getGhostedMap(0));
967 }
968 }
969
970 return ghostedRangeSpace_;
971}
972
973template <typename Traits,typename LocalOrdinalT>
974Teuchos::RCP<Thyra::VectorBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
976{
977 Teuchos::RCP<Thyra::VectorBase<double> > vec =
978 Thyra::createMember<double>(*getGhostedThyraDomainSpace());
979 Thyra::assign(vec.ptr(),0.0);
980
981 return vec;
982}
983
984template <typename Traits,typename LocalOrdinalT>
985Teuchos::RCP<Thyra::VectorBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
987{
988 Teuchos::RCP<Thyra::VectorBase<double> > vec =
989 Thyra::createMember<double>(*getGhostedThyraRangeSpace());
990 Thyra::assign(vec.ptr(),0.0);
991
992 return vec;
993}
994
995template <typename Traits,typename LocalOrdinalT>
996Teuchos::RCP<Thyra::LinearOpBase<double> > BlockedEpetraLinearObjFactory<Traits,LocalOrdinalT>::
998{
999 // return a flat epetra matrix
1000 if(!rowDOFManagerContainer_->containsBlockedDOFManager() &&
1001 !colDOFManagerContainer_->containsBlockedDOFManager()) {
1002 return Thyra::nonconstEpetraLinearOp(getGhostedEpetraMatrix(0,0));
1003 }
1004
1005 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockedOp = Thyra::defaultBlockedLinearOp<double>();
1006
1007 // get the block dimension
1008 std::size_t rBlockDim = getBlockRowCount();
1009 std::size_t cBlockDim = getBlockColCount();
1010
1011 // this operator will be square
1012 blockedOp->beginBlockFill(rBlockDim,cBlockDim);
1013
1014 // loop over each block
1015 for(std::size_t i=0;i<rBlockDim;i++) {
1016 for(std::size_t j=0;j<cBlockDim;j++) {
1017 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
1018 // build (i,j) block matrix and add it to blocked operator
1019 Teuchos::RCP<Thyra::LinearOpBase<double> > block = Thyra::nonconstEpetraLinearOp(getGhostedEpetraMatrix(i,j));
1020 blockedOp->setNonconstBlock(i,j,block);
1021 }
1022 }
1023 }
1024
1025 // all done
1026 blockedOp->endBlockFill();
1027
1028 return blockedOp;
1029}
1030
1031template <typename Traits,typename LocalOrdinalT>
1033ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<double> > & in,
1034 const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const
1035{
1036 using Teuchos::RCP;
1037 using Teuchos::rcp_dynamic_cast;
1039 using Thyra::get_Epetra_Vector;
1040
1041 std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
1042
1043 // get product vectors
1044 RCP<const ProductVectorBase<double> > prod_in = Thyra::castOrCreateProductVectorBase(in);
1045 RCP<ProductVectorBase<double> > prod_out = Thyra::castOrCreateNonconstProductVectorBase(out);
1046
1047 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
1048 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
1049
1050 for(std::size_t i=0;i<blockDim;i++) {
1051 // first get each Epetra vector
1052 RCP<const Epetra_Vector> ep_in;
1053 RCP<Epetra_Vector> ep_out;
1054
1055 if(not col) {
1056 ep_in = get_Epetra_Vector(*getGhostedMap(i),prod_in->getVectorBlock(i));
1057 ep_out = get_Epetra_Vector(*getMap(i),prod_out->getNonconstVectorBlock(i));
1058 } else {
1059 ep_in = get_Epetra_Vector(*getGhostedColMap(i),prod_in->getVectorBlock(i));
1060 ep_out = get_Epetra_Vector(*getColMap(i),prod_out->getNonconstVectorBlock(i));
1061 }
1062
1063 // use Epetra to do global communication
1064 ghostToGlobalEpetraVector(i,*ep_in,*ep_out,col);
1065 }
1066}
1067
1068template <typename Traits,typename LocalOrdinalT>
1071{
1072 using Teuchos::RCP;
1073 using Teuchos::rcp_dynamic_cast;
1074 using Teuchos::dyn_cast;
1075 using Thyra::LinearOpBase;
1076 using Thyra::PhysicallyBlockedLinearOpBase;
1077 using Thyra::get_Epetra_Operator;
1078
1079 // get the block dimension
1080 std::size_t rBlockDim = getBlockRowCount();
1081 std::size_t cBlockDim = getBlockColCount();
1082
1083 // get product vectors
1084 const PhysicallyBlockedLinearOpBase<double> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<double> >(in);
1085 PhysicallyBlockedLinearOpBase<double> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<double> >(out);
1086
1087 TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) rBlockDim);
1088 TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) cBlockDim);
1089 TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) rBlockDim);
1090 TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) cBlockDim);
1091
1092 for(std::size_t i=0;i<rBlockDim;i++) {
1093 for(std::size_t j=0;j<cBlockDim;j++) {
1094 if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
1095 // extract the blocks
1096 RCP<const LinearOpBase<double> > th_in = prod_in.getBlock(i,j);
1097 RCP<LinearOpBase<double> > th_out = prod_out.getNonconstBlock(i,j);
1098
1099 // sanity check
1100 TEUCHOS_ASSERT(th_in!=Teuchos::null);
1101 TEUCHOS_ASSERT(th_out!=Teuchos::null);
1102
1103 // get the epetra version of the blocks
1104 RCP<const Epetra_CrsMatrix> ep_in = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*th_in),true);
1105 RCP<Epetra_CrsMatrix> ep_out = rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*th_out),true);
1106
1107 // use Epetra to do global communication
1108 ghostToGlobalEpetraMatrix(i,*ep_in,*ep_out);
1109 }
1110 }
1111 }
1112}
1113
1114template <typename Traits,typename LocalOrdinalT>
1116globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<double> > & in,
1117 const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const
1118{
1119 using Teuchos::RCP;
1120 using Teuchos::rcp_dynamic_cast;
1122 using Thyra::get_Epetra_Vector;
1123
1124 std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
1125
1126 // get product vectors
1127 RCP<const ProductVectorBase<double> > prod_in = Thyra::castOrCreateProductVectorBase(in);
1128 RCP<ProductVectorBase<double> > prod_out = Thyra::castOrCreateNonconstProductVectorBase(out);
1129
1130 TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
1131 TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
1132
1133 for(std::size_t i=0;i<blockDim;i++) {
1134 // first get each Epetra vector
1135 RCP<const Epetra_Vector> ep_in;
1136 RCP<Epetra_Vector> ep_out;
1137
1138 if(not col) {
1139 ep_in = get_Epetra_Vector(*getMap(i),prod_in->getVectorBlock(i));
1140 ep_out = get_Epetra_Vector(*getGhostedMap(i),prod_out->getNonconstVectorBlock(i));
1141 }
1142 else {
1143 ep_in = get_Epetra_Vector(*getColMap(i),prod_in->getVectorBlock(i));
1144 ep_out = get_Epetra_Vector(*getGhostedColMap(i),prod_out->getNonconstVectorBlock(i));
1145 }
1146
1147 // use Epetra to do global communication
1148 globalToGhostEpetraVector(i,*ep_in,*ep_out,col);
1149 }
1150}
1151
1152// Epetra methods
1154
1155template <typename Traits,typename LocalOrdinalT>
1157ghostToGlobalEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const
1158{
1159 using Teuchos::RCP;
1160
1161 // do the global distribution
1162 RCP<Epetra_Export> exporter = col ? getGhostedColExport(i) : getGhostedExport(i);
1163 out.PutScalar(0.0);
1164 int err = out.Export(in,*exporter,Add);
1165 TEUCHOS_ASSERT_EQUALITY(err,0);
1166}
1167
1168template <typename Traits,typename LocalOrdinalT>
1170ghostToGlobalEpetraMatrix(int blockRow,const Epetra_CrsMatrix & in,Epetra_CrsMatrix & out) const
1171{
1172 using Teuchos::RCP;
1173
1174 // do the global distribution
1175 RCP<Epetra_Export> exporter = getGhostedExport(blockRow);
1176 out.PutScalar(0.0);
1177 int err = out.Export(in,*exporter,Add);
1178 TEUCHOS_ASSERT_EQUALITY(err,0);
1179}
1180
1181template <typename Traits,typename LocalOrdinalT>
1183globalToGhostEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const
1184{
1185 using Teuchos::RCP;
1186
1187 // do the global distribution
1188 RCP<Epetra_Import> importer = col ? getGhostedColImport(i) : getGhostedImport(i);
1189 out.PutScalar(0.0);
1190 int err = out.Import(in,*importer,Insert);
1191 TEUCHOS_ASSERT_EQUALITY(err,0);
1192}
1193
1194// get the map from the matrix
1195template <typename Traits,typename LocalOrdinalT>
1197getMap(int i) const
1198{
1199 if(maps_[i]==Teuchos::null)
1200 maps_[i] = buildMap(i);
1201
1202 return maps_[i];
1203}
1204
1205// get the map from the matrix
1206template <typename Traits,typename LocalOrdinalT>
1208getColMap(int i) const
1209{
1210 if(not useColGidProviders_)
1211 return getMap(i);
1212
1213 if(colMaps_[i]==Teuchos::null)
1214 colMaps_[i] = buildColMap(i);
1215
1216 return colMaps_[i];
1217}
1218
1220//
1221// getGhostedMap()
1222//
1224template<typename Traits, typename LocalOrdinalT>
1225const Teuchos::RCP<Epetra_Map>
1228 int i) const
1229{
1230 if (ghostedMaps_[i].is_null())
1231 ghostedMaps_[i] = buildGhostedMap(i);
1232 return ghostedMaps_[i];
1233} // end of getGhostedMap()
1234
1236//
1237// getGhostedMap2()
1238//
1240template<typename Traits, typename LocalOrdinalT>
1241const Teuchos::RCP<Epetra_Map>
1244 int i) const
1245{
1246 if (ghostedMaps2_[i].is_null())
1247 ghostedMaps2_[i] = buildGhostedMap2(i);
1248 return ghostedMaps2_[i];
1249} // end of getGhostedMap2()
1250
1252//
1253// getGhostedColMap()
1254//
1256template<typename Traits, typename LocalOrdinalT>
1257const Teuchos::RCP<Epetra_Map>
1260 int i) const
1261{
1262 if (not useColGidProviders_)
1263 return getGhostedMap(i);
1264 if (colGhostedMaps_[i].is_null())
1265 colGhostedMaps_[i] = buildColGhostedMap(i);
1266 return colGhostedMaps_[i];
1267} // end of getGhostedColMap()
1268
1270//
1271// getGhostedColMap2()
1272//
1274template<typename Traits, typename LocalOrdinalT>
1275const Teuchos::RCP<Epetra_Map>
1278 int i) const
1279{
1280 if (not useColGidProviders_)
1281 return getGhostedMap2(i);
1282 if (colGhostedMaps2_[i].is_null())
1283 colGhostedMaps2_[i] = buildColGhostedMap2(i);
1284 return colGhostedMaps2_[i];
1285} // end of getGhostedColMap2()
1286
1287// get the graph of the crs matrix
1288template <typename Traits,typename LocalOrdinalT>
1290getGraph(int i,int j) const
1291{
1292 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> GraphMap;
1293
1294 GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
1295 Teuchos::RCP<Epetra_CrsGraph> graph;
1296 if(itr==graphs_.end()) {
1297 graph = buildGraph(i,j);
1298 graphs_[std::make_pair(i,j)] = graph;
1299 }
1300 else
1301 graph = itr->second;
1302
1303 TEUCHOS_ASSERT(graph!=Teuchos::null);
1304 return graph;
1305}
1306
1307template <typename Traits,typename LocalOrdinalT>
1309getGhostedGraph(int i,int j) const
1310{
1311 typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> GraphMap;
1312
1313 GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
1314 Teuchos::RCP<Epetra_CrsGraph> ghostedGraph;
1315 if(itr==ghostedGraphs_.end()) {
1316 ghostedGraph = buildGhostedGraph(i,j,true);
1317 ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
1318 }
1319 else
1320 ghostedGraph = itr->second;
1321
1322 TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
1323 return ghostedGraph;
1324}
1325
1327//
1328// getGhostedImport()
1329//
1331template<typename Traits, typename LocalOrdinalT>
1332const Teuchos::RCP<Epetra_Import>
1335 int i) const
1336{
1337 using Teuchos::rcp;
1338 if (importers_[i].is_null())
1339 importers_[i] = rcp(new Epetra_Import(*getGhostedMap(i), *getMap(i)));
1340 return importers_[i];
1341} // end of getGhostedImport()
1342
1344//
1345// getGhostedImport2()
1346//
1348template<typename Traits, typename LocalOrdinalT>
1349const Teuchos::RCP<Epetra_Import>
1352 int i) const
1353{
1354 using Teuchos::rcp;
1355 if (importers2_[i].is_null())
1356 importers2_[i] = rcp(new Epetra_Import(*getGhostedMap2(i), *getMap(i)));
1357 return importers2_[i];
1358} // end of getGhostedImport2()
1359
1361//
1362// getGhostedColImport()
1363//
1365template<typename Traits, typename LocalOrdinalT>
1366const Teuchos::RCP<Epetra_Import>
1369 int i) const
1370{
1371 using Teuchos::rcp;
1372 if (not useColGidProviders_)
1373 return getGhostedImport(i);
1374 if (colImporters_[i].is_null())
1375 colImporters_[i] =
1376 rcp(new Epetra_Import(*getGhostedColMap(i), *getColMap(i)));
1377 return colImporters_[i];
1378} // end of getGhostedColImport()
1379
1381//
1382// getGhostedColImport2()
1383//
1385template<typename Traits, typename LocalOrdinalT>
1386const Teuchos::RCP<Epetra_Import>
1389 int i) const
1390{
1391 using Teuchos::rcp;
1392 if (not useColGidProviders_)
1393 return getGhostedImport2(i);
1394 if (colImporters2_[i].is_null())
1395 colImporters2_[i] =
1396 rcp(new Epetra_Import(*getGhostedColMap2(i), *getColMap(i)));
1397 return colImporters2_[i];
1398} // end of getGhostedColImport2()
1399
1401//
1402// getGhostedExport()
1403//
1405template<typename Traits, typename LocalOrdinalT>
1406const Teuchos::RCP<Epetra_Export>
1409 int i) const
1410{
1411 using Teuchos::rcp;
1412 if (exporters_[i].is_null())
1413 exporters_[i] = rcp(new Epetra_Export(*getGhostedMap(i), *getMap(i)));
1414 return exporters_[i];
1415} // end of getGhostedExport()
1416
1418//
1419// getGhostedExport2()
1420//
1422template<typename Traits, typename LocalOrdinalT>
1423const Teuchos::RCP<Epetra_Export>
1426 int i) const
1427{
1428 using Teuchos::rcp;
1429 if (exporters_[i].is_null())
1430 exporters_[i] = rcp(new Epetra_Export(*getGhostedMap2(i), *getMap(i)));
1431 return exporters_[i];
1432} // end of getGhostedExport2()
1433
1435//
1436// getGhostedColExport()
1437//
1439template<typename Traits, typename LocalOrdinalT>
1440const Teuchos::RCP<Epetra_Export>
1443 int i) const
1444{
1445 using Teuchos::rcp;
1446 if (not useColGidProviders_)
1447 return getGhostedExport(i);
1448 if (colExporters_[i].is_null())
1449 colExporters_[i] = rcp(new Epetra_Export(*getGhostedColMap(i),
1450 *getColMap(i)));
1451 return colExporters_[i];
1452} // end of getGhostedColExport()
1453
1455//
1456// getGhostedColExport2()
1457//
1459template<typename Traits, typename LocalOrdinalT>
1460const Teuchos::RCP<Epetra_Export>
1463 int i) const
1464{
1465 using Teuchos::rcp;
1466 if (not useColGidProviders_)
1467 return getGhostedExport2(i);
1468 if (colExporters_[i].is_null())
1469 colExporters_[i] = rcp(new Epetra_Export(*getGhostedColMap2(i),
1470 *getColMap(i)));
1471 return colExporters_[i];
1472} // end of getGhostedColExport2()
1473
1474template <typename Traits,typename LocalOrdinalT>
1476buildMap(int i) const
1477{
1478 std::vector<int> indices;
1479
1480 // get the global indices
1481 getGlobalIndexer(i)->getOwnedIndicesAsInt(indices);
1482
1483 return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1484}
1485
1486template <typename Traits,typename LocalOrdinalT>
1488buildColMap(int i) const
1489{
1490 if(not useColGidProviders_)
1491 return buildMap(i);
1492
1493 std::vector<int> indices;
1494
1495 // get the global indices
1496 getColGlobalIndexer(i)->getOwnedIndicesAsInt(indices);
1497
1498 return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1499}
1500
1502//
1503// buildGhostedMap()
1504//
1506template<typename Traits, typename LocalOrdinalT>
1507const Teuchos::RCP<Epetra_Map>
1510 int i) const
1511{
1512 using std::vector;
1513 using Teuchos::rcp;
1514 vector<int> indices;
1515 getGlobalIndexer(i)->getOwnedAndGhostedIndicesAsInt(indices);
1516 return rcp(new Epetra_Map(-1, indices.size(), &indices[0], 0, *eComm_));
1517} // end of buildGhostedMap()
1518
1520//
1521// buildGhostedMap2()
1522//
1524template<typename Traits, typename LocalOrdinalT>
1525const Teuchos::RCP<Epetra_Map>
1528 int i) const
1529{
1530 using std::vector;
1531 using Teuchos::rcp;
1532 vector<int> indices;
1533 getGlobalIndexer(i)->getGhostedIndicesAsInt(indices);
1534 return rcp(new Epetra_Map(-1, indices.size(), &indices[0], 0, *eComm_));
1535} // end of buildGhostedMap2()
1536
1538//
1539// buildColGhostedMap()
1540//
1542template<typename Traits, typename LocalOrdinalT>
1543const Teuchos::RCP<Epetra_Map>
1546 int i) const
1547{
1548 using std::vector;
1549 using Teuchos::rcp;
1550 if (not useColGidProviders_)
1551 return buildGhostedMap(i);
1552 vector<int> indices;
1553 getColGlobalIndexer(i)->getOwnedAndGhostedIndicesAsInt(indices);
1554 return rcp(new Epetra_Map(-1, indices.size(), &indices[0], 0, *eComm_));
1555} // end of buildColGhostedMap()
1556
1558//
1559// buildColGhostedMap2()
1560//
1562template<typename Traits, typename LocalOrdinalT>
1563const Teuchos::RCP<Epetra_Map>
1566 int i) const
1567{
1568 using std::vector;
1569 using Teuchos::rcp;
1570 if (not useColGidProviders_)
1571 return buildGhostedMap2(i);
1572 vector<int> indices;
1573 getColGlobalIndexer(i)->getGhostedIndicesAsInt(indices);
1574 return rcp(new Epetra_Map(-1, indices.size(), &indices[0], 0, *eComm_));
1575} // end of buildColGhostedMap2()
1576
1577// get the graph of the crs matrix
1578template <typename Traits,typename LocalOrdinalT>
1580buildGraph(int i,int j) const
1581{
1582 using Teuchos::RCP;
1583 using Teuchos::rcp;
1584
1585 // build the map and allocate the space for the graph and
1586 // grab the ghosted graph
1587 RCP<Epetra_Map> map_i = getMap(i);
1588 RCP<Epetra_Map> map_j = getColMap(j);
1589
1590 TEUCHOS_ASSERT(map_i!=Teuchos::null);
1591 TEUCHOS_ASSERT(map_j!=Teuchos::null);
1592
1593 RCP<Epetra_CrsGraph> graph = rcp(new Epetra_CrsGraph(Copy,*map_i,0));
1594 RCP<Epetra_CrsGraph> oGraph = buildFilteredGhostedGraph(i,j);
1595 // this is the only place buildFilteredGhostedGraph is called. That is because
1596 // only the unghosted graph should reflect any of the filtering.
1597
1598 // perform the communication to finish building graph
1599 RCP<Epetra_Export> exporter = getGhostedExport(i);
1600 int err = graph->Export( *oGraph, *exporter, Insert );
1601 TEUCHOS_ASSERT_EQUALITY(err,0);
1602 graph->FillComplete(*map_j,*map_i);
1603 graph->OptimizeStorage();
1604
1605 return graph;
1606}
1607
1608template <typename Traits,typename LocalOrdinalT>
1610buildGhostedGraph(int i,int j,bool optimizeStorage) const
1611{
1612 // build the map and allocate the space for the graph
1613 Teuchos::RCP<Epetra_Map> rowMap = getGhostedMap(i);
1614 Teuchos::RCP<Epetra_Map> colMap = getGhostedColMap(j);
1615 Teuchos::RCP<Epetra_CrsGraph> graph = Teuchos::rcp(new Epetra_CrsGraph(Copy,*rowMap,*colMap,0));
1616
1617 std::vector<std::string> elementBlockIds;
1618
1619 Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
1620
1621 rowProvider = getGlobalIndexer(i);
1622 colProvider = getColGlobalIndexer(j);
1623
1624 rowProvider->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
1625 // same element blocks
1626
1627 const Teuchos::RCP<const ConnManager> conn_mgr = colProvider->getConnManager();
1628 const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors();
1629
1630 // graph information about the mesh
1631 std::vector<std::string>::const_iterator blockItr;
1632 for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1633 std::string blockId = *blockItr;
1634
1635 // grab elements for this block
1636 const std::vector<LocalOrdinalT> & elements = rowProvider->getElementBlock(blockId); // each sub provider "should" have the
1637 // same elements in each element block
1638
1639 // get information about number of indicies
1640 std::vector<int> row_gids;
1641 std::vector<int> col_gids;
1642
1643 // loop over the elemnts
1644 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1645 rowProvider->getElementGIDsAsInt(elements[elmt],row_gids);
1646 colProvider->getElementGIDsAsInt(elements[elmt],col_gids);
1647
1648 if (han) {
1649 const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[elmt]);
1650 for (typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
1651 eit != aes.end(); ++eit) {
1652 std::vector<int> other_col_gids;
1653 colProvider->getElementGIDsAsInt(*eit, other_col_gids);
1654 col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
1655 }
1656 }
1657
1658 for(std::size_t row=0;row<row_gids.size();row++)
1659 graph->InsertGlobalIndices(row_gids[row],col_gids.size(),&col_gids[0]);
1660 }
1661 }
1662
1663 // finish filling the graph: Make sure the colmap and row maps coincide to
1664 // minimize calls to LID lookups
1665 graph->FillComplete(*colMap,*rowMap);
1666 if(optimizeStorage)
1667 graph->OptimizeStorage();
1668
1669 return graph;
1670}
1671
1672template <typename Traits,typename LocalOrdinalT>
1674buildFilteredGhostedGraph(int i,int j) const
1675{
1676 using Teuchos::RCP;
1677 using Teuchos::rcp;
1678 using Teuchos::rcp_dynamic_cast;
1679
1680 // figure out if the domain is filtered
1681 RCP<const Filtered_GlobalIndexer> filtered_ugi
1682 = rcp_dynamic_cast<const Filtered_GlobalIndexer>(getColGlobalIndexer(j));
1683
1684 // domain is unfiltered, a filtered graph is just the original ghosted graph
1685 if(filtered_ugi==Teuchos::null)
1686 return buildGhostedGraph(i,j,true);
1687
1688 // get all local indices that are active (i.e. unfiltered)
1689 std::vector<int> ghostedActive;
1690 filtered_ugi->getOwnedAndGhostedNotFilteredIndicator(ghostedActive);
1691
1692 // This will build a new ghosted graph without optimized storage so entries can be removed.
1693 Teuchos::RCP<Epetra_CrsGraph> filteredGraph = buildGhostedGraph(i,j,false);
1694 // false implies that storage is not optimzied
1695
1696 // remove filtered column entries
1697 for(int k=0;k<filteredGraph->NumMyRows();++k) {
1698 std::vector<int> removedIndices;
1699 int numIndices = 0;
1700 int * indices = 0;
1701 TEUCHOS_ASSERT(filteredGraph->ExtractMyRowView(k,numIndices,indices)==0);
1702
1703 for(int m=0;m<numIndices;++m) {
1704 if(ghostedActive[indices[m]]==0)
1705 removedIndices.push_back(indices[m]);
1706 }
1707
1708 TEUCHOS_ASSERT(filteredGraph->RemoveMyIndices(k,Teuchos::as<int>(removedIndices.size()),&removedIndices[0])==0);
1709 }
1710
1711 // finish filling the graph
1712 Teuchos::RCP<Epetra_Map> rowMap = getGhostedMap(i);
1713 Teuchos::RCP<Epetra_Map> colMap = getGhostedColMap(j);
1714
1715 TEUCHOS_ASSERT(filteredGraph->FillComplete(*colMap,*rowMap)==0);
1716 TEUCHOS_ASSERT(filteredGraph->OptimizeStorage()==0);
1717
1718 return filteredGraph;
1719}
1720
1721template <typename Traits,typename LocalOrdinalT>
1723getEpetraMatrix(int i,int j) const
1724{
1725 Teuchos::RCP<Epetra_CrsGraph> eGraph = getGraph(i,j);
1726 Teuchos::RCP<Epetra_CrsMatrix> mat = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *eGraph));
1727 TEUCHOS_ASSERT(mat->Filled());
1728 return mat;
1729}
1730
1731template <typename Traits,typename LocalOrdinalT>
1733getGhostedEpetraMatrix(int i,int j) const
1734{
1735 Teuchos::RCP<Epetra_CrsGraph> eGraph = getGhostedGraph(i,j);
1736 Teuchos::RCP<Epetra_CrsMatrix> mat = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *eGraph));
1737 TEUCHOS_ASSERT(mat->Filled());
1738 return mat;
1739}
1740
1741template <typename Traits,typename LocalOrdinalT>
1743getEpetraComm() const
1744{
1745 return eComm_;
1746}
1747
1748template <typename Traits,typename LocalOrdinalT>
1750getBlockRowCount() const
1751{
1752 return rowDOFManagerContainer_->getFieldBlocks();
1753}
1754
1755template <typename Traits,typename LocalOrdinalT>
1757getBlockColCount() const
1758{
1759 return colDOFManagerContainer_->getFieldBlocks();
1760}
1761
1762}
1763
1764#endif // __Panzer_BlockedEpetraLinearObjFactory_impl_hpp__
Insert
Add
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we're performing.
int PutScalar(double ScalarConstant)
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int MyLength() const
int PutScalar(double ScalarConstant)
virtual const Teuchos::RCP< Epetra_Map > buildGhostedMap2(int i) const
Build the i-th ghosted map from the ghosted indices of the i-th global indexer.
void ghostToGlobalEpetraMatrix(int blockRow, const Epetra_CrsMatrix &in, Epetra_CrsMatrix &out) const
void globalToGhostEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
void initializeGhostedContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< Thyra::LinearOpBase< double > > getThyraMatrix() const
Get a Thyra operator.
virtual const Teuchos::RCP< Epetra_Export > getGhostedColExport(int j) const
get exporter for converting an overalapped object to a "normal" object
virtual const Teuchos::RCP< Epetra_Map > buildMap(int i) const
Build the i-th owned map from the owned indices of the i-th global indexer.
virtual const Teuchos::RCP< Epetra_Map > getGhostedColMap(int i) const
get the ghosted map from the matrix
Teuchos::RCP< const panzer::BlockedDOFManager > getGlobalIndexer() const
virtual const Teuchos::RCP< Epetra_CrsGraph > buildFilteredGhostedGraph(int i, int j) const
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
virtual const Teuchos::RCP< const Epetra_Comm > getEpetraComm() const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraDomainSpace2() const
Get or create the ghosted Thyra domain space.
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
virtual const Teuchos::RCP< Epetra_Map > buildColGhostedMap2(int i) const
Build the i-th ghosted column map from the ghosted indices of the i-th (column) global indexer.
void makeRoomForBlocks(std::size_t blockCnt, std::size_t colBlockCnt=0)
Allocate the space in the std::vector objects so we can fill with appropriate Epetra data.
virtual void readVector(const std::string &identifier, LinearObjContainer &loc, int id) const
virtual const Teuchos::RCP< Epetra_Map > getGhostedMap(int i) const
get the ghosted map from the matrix
virtual const Teuchos::RCP< Epetra_Import > getGhostedColImport2(int i) const
Get or create the i-th ghosted column importer corresponding to the i-th ghosted column map.
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
Teuchos::RCP< const DOFManagerContainer > colDOFManagerContainer_
Teuchos::RCP< const DOFManagerContainer > rowDOFManagerContainer_
Teuchos::RCP< Thyra::VectorBase< double > > getThyraRangeVector() const
Get a range vector.
virtual const Teuchos::RCP< Epetra_CrsGraph > getGraph(int i, int j) const
get the graph of the crs matrix
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
virtual const Teuchos::RCP< Epetra_Map > buildGhostedMap(int i) const
Build the i-th ghosted map from the owned and ghosted indices of the i-th global indexer.
void ghostToGlobalEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraRangeSpace() const
Get the range vector space (f)
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< double > &in, Thyra::LinearOpBase< double > &out) const
Teuchos::RCP< const Teuchos::OpaqueWrapper< MPI_Comm > > rawMpiComm_
virtual const Teuchos::RCP< Epetra_Import > getGhostedImport2(int i) const
Get or create the i-th ghosted importer corresponding to the i-th ghosted map.
void initializeContainer_internal(int mem, ThyraObjContainer< double > &loc) const
virtual const Teuchos::RCP< Epetra_Map > getMap(int i) const
get the map from the matrix
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraDomainVector() const
Get a domain vector.
virtual const Teuchos::RCP< Epetra_Map > buildColMap(int i) const
Build the i-th owned column map from the owned indices of the i-th (column) global indexer.
virtual const Teuchos::RCP< Epetra_Export > getGhostedExport2(int i) const
Get or create the i-th ghosted exporter corresponding to the i-th ghosted map.
virtual const Teuchos::RCP< Epetra_Import > getGhostedImport(int i) const
get importer for converting an overalapped object to a "normal" object
Teuchos::RCP< Thyra::VectorBase< double > > getThyraDomainVector() const
Get a domain vector.
virtual const Teuchos::RCP< Epetra_CrsGraph > buildGhostedGraph(int i, int j, bool optimizeStorage) const
virtual const Teuchos::RCP< Epetra_Map > buildColGhostedMap(int i) const
Build the i-th ghosted column map from the owned and ghosted indices of the i-th (column) global inde...
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
BlockedEpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const GlobalIndexer > &gidProvider, bool useDiscreteAdjoint=false)
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraRangeVector() const
Get a range vector.
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
void initializeGhostedContainer_internal(int mem, ThyraObjContainer< double > &loc) const
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildReadOnlyDomainContainer() const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
Teuchos::RCP< Epetra_CrsMatrix > getGhostedEpetraMatrix(int i, int j) const
virtual const Teuchos::RCP< Epetra_Import > getGhostedColImport(int i) const
get importer for converting an overalapped object to a "normal" object
virtual const Teuchos::RCP< Epetra_CrsGraph > getGhostedGraph(int i, int j) const
get the ghosted graph of the crs matrix
virtual const Teuchos::RCP< Epetra_Export > getGhostedColExport2(int i) const
Get or create the i-th ghosted column exporter corresponding to the i-th ghosted column map.
virtual const Teuchos::RCP< Epetra_CrsGraph > buildGraph(int i, int j) const
Teuchos::RCP< Thyra::LinearOpBase< double > > getGhostedThyraMatrix() const
Get a Thyra operator.
virtual const Teuchos::RCP< Epetra_Map > getColMap(int i) const
get the map from the matrix
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
virtual const Teuchos::RCP< Epetra_Map > getGhostedMap2(int i) const
Get or create the i-th ghosted map.
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
virtual const Teuchos::RCP< Epetra_Map > getGhostedColMap2(int i) const
Get or create the i-th ghosted column map.
virtual Teuchos::RCP< WriteVector_GlobalEvaluationData > buildWriteDomainContainer() const
Teuchos::RCP< const GlobalIndexer > getColGlobalIndexer(int i) const
virtual const Teuchos::RCP< Epetra_Export > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< Epetra_CrsMatrix > getEpetraMatrix(int i, int j) const
virtual void writeVector(const std::string &identifier, const LinearObjContainer &loc, int id) const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors.
This class encapsulates the needs of a gather operation to do a // halo exchange for blocked vectors....
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
const Teuchos::RCP< Epetra_Vector > get_dxdt() const
const Teuchos::RCP< Epetra_Vector > get_f() const
const Teuchos::RCP< Epetra_Vector > get_x() const
This class provides a boundary exchange communication mechanism for vectors.
This class provides a boundary exchange communication mechanism for vectors.
void buildGatherScatterEvaluators(const BuilderT &builder)
virtual void set_f_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
virtual void set_A_th(const Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > &in)=0
virtual Teuchos::RCP< Thyra::VectorBase< ScalarT > > get_f_th() const =0
virtual void set_dxdt_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
virtual void set_x_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
int MatrixMarketFileToVector(const char *filename, const Epetra_BlockMap &map, Epetra_Vector *&A)
int VectorToMatrixMarketFile(const char *filename, const Epetra_Vector &A, const char *matrixName=0, const char *matrixDescription=0, bool writeHeader=true)