Panzer Version of the Day
Loading...
Searching...
No Matches
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
48#ifdef PANZER_HAVE_EPETRA_STACK
49#include "Panzer_EpetraVector_Write_GlobalEvaluationData.hpp" // JMG: Remove this eventually.
50#endif
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
63namespace panzer {
64
65using Teuchos::RCP;
66
67// ************************************************************
68// class TpetraLinearObjFactory
69// ************************************************************
70
71template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
73TpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::Comm<int> > & comm,
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.
82}
83
84template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
86TpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::Comm<int> > & comm,
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.
96}
97
98template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
102
103// LinearObjectFactory functions
105
106template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
107Teuchos::RCP<LinearObjContainer>
110{
111 Teuchos::RCP<ContainerType> container = Teuchos::rcp(new ContainerType(getColMap(),getMap()));
112
113 return container;
114}
115
116template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
117Teuchos::RCP<LinearObjContainer>
120{
121 Teuchos::RCP<ContainerType> container = Teuchos::rcp(new ContainerType(getGhostedMap(),getGhostedMap()));
122
123 return container;
124}
125
126template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
127void
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);
136 ContainerType & t_out = Teuchos::dyn_cast<ContainerType>(out);
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
150template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
151void
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);
161 ContainerType & t_out = Teuchos::dyn_cast<ContainerType>(out);
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
175template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
176void
178ghostToGlobalTpetraVector(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
189template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
190void
192ghostToGlobalTpetraMatrix(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
206template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
207void
209globalToGhostTpetraVector(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
220template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
221void
224 const LinearObjContainer & globalBCRows,
225 LinearObjContainer & ghostedObjs,
226 bool zeroVectorRows, bool adjustX) const
227{
228 typedef Teuchos::ArrayRCP<const double>::Ordinal Ordinal;
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
298template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
299void
301applyDirichletBCs(const LinearObjContainer & /* counter */,
302 LinearObjContainer & /* result */) const
303{
304 TEUCHOS_ASSERT(false); // not yet implemented
305}
306
308//
309// buildReadOnlyDomainContainer()
310//
312template<typename Traits, typename ScalarT, typename LocalOrdinalT,
313 typename GlobalOrdinalT, typename NodeT>
314Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData>
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
328//
329// buildWriteDomainContainer()
330//
332template<typename Traits, typename ScalarT, typename LocalOrdinalT,
333 typename GlobalOrdinalT, typename NodeT>
334Teuchos::RCP<WriteVector_GlobalEvaluationData>
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
347template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
349getComm() const
350{
351 return *Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(getTeuchosComm());
352}
353
355template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
356Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
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
371template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
372Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> >
374getThyraRangeSpace() const
375{
376 if(rangeSpace_==Teuchos::null)
377 rangeSpace_ = Thyra::tpetraVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap());
378
379 return rangeSpace_;
380}
381
383template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
384Teuchos::RCP<Thyra::LinearOpBase<ScalarT> >
386getThyraMatrix() const
387{
388 return Thyra::tpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getThyraRangeSpace(),getThyraDomainSpace(),getTpetraMatrix());
389}
390
391// Functions for initalizing a container
393
394template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
395void
397initializeContainer(int mem,LinearObjContainer & loc) const
398{
399 ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
400 initializeContainer(mem,tloc);
401}
402
403template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
404void
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
425template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
426void
429{
430 ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
431 initializeGhostedContainer(mem,tloc);
432}
433
434template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
435void
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
464template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
465const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
467getMap() const
468{
469 if(map_==Teuchos::null) map_ = buildMap();
470
471 return map_;
472}
473
474// get the map from the matrix
475template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
476const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
478getColMap() const
479{
480 if(cMap_==Teuchos::null) cMap_ = buildColMap();
481
482 return cMap_;
483}
484
485template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
486const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
488getGhostedMap() const
489{
490 if(ghostedMap_==Teuchos::null) ghostedMap_ = buildGhostedMap();
491
492 return ghostedMap_;
493}
494
495template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
496const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
498getGhostedColMap() const
499{
500 if(cGhostedMap_==Teuchos::null) cGhostedMap_ = buildGhostedColMap();
501
502 return cGhostedMap_;
503}
504
505// get the graph of the crs matrix
506template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
507const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
509getGraph() const
510{
511 if(graph_==Teuchos::null) graph_ = buildGraph();
512
513 return graph_;
514}
515
516template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
517const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
519getGhostedGraph() const
520{
521 if(ghostedGraph_==Teuchos::null) ghostedGraph_ = buildGhostedGraph();
522
523 return ghostedGraph_;
524}
525
526template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
527const Teuchos::RCP<Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
529getGhostedImport() const
530{
531 if(ghostedImporter_==Teuchos::null)
532 ghostedImporter_ = Teuchos::rcp(new ImportType(getMap(),getGhostedMap()));
533
534 return ghostedImporter_;
535}
536
537template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
538const Teuchos::RCP<Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,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
551template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
552const Teuchos::RCP<Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
554getGhostedExport() const
555{
556 if(ghostedExporter_==Teuchos::null)
557 ghostedExporter_ = Teuchos::rcp(new ExportType(getGhostedMap(),getMap()));
558
559 return ghostedExporter_;
560}
561
562template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
563const Teuchos::RCP<Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,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
579template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
580const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
582buildMap() const
583{
584 std::vector<GlobalOrdinalT> indices;
585
586 // get the global indices
587 gidProvider_->getOwnedIndices(indices);
588
589 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
590}
591
592template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
593const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
595buildColMap() 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
605 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
606}
607
608// build the ghosted map
609template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
610const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
612buildGhostedMap() const
613{
614 std::vector<GlobalOrdinalT> indices;
615
616 // get the global indices
617 gidProvider_->getOwnedAndGhostedIndices(indices);
618
619 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
620}
621
622// build the ghosted map
623template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
624const Teuchos::RCP<Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
626buildGhostedColMap() const
627{
628 if(!hasColProvider_)
629 return buildGhostedMap();
630
631 std::vector<GlobalOrdinalT> indices;
632
633 // get the global indices
634 colGidProvider_->getOwnedAndGhostedIndices(indices);
635
636 return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
637}
638
639// get the graph of the crs matrix
640template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
641const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
643buildGraph() 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
663template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
664const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
666buildGhostedGraph() const
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
675 const Teuchos::RCP<const GlobalIndexer>
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);
718 Teuchos::RCP<CrsGraphType> graph = Teuchos::rcp(new CrsGraphType(rMap,cMap,
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
758template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
759Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
762{
763 Teuchos::RCP<const MapType> tMap = getGhostedMap();
764 return Teuchos::rcp(new VectorType(tMap));
765}
766
767template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
768Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
771{
772 Teuchos::RCP<const MapType> tMap = getGhostedColMap();
773 return Teuchos::rcp(new VectorType(tMap));
774}
775
776template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
777Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
779getTpetraVector() const
780{
781 Teuchos::RCP<const MapType> tMap = getMap();
782 return Teuchos::rcp(new VectorType(tMap));
783}
784
785template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
786Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
788getTpetraColVector() const
789{
790 Teuchos::RCP<const MapType> tMap = getColMap();
791 return Teuchos::rcp(new VectorType(tMap));
792}
793
794template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
795Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
797getTpetraMatrix() const
798{
799 Teuchos::RCP<CrsGraphType> tGraph = getGraph();
800 Teuchos::RCP<CrsMatrixType> tMat = Teuchos::rcp(new CrsMatrixType(tGraph));
801 tMat->fillComplete(tMat->getDomainMap(),tMat->getRangeMap());
802
803 return tMat;
804}
805
806template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
807Teuchos::RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
810{
811 Teuchos::RCP<CrsGraphType> tGraph = getGhostedGraph();
812 Teuchos::RCP<CrsMatrixType> tMat = Teuchos::rcp(new CrsMatrixType(tGraph));
813 tMat->fillComplete(tMat->getDomainMap(),tMat->getRangeMap());
814
815 return tMat;
816}
817
818template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
819const Teuchos::RCP<const Teuchos::Comm<int> >
825
826template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
828beginFill(LinearObjContainer & loc) const
829{
830 ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
831 Teuchos::RCP<CrsMatrixType> A = tloc.get_A();
832 if(A!=Teuchos::null)
833 A->resumeFill();
834}
835
836template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
838endFill(LinearObjContainer & loc) const
839{
840 ContainerType & tloc = Teuchos::dyn_cast<ContainerType>(loc);
841 Teuchos::RCP<CrsMatrixType> A = tloc.get_A();
842 if(A!=Teuchos::null)
843 A->fillComplete(A->getDomainMap(),A->getRangeMap());
844}
845
846}
847
848#endif // __Panzer_TpetraLinearObjFactory_impl_hpp__
This class provides a boundary exchange communication mechanism for vectors.
void buildGatherScatterEvaluators(const BuilderT &builder)
void set_dxdt(const Teuchos::RCP< VectorType > &in)
void set_A(const Teuchos::RCP< CrsMatrixType > &in)
void set_x(const Teuchos::RCP< VectorType > &in)
const Teuchos::RCP< VectorType > get_x() const
const Teuchos::RCP< VectorType > get_dxdt() const
const Teuchos::RCP< VectorType > get_f() const
void set_f(const Teuchos::RCP< VectorType > &in)
const Teuchos::RCP< CrsMatrixType > get_A() const
void ghostToGlobalTpetraMatrix(const Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out) const
virtual const Teuchos::RCP< Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedExport() const
get exporter for converting an overalapped object to a "normal" object
void initializeContainer(int, LinearObjContainer &loc) const
virtual void applyDirichletBCs(const LinearObjContainer &counter, LinearObjContainer &result) const
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const
virtual const Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const
get exporter for converting an overalapped object to a "normal" object
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraColVector() const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getColMap() const
TpetraLinearObjFactory(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< const GlobalIndexer > &gidProvider)
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedColMap() const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColMap() const
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGraph() const
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedGraph() const
Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > MapType
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraVector() const
virtual const Teuchos::RCP< Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColExport() const
virtual void beginFill(LinearObjContainer &loc) const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedMap() const
get the ghosted map from the matrix
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const
Tpetra::Export< LocalOrdinalT, GlobalOrdinalT, NodeT > ExportType
virtual Teuchos::RCP< ReadOnlyVector_GlobalEvaluationData > buildReadOnlyDomainContainer() const
virtual Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const
Get a matrix operator.
Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > CrsGraphType
Teuchos::RCP< Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedTpetraMatrix() const
void initializeGhostedContainer(int, LinearObjContainer &loc) const
virtual const Teuchos::RCP< Tpetra::Map< LocalOrdinalT, GlobalOrdinalT, NodeT > > buildGhostedMap() const
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 > > getGhostedTpetraVector() const
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 > > buildMap() const
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const
virtual void endFill(LinearObjContainer &loc) const
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const
Get the domain space.
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const
Teuchos::RCP< const GlobalIndexer > colGidProvider_
Teuchos::RCP< Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraColVector() const
virtual const Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedGraph() const
get the ghosted graph of the crs matrix
void globalToGhostTpetraVector(const Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &in, Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > &out, bool col) const
Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > CrsMatrixType
virtual const Teuchos::RCP< Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedColImport() const
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraRangeSpace() const
Get the range space.
virtual Teuchos::MpiComm< int > getComm() 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 > > buildColMap() const
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const
virtual const Teuchos::RCP< Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > > getGhostedImport() const
get importer for converting an overalapped object to a "normal" object
Tpetra::Vector< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > VectorType
Tpetra::Import< LocalOrdinalT, GlobalOrdinalT, NodeT > ImportType
Teuchos::RCP< Tpetra::CrsMatrix< ScalarT, LocalOrdinalT, GlobalOrdinalT, NodeT > > getTpetraMatrix() const