Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_EpetraLinearObjFactory_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_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
44 #define PANZER_EPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
45 
46 #include "Panzer_UniqueGlobalIndexer.hpp"
47 #include "Panzer_Filtered_UniqueGlobalIndexer.hpp"
48 #include "Panzer_ConnManager.hpp"
49 #include "Panzer_ThyraObjContainer.hpp"
50 #include "Panzer_EpetraVector_ReadOnly_GlobalEvaluationData.hpp"
51 
52 #include "Thyra_SpmdVectorBase.hpp"
53 
54 #include "Epetra_MultiVector.h"
55 #include "Epetra_Vector.h"
56 #include "Epetra_CrsMatrix.h"
57 #include "Epetra_MpiComm.h"
58 
59 #include "EpetraExt_VectorOut.h"
60 #include "EpetraExt_VectorIn.h"
61 
62 namespace panzer {
63 
64 using Teuchos::RCP;
65 
66 // ************************************************************
67 // class EpetraLinearObjFactory
68 // ************************************************************
69 
70 template <typename Traits,typename LocalOrdinalT>
71 EpetraLinearObjFactory<Traits,LocalOrdinalT>::EpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
72  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,int> > & gidProvider,
73  bool useDiscreteAdjoint)
74  : comm_(Teuchos::null), gidProvider_(gidProvider), rawMpiComm_(comm->getRawMpiComm()), useDiscreteAdjoint_(useDiscreteAdjoint)
75 {
76  comm_ = Teuchos::rcp(new Epetra_MpiComm(*rawMpiComm_));
77  hasColProvider_ = colGidProvider_!=Teuchos::null;
78 
79  // build and register the gather/scatter evaluators with
80  // the base class.
81  this->buildGatherScatterEvaluators(*this);
82 }
83 
84 template <typename Traits,typename LocalOrdinalT>
85 EpetraLinearObjFactory<Traits,LocalOrdinalT>::EpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
86  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,int> > & gidProvider,
87  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,int> > & colGidProvider,
88  bool useDiscreteAdjoint)
89  : comm_(Teuchos::null), gidProvider_(gidProvider), colGidProvider_(colGidProvider), rawMpiComm_(comm->getRawMpiComm()), useDiscreteAdjoint_(useDiscreteAdjoint)
90 {
91  comm_ = Teuchos::rcp(new Epetra_MpiComm(*rawMpiComm_));
92  hasColProvider_ = colGidProvider_!=Teuchos::null;
93 
94  // build and register the gather/scatter evaluators with
95  // the base class.
96  this->buildGatherScatterEvaluators(*this);
97 }
98 
99 template <typename Traits,typename LocalOrdinalT>
100 EpetraLinearObjFactory<Traits,LocalOrdinalT>::EpetraLinearObjFactory(const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,int> > & gidProvider)
101  : comm_(Teuchos::null), gidProvider_(gidProvider), useDiscreteAdjoint_(false)
102 {
103  hasColProvider_ = colGidProvider_!=Teuchos::null;
104 
105  // build and register the gather/scatter evaluators with
106  // the base class.
107  this->buildGatherScatterEvaluators(*this);
108 }
109 
110 template <typename Traits,typename LocalOrdinalT>
111 EpetraLinearObjFactory<Traits,LocalOrdinalT>::~EpetraLinearObjFactory()
112 { }
113 
114 // LinearObjectFactory functions
116 
117 template <typename Traits,typename LocalOrdinalT>
118 void
119 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
120 readVector(const std::string & identifier,LinearObjContainer & loc,int id) const
121 {
122  using Teuchos::RCP;
123  using Teuchos::rcp_dynamic_cast;
124  using Teuchos::dyn_cast;
125 
126  EpetraLinearObjContainer & eloc = dyn_cast<EpetraLinearObjContainer>(loc);
127 
128  // extract the vector from linear object container
129  RCP<Thyra::VectorBase<double> > vec;
130  switch(id) {
131  case LinearObjContainer::X:
132  vec = eloc.get_x_th();
133  break;
134  case LinearObjContainer::DxDt:
135  vec = eloc.get_dxdt_th();
136  break;
137  case LinearObjContainer::F:
138  vec = eloc.get_f_th();
139  break;
140  default:
141  TEUCHOS_ASSERT(false);
142  break;
143  };
144 
145  // convert to Epetra then read in from a file
146  RCP<Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(),vec);
147 
148  // build the file name from the identifier
149  std::stringstream ss;
150  ss << identifier << ".mm";
151 
152  // read in vector (wow the MM to Vector is a poorly designed interface!)
153  Epetra_Vector * ptr_ex = 0;
154  TEUCHOS_ASSERT(0==EpetraExt::MatrixMarketFileToVector(ss.str().c_str(),*getMap(),ptr_ex));
155 
156  *ex = *ptr_ex;
157  delete ptr_ex;
158 }
159 
160 template <typename Traits,typename LocalOrdinalT>
161 void
162 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
163 writeVector(const std::string & identifier,const LinearObjContainer & loc,int id) const
164 {
165  using Teuchos::RCP;
166  using Teuchos::rcp_dynamic_cast;
167  using Teuchos::dyn_cast;
168 
169  const EpetraLinearObjContainer & eloc = dyn_cast<const EpetraLinearObjContainer>(loc);
170 
171  // extract the vector from linear object container
172  RCP<const Thyra::VectorBase<double> > vec;
173  switch(id) {
174  case LinearObjContainer::X:
175  vec = eloc.get_x_th();
176  break;
177  case LinearObjContainer::DxDt:
178  vec = eloc.get_dxdt_th();
179  break;
180  case LinearObjContainer::F:
181  vec = eloc.get_f_th();
182  break;
183  default:
184  TEUCHOS_ASSERT(false);
185  break;
186  };
187 
188  // convert to Epetra then write out to a file
189  RCP<const Epetra_Vector> ex = Thyra::get_Epetra_Vector(*getMap(),vec);
190 
191  // build the file name from the identifier
192  std::stringstream ss;
193  ss << identifier << ".mm";
194 
195  // write out vector
196  TEUCHOS_ASSERT(0==EpetraExt::VectorToMatrixMarketFile(ss.str().c_str(),*ex));
197 }
198 
199 template <typename Traits,typename LocalOrdinalT>
200 Teuchos::RCP<LinearObjContainer> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildLinearObjContainer() const
201 {
202  Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(new EpetraLinearObjContainer(getColMap(),getMap()));
203 
204  return container;
205 }
206 
207 template <typename Traits,typename LocalOrdinalT>
208 Teuchos::RCP<LinearObjContainer> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedLinearObjContainer() const
209 {
210  Teuchos::RCP<EpetraLinearObjContainer> container = Teuchos::rcp(new EpetraLinearObjContainer(getGhostedColMap(),getGhostedMap()));
211 
212  return container;
213 }
214 
215 template <typename Traits,typename LocalOrdinalT>
216 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::globalToGhostContainer(const LinearObjContainer & in,
217  LinearObjContainer & out,int mem) const
218 {
219  using Teuchos::is_null;
220 
221  typedef LinearObjContainer LOC;
222  const EpetraLinearObjContainer & e_in = Teuchos::dyn_cast<const EpetraLinearObjContainer>(in);
223  EpetraLinearObjContainer & e_out = Teuchos::dyn_cast<EpetraLinearObjContainer>(out);
224 
225  // Operations occur if the GLOBAL container has the correct targets!
226  // Users set the GLOBAL continer arguments
227  if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
228  globalToGhostEpetraVector(*e_in.get_x(),*e_out.get_x(),true);
229 
230  if ( !is_null(e_in.get_dxdt()) && !is_null(e_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
231  globalToGhostEpetraVector(*e_in.get_dxdt(),*e_out.get_dxdt(),true);
232 
233  if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
234  globalToGhostEpetraVector(*e_in.get_f(),*e_out.get_f(),false);
235 }
236 
237 template <typename Traits,typename LocalOrdinalT>
238 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::ghostToGlobalContainer(const LinearObjContainer & in,
239  LinearObjContainer & out,int mem) const
240 {
241  using Teuchos::is_null;
242 
243  typedef LinearObjContainer LOC;
244  const EpetraLinearObjContainer & e_in = Teuchos::dyn_cast<const EpetraLinearObjContainer>(in);
245  EpetraLinearObjContainer & e_out = Teuchos::dyn_cast<EpetraLinearObjContainer>(out);
246 
247  // Operations occur if the GLOBAL container has the correct targets!
248  // Users set the GLOBAL continer arguments
249  if ( !is_null(e_in.get_x()) && !is_null(e_out.get_x()) && ((mem & LOC::X)==LOC::X))
250  ghostToGlobalEpetraVector(*e_in.get_x(),*e_out.get_x(),true);
251 
252  if ( !is_null(e_in.get_f()) && !is_null(e_out.get_f()) && ((mem & LOC::F)==LOC::F))
253  ghostToGlobalEpetraVector(*e_in.get_f(),*e_out.get_f(),false);
254 
255  if ( !is_null(e_in.get_A()) && !is_null(e_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
256  ghostToGlobalEpetraMatrix(*e_in.get_A(),*e_out.get_A());
257 }
258 
259 template <typename Traits,typename LocalOrdinalT>
260 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::ghostToGlobalEpetraVector(const Epetra_Vector & in,Epetra_Vector & out,bool col) const
261 {
262  using Teuchos::RCP;
263 
264  // do the global distribution
265  RCP<Epetra_Export> exporter = col ? getGhostedColExport() : getGhostedExport();
266  TEUCHOS_ASSERT(out.PutScalar(0.0)==0);
267 
268  int errCode = out.Export(in,*exporter,Add);
269  TEUCHOS_TEST_FOR_EXCEPTION(errCode!=0,std::runtime_error,
270  "Epetra_Vector::Export returned an error code of " << errCode << "!");
271 }
272 
273 template <typename Traits,typename LocalOrdinalT>
274 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::ghostToGlobalEpetraMatrix(const Epetra_CrsMatrix & in,Epetra_CrsMatrix & out) const
275 {
276  using Teuchos::RCP;
277 
278  // do the global distribution
279  RCP<Epetra_Export> exporter = getGhostedExport();
280  out.PutScalar(0.0);
281  out.Export(in,*exporter,Add);
282 }
283 
284 template <typename Traits,typename LocalOrdinalT>
285 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::globalToGhostEpetraVector(const Epetra_Vector & in,Epetra_Vector & out,bool col) const
286 {
287  using Teuchos::RCP;
288 
289  // do the global distribution
290  RCP<Epetra_Import> importer = col ? getGhostedColImport() : getGhostedImport();
291  out.PutScalar(0.0);
292  out.Import(in,*importer,Insert);
293  // NOTE: These lines cause the response_residual test to fail!
294  // int retval = out.Import(in,*importer,Insert);
295  // TEUCHOS_TEST_FOR_EXCEPTION(0!=retval,std::logic_error,"panzer::EpetraLOF::globalToGhostEpetraVector "
296  // "import call failed with return value retval = " << retval);
297 }
298 
299 template <typename Traits,typename LocalOrdinalT>
300 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::
301 adjustForDirichletConditions(const LinearObjContainer & localBCRows,
302  const LinearObjContainer & globalBCRows,
303  LinearObjContainer & ghostedObjs,
304  bool zeroVectorRows, bool adjustX) const
305 
306 {
307  const EpetraLinearObjContainer & e_localBCRows = Teuchos::dyn_cast<const EpetraLinearObjContainer>(localBCRows);
308  const EpetraLinearObjContainer & e_globalBCRows = Teuchos::dyn_cast<const EpetraLinearObjContainer>(globalBCRows);
309  EpetraLinearObjContainer & e_ghosted = Teuchos::dyn_cast<EpetraLinearObjContainer>(ghostedObjs);
310 
311  TEUCHOS_ASSERT(!Teuchos::is_null(e_localBCRows.get_f()));
312  TEUCHOS_ASSERT(!Teuchos::is_null(e_globalBCRows.get_f()));
313 
314  // pull out jacobian and vector
315  Teuchos::RCP<Epetra_CrsMatrix> A = e_ghosted.get_A();
316  Teuchos::RCP<Epetra_Vector> f = e_ghosted.get_f();
317  if(adjustX) f = e_ghosted.get_x();
318 
319  const Epetra_Vector & local_bcs = *(e_localBCRows.get_f());
320  const Epetra_Vector & global_bcs = *(e_globalBCRows.get_f());
321 
322  TEUCHOS_ASSERT(local_bcs.MyLength()==global_bcs.MyLength());
323  for(int i=0;i<local_bcs.MyLength();i++) {
324  if(global_bcs[i]==0.0)
325  continue;
326 
327  int numEntries = 0;
328  double * values = 0;
329  int * indices = 0;
330 
331  if(local_bcs[i]==0.0 || zeroVectorRows) {
332  // this boundary condition was NOT set by this processor
333  // or the user requrested that every row be zeroed
334 
335  // if they exist put 0.0 in each entry
336  if(!Teuchos::is_null(f))
337  (*f)[i] = 0.0;
338  if(!Teuchos::is_null(A)) {
339  A->ExtractMyRowView(i,numEntries,values,indices);
340  for(int c=0;c<numEntries;c++)
341  values[c] = 0.0;
342  }
343  }
344  else {
345  // this boundary condition was set by this processor
346 
347  double scaleFactor = global_bcs[i];
348 
349  // if they exist scale linear objects by scale factor
350  if(!Teuchos::is_null(f))
351  (*f)[i] /= scaleFactor;
352  if(!Teuchos::is_null(A)) {
353  A->ExtractMyRowView(i,numEntries,values,indices);
354  for(int c=0;c<numEntries;c++)
355  values[c] /= scaleFactor;
356  }
357  }
358  }
359 }
360 
361 template <typename Traits,typename LocalOrdinalT>
362 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::
363 applyDirichletBCs(const LinearObjContainer & counter,
364  LinearObjContainer & result) const
365 {
366  using Teuchos::RCP;
367  using Teuchos::rcp_dynamic_cast;
368  using Teuchos::dyn_cast;
369 
370  const ThyraObjContainer<double> & th_counter = dyn_cast<const ThyraObjContainer<double> >(counter);
371  ThyraObjContainer<double> & th_result = dyn_cast<ThyraObjContainer<double> >(result);
372 
373  RCP<const Thyra::VectorBase<double> > count = th_counter.get_f_th();
374  RCP<const Thyra::VectorBase<double> > f_in = th_counter.get_f_th();
375  RCP<Thyra::VectorBase<double> > f_out = th_result.get_f_th();
376 
377  Teuchos::ArrayRCP<const double> count_array,f_in_array;
378  Teuchos::ArrayRCP<double> f_out_array;
379 
380  rcp_dynamic_cast<const Thyra::SpmdVectorBase<double> >(count,true)->getLocalData(Teuchos::ptrFromRef(count_array));
381  rcp_dynamic_cast<const Thyra::SpmdVectorBase<double> >(f_in,true)->getLocalData(Teuchos::ptrFromRef(f_in_array));
382  rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(f_out,true)->getNonconstLocalData(Teuchos::ptrFromRef(f_out_array));
383 
384  TEUCHOS_ASSERT(count_array.size()==f_in_array.size());
385  TEUCHOS_ASSERT(count_array.size()==f_out_array.size());
386 
387  for(Teuchos::ArrayRCP<double>::size_type i=0;i<count_array.size();i++) {
388  if(count_array[i]!=0.0)
389  f_out_array[i] = f_in_array[i];
390  }
391 }
392 
393 template <typename Traits,typename LocalOrdinalT>
394 Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
395 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
396 buildDomainContainer() const
397 {
398  Teuchos::RCP<EpetraVector_ReadOnly_GlobalEvaluationData> vec_ged
399  = Teuchos::rcp(new EpetraVector_ReadOnly_GlobalEvaluationData);
400  vec_ged->initialize(getGhostedImport(),getGhostedColMap(),getColMap());
401 
402  return vec_ged;
403 }
404 
405 template <typename Traits,typename LocalOrdinalT>
406 Teuchos::MpiComm<int> EpetraLinearObjFactory<Traits,LocalOrdinalT>::
407 getComm() const
408 {
409  return Teuchos::MpiComm<int>(Teuchos::opaqueWrapper(dynamic_cast<const Epetra_MpiComm &>(*getEpetraComm()).Comm()));
410 }
411 
412 template <typename Traits,typename LocalOrdinalT>
413 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
414 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
415 getThyraDomainSpace() const
416 {
417  if(domainSpace_ == Teuchos::null) {
418  // in the first case, range is domain,
419  // in the second domain must be constructed
420  if(!hasColProvider_)
421  domainSpace_ = getThyraRangeSpace();
422  else
423  domainSpace_ = Thyra::create_VectorSpace(getColMap());
424  }
425 
426  return domainSpace_;
427 }
428 
429 template <typename Traits,typename LocalOrdinalT>
430 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
431 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
432 getThyraRangeSpace() const
433 {
434  if(rangeSpace_ == Teuchos::null)
435  rangeSpace_ = Thyra::create_VectorSpace(getMap());
436 
437  return rangeSpace_;
438 }
439 
440 template <typename Traits,typename LocalOrdinalT>
441 Teuchos::RCP<Thyra::LinearOpBase<double> >
442 EpetraLinearObjFactory<Traits,LocalOrdinalT>::
443 getThyraMatrix() const
444 {
445  return Thyra::nonconstEpetraLinearOp(getEpetraMatrix());
446 }
447 
448 // Functions for initalizing a container
450 
451 template <typename Traits,typename LocalOrdinalT>
452 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeContainer(int mem,LinearObjContainer & loc) const
453 {
454  EpetraLinearObjContainer & eloc = Teuchos::dyn_cast<EpetraLinearObjContainer>(loc);
455  initializeContainer(mem,eloc);
456 }
457 
458 template <typename Traits,typename LocalOrdinalT>
459 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeContainer(int mem,EpetraLinearObjContainer & loc) const
460 {
461  typedef EpetraLinearObjContainer ELOC;
462 
463  loc.clear();
464 
465  if((mem & ELOC::X) == ELOC::X)
466  loc.set_x(getEpetraColVector());
467 
468  if((mem & ELOC::DxDt) == ELOC::DxDt)
469  loc.set_dxdt(getEpetraColVector());
470 
471  if((mem & ELOC::F) == ELOC::F)
472  loc.set_f(getEpetraVector());
473 
474  if((mem & ELOC::Mat) == ELOC::Mat)
475  loc.set_A(getEpetraMatrix());
476 }
477 
478 template <typename Traits,typename LocalOrdinalT>
479 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeGhostedContainer(int mem,LinearObjContainer & loc) const
480 {
481  EpetraLinearObjContainer & eloc = Teuchos::dyn_cast<EpetraLinearObjContainer>(loc);
482  initializeGhostedContainer(mem,eloc);
483 }
484 
485 template <typename Traits,typename LocalOrdinalT>
486 void EpetraLinearObjFactory<Traits,LocalOrdinalT>::initializeGhostedContainer(int mem,EpetraLinearObjContainer & loc) const
487 {
488  typedef EpetraLinearObjContainer ELOC;
489 
490  loc.clear();
491 
492  if((mem & ELOC::X) == ELOC::X)
493  loc.set_x(getGhostedEpetraColVector());
494 
495  if((mem & ELOC::DxDt) == ELOC::DxDt)
496  loc.set_dxdt(getGhostedEpetraColVector());
497 
498  if((mem & ELOC::F) == ELOC::F) {
499  loc.set_f(getGhostedEpetraVector());
500  loc.setRequiresDirichletAdjustment(true);
501  }
502 
503  if((mem & ELOC::Mat) == ELOC::Mat) {
504  loc.set_A(getGhostedEpetraMatrix());
505  loc.setRequiresDirichletAdjustment(true);
506  }
507 }
508 
509 // "Get" functions
511 
512 // get the map from the matrix
513 template <typename Traits,typename LocalOrdinalT>
514 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getMap() const
515 {
516  if(map_==Teuchos::null) map_ = buildMap();
517 
518  return map_;
519 }
520 
521 // get the map from the matrix
522 template <typename Traits,typename LocalOrdinalT>
523 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getColMap() const
524 {
525  if(cMap_==Teuchos::null) cMap_ = buildColMap();
526 
527  return cMap_;
528 }
529 
530 template <typename Traits,typename LocalOrdinalT>
531 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedMap() const
532 {
533  if(ghostedMap_==Teuchos::null) ghostedMap_ = buildGhostedMap();
534 
535  return ghostedMap_;
536 }
537 
538 template <typename Traits,typename LocalOrdinalT>
539 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedColMap() const
540 {
541  if(cGhostedMap_==Teuchos::null) cGhostedMap_ = buildGhostedColMap();
542 
543  return cGhostedMap_;
544 }
545 
546 // get the graph of the crs matrix
547 template <typename Traits,typename LocalOrdinalT>
548 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGraph() const
549 {
550  if(graph_==Teuchos::null) graph_ = buildGraph();
551 
552  return graph_;
553 }
554 
555 template <typename Traits,typename LocalOrdinalT>
556 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedGraph() const
557 {
558  if(ghostedGraph_==Teuchos::null) ghostedGraph_ = buildGhostedGraph(true);
559 
560  return ghostedGraph_;
561 }
562 
563 template <typename Traits,typename LocalOrdinalT>
564 const Teuchos::RCP<Epetra_Import> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedImport() const
565 {
566  if(importer_==Teuchos::null)
567  importer_ = Teuchos::rcp(new Epetra_Import(*getGhostedMap(),*getMap()));
568 
569  return importer_;
570 }
571 
572 template <typename Traits,typename LocalOrdinalT>
573 const Teuchos::RCP<Epetra_Import> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedColImport() const
574 {
575  if(!hasColProvider_)
576  colImporter_ = getGhostedImport(); // they are the same in this case
577 
578  if(colImporter_==Teuchos::null)
579  colImporter_ = Teuchos::rcp(new Epetra_Import(*getGhostedColMap(),*getColMap()));
580 
581  return colImporter_;
582 }
583 
584 template <typename Traits,typename LocalOrdinalT>
585 const Teuchos::RCP<Epetra_Export> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedExport() const
586 {
587  if(exporter_==Teuchos::null)
588  exporter_ = Teuchos::rcp(new Epetra_Export(*getGhostedMap(),*getMap()));
589 
590  return exporter_;
591 }
592 
593 template <typename Traits,typename LocalOrdinalT>
594 const Teuchos::RCP<Epetra_Export> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedColExport() const
595 {
596  if(!hasColProvider_)
597  colExporter_ = getGhostedExport(); // they are the same in this case
598 
599  if(colExporter_==Teuchos::null)
600  colExporter_ = Teuchos::rcp(new Epetra_Export(*getGhostedColMap(),*getColMap()));
601 
602  return colExporter_;
603 }
604 
605 // "Build" functions
607 
608 template <typename Traits,typename LocalOrdinalT>
609 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildMap() const
610 {
611  Teuchos::RCP<Epetra_Map> map; // result
612  std::vector<int> indices;
613 
614  // get the global indices
615  gidProvider_->getOwnedIndices(indices);
616 
617  map = Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
618 
619  return map;
620 }
621 
622 template <typename Traits,typename LocalOrdinalT>
623 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildColMap() const
624 {
625  if(!hasColProvider_)
626  return buildMap();
627 
628  std::vector<int> indices;
629 
630  // get the global indices
631  colGidProvider_->getOwnedIndices(indices);
632 
633  return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
634 }
635 
636 // build the ghosted map
637 template <typename Traits,typename LocalOrdinalT>
638 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedMap() const
639 {
640  std::vector<int> indices;
641 
642  // get the global indices
643  gidProvider_->getOwnedAndGhostedIndices(indices);
644 
645  return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
646 }
647 
648 // build the ghosted map
649 template <typename Traits,typename LocalOrdinalT>
650 const Teuchos::RCP<Epetra_Map> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedColMap() const
651 {
652  if(!hasColProvider_)
653  return buildGhostedMap();
654 
655  std::vector<int> indices;
656 
657  // get the global indices
658  colGidProvider_->getOwnedAndGhostedIndices(indices);
659 
660  return Teuchos::rcp(new Epetra_Map(-1,indices.size(),&indices[0],0,*comm_));
661 }
662 
663 // get the graph of the crs matrix
664 template <typename Traits,typename LocalOrdinalT>
665 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGraph() const
666 {
667  using Teuchos::RCP;
668  using Teuchos::rcp;
669 
670  // build the map and allocate the space for the graph and
671  // grab the ghosted graph
672  RCP<Epetra_Map> rMap = getMap();
673  RCP<Epetra_Map> cMap = getColMap();
674  RCP<Epetra_CrsGraph> graph = rcp(new Epetra_CrsGraph(Copy,*rMap,0));
675  // RCP<Epetra_CrsGraph> oGraph = getGhostedGraph();
676  RCP<Epetra_CrsGraph> oGraph = buildFilteredGhostedGraph();
677  // this method calls getGhosted graph if no filtering is enabled
678  // otherwise a filtered matrix is constructed
679 
680  // perform the communication to finish building graph
681  RCP<Epetra_Export> exporter = getGhostedExport();
682  graph->Export( *oGraph, *exporter, Insert );
683  graph->FillComplete(*cMap,*rMap);
684  graph->OptimizeStorage();
685  return graph;
686 }
687 
688 template <typename Traits,typename LocalOrdinalT>
689 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildGhostedGraph(bool optimizeStorage) const
690 {
691  // build the map and allocate the space for the graph
692  Teuchos::RCP<Epetra_Map> rMap = getGhostedMap();
693  Teuchos::RCP<Epetra_Map> cMap = getGhostedColMap();
694  Teuchos::RCP<Epetra_CrsGraph> graph = Teuchos::rcp(new Epetra_CrsGraph(Copy,*rMap,*cMap,0));
695 
696  std::vector<std::string> elementBlockIds;
697  gidProvider_->getElementBlockIds(elementBlockIds);
698 
699  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,int> >
700  colGidProvider = hasColProvider_ ? colGidProvider_ : gidProvider_;
701  const Teuchos::RCP<const ConnManagerBase<LocalOrdinalT> > conn_mgr = colGidProvider->getConnManagerBase();
702  const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors();
703 
704  // graph information about the mesh
705  std::vector<std::string>::const_iterator blockItr;
706  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
707  std::string blockId = *blockItr;
708 
709  // grab elements for this block
710  const std::vector<LocalOrdinalT> & elements = gidProvider_->getElementBlock(blockId);
711 
712  // get information about number of indicies
713  std::vector<int> gids;
714  std::vector<int> col_gids;
715 
716  // loop over the elemnts
717  for(std::size_t i=0;i<elements.size();i++) {
718  gidProvider_->getElementGIDs(elements[i],gids);
719 
720  colGidProvider->getElementGIDs(elements[i],col_gids);
721  if (han) {
722  const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[i]);
723  for (typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
724  eit != aes.end(); ++eit) {
725  std::vector<int> other_col_gids;
726  colGidProvider->getElementGIDs(*eit, other_col_gids);
727  col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
728  }
729  }
730 
731  for(std::size_t j=0;j<gids.size();j++)
732  graph->InsertGlobalIndices(gids[j],col_gids.size(),&col_gids[0]);
733  }
734  }
735 
736  // finish filling the graph
737  graph->FillComplete(*cMap,*rMap);
738  if(optimizeStorage) {
739  graph->OptimizeStorage();
740  }
741 
742  return graph;
743 }
744 
745 template <typename Traits,typename LocalOrdinalT>
746 const Teuchos::RCP<Epetra_CrsGraph> EpetraLinearObjFactory<Traits,LocalOrdinalT>::buildFilteredGhostedGraph() const
747 {
748  using Teuchos::rcp;
749  using Teuchos::rcp_dynamic_cast;
750 
751  // figure out if the domain is filtered
752  RCP<const Filtered_UniqueGlobalIndexer<LocalOrdinalT,int> > filtered_ugi
753  = rcp_dynamic_cast<const Filtered_UniqueGlobalIndexer<LocalOrdinalT,int> >(getDomainGlobalIndexer());
754 
755  // domain is unfiltered, a filtered graph is just the original ghosted graph
756  if(filtered_ugi==Teuchos::null)
757  return getGhostedGraph();
758 
759  // get all local indices that are active (i.e. unfiltered)
760  std::vector<int> ghostedActive;
761  filtered_ugi->getOwnedAndGhostedNotFilteredIndicator(ghostedActive);
762 
763  // not sure how the deep copy works, so we will do this instead.
764  // This will build a new ghosted graph // without optimized storage so entries can be removed.
765  Teuchos::RCP<Epetra_CrsGraph> filteredGraph = buildGhostedGraph(false);
766  // false implies that storage is not optimzied
767 
768  // remove filtered column entries
769  for(int i=0;i<filteredGraph->NumMyRows();i++) {
770  std::vector<int> removedIndices;
771  int numIndices = 0;
772  int * indices = 0;
773  TEUCHOS_ASSERT(filteredGraph->ExtractMyRowView(i,numIndices,indices)==0);
774 
775  for(int j=0;j<numIndices;j++) {
776  if(ghostedActive[indices[j]]==0)
777  removedIndices.push_back(indices[j]);
778  }
779 
780  TEUCHOS_ASSERT(filteredGraph->RemoveMyIndices(i,Teuchos::as<int>(removedIndices.size()),&removedIndices[0])==0);
781  }
782 
783  // finish filling the graph
784  Teuchos::RCP<Epetra_Map> rMap = getGhostedMap();
785  Teuchos::RCP<Epetra_Map> cMap = getGhostedColMap();
786 
787  TEUCHOS_ASSERT(filteredGraph->FillComplete(*cMap,*rMap)==0);
788  TEUCHOS_ASSERT(filteredGraph->OptimizeStorage()==0);
789 
790  return filteredGraph;
791 }
792 
793 template <typename Traits,typename LocalOrdinalT>
794 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedEpetraColVector() const
795 {
796  Teuchos::RCP<const Epetra_Map> eMap = getGhostedColMap();
797  return Teuchos::rcp(new Epetra_Vector(*eMap));
798 }
799 
800 template <typename Traits,typename LocalOrdinalT>
801 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedEpetraVector() const
802 {
803  Teuchos::RCP<const Epetra_Map> eMap = getGhostedMap();
804  return Teuchos::rcp(new Epetra_Vector(*eMap));
805 }
806 
807 template <typename Traits,typename LocalOrdinalT>
808 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraColVector() const
809 {
810  Teuchos::RCP<const Epetra_Map> eMap = getColMap();
811  return Teuchos::rcp(new Epetra_Vector(*eMap));
812 }
813 
814 template <typename Traits,typename LocalOrdinalT>
815 Teuchos::RCP<Epetra_Vector> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraVector() const
816 {
817  Teuchos::RCP<const Epetra_Map> eMap = getMap();
818  return Teuchos::rcp(new Epetra_Vector(*eMap));
819 }
820 
821 template <typename Traits,typename LocalOrdinalT>
822 Teuchos::RCP<Epetra_CrsMatrix> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraMatrix() const
823 {
824  Teuchos::RCP<Epetra_CrsGraph> eGraph = getGraph();
825  return Teuchos::rcp(new Epetra_CrsMatrix(Copy, *eGraph));
826 }
827 
828 template <typename Traits,typename LocalOrdinalT>
829 Teuchos::RCP<Epetra_CrsMatrix> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getGhostedEpetraMatrix() const
830 {
831  Teuchos::RCP<Epetra_CrsGraph> eGraph = getGhostedGraph();
832  return Teuchos::rcp(new Epetra_CrsMatrix(Copy, *eGraph));
833 }
834 
835 template <typename Traits,typename LocalOrdinalT>
836 const Teuchos::RCP<const Epetra_Comm> EpetraLinearObjFactory<Traits,LocalOrdinalT>::getEpetraComm() const
837 {
838  return comm_;
839 }
840 
841 }
842 
843 #endif
Panzer_EpetraVector_ReadOnly_GlobalEvaluationData.hpp
Epetra_Map
Epetra_CrsGraph::FillComplete
int FillComplete()
is_null
bool is_null(const boost::shared_ptr< T > &p)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::dyn_cast
T_To & dyn_cast(T_From &from)
Epetra_CrsGraph::NumMyRows
int NumMyRows() const
Epetra_Export
Epetra_Vector
Epetra_MpiComm
Teuchos::ArrayRCP::size
size_type size() const
Epetra_Import
Epetra_CrsMatrix::ExtractMyRowView
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Epetra_CrsMatrix::PutScalar
int PutScalar(double ScalarConstant)
panzer::EpetraLinearObjFactory
BlockedEpetraLinearObjFactory< Traits, LocalOrdinalT > EpetraLinearObjFactory
Definition: Panzer_EpetraLinearObjFactory_decl.hpp:341
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Epetra_CrsGraph::OptimizeStorage
int OptimizeStorage()
panzer::values
std::vector< ScalarT > values
Definition: Panzer_ConstantFlux_decl.hpp:57
Epetra_CrsGraph::InsertGlobalIndices
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Epetra_CrsGraph::ExtractMyRowView
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Panzer_ConnManager.hpp
EpetraExt_VectorIn.h
Panzer_UniqueGlobalIndexer.hpp
EpetraExt::VectorToMatrixMarketFile
int VectorToMatrixMarketFile(const char *filename, const Epetra_Vector &A, const char *matrixName=0, const char *matrixDescription=0, bool writeHeader=true)
Panzer_Filtered_UniqueGlobalIndexer.hpp
Epetra_CrsGraph::RemoveMyIndices
int RemoveMyIndices(int LocalRow, int NumIndices, int *Indices)
panzer::comm
Teuchos::RCP< const Teuchos::Comm< int > > comm
Definition: Panzer_GlobalStatistics_decl.hpp:82
Epetra_CrsMatrix
EpetraExt::MatrixMarketFileToVector
int MatrixMarketFileToVector(const char *filename, const Epetra_BlockMap &map, Epetra_Vector *&A)
Panzer_ThyraObjContainer.hpp
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
EpetraExt_VectorOut.h
Epetra_CrsGraph
Teuchos::RCP::is_null
bool is_null() const
Generated on Wed Oct 25 2017 14:35:47 for Panzer by   doxygen 1.8.6