Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
64 #include "Panzer_GlobalIndexer.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 
76 using Teuchos::RCP;
77 
78 namespace panzer {
79 
80 // ************************************************************
81 // class BlockedEpetraLinearObjFactory
82 // ************************************************************
83 
84 template <typename Traits,typename LocalOrdinalT>
87  const Teuchos::RCP<const GlobalIndexer> & gidProvider,
88  bool useDiscreteAdjoint)
89  : useColGidProviders_(false), eComm_(Teuchos::null)
90  , rawMpiComm_(comm->getRawMpiComm())
91  , useDiscreteAdjoint_(useDiscreteAdjoint)
92 {
95 
97 
98  makeRoomForBlocks(rowDOFManagerContainer_->getFieldBlocks());
99 
100  // build and register the gather/scatter evaluators with
101  // the base class.
102  this->buildGatherScatterEvaluators(*this);
103 
105 }
106 
107 template <typename Traits,typename LocalOrdinalT>
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 {
119 
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 
131 }
132 
133 template <typename Traits,typename LocalOrdinalT>
135 { }
136 
137 // LinearObjectFactory functions
139 
140 template <typename Traits,typename LocalOrdinalT>
141 void
143 readVector(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 
151 
152  // extract the vector from linear object container
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 
190 template <typename Traits,typename LocalOrdinalT>
191 void
193 writeVector(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 
201 
202  // extract the vector from linear object container
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 
236 template <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 
252  container->setMapsForBlocks(blockMaps);
253 
254  return container;
255 }
256 
257 template <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 
273  container->setMapsForBlocks(blockMaps);
274 
275  return container;
276 }
277 
278 template <typename Traits,typename LocalOrdinalT>
280  LinearObjContainer & out,int mem) const
281 {
282  using Teuchos::is_null;
283 
284  typedef LinearObjContainer LOC;
285  typedef BlockedEpetraLinearObjContainer BLOC;
286 
287  if( !rowDOFManagerContainer_->containsBlockedDOFManager()
288  && !colDOFManagerContainer_->containsBlockedDOFManager()) {
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 
320 template <typename Traits,typename LocalOrdinalT>
322  LinearObjContainer & out,int mem) const
323 {
324  using Teuchos::is_null;
325 
326  typedef LinearObjContainer LOC;
327  typedef BlockedEpetraLinearObjContainer BLOC;
328 
329  if( !rowDOFManagerContainer_->containsBlockedDOFManager()
330  && !colDOFManagerContainer_->containsBlockedDOFManager()) {
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 
362 template <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;
373  using Thyra::LinearOpBase;
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
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 
450 template <typename Traits,typename LocalOrdinalT>
453  const Epetra_Vector & global_bcs,
454  const Teuchos::Ptr<Epetra_Vector> & f,
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 
499 template <typename Traits,typename LocalOrdinalT>
502  LinearObjContainer & result) const
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);
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 //
542 template<typename Traits, typename LocalOrdinalT>
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 //
584 template<typename Traits, typename LocalOrdinalT>
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 
621 template <typename Traits,typename LocalOrdinalT>
623 getComm() const
624 {
625  return *tComm_;
626 }
627 
628 template <typename Traits,typename LocalOrdinalT>
631 {
632  typedef ThyraObjContainer<double> TOC;
633 
634  TOC & toc = Teuchos::dyn_cast<TOC>(loc);
635  initializeContainer_internal(mem,toc);
636 }
637 
638 template <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()) {
649  typedef BlockedEpetraLinearObjContainer BLOC;
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 {
661 
662  if((mem & LOC::F) == LOC::F)
664 
665  if((mem & LOC::Mat) == LOC::Mat)
667  }
668 }
669 
670 // Generic methods
672 
673 template <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 
694 template <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 
715 template <typename Traits,typename LocalOrdinalT>
717 addExcludedPair(int rowBlock,int colBlock)
718 {
719  excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
720 }
721 
722 template <typename Traits,typename LocalOrdinalT>
724 addExcludedPairs(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 
730 template <typename Traits,typename LocalOrdinalT>
732 getGlobalIndexer(int i) const
733 {
734  return rowDOFManagerContainer_->getFieldDOFManagers()[i];
735 }
736 
737 template <typename Traits,typename LocalOrdinalT>
740 {
741  return colDOFManagerContainer_->getFieldDOFManagers()[i];
742 }
743 
745 //
746 // makeRoomForBlocks()
747 //
749 template<typename Traits, typename LocalOrdinalT>
750 void
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 
776 template <typename Traits,typename 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 
799 template <typename Traits,typename LocalOrdinalT>
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 
822 template <typename Traits,typename LocalOrdinalT>
825 {
827  Thyra::createMember<double>(*getThyraDomainSpace());
828  Thyra::assign(vec.ptr(),0.0);
829 
830  return vec;
831 }
832 
833 template <typename Traits,typename LocalOrdinalT>
836 {
838  Thyra::createMember<double>(*getThyraRangeSpace());
839  Thyra::assign(vec.ptr(),0.0);
840 
841  return vec;
842 }
843 
844 template <typename Traits,typename LocalOrdinalT>
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 //
885 template<typename Traits, typename LocalOrdinalT>
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 //
920 template<typename Traits, typename LocalOrdinalT>
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 
950 template <typename Traits,typename 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 
973 template <typename Traits,typename LocalOrdinalT>
976 {
978  Thyra::createMember<double>(*getGhostedThyraDomainSpace());
979  Thyra::assign(vec.ptr(),0.0);
980 
981  return vec;
982 }
983 
984 template <typename Traits,typename LocalOrdinalT>
987 {
989  Thyra::createMember<double>(*getGhostedThyraRangeSpace());
990  Thyra::assign(vec.ptr(),0.0);
991 
992  return vec;
993 }
994 
995 template <typename Traits,typename 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 
1031 template <typename Traits,typename LocalOrdinalT>
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
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 
1068 template <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 
1114 template <typename Traits,typename LocalOrdinalT>
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 
1155 template <typename Traits,typename LocalOrdinalT>
1157 ghostToGlobalEpetraVector(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 
1168 template <typename Traits,typename LocalOrdinalT>
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 
1181 template <typename Traits,typename LocalOrdinalT>
1183 globalToGhostEpetraVector(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
1195 template <typename Traits,typename LocalOrdinalT>
1197 getMap(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
1206 template <typename Traits,typename LocalOrdinalT>
1208 getColMap(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 //
1224 template<typename Traits, typename LocalOrdinalT>
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 //
1240 template<typename Traits, typename LocalOrdinalT>
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 //
1256 template<typename Traits, typename LocalOrdinalT>
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 //
1274 template<typename Traits, typename LocalOrdinalT>
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
1288 template <typename Traits,typename LocalOrdinalT>
1290 getGraph(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 
1307 template <typename Traits,typename LocalOrdinalT>
1309 getGhostedGraph(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 //
1331 template<typename Traits, typename LocalOrdinalT>
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 //
1348 template<typename Traits, typename LocalOrdinalT>
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 //
1365 template<typename Traits, typename LocalOrdinalT>
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 //
1385 template<typename Traits, typename LocalOrdinalT>
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 //
1405 template<typename Traits, typename LocalOrdinalT>
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 //
1422 template<typename Traits, typename LocalOrdinalT>
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 //
1439 template<typename Traits, typename LocalOrdinalT>
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 //
1459 template<typename Traits, typename LocalOrdinalT>
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 
1474 template <typename Traits,typename LocalOrdinalT>
1476 buildMap(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 
1486 template <typename Traits,typename LocalOrdinalT>
1488 buildColMap(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 //
1506 template<typename Traits, typename LocalOrdinalT>
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 //
1524 template<typename Traits, typename LocalOrdinalT>
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 //
1542 template<typename Traits, typename LocalOrdinalT>
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 //
1562 template<typename Traits, typename LocalOrdinalT>
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
1578 template <typename Traits,typename LocalOrdinalT>
1580 buildGraph(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 
1608 template <typename Traits,typename LocalOrdinalT>
1610 buildGhostedGraph(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 
1672 template <typename Traits,typename LocalOrdinalT>
1674 buildFilteredGhostedGraph(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 
1721 template <typename Traits,typename LocalOrdinalT>
1723 getEpetraMatrix(int i,int j) const
1724 {
1725  Teuchos::RCP<Epetra_CrsGraph> eGraph = getGraph(i,j);
1727  TEUCHOS_ASSERT(mat->Filled());
1728  return mat;
1729 }
1730 
1731 template <typename Traits,typename LocalOrdinalT>
1733 getGhostedEpetraMatrix(int i,int j) const
1734 {
1735  Teuchos::RCP<Epetra_CrsGraph> eGraph = getGhostedGraph(i,j);
1737  TEUCHOS_ASSERT(mat->Filled());
1738  return mat;
1739 }
1740 
1741 template <typename Traits,typename LocalOrdinalT>
1744 {
1745  return eComm_;
1746 }
1747 
1748 template <typename Traits,typename LocalOrdinalT>
1751 {
1752  return rowDOFManagerContainer_->getFieldBlocks();
1753 }
1754 
1755 template <typename Traits,typename LocalOrdinalT>
1758 {
1759  return colDOFManagerContainer_->getFieldBlocks();
1760 }
1761 
1762 }
1763 
1764 #endif // __Panzer_BlockedEpetraLinearObjFactory_impl_hpp__
RCP< const T > getConst() const
void getOwnedAndGhostedNotFilteredIndicator(std::vector< int > &indicator) const
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...
Teuchos::RCP< Epetra_CrsMatrix > getEpetraMatrix(int i, int j) const
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraDomainVector() const
Get a domain vector.
bool is_null(const boost::shared_ptr< T > &p)
virtual const Teuchos::RCP< Epetra_Export > getGhostedColExport(int j) const
get exporter for converting an overalapped object to a &quot;normal&quot; object
virtual const Teuchos::RCP< Epetra_Map > getGhostedColMap(int i) const
get the ghosted map from the matrix
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
const Teuchos::RCP< Epetra_Vector > get_x() const
virtual void writeVector(const std::string &identifier, const LinearObjContainer &loc, int id) const
void initializeContainer_internal(int mem, ThyraObjContainer< double > &loc) const
bool is_null(const std::shared_ptr< T > &p)
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...
virtual void set_x_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
const std::vector< Teuchos::RCP< GlobalIndexer > > & getFieldDOFManagers() const
T_To & dyn_cast(T_From &from)
Teuchos::RCP< Thyra::VectorBase< double > > getThyraDomainVector() const
Get a domain vector.
int NumMyRows() const
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
void initializeContainer(int, LinearObjContainer &loc) const
virtual const std::vector< panzer::LocalOrdinal > & getElementBlock(const std::string &blockId) const =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
Teuchos::RCP< const GlobalIndexer > getColGlobalIndexer(int i) 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.
virtual void getElementBlockIds(std::vector< std::string > &elementBlockIds) const =0
size_type size() const
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
This class encapsulates the needs of a gather operation to do a // halo exchange for blocked vectors...
virtual const Teuchos::RCP< const Epetra_Comm > getEpetraComm() const
get exporter for converting an overalapped object to a &quot;normal&quot; object
Teuchos::RCP< Thyra::LinearOpBase< double > > getThyraMatrix() const
Get a Thyra operator.
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< Thyra::LinearOpBase< double > > getGhostedThyraMatrix() const
Get a Thyra operator.
virtual void getElementGIDsAsInt(panzer::LocalOrdinal localElmtId, std::vector< int > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
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_CrsGraph > buildGraph(int i, int j) const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraRangeVector() const
Get a range vector.
int PutScalar(double ScalarConstant)
virtual const std::vector< LocalOrdinal > & getAssociatedNeighbors(const LocalOrdinal &el) const =0
Teuchos::RCP< Epetra_CrsMatrix > getGhostedEpetraMatrix(int i, int j) const
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we&#39;re performing.
virtual const Teuchos::RCP< Epetra_CrsGraph > buildFilteredGhostedGraph(int i, int j) const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
virtual const Teuchos::RCP< Epetra_Map > getMap(int i) const
get the map from the matrix
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.
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
Teuchos::RCP< const panzer::BlockedDOFManager > getGlobalIndexer() const
Teuchos::RCP< const DOFManagerContainer > colDOFManagerContainer_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual const Teuchos::RCP< Epetra_Map > getGhostedMap2(int i) const
Get or create the i-th ghosted map.
void initializeGhostedContainer_internal(int mem, ThyraObjContainer< double > &loc) const
Teuchos::RCP< const Teuchos::OpaqueWrapper< MPI_Comm > > rawMpiComm_
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 const Teuchos::RCP< Epetra_CrsGraph > buildGhostedGraph(int i, int j, bool optimizeStorage) const
BlockedEpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const GlobalIndexer > &gidProvider, bool useDiscreteAdjoint=false)
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
virtual void set_f_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
This class provides a boundary exchange communication mechanism for vectors.
virtual const Teuchos::RCP< Epetra_Import > getGhostedImport(int i) const
get importer for converting an overalapped object to a &quot;normal&quot; object
Ptr< T > ptr() const
void ghostToGlobalEpetraMatrix(int blockRow, const Epetra_CrsMatrix &in, Epetra_CrsMatrix &out) const
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 Teuchos::RCP< const ConnManager > getConnManager() const =0
Returns the connection manager currently being used.
virtual const Teuchos::RCP< Epetra_Import > getGhostedColImport(int i) const
get importer for converting an overalapped object to a &quot;normal&quot; object
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
const Teuchos::RCP< Epetra_Vector > get_f() const
void ghostToGlobalEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
virtual const Teuchos::RCP< Epetra_Map > getColMap(int i) const
get the map from the matrix
This class provides a boundary exchange communication mechanism for vectors.
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...
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
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.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
virtual const Teuchos::RCP< Epetra_CrsGraph > getGraph(int i, int j) const
get the graph of the crs matrix
void initializeGhostedContainer(int, LinearObjContainer &loc) const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< double > &in, Thyra::LinearOpBase< double > &out) const
virtual const Teuchos::RCP< Epetra_Export > getGhostedExport(int j) const
get exporter for converting an overalapped object to a &quot;normal&quot; object
virtual const Teuchos::RCP< Epetra_Map > getGhostedColMap2(int i) const
Get or create the i-th ghosted column map.
void setMapsForBlocks(const std::vector< Teuchos::RCP< const Epetra_Map > > &blockMaps)
int VectorToMatrixMarketFile(const char *filename, const Epetra_Vector &A, const char *matrixName=0, const char *matrixDescription=0, bool writeHeader=true)
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...
int RemoveMyIndices(int LocalRow, int NumIndices, int *Indices)
Teuchos::RCP< Thyra::VectorBase< double > > getThyraRangeVector() const
Get a range vector.
virtual void set_A_th(const Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > &in)=0
bool Filled() const
virtual const Teuchos::RCP< Epetra_Map > getGhostedMap(int i) const
get the ghosted map from the matrix
int MatrixMarketFileToVector(const char *filename, const Epetra_BlockMap &map, Epetra_Vector *&A)
virtual void readVector(const std::string &identifier, LinearObjContainer &loc, int id) const
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildReadOnlyDomainContainer() const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
virtual bool hasAssociatedNeighbors() const =0
void buildGatherScatterEvaluators(const BuilderT &builder)
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
void globalToGhostEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
virtual Teuchos::RCP< WriteVector_GlobalEvaluationData > buildWriteDomainContainer() const
const Teuchos::RCP< Epetra_Vector > get_dxdt() const
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 Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< const DOFManagerContainer > rowDOFManagerContainer_
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraDomainSpace2() const
Get or create the ghosted Thyra domain space.
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.
bool is_null() const