Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_TpetraLinearObjFactory_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_TpetraLinearObjFactory_impl_hpp__
44 #define __Panzer_TpetraLinearObjFactory_impl_hpp__
45 
46 // Panzer
47 #include "Panzer_ConnManager.hpp"
48 #ifdef PANZER_HAVE_EPETRA_STACK
49 #include "Panzer_EpetraVector_Write_GlobalEvaluationData.hpp" // JMG: Remove this eventually.
50 #endif
52 #include "Panzer_GlobalIndexer.hpp"
53 
54 // Thyra
55 #include "Thyra_TpetraVectorSpace.hpp"
56 #include "Thyra_TpetraLinearOp.hpp"
57 
58 // Tpetra
59 #include "Tpetra_MultiVector.hpp"
60 #include "Tpetra_Vector.hpp"
61 #include "Tpetra_CrsMatrix.hpp"
62 
63 namespace panzer {
64 
65 using Teuchos::RCP;
66 
67 // ************************************************************
68 // class TpetraLinearObjFactory
69 // ************************************************************
70 
71 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
74  const Teuchos::RCP<const GlobalIndexer> & gidProvider)
75  : comm_(comm), gidProvider_(gidProvider)
76 {
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 ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
87  const Teuchos::RCP<const GlobalIndexer> & gidProvider,
88  const Teuchos::RCP<const GlobalIndexer> & colGidProvider)
89  : comm_(comm), gidProvider_(gidProvider), colGidProvider_(colGidProvider)
90 {
91  hasColProvider_ = colGidProvider_!=Teuchos::null;
92 
93  // build and register the gather/scatter evaluators with
94  // the base class.
95  this->buildGatherScatterEvaluators(*this);
96 }
97 
98 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
101 { }
102 
103 // LinearObjectFactory functions
105 
106 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
110 {
111  Teuchos::RCP<ContainerType> container = Teuchos::rcp(new ContainerType(getColMap(),getMap()));
112 
113  return container;
114 }
115 
116 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
120 {
121  Teuchos::RCP<ContainerType> container = Teuchos::rcp(new ContainerType(getGhostedMap(),getGhostedMap()));
122 
123  return container;
124 }
125 
126 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
127 void
130  LinearObjContainer & out,int mem) const
131 {
132  using Teuchos::is_null;
133  typedef LinearObjContainer LOC;
134 
135  const ContainerType & t_in = Teuchos::dyn_cast<const ContainerType>(in);
137 
138  // Operations occur if the GLOBAL container has the correct targets!
139  // Users set the GLOBAL continer arguments
140  if ( !is_null(t_in.get_x()) && !is_null(t_out.get_x()) && ((mem & LOC::X)==LOC::X))
141  globalToGhostTpetraVector(*t_in.get_x(),*t_out.get_x(),true);
142 
143  if ( !is_null(t_in.get_dxdt()) && !is_null(t_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
144  globalToGhostTpetraVector(*t_in.get_dxdt(),*t_out.get_dxdt(),true);
145 
146  if ( !is_null(t_in.get_f()) && !is_null(t_out.get_f()) && ((mem & LOC::F)==LOC::F))
147  globalToGhostTpetraVector(*t_in.get_f(),*t_out.get_f(),false);
148 }
149 
150 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
151 void
154  LinearObjContainer & out,int mem) const
155 {
156  using Teuchos::is_null;
157 
158  typedef LinearObjContainer LOC;
159 
160  const ContainerType & t_in = Teuchos::dyn_cast<const ContainerType>(in);
162 
163  // Operations occur if the GLOBAL container has the correct targets!
164  // Users set the GLOBAL continer arguments
165  if ( !is_null(t_in.get_x()) && !is_null(t_out.get_x()) && ((mem & LOC::X)==LOC::X))
166  ghostToGlobalTpetraVector(*t_in.get_x(),*t_out.get_x(),true);
167 
168  if ( !is_null(t_in.get_f()) && !is_null(t_out.get_f()) && ((mem & LOC::F)==LOC::F))
169  ghostToGlobalTpetraVector(*t_in.get_f(),*t_out.get_f(),false);
170 
171  if ( !is_null(t_in.get_A()) && !is_null(t_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
172  ghostToGlobalTpetraMatrix(*t_in.get_A(),*t_out.get_A());
173 }
174 
175 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
176 void
178 ghostToGlobalTpetraVector(const Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & in,
179  Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & out, bool col) const
180 {
181  using Teuchos::RCP;
182 
183  // do the global distribution
184  RCP<ExportType> exporter = col ? getGhostedColExport() : getGhostedExport();
185  out.putScalar(0.0);
186  out.doExport(in,*exporter,Tpetra::ADD);
187 }
188 
189 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
190 void
192 ghostToGlobalTpetraMatrix(const Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & in,
193  Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & out) const
194 {
195  using Teuchos::RCP;
196 
197  // do the global distribution
198  RCP<ExportType> exporter = getGhostedExport();
199 
200  out.resumeFill();
201  out.setAllToScalar(0.0);
202  out.doExport(in,*exporter,Tpetra::ADD);
203  out.fillComplete(out.getDomainMap(),out.getRangeMap());
204 }
205 
206 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
207 void
209 globalToGhostTpetraVector(const Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & in,
210  Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & out, bool col) const
211 {
212  using Teuchos::RCP;
213 
214  // do the global distribution
215  RCP<ImportType> importer = col ? getGhostedColImport() : getGhostedImport();
216  out.putScalar(0.0);
217  out.doImport(in,*importer,Tpetra::INSERT);
218 }
219 
220 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
221 void
224  const LinearObjContainer & globalBCRows,
225  LinearObjContainer & ghostedObjs,
226  bool zeroVectorRows, bool adjustX) const
227 {
229 
230  const ContainerType & t_localBCRows = Teuchos::dyn_cast<const ContainerType>(localBCRows);
231  const ContainerType & t_globalBCRows = Teuchos::dyn_cast<const ContainerType>(globalBCRows);
232  ContainerType & t_ghosted = Teuchos::dyn_cast<ContainerType>(ghostedObjs);
233 
234  TEUCHOS_ASSERT(!Teuchos::is_null(t_localBCRows.get_f()));
235  TEUCHOS_ASSERT(!Teuchos::is_null(t_globalBCRows.get_f()));
236 
237  // pull out jacobian and vector
238  Teuchos::RCP<CrsMatrixType> A = t_ghosted.get_A();
239  Teuchos::RCP<VectorType> f = t_ghosted.get_f();
240  if(adjustX) f = t_ghosted.get_x();
241  Teuchos::ArrayRCP<double> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
242 
243  const VectorType & local_bcs = *(t_localBCRows.get_f());
244  const VectorType & global_bcs = *(t_globalBCRows.get_f());
245  Teuchos::ArrayRCP<const double> local_bcs_array = local_bcs.get1dView();
246  Teuchos::ArrayRCP<const double> global_bcs_array = global_bcs.get1dView();
247 
248  TEUCHOS_ASSERT(local_bcs_array.size()==global_bcs_array.size());
249  for(Ordinal i=0;i<local_bcs_array.size();i++) {
250  if(global_bcs_array[i]==0.0)
251  continue;
252 
253  if(local_bcs_array[i]==0.0 || zeroVectorRows) {
254  // this boundary condition was NOT set by this processor
255 
256  // if they exist put 0.0 in each entry
257  if(!Teuchos::is_null(f))
258  f_array[i] = 0.0;
259  if(!Teuchos::is_null(A)) {
260  std::size_t numEntries = 0;
261  std::size_t sz = A->getNumEntriesInLocalRow(i);
262  typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
263  typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
264 
265  A->getLocalRowCopy(i,indices,values,numEntries);
266 
267  for(std::size_t c=0;c<numEntries;c++)
268  values(c) = 0.0;
269 
270  A->replaceLocalValues(i,indices,values);
271  }
272  }
273  else {
274  // this boundary condition was set by this processor
275 
276  double scaleFactor = global_bcs_array[i];
277 
278  // if they exist scale linear objects by scale factor
279  if(!Teuchos::is_null(f))
280  f_array[i] /= scaleFactor;
281  if(!Teuchos::is_null(A)) {
282  std::size_t numEntries = 0;
283  std::size_t sz = A->getNumEntriesInLocalRow(i);
284  typename CrsMatrixType::nonconst_local_inds_host_view_type indices("indices", sz);
285  typename CrsMatrixType::nonconst_values_host_view_type values("values", sz);
286 
287  A->getLocalRowCopy(i,indices,values,numEntries);
288 
289  for(std::size_t c=0;c<numEntries;c++)
290  values(c) /= scaleFactor;
291 
292  A->replaceLocalValues(i,indices,values);
293  }
294  }
295  }
296 }
297 
298 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
299 void
301 applyDirichletBCs(const LinearObjContainer & /* counter */,
302  LinearObjContainer & /* result */) const
303 {
304  TEUCHOS_ASSERT(false); // not yet implemented
305 }
306 
308 //
309 // buildReadOnlyDomainContainer()
310 //
312 template<typename Traits, typename ScalarT, typename LocalOrdinalT,
313  typename GlobalOrdinalT, typename NodeT>
317 {
318  using Teuchos::rcp;
319  using TVROGED = TpetraVector_ReadOnly_GlobalEvaluationData<ScalarT,
320  LocalOrdinalT, GlobalOrdinalT, NodeT>;
321  auto ged = rcp(new TVROGED);
322  ged->initialize(getGhostedImport(), getGhostedColMap(), getColMap());
323  return ged;
324 } // end of buildReadOnlyDomainContainer()
325 
326 #ifdef PANZER_HAVE_EPETRA_STACK
327 //
329 // buildWriteDomainContainer()
330 //
332 template<typename Traits, typename ScalarT, typename LocalOrdinalT,
333  typename GlobalOrdinalT, typename NodeT>
337 {
338  using std::logic_error;
339  using Teuchos::rcp;
341  auto ged = rcp(new EVWGED);
342  TEUCHOS_TEST_FOR_EXCEPTION(true, logic_error, "NOT IMPLEMENTED YET")
343  return ged;
344 } // end of buildWriteDomainContainer()
345 #endif // PANZER_HAVE_EPETRA_STACK
346 
347 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
349 getComm() const
350 {
351  return *Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(getTeuchosComm());
352 }
353 
355 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
359 {
360  if(domainSpace_==Teuchos::null) {
361  if(!hasColProvider_)
362  domainSpace_ = Thyra::tpetraVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap());
363  else
364  domainSpace_ = Thyra::tpetraVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getColMap());
365  }
366 
367  return domainSpace_;
368 }
369 
371 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
375 {
376  if(rangeSpace_==Teuchos::null)
377  rangeSpace_ = Thyra::tpetraVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap());
378 
379  return rangeSpace_;
380 }
381 
383 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
387 {
388  return Thyra::tpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getThyraRangeSpace(),getThyraDomainSpace(),getTpetraMatrix());
389 }
390 
391 // Functions for initalizing a container
393 
394 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
395 void
398 {
400  initializeContainer(mem,tloc);
401 }
402 
403 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
404 void
407 {
408  typedef LinearObjContainer LOC;
409 
410  loc.clear();
411 
412  if((mem & LOC::X) == LOC::X)
413  loc.set_x(getTpetraColVector());
414 
415  if((mem & LOC::DxDt) == LOC::DxDt)
416  loc.set_dxdt(getTpetraColVector());
417 
418  if((mem & LOC::F) == LOC::F)
419  loc.set_f(getTpetraVector());
420 
421  if((mem & LOC::Mat) == LOC::Mat)
422  loc.set_A(getTpetraMatrix());
423 }
424 
425 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
426 void
429 {
431  initializeGhostedContainer(mem,tloc);
432 }
433 
434 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
435 void
438 {
439  typedef LinearObjContainer LOC;
440 
441  loc.clear();
442 
443  if((mem & LOC::X) == LOC::X)
444  loc.set_x(getGhostedTpetraColVector());
445 
446  if((mem & LOC::DxDt) == LOC::DxDt)
447  loc.set_dxdt(getGhostedTpetraColVector());
448 
449  if((mem & LOC::F) == LOC::F) {
450  loc.set_f(getGhostedTpetraVector());
452  }
453 
454  if((mem & LOC::Mat) == LOC::Mat) {
455  loc.set_A(getGhostedTpetraMatrix());
457  }
458 }
459 
460 // "Get" functions
462 
463 // get the map from the matrix
464 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
467 getMap() const
468 {
469  if(map_==Teuchos::null) map_ = buildMap();
470 
471  return map_;
472 }
473 
474 // get the map from the matrix
475 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
478 getColMap() const
479 {
480  if(cMap_==Teuchos::null) cMap_ = buildColMap();
481 
482  return cMap_;
483 }
484 
485 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
489 {
490  if(ghostedMap_==Teuchos::null) ghostedMap_ = buildGhostedMap();
491 
492  return ghostedMap_;
493 }
494 
495 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
499 {
500  if(cGhostedMap_==Teuchos::null) cGhostedMap_ = buildGhostedColMap();
501 
502  return cGhostedMap_;
503 }
504 
505 // get the graph of the crs matrix
506 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
509 getGraph() const
510 {
511  if(graph_==Teuchos::null) graph_ = buildGraph();
512 
513  return graph_;
514 }
515 
516 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
520 {
521  if(ghostedGraph_==Teuchos::null) ghostedGraph_ = buildGhostedGraph();
522 
523  return ghostedGraph_;
524 }
525 
526 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
530 {
531  if(ghostedImporter_==Teuchos::null)
532  ghostedImporter_ = Teuchos::rcp(new ImportType(getMap(),getGhostedMap()));
533 
534  return ghostedImporter_;
535 }
536 
537 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
541 {
542  if(!hasColProvider_)
543  ghostedColImporter_ = getGhostedImport(); // they are the same in this case
544 
545  if(ghostedColImporter_==Teuchos::null)
546  ghostedColImporter_ = Teuchos::rcp(new ImportType(getColMap(),getGhostedColMap()));
547 
548  return ghostedColImporter_;
549 }
550 
551 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
555 {
556  if(ghostedExporter_==Teuchos::null)
557  ghostedExporter_ = Teuchos::rcp(new ExportType(getGhostedMap(),getMap()));
558 
559  return ghostedExporter_;
560 }
561 
562 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
566 {
567  if(!hasColProvider_)
568  ghostedColExporter_ = getGhostedExport(); // they are the same in this case
569 
570  if(ghostedColExporter_==Teuchos::null)
571  ghostedColExporter_ = Teuchos::rcp(new ExportType(getGhostedColMap(),getColMap()));
572 
573  return ghostedColExporter_;
574 }
575 
576 // "Build" functions
578 
579 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
582 buildMap() const
583 {
584  std::vector<GlobalOrdinalT> indices;
585 
586  // get the global indices
587  gidProvider_->getOwnedIndices(indices);
588 
590 }
591 
592 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
595 buildColMap() const
596 {
597  if(!hasColProvider_)
598  return buildMap();
599 
600  std::vector<GlobalOrdinalT> indices;
601 
602  // get the global indices
603  colGidProvider_->getOwnedIndices(indices);
604 
606 }
607 
608 // build the ghosted map
609 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
613 {
614  std::vector<GlobalOrdinalT> indices;
615 
616  // get the global indices
617  gidProvider_->getOwnedAndGhostedIndices(indices);
618 
620 }
621 
622 // build the ghosted map
623 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
627 {
628  if(!hasColProvider_)
629  return buildGhostedMap();
630 
631  std::vector<GlobalOrdinalT> indices;
632 
633  // get the global indices
634  colGidProvider_->getOwnedAndGhostedIndices(indices);
635 
637 }
638 
639 // get the graph of the crs matrix
640 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
643 buildGraph() const
644 {
645  using Teuchos::RCP;
646  using Teuchos::rcp;
647 
648  // build the map and allocate the space for the graph and
649  // grab the ghosted graph
650  RCP<MapType> rMap = getMap();
651  RCP<MapType> cMap = getColMap();
652  RCP<CrsGraphType> graph = rcp(new CrsGraphType(rMap,0));
653  RCP<CrsGraphType> oGraph = getGhostedGraph();
654 
655  // perform the communication to finish building graph
656  RCP<ExportType> exporter = getGhostedExport();
657  graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
658  graph->fillComplete(cMap,rMap);
659 
660  return graph;
661 }
662 
663 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
667 {
668  // build the map and allocate the space for the graph
669  Teuchos::RCP<MapType> rMap = getGhostedMap();
670  Teuchos::RCP<MapType> cMap = getGhostedColMap();
671 
672  std::vector<std::string> elementBlockIds;
673  gidProvider_->getElementBlockIds(elementBlockIds);
674 
676  colGidProvider = hasColProvider_ ? colGidProvider_ : gidProvider_;
677  const Teuchos::RCP<const ConnManager> conn_mgr = colGidProvider->getConnManager();
678  const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors();
679 
680  // graph information about the mesh
681  // Count number of entries per graph row; needed for graph constructor
682  std::vector<size_t> nEntriesPerRow(rMap->getLocalNumElements(), 0);
683 
684  std::vector<std::string>::const_iterator blockItr;
685  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
686  std::string blockId = *blockItr;
687 
688  // grab elements for this block
689  const std::vector<LocalOrdinalT> & elements = gidProvider_->getElementBlock(blockId);
690 
691  // get information about number of indicies
692  std::vector<GlobalOrdinalT> gids;
693  std::vector<GlobalOrdinalT> col_gids;
694 
695  // loop over the elemnts
696  for(std::size_t i=0;i<elements.size();i++) {
697  gidProvider_->getElementGIDs(elements[i],gids);
698 
699  colGidProvider->getElementGIDs(elements[i],col_gids);
700  if (han) {
701  const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[i]);
702  for (typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
703  eit != aes.end(); ++eit) {
704  std::vector<GlobalOrdinalT> other_col_gids;
705  colGidProvider->getElementGIDs(*eit, other_col_gids);
706  col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
707  }
708  }
709 
710  for(std::size_t j=0;j<gids.size();j++){
711  LocalOrdinalT lid = rMap->getLocalElement(gids[j]);
712  nEntriesPerRow[lid] += col_gids.size();
713  }
714  }
715  }
716 
717  Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
719  nEntriesPerRowView));
720 
721  // Now insert entries into the graph
722  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
723  std::string blockId = *blockItr;
724 
725  // grab elements for this block
726  const std::vector<LocalOrdinalT> & elements = gidProvider_->getElementBlock(blockId);
727 
728  // get information about number of indicies
729  std::vector<GlobalOrdinalT> gids;
730  std::vector<GlobalOrdinalT> col_gids;
731 
732  // loop over the elemnts
733  for(std::size_t i=0;i<elements.size();i++) {
734  gidProvider_->getElementGIDs(elements[i],gids);
735 
736  colGidProvider->getElementGIDs(elements[i],col_gids);
737  if (han) {
738  const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[i]);
739  for (typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
740  eit != aes.end(); ++eit) {
741  std::vector<GlobalOrdinalT> other_col_gids;
742  colGidProvider->getElementGIDs(*eit, other_col_gids);
743  col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
744  }
745  }
746 
747  for(std::size_t j=0;j<gids.size();j++)
748  graph->insertGlobalIndices(gids[j],col_gids);
749  }
750  }
751 
752  // finish filling the graph
753  graph->fillComplete(cMap,rMap);
754 
755  return graph;
756 }
757 
758 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
762 {
763  Teuchos::RCP<const MapType> tMap = getGhostedMap();
764  return Teuchos::rcp(new VectorType(tMap));
765 }
766 
767 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
771 {
772  Teuchos::RCP<const MapType> tMap = getGhostedColMap();
773  return Teuchos::rcp(new VectorType(tMap));
774 }
775 
776 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
780 {
781  Teuchos::RCP<const MapType> tMap = getMap();
782  return Teuchos::rcp(new VectorType(tMap));
783 }
784 
785 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
789 {
790  Teuchos::RCP<const MapType> tMap = getColMap();
791  return Teuchos::rcp(new VectorType(tMap));
792 }
793 
794 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
798 {
799  Teuchos::RCP<CrsGraphType> tGraph = getGraph();
801  tMat->fillComplete(tMat->getDomainMap(),tMat->getRangeMap());
802 
803  return tMat;
804 }
805 
806 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
810 {
811  Teuchos::RCP<CrsGraphType> tGraph = getGhostedGraph();
813  tMat->fillComplete(tMat->getDomainMap(),tMat->getRangeMap());
814 
815  return tMat;
816 }
817 
818 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
822 {
823  return comm_;
824 }
825 
826 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
829 {
832  if(A!=Teuchos::null)
833  A->resumeFill();
834 }
835 
836 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
839 {
842  if(A!=Teuchos::null)
843  A->fillComplete(A->getDomainMap(),A->getRangeMap());
844 }
845 
846 }
847 
848 #endif // __Panzer_TpetraLinearObjFactory_impl_hpp__
virtual const Teuchos::RCP< Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedExport() const
get exporter for converting an overalapped object to a &quot;normal&quot; object
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraVector() const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColMap() const
const Teuchos::RCP< CrsMatrixType > get_A() const
virtual Teuchos::MpiComm< int > getComm() const
bool is_null(const boost::shared_ptr< T > &p)
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
void initializeGhostedContainer(int, LinearObjContainer &loc) const
const Teuchos::RCP< VectorType > get_x() const
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildReadOnlyDomainContainer() const
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraVector() const
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedGraph() const
bool is_null(const std::shared_ptr< T > &p)
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain space.
Teuchos::RCP< Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraMatrix() const
virtual const Teuchos::RCP< Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedImport() const
get importer for converting an overalapped object to a &quot;normal&quot; object
const Teuchos::RCP< VectorType > get_dxdt() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void initializeContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraColVector() const
T_To & dyn_cast(T_From &from)
virtual void endFill(LinearObjContainer &loc) const
size_type size() const
virtual const Teuchos::RCP< Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColImport() const
virtual void beginFill(LinearObjContainer &loc) const
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
virtual const Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
get exporter for converting an overalapped object to a &quot;normal&quot; object
Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > ImportType
virtual const Teuchos::RCP< Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColExport() const
Teuchos::RCP< Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraMatrix() const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getColMap() const
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
virtual const std::vector< LocalOrdinal > & getAssociatedNeighbors(const LocalOrdinal &el) const =0
void globalToGhostTpetraVector(const Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out, bool col) const
void set_x(const Teuchos::RCP< VectorType > &in)
void ghostToGlobalTpetraMatrix(const Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out) const
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
Teuchos_Ordinal Ordinal
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedMap() const
get the ghosted map from the matrix
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildColMap() const
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGraph() const
virtual Teuchos::RCP< const ConnManager > getConnManager() const =0
Returns the connection manager currently being used.
void set_dxdt(const Teuchos::RCP< VectorType > &in)
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraColVector() const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getMap() const
get the map from the matrix
This class provides a boundary exchange communication mechanism for vectors.
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGraph() const
get the graph of the crs matrix
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
TpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< const GlobalIndexer > &gidProvider)
void set_f(const Teuchos::RCP< VectorType > &in)
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildMap() const
Teuchos::RCP< const GlobalIndexer > colGidProvider_
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range space.
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedColMap() const
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
virtual void getElementGIDs(panzer::LocalOrdinal localElmtId, std::vector< panzer::GlobalOrdinal > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
virtual bool hasAssociatedNeighbors() const =0
void buildGatherScatterEvaluators(const BuilderT &builder)
virtual Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a matrix operator.
#define TEUCHOS_ASSERT(assertion_test)
void ghostToGlobalTpetraVector(const Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out, bool col) const
const Teuchos::RCP< VectorType > get_f() const
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedGraph() const
get the ghosted graph of the crs matrix
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedMap() const
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
bool is_null() const