Panzer  Version of the Day
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_TPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
44 #define PANZER_TPETRA_LINEAR_OBJ_FACTORY_IMPL_HPP
45 
48 #include "Panzer_ConnManager.hpp"
49 
50 #include "Tpetra_MultiVector.hpp"
51 #include "Tpetra_Vector.hpp"
52 #include "Tpetra_CrsMatrix.hpp"
53 
54 #include "Thyra_TpetraVectorSpace.hpp"
55 #include "Thyra_TpetraLinearOp.hpp"
56 
57 using Teuchos::RCP;
58 
59 namespace panzer {
60 
61 // ************************************************************
62 // class TpetraLinearObjFactory
63 // ************************************************************
64 
65 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
67 TpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::Comm<int> > & comm,
68  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> > & gidProvider)
69  : comm_(comm), gidProvider_(gidProvider)
70 {
71  hasColProvider_ = colGidProvider_!=Teuchos::null;
72 
73  // build and register the gather/scatter evaluators with
74  // the base class.
75  this->buildGatherScatterEvaluators(*this);
76 }
77 
78 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
80 TpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::Comm<int> > & comm,
81  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> > & gidProvider,
82  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> > & colGidProvider)
83  : comm_(comm), gidProvider_(gidProvider), colGidProvider_(colGidProvider)
84 {
85  hasColProvider_ = colGidProvider_!=Teuchos::null;
86 
87  // build and register the gather/scatter evaluators with
88  // the base class.
89  this->buildGatherScatterEvaluators(*this);
90 }
91 
92 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
95 { }
96 
97 // LinearObjectFactory functions
99 
100 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
101 Teuchos::RCP<LinearObjContainer>
104 {
105  Teuchos::RCP<ContainerType> container = Teuchos::rcp(new ContainerType(getColMap(),getMap()));
106 
107  return container;
108 }
109 
110 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
111 Teuchos::RCP<LinearObjContainer>
114 {
115  Teuchos::RCP<ContainerType> container = Teuchos::rcp(new ContainerType(getGhostedMap(),getGhostedMap()));
116 
117  return container;
118 }
119 
120 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
121 void
124  LinearObjContainer & out,int mem) const
125 {
126  using Teuchos::is_null;
127  typedef LinearObjContainer LOC;
128 
129  const ContainerType & t_in = Teuchos::dyn_cast<const ContainerType>(in);
130  ContainerType & t_out = Teuchos::dyn_cast<ContainerType>(out);
131 
132  // Operations occur if the GLOBAL container has the correct targets!
133  // Users set the GLOBAL continer arguments
134  if ( !is_null(t_in.get_x()) && !is_null(t_out.get_x()) && ((mem & LOC::X)==LOC::X))
135  globalToGhostTpetraVector(*t_in.get_x(),*t_out.get_x(),true);
136 
137  if ( !is_null(t_in.get_dxdt()) && !is_null(t_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
138  globalToGhostTpetraVector(*t_in.get_dxdt(),*t_out.get_dxdt(),true);
139 
140  if ( !is_null(t_in.get_f()) && !is_null(t_out.get_f()) && ((mem & LOC::F)==LOC::F))
141  globalToGhostTpetraVector(*t_in.get_f(),*t_out.get_f(),false);
142 }
143 
144 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
145 void
148  LinearObjContainer & out,int mem) const
149 {
150  using Teuchos::is_null;
151 
152  typedef LinearObjContainer LOC;
153 
154  const ContainerType & t_in = Teuchos::dyn_cast<const ContainerType>(in);
155  ContainerType & t_out = Teuchos::dyn_cast<ContainerType>(out);
156 
157  // Operations occur if the GLOBAL container has the correct targets!
158  // Users set the GLOBAL continer arguments
159  if ( !is_null(t_in.get_x()) && !is_null(t_out.get_x()) && ((mem & LOC::X)==LOC::X))
160  ghostToGlobalTpetraVector(*t_in.get_x(),*t_out.get_x(),true);
161 
162  if ( !is_null(t_in.get_f()) && !is_null(t_out.get_f()) && ((mem & LOC::F)==LOC::F))
163  ghostToGlobalTpetraVector(*t_in.get_f(),*t_out.get_f(),false);
164 
165  if ( !is_null(t_in.get_A()) && !is_null(t_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
166  ghostToGlobalTpetraMatrix(*t_in.get_A(),*t_out.get_A());
167 }
168 
169 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
170 void
172 ghostToGlobalTpetraVector(const Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & in,
173  Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & out, bool col) const
174 {
175  using Teuchos::RCP;
176 
177  // do the global distribution
178  RCP<ExportType> exporter = col ? getGhostedColExport() : getGhostedExport();
179  out.putScalar(0.0);
180  out.doExport(in,*exporter,Tpetra::ADD);
181 }
182 
183 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
184 void
186 ghostToGlobalTpetraMatrix(const Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & in,
187  Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & out) const
188 {
189  using Teuchos::RCP;
190 
191  // do the global distribution
192  RCP<ExportType> exporter = getGhostedExport();
193 
194  out.resumeFill();
195  out.setAllToScalar(0.0);
196  out.doExport(in,*exporter,Tpetra::ADD);
197  out.fillComplete(out.getDomainMap(),out.getRangeMap());
198 }
199 
200 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
201 void
203 globalToGhostTpetraVector(const Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & in,
204  Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> & out, bool col) const
205 {
206  using Teuchos::RCP;
207 
208  // do the global distribution
209  RCP<ImportType> importer = col ? getGhostedColImport() : getGhostedImport();
210  out.putScalar(0.0);
211  out.doImport(in,*importer,Tpetra::INSERT);
212 }
213 
214 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
215 void
218  const LinearObjContainer & globalBCRows,
219  LinearObjContainer & ghostedObjs,
220  bool zeroVectorRows, bool adjustX) const
221 {
222  typedef Teuchos::ArrayRCP<const double>::Ordinal Ordinal;
223 
224  const ContainerType & t_localBCRows = Teuchos::dyn_cast<const ContainerType>(localBCRows);
225  const ContainerType & t_globalBCRows = Teuchos::dyn_cast<const ContainerType>(globalBCRows);
226  ContainerType & t_ghosted = Teuchos::dyn_cast<ContainerType>(ghostedObjs);
227 
228  TEUCHOS_ASSERT(!Teuchos::is_null(t_localBCRows.get_f()));
229  TEUCHOS_ASSERT(!Teuchos::is_null(t_globalBCRows.get_f()));
230 
231  // pull out jacobian and vector
232  Teuchos::RCP<CrsMatrixType> A = t_ghosted.get_A();
233  Teuchos::RCP<VectorType> f = t_ghosted.get_f();
234  if(adjustX) f = t_ghosted.get_x();
235  Teuchos::ArrayRCP<double> f_array = f!=Teuchos::null ? f->get1dViewNonConst() : Teuchos::null;
236 
237  const VectorType & local_bcs = *(t_localBCRows.get_f());
238  const VectorType & global_bcs = *(t_globalBCRows.get_f());
239  Teuchos::ArrayRCP<const double> local_bcs_array = local_bcs.get1dView();
240  Teuchos::ArrayRCP<const double> global_bcs_array = global_bcs.get1dView();
241 
242  TEUCHOS_ASSERT(local_bcs_array.size()==global_bcs_array.size());
243  for(Ordinal i=0;i<local_bcs_array.size();i++) {
244  if(global_bcs_array[i]==0.0)
245  continue;
246 
247  if(local_bcs_array[i]==0.0 || zeroVectorRows) {
248  // this boundary condition was NOT set by this processor
249 
250  // if they exist put 0.0 in each entry
251  if(!Teuchos::is_null(f))
252  f_array[i] = 0.0;
253  if(!Teuchos::is_null(A)) {
254  std::size_t numEntries = 0;
255  std::size_t sz = A->getNumEntriesInLocalRow(i);
256  Teuchos::Array<LocalOrdinalT> indices(sz);
257  Teuchos::Array<double> values(sz);
258 
259  A->getLocalRowCopy(i,indices,values,numEntries);
260 
261  for(std::size_t c=0;c<numEntries;c++)
262  values[c] = 0.0;
263 
264  A->replaceLocalValues(i,indices,values);
265  }
266  }
267  else {
268  // this boundary condition was set by this processor
269 
270  double scaleFactor = global_bcs_array[i];
271 
272  // if they exist scale linear objects by scale factor
273  if(!Teuchos::is_null(f))
274  f_array[i] /= scaleFactor;
275  if(!Teuchos::is_null(A)) {
276  std::size_t numEntries = 0;
277  std::size_t sz = A->getNumEntriesInLocalRow(i);
278  Teuchos::Array<LocalOrdinalT> indices(sz);
279  Teuchos::Array<double> values(sz);
280 
281  A->getLocalRowCopy(i,indices,values,numEntries);
282 
283  for(std::size_t c=0;c<numEntries;c++)
284  values[c] /= scaleFactor;
285 
286  A->replaceLocalValues(i,indices,values);
287  }
288  }
289  }
290 }
291 
292 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
293 void
296  LinearObjContainer & result) const
297 {
298  TEUCHOS_ASSERT(false); // not yet implemented
299 }
300 
301 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
302 Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
305 {
306  // TEUCHOS_ASSERT(false);
307  //return Teuchos::null;
308  Teuchos::RCP<TpetraVector_ReadOnly_GlobalEvaluationData<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > vec_ged
310  vec_ged->initialize(getGhostedImport(),getGhostedColMap(),getColMap());
311 
312  return vec_ged;
313 }
314 
315 
316 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
318 getComm() const
319 {
320  return *Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(getTeuchosComm());
321 }
322 
324 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
325 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
328 {
329  if(domainSpace_==Teuchos::null) {
330  if(!hasColProvider_)
331  domainSpace_ = Thyra::tpetraVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap());
332  else
333  domainSpace_ = Thyra::tpetraVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getColMap());
334  }
335 
336  return domainSpace_;
337 }
338 
340 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
341 Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
344 {
345  if(rangeSpace_==Teuchos::null)
346  rangeSpace_ = Thyra::tpetraVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap());
347 
348  return rangeSpace_;
349 }
350 
352 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
353 Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
356 {
357  return Thyra::tpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getThyraRangeSpace(),getThyraDomainSpace(),getTpetraMatrix());
358 }
359 
360 // Functions for initalizing a container
362 
363 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
364 void
367 {
368  ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
369  initializeContainer(mem,tloc);
370 }
371 
372 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
373 void
376 {
377  typedef LinearObjContainer LOC;
378 
379  loc.clear();
380 
381  if((mem & LOC::X) == LOC::X)
382  loc.set_x(getTpetraColVector());
383 
384  if((mem & LOC::DxDt) == LOC::DxDt)
385  loc.set_dxdt(getTpetraColVector());
386 
387  if((mem & LOC::F) == LOC::F)
388  loc.set_f(getTpetraVector());
389 
390  if((mem & LOC::Mat) == LOC::Mat)
391  loc.set_A(getTpetraMatrix());
392 }
393 
394 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
395 void
398 {
399  ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
400  initializeGhostedContainer(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(getGhostedTpetraColVector());
414 
415  if((mem & LOC::DxDt) == LOC::DxDt)
416  loc.set_dxdt(getGhostedTpetraColVector());
417 
418  if((mem & LOC::F) == LOC::F) {
419  loc.set_f(getGhostedTpetraVector());
421  }
422 
423  if((mem & LOC::Mat) == LOC::Mat) {
424  loc.set_A(getGhostedTpetraMatrix());
426  }
427 }
428 
429 // "Get" functions
431 
432 // get the map from the matrix
433 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
434 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
436 getMap() const
437 {
438  if(map_==Teuchos::null) map_ = buildMap();
439 
440  return map_;
441 }
442 
443 // get the map from the matrix
444 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
445 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
447 getColMap() const
448 {
449  if(cMap_==Teuchos::null) cMap_ = buildColMap();
450 
451  return cMap_;
452 }
453 
454 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
455 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
458 {
459  if(ghostedMap_==Teuchos::null) ghostedMap_ = buildGhostedMap();
460 
461  return ghostedMap_;
462 }
463 
464 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
465 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
468 {
469  if(cGhostedMap_==Teuchos::null) cGhostedMap_ = buildGhostedColMap();
470 
471  return cGhostedMap_;
472 }
473 
474 // get the graph of the crs matrix
475 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
476 const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
478 getGraph() const
479 {
480  if(graph_==Teuchos::null) graph_ = buildGraph();
481 
482  return graph_;
483 }
484 
485 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
486 const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
489 {
490  if(ghostedGraph_==Teuchos::null) ghostedGraph_ = buildGhostedGraph();
491 
492  return ghostedGraph_;
493 }
494 
495 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
496 const Teuchos::RCP<Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
499 {
500  if(ghostedImporter_==Teuchos::null)
501  ghostedImporter_ = Teuchos::rcp(new ImportType(getMap(),getGhostedMap()));
502 
503  return ghostedImporter_;
504 }
505 
506 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
507 const Teuchos::RCP<Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
510 {
511  if(!hasColProvider_)
512  ghostedColImporter_ = getGhostedImport(); // they are the same in this case
513 
514  if(ghostedColImporter_==Teuchos::null)
515  ghostedColImporter_ = Teuchos::rcp(new ImportType(getColMap(),getGhostedColMap()));
516 
517  return ghostedColImporter_;
518 }
519 
520 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
521 const Teuchos::RCP<Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
524 {
525  if(ghostedExporter_==Teuchos::null)
526  ghostedExporter_ = Teuchos::rcp(new ExportType(getGhostedMap(),getMap()));
527 
528  return ghostedExporter_;
529 }
530 
531 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
532 const Teuchos::RCP<Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
535 {
536  if(!hasColProvider_)
537  ghostedColExporter_ = getGhostedExport(); // they are the same in this case
538 
539  if(ghostedColExporter_==Teuchos::null)
540  ghostedColExporter_ = Teuchos::rcp(new ExportType(getGhostedColMap(),getColMap()));
541 
542  return ghostedColExporter_;
543 }
544 
545 // "Build" functions
547 
548 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
549 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
551 buildMap() const
552 {
553  std::vector<GlobalOrdinalT> indices;
554 
555  // get the global indices
556  gidProvider_->getOwnedIndices(indices);
557 
558  return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
559 }
560 
561 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
562 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
564 buildColMap() const
565 {
566  if(!hasColProvider_)
567  return buildMap();
568 
569  std::vector<GlobalOrdinalT> indices;
570 
571  // get the global indices
572  colGidProvider_->getOwnedIndices(indices);
573 
574  return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
575 }
576 
577 // build the ghosted map
578 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
579 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
582 {
583  std::vector<GlobalOrdinalT> indices;
584 
585  // get the global indices
586  gidProvider_->getOwnedAndGhostedIndices(indices);
587 
588  return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
589 }
590 
591 // build the ghosted map
592 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
593 const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
596 {
597  if(!hasColProvider_)
598  return buildGhostedMap();
599 
600  std::vector<GlobalOrdinalT> indices;
601 
602  // get the global indices
603  colGidProvider_->getOwnedAndGhostedIndices(indices);
604 
605  return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
606 }
607 
608 // get the graph of the crs matrix
609 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
610 const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
612 buildGraph() const
613 {
614  using Teuchos::RCP;
615  using Teuchos::rcp;
616 
617  // build the map and allocate the space for the graph and
618  // grab the ghosted graph
619  RCP<MapType> rMap = getMap();
620  RCP<MapType> cMap = getColMap();
621  RCP<CrsGraphType> graph = rcp(new CrsGraphType(rMap,0));
622  RCP<CrsGraphType> oGraph = getGhostedGraph();
623 
624  // perform the communication to finish building graph
625  RCP<ExportType> exporter = getGhostedExport();
626  graph->doExport( *oGraph, *exporter, Tpetra::INSERT );
627  graph->fillComplete(cMap,rMap);
628 
629  return graph;
630 }
631 
632 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
633 const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
636 {
637  // build the map and allocate the space for the graph
638  Teuchos::RCP<MapType> rMap = getGhostedMap();
639  Teuchos::RCP<MapType> cMap = getGhostedColMap();
640  Teuchos::RCP<CrsGraphType> graph = Teuchos::rcp(new CrsGraphType(rMap,cMap,0));
641 
642  std::vector<std::string> elementBlockIds;
643  gidProvider_->getElementBlockIds(elementBlockIds);
644 
645  const Teuchos::RCP<const UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> >
646  colGidProvider = hasColProvider_ ? colGidProvider_ : gidProvider_;
647  const Teuchos::RCP<const ConnManagerBase<LocalOrdinalT> > conn_mgr = colGidProvider->getConnManagerBase();
648  const bool han = conn_mgr.is_null() ? false : conn_mgr->hasAssociatedNeighbors();
649 
650  // graph information about the mesh
651  std::vector<std::string>::const_iterator blockItr;
652  for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
653  std::string blockId = *blockItr;
654 
655  // grab elements for this block
656  const std::vector<LocalOrdinalT> & elements = gidProvider_->getElementBlock(blockId);
657 
658  // get information about number of indicies
659  std::vector<GlobalOrdinalT> gids;
660  std::vector<GlobalOrdinalT> col_gids;
661 
662  // loop over the elemnts
663  for(std::size_t i=0;i<elements.size();i++) {
664  gidProvider_->getElementGIDs(elements[i],gids);
665 
666  colGidProvider->getElementGIDs(elements[i],col_gids);
667  if (han) {
668  const std::vector<LocalOrdinalT>& aes = conn_mgr->getAssociatedNeighbors(elements[i]);
669  for (typename std::vector<LocalOrdinalT>::const_iterator eit = aes.begin();
670  eit != aes.end(); ++eit) {
671  std::vector<GlobalOrdinalT> other_col_gids;
672  colGidProvider->getElementGIDs(*eit, other_col_gids);
673  col_gids.insert(col_gids.end(), other_col_gids.begin(), other_col_gids.end());
674  }
675  }
676 
677  for(std::size_t j=0;j<gids.size();j++)
678  graph->insertGlobalIndices(gids[j],col_gids);
679  }
680  }
681 
682  // finish filling the graph
683  graph->fillComplete(cMap,rMap);
684 
685  return graph;
686 }
687 
688 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
689 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
692 {
693  Teuchos::RCP<const MapType> tMap = getGhostedMap();
694  return Teuchos::rcp(new VectorType(tMap));
695 }
696 
697 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
698 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
701 {
702  Teuchos::RCP<const MapType> tMap = getGhostedColMap();
703  return Teuchos::rcp(new VectorType(tMap));
704 }
705 
706 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
707 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
710 {
711  Teuchos::RCP<const MapType> tMap = getMap();
712  return Teuchos::rcp(new VectorType(tMap));
713 }
714 
715 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
716 Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
719 {
720  Teuchos::RCP<const MapType> tMap = getColMap();
721  return Teuchos::rcp(new VectorType(tMap));
722 }
723 
724 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
725 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
728 {
729  Teuchos::RCP<CrsGraphType> tGraph = getGraph();
730  Teuchos::RCP<CrsMatrixType> tMat = Teuchos::rcp(new CrsMatrixType(tGraph));
731  tMat->fillComplete(tMat->getDomainMap(),tMat->getRangeMap());
732 
733  return tMat;
734 }
735 
736 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
737 Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
740 {
741  Teuchos::RCP<CrsGraphType> tGraph = getGhostedGraph();
742  Teuchos::RCP<CrsMatrixType> tMat = Teuchos::rcp(new CrsMatrixType(tGraph));
743  tMat->fillComplete(tMat->getDomainMap(),tMat->getRangeMap());
744 
745  return tMat;
746 }
747 
748 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
749 const Teuchos::RCP<const Teuchos::Comm<int> >
752 {
753  return comm_;
754 }
755 
756 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
759 {
760  ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
761  Teuchos::RCP<CrsMatrixType> A = tloc.get_A();
762  if(A!=Teuchos::null)
763  A->resumeFill();
764 }
765 
766 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
769 {
770  ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
771  Teuchos::RCP<CrsMatrixType> A = tloc.get_A();
772  if(A!=Teuchos::null)
773  A->fillComplete(A->getDomainMap(),A->getRangeMap());
774 }
775 
776 }
777 
778 #endif
virtual const Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
get exporter for converting an overalapped object to a "normal" object
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
void globalToGhostTpetraVector(const Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out, bool col) const
void initializeContainer(int, LinearObjContainer &loc) const
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range space.
void ghostToGlobalTpetraMatrix(const Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out) const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedMap() const
get the ghosted map from the matrix
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
virtual void beginFill(LinearObjContainer &loc) const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getMap() const
get the map from the matrix
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedMap() const
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedGraph() const
get the ghosted graph of the crs matrix
Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > ImportType
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedGraph() const
void set_x(const Teuchos::RCP< VectorType > &in)
virtual const Teuchos::RCP< Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColExport() const
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
virtual const Teuchos::RCP< Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColImport() const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildMap() const
TpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< const UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > &gidProvider)
std::vector< ScalarT > values
void ghostToGlobalTpetraVector(const Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out, bool col) const
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraColVector() const
virtual const Teuchos::RCP< Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedImport() const
get importer for converting an overalapped object to a "normal" object
Teuchos::RCP< const UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > colGidProvider_
void set_dxdt(const Teuchos::RCP< VectorType > &in)
const Teuchos::RCP< CrsMatrixType > get_A() const
Teuchos::RCP< Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraMatrix() const
const Teuchos::RCP< VectorType > get_x() const
Teuchos::RCP< Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraMatrix() const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedColMap() const
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGraph() const
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
virtual Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a matrix operator.
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraVector() const
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraVector() const
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
const Teuchos::RCP< VectorType > get_dxdt() const
void set_f(const Teuchos::RCP< VectorType > &in)
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColMap() const
virtual void endFill(LinearObjContainer &loc) const
Teuchos::RCP< const Teuchos::Comm< int > > comm
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildDomainContainer() const
virtual Teuchos::MpiComm< int > getComm() const
const Teuchos::RCP< VectorType > get_f() const
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
virtual const Teuchos::RCP< Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedExport() const
get exporter for converting an overalapped object to a "normal" object
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain space.
void buildGatherScatterEvaluators(const BuilderT &builder)
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGraph() const
get the graph of the crs matrix
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getColMap() const
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildColMap() const
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
void initializeGhostedContainer(int, LinearObjContainer &loc) const
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraColVector() const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const