Panzer  Version of the Day
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_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
44 #define PANZER_BLOCKED_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
45 
48 #include "Panzer_HashUtils.hpp"
49 
50 #include "Epetra_MultiVector.h"
51 #include "Epetra_Vector.h"
52 #include "Epetra_CrsMatrix.h"
53 #include "Epetra_MpiComm.h"
54 
55 #include "EpetraExt_VectorOut.h"
56 #include "EpetraExt_VectorIn.h"
57 
58 #include "Thyra_EpetraThyraWrappers.hpp"
59 #include "Thyra_DefaultProductVectorSpace.hpp"
60 #include "Thyra_DefaultProductVector.hpp"
61 #include "Thyra_DefaultBlockedLinearOp.hpp"
62 #include "Thyra_EpetraLinearOp.hpp"
63 #include "Thyra_SpmdVectorBase.hpp"
64 #include "Thyra_get_Epetra_Operator.hpp"
65 #include "Thyra_VectorStdOps.hpp"
66 
69 
70 using Teuchos::RCP;
71 
72 namespace panzer {
73 
74 // ************************************************************
75 // class BlockedEpetraLinearObjFactory
76 // ************************************************************
77 
78 template <typename Traits,typename LocalOrdinalT>
80 BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
82  bool useDiscreteAdjoint)
83  : useColGidProviders_(false), eComm_(Teuchos::null)
84  , rawMpiComm_(comm->getRawMpiComm())
85  , useDiscreteAdjoint_(useDiscreteAdjoint)
86 {
89 
91 
92  makeRoomForBlocks(rowDOFManagerContainer_->getFieldBlocks());
93 
94  // build and register the gather/scatter evaluators with
95  // the base class.
96  this->buildGatherScatterEvaluators(*this);
97 
98  tComm_ = Teuchos::rcp(new Teuchos::MpiComm<int>(rawMpiComm_));
99 }
100 
101 template <typename Traits,typename LocalOrdinalT>
103 BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
105  const Teuchos::RCP<const UniqueGlobalIndexerBase> & colGidProvider,
106  bool useDiscreteAdjoint)
107  : eComm_(Teuchos::null)
108  , rawMpiComm_(comm->getRawMpiComm())
109  , useDiscreteAdjoint_(useDiscreteAdjoint)
110 {
113 
115 
116  useColGidProviders_ = true;
117 
118  makeRoomForBlocks(rowDOFManagerContainer_->getFieldBlocks(),colDOFManagerContainer_->getFieldBlocks());
119 
120  // build and register the gather/scatter evaluators with
121  // the base class.
122  this->buildGatherScatterEvaluators(*this);
123 
124  tComm_ = Teuchos::rcp(new Teuchos::MpiComm<int>(rawMpiComm_));
125 }
126 
127 template <typename Traits,typename LocalOrdinalT>
129 { }
130 
131 // LinearObjectFactory functions
133 
134 template <typename Traits,typename LocalOrdinalT>
135 void
137 readVector(const std::string & identifier,LinearObjContainer & loc,int id) const
138 {
139  using Teuchos::RCP;
140  using Teuchos::rcp_dynamic_cast;
141  using Teuchos::dyn_cast;
143 
145 
146  // extract the vector from linear object container
148  switch(id) {
150  vec = eloc.get_x();
151  break;
153  vec = eloc.get_dxdt();
154  break;
156  vec = eloc.get_f();
157  break;
158  default:
159  TEUCHOS_ASSERT(false);
160  break;
161  };
162 
163  int blockRows = this->getBlockRowCount();
164  RCP<ProductVectorBase<double> > b_vec = Thyra::nonconstProductVectorBase(vec);
165 
166  // convert to Epetra then write out each vector to file
167  for(int i=0;i<blockRows;i++) {
168  RCP<Thyra::VectorBase<double> > x = b_vec->getNonconstVectorBlock(i);
169  RCP<Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(i),x);
170 
171  // build the file name from the identifier
172  std::stringstream ss;
173  ss << identifier << "-" << i << ".mm";
174 
175  // read in vector (wow the MM to Vector is a poorly designed interface!)
176  Epetra_Vector * ptr_ex = 0;
177  TEUCHOS_ASSERT(0==EpetraExt::MatrixMarketFileToVector(ss.str().c_str(),*getMap(i),ptr_ex));
178 
179  *ex = *ptr_ex;
180  delete ptr_ex;
181  }
182 }
183 
184 template <typename Traits,typename LocalOrdinalT>
185 void
187 writeVector(const std::string & identifier,const LinearObjContainer & loc,int id) const
188 {
189  using Teuchos::RCP;
190  using Teuchos::rcp_dynamic_cast;
191  using Teuchos::dyn_cast;
193 
195 
196  // extract the vector from linear object container
198  switch(id) {
200  vec = eloc.get_x();
201  break;
203  vec = eloc.get_dxdt();
204  break;
206  vec = eloc.get_f();
207  break;
208  default:
209  TEUCHOS_ASSERT(false);
210  break;
211  };
212 
213  int blockRows = this->getBlockRowCount();
214  RCP<const ProductVectorBase<double> > b_vec = Thyra::productVectorBase(vec);
215 
216  // convert to Epetra then write out each vector to file
217  for(int i=0;i<blockRows;i++) {
218  RCP<const Thyra::VectorBase<double> > x = b_vec->getVectorBlock(i);
219  RCP<const Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(i),x);
220 
221  // build the file name from the identifier
222  std::stringstream ss;
223  ss << identifier << "-" << i << ".mm";
224 
225  // write out vector
226  TEUCHOS_ASSERT(0==EpetraExt::VectorToMatrixMarketFile(ss.str().c_str(),*ex));
227  }
228 }
229 
230 template <typename Traits,typename LocalOrdinalT>
232 {
233  // if a "single field" DOFManager is used
234  if(!rowDOFManagerContainer_->containsBlockedDOFManager() && !colDOFManagerContainer_->containsBlockedDOFManager()) {
235  Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(new EpetraLinearObjContainer(getColMap(0),getMap(0)));
236 
237  return container;
238  }
239 
240  std::vector<Teuchos::RCP<const Epetra_Map> > blockMaps;
241  std::size_t blockDim = getBlockRowCount();
242  for(std::size_t i=0;i<blockDim;i++)
243  blockMaps.push_back(getMap(i));
244 
246  container->setMapsForBlocks(blockMaps);
247 
248  return container;
249 }
250 
251 template <typename Traits,typename LocalOrdinalT>
253 {
254  // if a "single field" DOFManager is used
255  if(!rowDOFManagerContainer_->containsBlockedDOFManager() && !colDOFManagerContainer_->containsBlockedDOFManager()) {
256  Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(new EpetraLinearObjContainer(getGhostedColMap(0),getGhostedMap(0)));
257 
258  return container;
259  }
260 
261  std::vector<Teuchos::RCP<const Epetra_Map> > blockMaps;
262  std::size_t blockDim = getBlockRowCount();
263  for(std::size_t i=0;i<blockDim;i++)
264  blockMaps.push_back(getGhostedMap(i));
265 
267  container->setMapsForBlocks(blockMaps);
268 
269  return container;
270 }
271 
272 template <typename Traits,typename LocalOrdinalT>
274  LinearObjContainer & out,int mem) const
275 {
276  using Teuchos::is_null;
277 
278  typedef LinearObjContainer LOC;
279  typedef BlockedEpetraLinearObjContainer BLOC;
280 
281  if( !rowDOFManagerContainer_->containsBlockedDOFManager()
282  && !colDOFManagerContainer_->containsBlockedDOFManager()) {
285 
286  // Operations occur if the GLOBAL container has the correct targets!
287  // Users set the GLOBAL continer arguments
288  if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
289  globalToGhostEpetraVector(0,*e_in.get_x(),*e_out.get_x(),true);
290 
291  if ( !is_null(e_in.get_dxdt()) && !is_null(e_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
292  globalToGhostEpetraVector(0,*e_in.get_dxdt(),*e_out.get_dxdt(),true);
293 
294  if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
295  globalToGhostEpetraVector(0,*e_in.get_f(),*e_out.get_f(),false);
296  }
297  else {
298  const BLOC & b_in = Teuchos::dyn_cast<const BLOC>(in);
299  BLOC & b_out = Teuchos::dyn_cast<BLOC>(out);
300 
301  // Operations occur if the GLOBAL container has the correct targets!
302  // Users set the GLOBAL continer arguments
303  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
304  globalToGhostThyraVector(b_in.get_x(),b_out.get_x(),true);
305 
306  if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
307  globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt(),true);
308 
309  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
310  globalToGhostThyraVector(b_in.get_f(),b_out.get_f(),false);
311  }
312 }
313 
314 template <typename Traits,typename LocalOrdinalT>
316  LinearObjContainer & out,int mem) const
317 {
318  using Teuchos::is_null;
319 
320  typedef LinearObjContainer LOC;
321  typedef BlockedEpetraLinearObjContainer BLOC;
322 
323  if( !rowDOFManagerContainer_->containsBlockedDOFManager()
324  && !colDOFManagerContainer_->containsBlockedDOFManager()) {
327 
328  // Operations occur if the GLOBAL container has the correct targets!
329  // Users set the GLOBAL continer arguments
330  if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
331  ghostToGlobalEpetraVector(0,*e_in.get_x(),*e_out.get_x(),true);
332 
333  if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
334  ghostToGlobalEpetraVector(0,*e_in.get_f(),*e_out.get_f(),false);
335 
336  if ( !is_null(e_in.get_A()) && !is_null(e_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
337  ghostToGlobalEpetraMatrix(0,*e_in.get_A(),*e_out.get_A());
338  }
339  else {
340  const BLOC & b_in = Teuchos::dyn_cast<const BLOC>(in);
341  BLOC & b_out = Teuchos::dyn_cast<BLOC>(out);
342 
343  // Operations occur if the GLOBAL container has the correct targets!
344  // Users set the GLOBAL continer arguments
345  if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
346  ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x(),true);
347 
348  if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
349  ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f(),false);
350 
351  if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
352  ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
353  }
354 }
355 
356 template <typename Traits,typename LocalOrdinalT>
359  const LinearObjContainer & globalBCRows,
360  LinearObjContainer & ghostedObjs,
361  bool zeroVectorRows, bool adjustX) const
362 {
363  typedef ThyraObjContainer<double> TOC;
364 
365  using Teuchos::RCP;
366  using Teuchos::rcp_dynamic_cast;
367  using Thyra::LinearOpBase;
368  using Thyra::PhysicallyBlockedLinearOpBase;
369  using Thyra::VectorBase;
371  using Thyra::get_Epetra_Vector;
372  using Thyra::get_Epetra_Operator;
373 
374  int rBlockDim = getBlockRowCount();
375  int cBlockDim = getBlockColCount();
376 
377  // first cast to block LOCs
378  const TOC & b_localBCRows = Teuchos::dyn_cast<const TOC>(localBCRows);
379  const TOC & b_globalBCRows = Teuchos::dyn_cast<const TOC>(globalBCRows);
380  TOC & b_ghosted = Teuchos::dyn_cast<TOC>(ghostedObjs);
381 
382  TEUCHOS_ASSERT(b_localBCRows.get_f_th()!=Teuchos::null);
383  TEUCHOS_ASSERT(b_globalBCRows.get_f_th()!=Teuchos::null);
384 
385  // cast each component as needed to their product form
386  RCP<PhysicallyBlockedLinearOpBase<double> > A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<double> >(b_ghosted.get_A_th());
387  if(A==Teuchos::null && b_ghosted.get_A_th()!=Teuchos::null) {
388  // assume it isn't physically blocked, for convenience physically block it
389  A = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<double> >(Thyra::nonconstBlock1x1(b_ghosted.get_A_th()));
390  }
391 
392  RCP<ProductVectorBase<double> > f = b_ghosted.get_f_th()==Teuchos::null
393  ? Teuchos::null
394  : Thyra::castOrCreateNonconstProductVectorBase(b_ghosted.get_f_th());
395  RCP<ProductVectorBase<double> > local_bcs = b_localBCRows.get_f_th()==Teuchos::null
396  ? Teuchos::null
397  : Thyra::castOrCreateNonconstProductVectorBase(b_localBCRows.get_f_th());
398  RCP<ProductVectorBase<double> > global_bcs = b_globalBCRows.get_f_th()==Teuchos::null
399  ? Teuchos::null
400  : Thyra::castOrCreateNonconstProductVectorBase(b_globalBCRows.get_f_th());
401 
402  if(adjustX) f = Thyra::castOrCreateNonconstProductVectorBase(b_ghosted.get_x_th());
403 
404  // sanity check!
405  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==rBlockDim);
406  if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==cBlockDim);
407  if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==rBlockDim);
408  TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==rBlockDim);
409  TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==rBlockDim);
410 
411  for(int i=0;i<rBlockDim;i++) {
412  // grab epetra vector
413  RCP<const Epetra_Vector> e_local_bcs = get_Epetra_Vector(*getGhostedMap(i),local_bcs->getVectorBlock(i));
414  RCP<const Epetra_Vector> e_global_bcs = get_Epetra_Vector(*getGhostedMap(i),global_bcs->getVectorBlock(i));
415 
416  // pull out epetra values
417  RCP<VectorBase<double> > th_f = (f==Teuchos::null) ? Teuchos::null : f->getNonconstVectorBlock(i);
418  RCP<Epetra_Vector> e_f;
419  if(th_f==Teuchos::null)
420  e_f = Teuchos::null;
421  else
422  e_f = get_Epetra_Vector(*getGhostedMap(i),th_f);
423 
424  for(int j=0;j<cBlockDim;j++) {
425 
426  // pull out epetra values
427  RCP<LinearOpBase<double> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
428 
429  // don't do anyting if opertor is null
431  if(th_A==Teuchos::null)
432  e_A = Teuchos::null;
433  else
434  e_A = rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*th_A),true);
435 
436  // adjust Block operator
437  adjustForDirichletConditions(*e_local_bcs,*e_global_bcs,e_f.ptr(),e_A.ptr(),zeroVectorRows);
438 
439  e_f = Teuchos::null; // this is so we only adjust it once on the first pass
440  }
441  }
442 }
443 
444 template <typename Traits,typename LocalOrdinalT>
447  const Epetra_Vector & global_bcs,
448  const Teuchos::Ptr<Epetra_Vector> & f,
450  bool zeroVectorRows) const
451 {
452  if(f==Teuchos::null && A==Teuchos::null)
453  return;
454 
455  TEUCHOS_ASSERT(local_bcs.MyLength()==global_bcs.MyLength());
456  for(int i=0;i<local_bcs.MyLength();i++) {
457  if(global_bcs[i]==0.0)
458  continue;
459 
460  int numEntries = 0;
461  double * values = 0;
462  int * indices = 0;
463 
464  if(local_bcs[i]==0.0 || zeroVectorRows) {
465  // this boundary condition was NOT set by this processor
466 
467  // if they exist put 0.0 in each entry
468  if(!Teuchos::is_null(f))
469  (*f)[i] = 0.0;
470  if(!Teuchos::is_null(A)) {
471  A->ExtractMyRowView(i,numEntries,values,indices);
472  for(int c=0;c<numEntries;c++)
473  values[c] = 0.0;
474  }
475  }
476  else {
477  // this boundary condition was set by this processor
478 
479  double scaleFactor = global_bcs[i];
480 
481  // if they exist scale linear objects by scale factor
482  if(!Teuchos::is_null(f))
483  (*f)[i] /= scaleFactor;
484  if(!Teuchos::is_null(A)) {
485  A->ExtractMyRowView(i,numEntries,values,indices);
486  for(int c=0;c<numEntries;c++)
487  values[c] /= scaleFactor;
488  }
489  }
490  }
491 }
492 
493 template <typename Traits,typename LocalOrdinalT>
496  LinearObjContainer & result) const
497 {
498  using Teuchos::RCP;
499  using Teuchos::rcp_dynamic_cast;
500  using Teuchos::dyn_cast;
501 
502  typedef Thyra::ProductVectorBase<double> PVector;
503 
504  const ThyraObjContainer<double> & th_counter = dyn_cast<const ThyraObjContainer<double> >(counter);
506 
507  RCP<const PVector> count = Thyra::castOrCreateProductVectorBase(th_counter.get_f_th().getConst());
508  RCP<const PVector> f_in = Thyra::castOrCreateProductVectorBase(th_counter.get_f_th().getConst());
509  RCP<PVector> f_out = Thyra::castOrCreateNonconstProductVectorBase(th_result.get_f_th());
510 
511  int rBlockDim = getBlockRowCount();
512  for(int i=0;i<rBlockDim;i++) {
513 
514  Teuchos::ArrayRCP<const double> count_array,f_in_array;
515  Teuchos::ArrayRCP<double> f_out_array;
516 
517  rcp_dynamic_cast<const Thyra::SpmdVectorBase<double> >(count->getVectorBlock(i),true)->getLocalData(Teuchos::ptrFromRef(count_array));
518  rcp_dynamic_cast<const Thyra::SpmdVectorBase<double> >(f_in->getVectorBlock(i),true)->getLocalData(Teuchos::ptrFromRef(f_in_array));
519  rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(f_out->getNonconstVectorBlock(i),true)->getNonconstLocalData(Teuchos::ptrFromRef(f_out_array));
520 
521  TEUCHOS_ASSERT(count_array.size()==f_in_array.size());
522  TEUCHOS_ASSERT(count_array.size()==f_out_array.size());
523 
524  for(Teuchos::ArrayRCP<double>::size_type i=0;i<count_array.size();i++) {
525  if(count_array[i]!=0.0)
526  f_out_array[i] = f_in_array[i];
527  }
528  }
529 }
530 
531 template <typename Traits,typename LocalOrdinalT>
535 {
536  using Teuchos::RCP;
537  using Teuchos::rcp;
538 
539  // if a "single field" DOFManager is used, return a single RO_GED
540  if(!colDOFManagerContainer_->containsBlockedDOFManager()) {
542  ged->initialize(getGhostedColImport(0),getGhostedColMap(0),getColMap(0));
543 
544  return ged;
545  }
546 
547  std::vector<RCP<ReadOnlyVector_GlobalEvaluationData> > gedBlocks;
548  for(int i=0;i<getBlockColCount();i++) {
550  vec_ged->initialize(getGhostedColImport(i),getGhostedColMap(i),getColMap(i));
551 
552  gedBlocks.push_back(vec_ged);
553  }
554 
556  ged->initialize(getGhostedThyraDomainSpace(),getThyraDomainSpace(),gedBlocks);
557 
558  return ged;
559 }
560 
561 
562 template <typename Traits,typename LocalOrdinalT>
564 getComm() const
565 {
566  return *tComm_;
567 }
568 
569 template <typename Traits,typename LocalOrdinalT>
572 {
573  typedef ThyraObjContainer<double> TOC;
574 
575  TOC & toc = Teuchos::dyn_cast<TOC>(loc);
576  initializeContainer_internal(mem,toc);
577 }
578 
579 template <typename Traits,typename LocalOrdinalT>
582 {
583  typedef LinearObjContainer LOC;
584  typedef ThyraObjContainer<double> TOC;
585 
586  TOC & toc = Teuchos::dyn_cast<TOC>(loc);
587  initializeGhostedContainer_internal(mem,toc);
588 
589  if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
590  typedef BlockedEpetraLinearObjContainer BLOC;
591 
592  BLOC & bloc = Teuchos::dyn_cast<BLOC>(loc);
593 
594  if((mem & LOC::F) == LOC::F)
595  bloc.setRequiresDirichletAdjustment(true);
596 
597  if((mem & LOC::Mat) == LOC::Mat)
598  bloc.setRequiresDirichletAdjustment(true);
599  }
600  else {
602 
603  if((mem & LOC::F) == LOC::F)
605 
606  if((mem & LOC::Mat) == LOC::Mat)
608  }
609 }
610 
611 // Generic methods
613 
614 template <typename Traits,typename LocalOrdinalT>
617 {
618  typedef LinearObjContainer LOC;
619 
620  loc.clear();
621 
622  if((mem & LOC::X) == LOC::X)
623  loc.set_x_th(getThyraDomainVector());
624 
625  if((mem & LOC::DxDt) == LOC::DxDt)
626  loc.set_dxdt_th(getThyraDomainVector());
627 
628  if((mem & LOC::F) == LOC::F)
629  loc.set_f_th(getThyraRangeVector());
630 
631  if((mem & LOC::Mat) == LOC::Mat)
632  loc.set_A_th(getThyraMatrix());
633 }
634 
635 template <typename Traits,typename LocalOrdinalT>
638 {
639  typedef LinearObjContainer LOC;
640 
641  loc.clear();
642 
643  if((mem & LOC::X) == LOC::X)
644  loc.set_x_th(getGhostedThyraDomainVector());
645 
646  if((mem & LOC::DxDt) == LOC::DxDt)
647  loc.set_dxdt_th(getGhostedThyraDomainVector());
648 
649  if((mem & LOC::F) == LOC::F)
650  loc.set_f_th(getGhostedThyraRangeVector());
651 
652  if((mem & LOC::Mat) == LOC::Mat)
653  loc.set_A_th(getGhostedThyraMatrix());
654 }
655 
656 template <typename Traits,typename LocalOrdinalT>
658 addExcludedPair(int rowBlock,int colBlock)
659 {
660  excludedPairs_.insert(std::make_pair(rowBlock,colBlock));
661 }
662 
663 template <typename Traits,typename LocalOrdinalT>
665 addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs)
666 {
667  for(std::size_t i=0;i<exPairs.size();i++)
668  excludedPairs_.insert(exPairs[i]);
669 }
670 
671 template <typename Traits,typename LocalOrdinalT>
673 getGlobalIndexer(int i) const
674 {
675  return rowDOFManagerContainer_->getFieldDOFManagers()[i];
676 }
677 
678 template <typename Traits,typename LocalOrdinalT>
681 {
682  return colDOFManagerContainer_->getFieldDOFManagers()[i];
683 }
684 
685 template <typename Traits,typename LocalOrdinalT>
687 makeRoomForBlocks(std::size_t blockCnt,std::size_t colBlockCnt)
688 {
689  maps_.resize(blockCnt);
690  ghostedMaps_.resize(blockCnt);
691  importers_.resize(blockCnt);
692  exporters_.resize(blockCnt);
693 
694  if(colBlockCnt>0) {
695  colMaps_.resize(colBlockCnt);
696  colGhostedMaps_.resize(colBlockCnt);
697  colImporters_.resize(colBlockCnt);
698  colExporters_.resize(colBlockCnt);
699  }
700 }
701 
702 // Thyra methods
704 
705 template <typename Traits,typename LocalOrdinalT>
708 {
709  if(domainSpace_==Teuchos::null) {
710  if(colDOFManagerContainer_->containsBlockedDOFManager()) {
711  // loop over all vectors and build the vector space
712  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
713  for(int i=0;i<getBlockColCount();i++)
714  vsArray.push_back(Thyra::create_VectorSpace(getColMap(i)));
715 
716  domainSpace_ = Thyra::productVectorSpace<double>(vsArray);
717  }
718  else {
719  // the domain space is not blocked (just an SPMD vector), build it from
720  // the zeroth column
721  domainSpace_ = Thyra::create_VectorSpace(getColMap(0));
722  }
723  }
724 
725  return domainSpace_;
726 }
727 
728 template <typename Traits,typename LocalOrdinalT>
731 {
732  if(rangeSpace_==Teuchos::null) {
733  if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
734  // loop over all vectors and build the vector space
735  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
736  for(int i=0;i<getBlockRowCount();i++)
737  vsArray.push_back(Thyra::create_VectorSpace(getMap(i)));
738 
739  rangeSpace_ = Thyra::productVectorSpace<double>(vsArray);
740  }
741  else {
742  // the range space is not blocked (just an SPMD vector), build it from
743  // the zeroth row
744  rangeSpace_ = Thyra::create_VectorSpace(getMap(0));
745  }
746  }
747 
748  return rangeSpace_;
749 }
750 
751 template <typename Traits,typename LocalOrdinalT>
754 {
756  Thyra::createMember<double>(*getThyraDomainSpace());
757  Thyra::assign(vec.ptr(),0.0);
758 
759  return vec;
760 }
761 
762 template <typename Traits,typename LocalOrdinalT>
765 {
767  Thyra::createMember<double>(*getThyraRangeSpace());
768  Thyra::assign(vec.ptr(),0.0);
769 
770  return vec;
771 }
772 
773 template <typename Traits,typename LocalOrdinalT>
776 {
777  // return a flat epetra matrix
778  if(!rowDOFManagerContainer_->containsBlockedDOFManager() &&
779  !colDOFManagerContainer_->containsBlockedDOFManager()) {
780  return Thyra::nonconstEpetraLinearOp(getEpetraMatrix(0,0));
781  }
782 
783  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockedOp = Thyra::defaultBlockedLinearOp<double>();
784 
785  // get the block dimension
786  std::size_t rBlockDim = getBlockRowCount();
787  std::size_t cBlockDim = getBlockColCount();
788 
789  // this operator will be square
790  blockedOp->beginBlockFill(rBlockDim,cBlockDim);
791 
792  // loop over each block
793  for(std::size_t i=0;i<rBlockDim;i++) {
794  for(std::size_t j=0;j<cBlockDim;j++) {
795  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
796  // build (i,j) block matrix and add it to blocked operator
797  Teuchos::RCP<Thyra::LinearOpBase<double> > block = Thyra::nonconstEpetraLinearOp(getEpetraMatrix(i,j));
798  blockedOp->setNonconstBlock(i,j,block);
799  }
800  }
801  }
802 
803  // all done
804  blockedOp->endBlockFill();
805 
806  return blockedOp;
807 }
808 
809 template <typename Traits,typename LocalOrdinalT>
812 {
813  if(ghostedDomainSpace_==Teuchos::null) {
814  if(colDOFManagerContainer_->containsBlockedDOFManager()) {
815  // loop over all vectors and build the vector space
816  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
817  for(int i=0;i<getBlockColCount();i++)
818  vsArray.push_back(Thyra::create_VectorSpace(getGhostedColMap(i)));
819 
820  ghostedDomainSpace_ = Thyra::productVectorSpace<double>(vsArray);
821  }
822  else {
823  // the domain space is not blocked (just an SPMD vector), build it from
824  // the zeroth column
825  ghostedDomainSpace_ = Thyra::create_VectorSpace(getGhostedColMap(0));
826  }
827  }
828 
829  return ghostedDomainSpace_;
830 }
831 
832 template <typename Traits,typename LocalOrdinalT>
835 {
836  if(ghostedRangeSpace_==Teuchos::null) {
837  if(rowDOFManagerContainer_->containsBlockedDOFManager()) {
838  // loop over all vectors and build the vector space
839  std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > vsArray;
840  for(int i=0;i<getBlockRowCount();i++)
841  vsArray.push_back(Thyra::create_VectorSpace(getGhostedMap(i)));
842 
843  ghostedRangeSpace_ = Thyra::productVectorSpace<double>(vsArray);
844  }
845  else {
846  // the range space is not blocked (just an SPMD vector), build it from
847  // the zeroth row
848  ghostedRangeSpace_ = Thyra::create_VectorSpace(getGhostedMap(0));
849  }
850  }
851 
852  return ghostedRangeSpace_;
853 }
854 
855 template <typename Traits,typename LocalOrdinalT>
858 {
860  Thyra::createMember<double>(*getGhostedThyraDomainSpace());
861  Thyra::assign(vec.ptr(),0.0);
862 
863  return vec;
864 }
865 
866 template <typename Traits,typename LocalOrdinalT>
869 {
871  Thyra::createMember<double>(*getGhostedThyraRangeSpace());
872  Thyra::assign(vec.ptr(),0.0);
873 
874  return vec;
875 }
876 
877 template <typename Traits,typename LocalOrdinalT>
880 {
881  // return a flat epetra matrix
882  if(!rowDOFManagerContainer_->containsBlockedDOFManager() &&
883  !colDOFManagerContainer_->containsBlockedDOFManager()) {
884  return Thyra::nonconstEpetraLinearOp(getGhostedEpetraMatrix(0,0));
885  }
886 
887  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockedOp = Thyra::defaultBlockedLinearOp<double>();
888 
889  // get the block dimension
890  std::size_t rBlockDim = getBlockRowCount();
891  std::size_t cBlockDim = getBlockColCount();
892 
893  // this operator will be square
894  blockedOp->beginBlockFill(rBlockDim,cBlockDim);
895 
896  // loop over each block
897  for(std::size_t i=0;i<rBlockDim;i++) {
898  for(std::size_t j=0;j<cBlockDim;j++) {
899  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
900  // build (i,j) block matrix and add it to blocked operator
901  Teuchos::RCP<Thyra::LinearOpBase<double> > block = Thyra::nonconstEpetraLinearOp(getGhostedEpetraMatrix(i,j));
902  blockedOp->setNonconstBlock(i,j,block);
903  }
904  }
905  }
906 
907  // all done
908  blockedOp->endBlockFill();
909 
910  return blockedOp;
911 }
912 
913 template <typename Traits,typename LocalOrdinalT>
916  const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const
917 {
918  using Teuchos::RCP;
919  using Teuchos::rcp_dynamic_cast;
921  using Thyra::get_Epetra_Vector;
922 
923  std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
924 
925  // get product vectors
926  RCP<const ProductVectorBase<double> > prod_in = Thyra::castOrCreateProductVectorBase(in);
927  RCP<ProductVectorBase<double> > prod_out = Thyra::castOrCreateNonconstProductVectorBase(out);
928 
929  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
930  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
931 
932  for(std::size_t i=0;i<blockDim;i++) {
933  // first get each Epetra vector
935  RCP<Epetra_Vector> ep_out;
936 
937  if(not col) {
938  ep_in = get_Epetra_Vector(*getGhostedMap(i),prod_in->getVectorBlock(i));
939  ep_out = get_Epetra_Vector(*getMap(i),prod_out->getNonconstVectorBlock(i));
940  } else {
941  ep_in = get_Epetra_Vector(*getGhostedColMap(i),prod_in->getVectorBlock(i));
942  ep_out = get_Epetra_Vector(*getColMap(i),prod_out->getNonconstVectorBlock(i));
943  }
944 
945  // use Epetra to do global communication
946  ghostToGlobalEpetraVector(i,*ep_in,*ep_out,col);
947  }
948 }
949 
950 template <typename Traits,typename LocalOrdinalT>
953 {
954  using Teuchos::RCP;
955  using Teuchos::rcp_dynamic_cast;
956  using Teuchos::dyn_cast;
957  using Thyra::LinearOpBase;
958  using Thyra::PhysicallyBlockedLinearOpBase;
959  using Thyra::get_Epetra_Operator;
960 
961  // get the block dimension
962  std::size_t rBlockDim = getBlockRowCount();
963  std::size_t cBlockDim = getBlockColCount();
964 
965  // get product vectors
966  const PhysicallyBlockedLinearOpBase<double> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<double> >(in);
967  PhysicallyBlockedLinearOpBase<double> & prod_out = dyn_cast<PhysicallyBlockedLinearOpBase<double> >(out);
968 
969  TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) rBlockDim);
970  TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) cBlockDim);
971  TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) rBlockDim);
972  TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) cBlockDim);
973 
974  for(std::size_t i=0;i<rBlockDim;i++) {
975  for(std::size_t j=0;j<cBlockDim;j++) {
976  if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
977  // extract the blocks
978  RCP<const LinearOpBase<double> > th_in = prod_in.getBlock(i,j);
979  RCP<LinearOpBase<double> > th_out = prod_out.getNonconstBlock(i,j);
980 
981  // sanity check
982  TEUCHOS_ASSERT(th_in!=Teuchos::null);
983  TEUCHOS_ASSERT(th_out!=Teuchos::null);
984 
985  // get the epetra version of the blocks
986  RCP<const Epetra_CrsMatrix> ep_in = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*th_in),true);
987  RCP<Epetra_CrsMatrix> ep_out = rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*th_out),true);
988 
989  // use Epetra to do global communication
990  ghostToGlobalEpetraMatrix(i,*ep_in,*ep_out);
991  }
992  }
993  }
994 }
995 
996 template <typename Traits,typename LocalOrdinalT>
999  const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const
1000 {
1001  using Teuchos::RCP;
1002  using Teuchos::rcp_dynamic_cast;
1004  using Thyra::get_Epetra_Vector;
1005 
1006  std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
1007 
1008  // get product vectors
1009  RCP<const ProductVectorBase<double> > prod_in = Thyra::castOrCreateProductVectorBase(in);
1010  RCP<ProductVectorBase<double> > prod_out = Thyra::castOrCreateNonconstProductVectorBase(out);
1011 
1012  TEUCHOS_ASSERT(prod_in->productSpace()->numBlocks()==(int) blockDim);
1013  TEUCHOS_ASSERT(prod_out->productSpace()->numBlocks()==(int) blockDim);
1014 
1015  for(std::size_t i=0;i<blockDim;i++) {
1016  // first get each Epetra vector
1017  RCP<const Epetra_Vector> ep_in;
1018  RCP<Epetra_Vector> ep_out;
1019 
1020  if(not col) {
1021  ep_in = get_Epetra_Vector(*getMap(i),prod_in->getVectorBlock(i));
1022  ep_out = get_Epetra_Vector(*getGhostedMap(i),prod_out->getNonconstVectorBlock(i));
1023  }
1024  else {
1025  ep_in = get_Epetra_Vector(*getColMap(i),prod_in->getVectorBlock(i));
1026  ep_out = get_Epetra_Vector(*getGhostedColMap(i),prod_out->getNonconstVectorBlock(i));
1027  }
1028 
1029  // use Epetra to do global communication
1030  globalToGhostEpetraVector(i,*ep_in,*ep_out,col);
1031  }
1032 }
1033 
1034 // Epetra methods
1036 
1037 template <typename Traits,typename LocalOrdinalT>
1039 ghostToGlobalEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const
1040 {
1041  using Teuchos::RCP;
1042 
1043  // do the global distribution
1044  RCP<Epetra_Export> exporter = col ? getGhostedColExport(i) : getGhostedExport(i);
1045  out.PutScalar(0.0);
1046  int err = out.Export(in,*exporter,Add);
1047  TEUCHOS_ASSERT_EQUALITY(err,0);
1048 }
1049 
1050 template <typename Traits,typename LocalOrdinalT>
1053 {
1054  using Teuchos::RCP;
1055 
1056  // do the global distribution
1057  RCP<Epetra_Export> exporter = getGhostedExport(blockRow);
1058  out.PutScalar(0.0);
1059  int err = out.Export(in,*exporter,Add);
1060  TEUCHOS_ASSERT_EQUALITY(err,0);
1061 }
1062 
1063 template <typename Traits,typename LocalOrdinalT>
1065 globalToGhostEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const
1066 {
1067  using Teuchos::RCP;
1068 
1069  // do the global distribution
1070  RCP<Epetra_Import> importer = col ? getGhostedColImport(i) : getGhostedImport(i);
1071  out.PutScalar(0.0);
1072  int err = out.Import(in,*importer,Insert);
1073  TEUCHOS_ASSERT_EQUALITY(err,0);
1074 }
1075 
1076 // get the map from the matrix
1077 template <typename Traits,typename LocalOrdinalT>
1079 getMap(int i) const
1080 {
1081  if(maps_[i]==Teuchos::null)
1082  maps_[i] = buildMap(i);
1083 
1084  return maps_[i];
1085 }
1086 
1087 // get the map from the matrix
1088 template <typename Traits,typename LocalOrdinalT>
1090 getColMap(int i) const
1091 {
1092  if(not useColGidProviders_)
1093  return getMap(i);
1094 
1095  if(colMaps_[i]==Teuchos::null)
1096  colMaps_[i] = buildColMap(i);
1097 
1098  return colMaps_[i];
1099 }
1100 
1101 template <typename Traits,typename LocalOrdinalT>
1103 getGhostedMap(int i) const
1104 {
1105  if(ghostedMaps_[i]==Teuchos::null)
1106  ghostedMaps_[i] = buildGhostedMap(i);
1107 
1108  return ghostedMaps_[i];
1109 }
1110 
1111 template <typename Traits,typename LocalOrdinalT>
1113 getGhostedColMap(int i) const
1114 {
1115  if(not useColGidProviders_)
1116  return getGhostedMap(i);
1117 
1118  if(colGhostedMaps_[i]==Teuchos::null)
1119  colGhostedMaps_[i] = buildColGhostedMap(i);
1120 
1121  return colGhostedMaps_[i];
1122 }
1123 
1124 // get the graph of the crs matrix
1125 template <typename Traits,typename LocalOrdinalT>
1127 getGraph(int i,int j) const
1128 {
1129  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> GraphMap;
1130 
1131  GraphMap::const_iterator itr = graphs_.find(std::make_pair(i,j));
1133  if(itr==graphs_.end()) {
1134  graph = buildGraph(i,j);
1135  graphs_[std::make_pair(i,j)] = graph;
1136  }
1137  else
1138  graph = itr->second;
1139 
1140  TEUCHOS_ASSERT(graph!=Teuchos::null);
1141  return graph;
1142 }
1143 
1144 template <typename Traits,typename LocalOrdinalT>
1146 getGhostedGraph(int i,int j) const
1147 {
1148  typedef std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> GraphMap;
1149 
1150  GraphMap::const_iterator itr = ghostedGraphs_.find(std::make_pair(i,j));
1151  Teuchos::RCP<Epetra_CrsGraph> ghostedGraph;
1152  if(itr==ghostedGraphs_.end()) {
1153  ghostedGraph = buildGhostedGraph(i,j,true);
1154  ghostedGraphs_[std::make_pair(i,j)] = ghostedGraph;
1155  }
1156  else
1157  ghostedGraph = itr->second;
1158 
1159  TEUCHOS_ASSERT(ghostedGraph!=Teuchos::null);
1160  return ghostedGraph;
1161 }
1162 
1163 template <typename Traits,typename LocalOrdinalT>
1165 getGhostedImport(int i) const
1166 {
1167  if(importers_[i]==Teuchos::null)
1168  importers_[i] = Teuchos::rcp(new Epetra_Import(*getGhostedMap(i),*getMap(i)));
1169 
1170  return importers_[i];
1171 }
1172 
1173 template <typename Traits,typename LocalOrdinalT>
1176 {
1177  if(not useColGidProviders_)
1178  return getGhostedImport(i);
1179 
1180  if(colImporters_[i]==Teuchos::null)
1181  colImporters_[i] = Teuchos::rcp(new Epetra_Import(*getGhostedColMap(i),*getColMap(i)));
1182 
1183  return colImporters_[i];
1184 }
1185 
1186 template <typename Traits,typename LocalOrdinalT>
1188 getGhostedExport(int i) const
1189 {
1190  if(exporters_[i]==Teuchos::null)
1191  exporters_[i] = Teuchos::rcp(new Epetra_Export(*getGhostedMap(i),*getMap(i)));
1192 
1193  return exporters_[i];
1194 }
1195 
1196 template <typename Traits,typename LocalOrdinalT>
1199 {
1200  if(not useColGidProviders_)
1201  return getGhostedExport(i);
1202 
1203  if(colExporters_[i]==Teuchos::null)
1204  colExporters_[i] = Teuchos::rcp(new Epetra_Export(*getGhostedColMap(i),*getColMap(i)));
1205 
1206  return colExporters_[i];
1207 }
1208 
1209 template <typename Traits,typename LocalOrdinalT>
1211 buildMap(int i) const
1212 {
1213  std::vector<int> indices;
1214 
1215  // get the global indices
1216  getGlobalIndexer(i)->getOwnedIndices(indices);
1217 
1218  return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1219 }
1220 
1221 template <typename Traits,typename LocalOrdinalT>
1223 buildColMap(int i) const
1224 {
1225  if(not useColGidProviders_)
1226  return buildMap(i);
1227 
1228  std::vector<int> indices;
1229 
1230  // get the global indices
1231  getColGlobalIndexer(i)->getOwnedIndices(indices);
1232 
1233  return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1234 }
1235 
1236 // build the ghosted map
1237 template <typename Traits,typename LocalOrdinalT>
1239 buildGhostedMap(int i) const
1240 {
1241  std::vector<int> indices;
1242 
1243  // get the global indices
1244  getGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
1245 
1246  return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1247 }
1248 
1249 // build the ghosted map
1250 template <typename Traits,typename LocalOrdinalT>
1253 {
1254  if(not useColGidProviders_)
1255  return buildGhostedMap(i);
1256 
1257  std::vector<int> indices;
1258 
1259  // get the global indices
1260  getColGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
1261 
1262  return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*eComm_));
1263 }
1264 
1265 // get the graph of the crs matrix
1266 template <typename Traits,typename LocalOrdinalT>
1268 buildGraph(int i,int j) const
1269 {
1270  using Teuchos::RCP;
1271  using Teuchos::rcp;
1272 
1273  // build the map and allocate the space for the graph and
1274  // grab the ghosted graph
1275  RCP<Epetra_Map> map_i = getMap(i);
1276  RCP<Epetra_Map> map_j = getColMap(j);
1277 
1278  TEUCHOS_ASSERT(map_i!=Teuchos::null);
1279  TEUCHOS_ASSERT(map_j!=Teuchos::null);
1280 
1281  RCP<Epetra_CrsGraph> graph = rcp(new Epetra_CrsGraph(Copy,*map_i,0));
1282  RCP<Epetra_CrsGraph> oGraph = buildFilteredGhostedGraph(i,j);
1283  // this is the only place buildFilteredGhostedGraph is called. That is because
1284  // only the unghosted graph should reflect any of the filtering.
1285 
1286  // perform the communication to finish building graph
1287  RCP<Epetra_Export> exporter = getGhostedExport(i);
1288  int err = graph->Export( *oGraph, *exporter, Insert );
1289  TEUCHOS_ASSERT_EQUALITY(err,0);
1290  graph->FillComplete(*map_j,*map_i);
1291  graph->OptimizeStorage();
1292 
1293  return graph;
1294 }
1295 
1296 template <typename Traits,typename LocalOrdinalT>
1298 buildGhostedGraph(int i,int j,bool optimizeStorage) const
1299 {
1300  // build the map and allocate the space for the graph
1301  Teuchos::RCP<Epetra_Map> rowMap = getGhostedMap(i);
1302  Teuchos::RCP<Epetra_Map> colMap = getGhostedColMap(j);
1303  Teuchos::RCP<Epetra_CrsGraph> graph = Teuchos::rcp(new Epetra_CrsGraph(Copy,*rowMap,*colMap,0));
1304 
1305  std::vector<std::string> elementBlockIds;
1306 
1308 
1309  rowProvider = getGlobalIndexer(i);
1310  colProvider = getColGlobalIndexer(j);
1311 
1312  rowProvider->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
1313  // same element blocks
1314 
1315  const Teuchos::RCP<const ConnManagerBase<LocalOrdinalT> > conn_mgr = colProvider->getConnManagerBase();
1316  const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors();
1317 
1318  // graph information about the mesh
1319  std::vector<std::string>::const_iterator blockItr;
1320  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
1321  std::string blockId = *blockItr;
1322 
1323  // grab elements for this block
1324  const std::vector<LocalOrdinalT> & elements = rowProvider->getElementBlock(blockId); // each sub provider "should" have the
1325  // same elements in each element block
1326 
1327  // get information about number of indicies
1328  std::vector<int> row_gids;
1329  std::vector<int> col_gids;
1330 
1331  // loop over the elemnts
1332  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
1333 
1334  rowProvider->getElementGIDs(elements[elmt],row_gids);
1335  colProvider->getElementGIDs(elements[elmt],col_gids);
1336 
1337  if (han) {
1338  const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[elmt]);
1339  for (typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
1340  eit != aes.end(); ++eit) {
1341  std::vector<int> other_col_gids;
1342  colProvider->getElementGIDs(*eit, other_col_gids);
1343  col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
1344  }
1345  }
1346 
1347  for(std::size_t row=0;row<row_gids.size();row++)
1348  graph->InsertGlobalIndices(row_gids[row],col_gids.size(),&col_gids[0]);
1349  }
1350  }
1351 
1352  // finish filling the graph: Make sure the colmap and row maps coincide to
1353  // minimize calls to LID lookups
1354  graph->FillComplete(*colMap,*rowMap);
1355  if(optimizeStorage)
1356  graph->OptimizeStorage();
1357 
1358  return graph;
1359 }
1360 
1361 template <typename Traits,typename LocalOrdinalT>
1363 buildFilteredGhostedGraph(int i,int j) const
1364 {
1365  using Teuchos::RCP;
1366  using Teuchos::rcp;
1367  using Teuchos::rcp_dynamic_cast;
1368 
1369  // figure out if the domain is filtered
1371  = rcp_dynamic_cast<const Filtered_UniqueGlobalIndexer<LocalOrdinalT,int> >(getColGlobalIndexer(j));
1372 
1373  // domain is unfiltered, a filtered graph is just the original ghosted graph
1374  if(filtered_ugi==Teuchos::null)
1375  return buildGhostedGraph(i,j,true);
1376 
1377  // get all local indices that are active (i.e. unfiltered)
1378  std::vector<int> ghostedActive;
1379  filtered_ugi->getOwnedAndGhostedNotFilteredIndicator(ghostedActive);
1380 
1381  // This will build a new ghosted graph without optimized storage so entries can be removed.
1382  Teuchos::RCP<Epetra_CrsGraph> filteredGraph = buildGhostedGraph(i,j,false);
1383  // false implies that storage is not optimzied
1384 
1385  // remove filtered column entries
1386  for(int i=0;i<filteredGraph->NumMyRows();i++) {
1387  std::vector<int> removedIndices;
1388  int numIndices = 0;
1389  int * indices = 0;
1390  TEUCHOS_ASSERT(filteredGraph->ExtractMyRowView(i,numIndices,indices)==0);
1391 
1392  for(int j=0;j<numIndices;j++) {
1393  if(ghostedActive[indices[j]]==0)
1394  removedIndices.push_back(indices[j]);
1395  }
1396 
1397  TEUCHOS_ASSERT(filteredGraph->RemoveMyIndices(i,Teuchos::as<int>(removedIndices.size()),&removedIndices[0])==0);
1398  }
1399 
1400  // finish filling the graph
1401  Teuchos::RCP<Epetra_Map> rowMap = getGhostedMap(i);
1402  Teuchos::RCP<Epetra_Map> colMap = getGhostedColMap(j);
1403 
1404  TEUCHOS_ASSERT(filteredGraph->FillComplete(*colMap,*rowMap)==0);
1405  TEUCHOS_ASSERT(filteredGraph->OptimizeStorage()==0);
1406 
1407  return filteredGraph;
1408 }
1409 
1410 template <typename Traits,typename LocalOrdinalT>
1412 getEpetraMatrix(int i,int j) const
1413 {
1414  Teuchos::RCP<Epetra_CrsGraph> eGraph = getGraph(i,j);
1416  TEUCHOS_ASSERT(mat->Filled());
1417  return mat;
1418 }
1419 
1420 template <typename Traits,typename LocalOrdinalT>
1422 getGhostedEpetraMatrix(int i,int j) const
1423 {
1424  Teuchos::RCP<Epetra_CrsGraph> eGraph = getGhostedGraph(i,j);
1426  TEUCHOS_ASSERT(mat->Filled());
1427  return mat;
1428 }
1429 
1430 template <typename Traits,typename LocalOrdinalT>
1433 {
1434  return eComm_;
1435 }
1436 
1437 template <typename Traits,typename LocalOrdinalT>
1440 {
1441  return rowDOFManagerContainer_->getFieldBlocks();
1442 }
1443 
1444 template <typename Traits,typename LocalOrdinalT>
1447 {
1448  return colDOFManagerContainer_->getFieldBlocks();
1449 }
1450 
1451 }
1452 
1453 #endif
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraRangeVector() const
Get a range vector.
void ghostToGlobalEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
bool is_null(const boost::shared_ptr< T > &p)
virtual const Teuchos::RCP< Epetra_Map > buildMap(int i) const
Teuchos::RCP< Epetra_CrsMatrix > getGhostedEpetraMatrix(int i, int j) const
virtual const Teuchos::RCP< Epetra_Export > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< Thyra::LinearOpBase< double > > getThyraMatrix() const
Get a Thyra operator.
int MyLength() const
const Teuchos::RCP< Epetra_CrsMatrix > get_A() const
size_type size() const
virtual void readVector(const std::string &identifier, LinearObjContainer &loc, int id) const
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
Teuchos::RCP< Thyra::VectorBase< double > > getThyraDomainVector() const
Get a domain vector.
bool is_null(const std::shared_ptr< T > &p)
virtual const Teuchos::RCP< Epetra_Map > buildColMap(int i) const
void globalToGhostEpetraVector(int i, const Epetra_Vector &in, Epetra_Vector &out, bool col) const
void globalToGhostThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
Teuchos::RCP< const UniqueGlobalIndexer< LocalOrdinalT, int > > getColGlobalIndexer(int i) const
int PutScalar(double ScalarConstant)
virtual void set_x_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
T_To & dyn_cast(T_From &from)
virtual const Teuchos::RCP< Epetra_CrsGraph > buildGhostedGraph(int i, int j, bool optimizeStorage) const
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildDomainContainer() const
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraRangeSpace() const
Get the range vector space (f)
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraRangeSpace() const
Get the range vector space (f)
Insert
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 const Teuchos::RCP< Epetra_Map > getMap(int i) const
get the map from the matrix
Teuchos::RCP< Thyra::VectorBase< double > > getThyraRangeVector() const
Get a range vector.
void ghostToGlobalEpetraMatrix(int blockRow, const Epetra_CrsMatrix &in, Epetra_CrsMatrix &out) const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual const Teuchos::RCP< Epetra_Import > getGhostedImport(int i) const
get importer for converting an overalapped object to a "normal" object
PHX::MDField< ScalarT > vector
void ghostToGlobalThyraVector(const Teuchos::RCP< const Thyra::VectorBase< double > > &in, const Teuchos::RCP< Thyra::VectorBase< double > > &out, bool col) const
bool is_null() const
void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase< double > &in, Thyra::LinearOpBase< double > &out) const
int PutScalar(double ScalarConstant)
int NumMyRows() const
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 > buildFilteredGhostedGraph(int i, int j) const
virtual const Teuchos::RCP< Epetra_Map > getColMap(int i) const
get the map from the matrix
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
Teuchos::RCP< const DOFManagerContainer > colDOFManagerContainer_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const Teuchos::OpaqueWrapper< MPI_Comm > > rawMpiComm_
void initialize(const Teuchos::RCP< const Thyra::VectorSpaceBase< double > > &ghostedSpace, const Teuchos::RCP< const Thyra::VectorSpaceBase< double > > &ownedSpace, const std::vector< Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > > &gedBlocks)
void addExcludedPairs(const std::vector< std::pair< int, int > > &exPairs)
exclude a vector of pairs from the matrix
const Teuchos::RCP< Epetra_Vector > get_x() const
virtual void set_f_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
Teuchos::RCP< Thyra::VectorBase< double > > getGhostedThyraDomainVector() const
Get a domain vector.
std::vector< ScalarT > values
virtual const Teuchos::RCP< const Epetra_Comm > getEpetraComm() const
get exporter for converting an overalapped object to a "normal" object
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...
void addExcludedPair(int rowBlock, int colBlock)
exclude a block pair from the matrix
void initializeGhostedContainer(int, LinearObjContainer &loc) const
virtual const Teuchos::RCP< Epetra_Map > getGhostedMap(int i) const
get the ghosted map from the matrix
Teuchos::RCP< Epetra_CrsMatrix > getEpetraMatrix(int i, int j) const
bool Filled() const
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
virtual const Teuchos::RCP< Epetra_CrsGraph > buildGraph(int i, int j) const
void setMapsForBlocks(const std::vector< Teuchos::RCP< const Epetra_Map > > &blockMaps)
void initializeContainer(int, LinearObjContainer &loc) const
virtual const Teuchos::RCP< Epetra_Map > buildColGhostedMap(int i) const
virtual const Teuchos::RCP< Epetra_Map > buildGhostedMap(int i) const
int RemoveMyIndices(int LocalRow, int NumIndices, int *Indices)
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void initializeGhostedContainer_internal(int mem, ThyraObjContainer< double > &loc) const
virtual void set_A_th(const Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > &in)=0
Teuchos::RCP< const Teuchos::Comm< int > > comm
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 > getGraph(int i, int j) const
get the graph of the crs matrix
Copy
virtual void writeVector(const std::string &identifier, const LinearObjContainer &loc, int id) const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Teuchos::RCP< Thyra::LinearOpBase< double > > getGhostedThyraMatrix() const
Get a Thyra operator.
virtual const Teuchos::RCP< Epetra_Map > getGhostedColMap(int i) const
get the ghosted map from the matrix
Teuchos::RCP< const panzer::BlockedDOFManager< int, int > > getGlobalIndexer() const
void buildGatherScatterEvaluators(const BuilderT &builder)
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
void initializeContainer_internal(int mem, ThyraObjContainer< double > &loc) const
RCP< const T > getConst() const
const Teuchos::RCP< Epetra_Vector > get_f() const
Add
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Ptr< T > ptr() const
Teuchos::RCP< const DOFManagerContainer > rowDOFManagerContainer_
virtual const Teuchos::RCP< Epetra_Export > getGhostedColExport(int j) const
get exporter for converting an overalapped object to a "normal" object
const Teuchos::RCP< Epetra_Vector > get_dxdt() const
BlockedEpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const UniqueGlobalIndexerBase > &gidProvider, bool useDiscreteAdjoint=false)
Teuchos::RCP< const Thyra::VectorSpaceBase< double > > getGhostedThyraDomainSpace() const
Get the domain vector space (x and dxdt)
void initialize(const Teuchos::RCP< const Epetra_Import > &importer, const Teuchos::RCP< const Epetra_Map > &ghostedMap, const Teuchos::RCP< const Epetra_Map > &ownedMap)