Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_BlockedTpetraLinearObjContainer_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Panzer: A partial differential equation assembly
4 // engine for strongly coupled complex multiphysics systems
5 //
6 // Copyright 2011 NTESS and the Panzer contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
11 #include "Thyra_VectorStdOps.hpp"
12 #include "Thyra_ProductVectorSpaceBase.hpp"
13 #include "Thyra_TpetraLinearOp.hpp"
14 
15 #include "Tpetra_CrsMatrix.hpp"
16 
17 namespace panzer {
18 
20 template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
21 bool BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
22 checkCompatibility() const
23 {
24  using Thyra::VectorSpaceBase;
25  using Teuchos::RCP;
26  using Teuchos::null;
27 
28  bool x_matches=false, f_matches=false, dxdt_matches=false;
29 
30  if(get_A()!=null) {
31  RCP<const VectorSpaceBase<ScalarT> > range = get_A()->range();
32  RCP<const VectorSpaceBase<ScalarT> > domain = get_A()->domain();
33 
34  if(get_x()!=null)
35  x_matches = range->isCompatible(*get_x()->space());
36  else
37  x_matches = true; // nothing to compare
38 
39  if(get_dxdt()!=null)
40  dxdt_matches = range->isCompatible(*get_dxdt()->space());
41  else
42  dxdt_matches = true; // nothing to compare
43 
44  if(get_f()!=null)
45  f_matches = range->isCompatible(*get_f()->space());
46  else
47  f_matches = true; // nothing to compare
48  }
49  else if(get_x()!=null && get_dxdt()!=null) {
50  f_matches = true; // nothing to compare f to
51  x_matches = get_x()->space()->isCompatible(*get_dxdt()->space()); // dxdt and x are in the same space
52  dxdt_matches = x_matches;
53  }
54  else {
55  f_matches = x_matches = dxdt_matches = true; // nothing to compare to
56  }
57 
58  return x_matches && dxdt_matches && f_matches;
59 }
60 
61 template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
62 void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
63 initialize()
64 {
65  using Thyra::LinearOpBase;
66  using Thyra::PhysicallyBlockedLinearOpBase;
67  using Thyra::ProductVectorSpaceBase;
68  using Teuchos::RCP;
69  using Teuchos::rcp_dynamic_cast;
70 
71  if(get_x()!=Teuchos::null) Thyra::assign<ScalarT>(x.ptr(),0.0);
72  if(get_dxdt()!=Teuchos::null) Thyra::assign<ScalarT>(get_dxdt().ptr(),0.0);
73  if(get_f()!=Teuchos::null) Thyra::assign<ScalarT>(get_f().ptr(),0.0);
74  if(get_A()!=Teuchos::null) {
75  RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
76  = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
77  RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
78  RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
79 
80  // loop over block entries
81  for(int i=0;i<range->numBlocks();i++) {
82  for(int j=0;j<domain->numBlocks();j++) {
83  RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
84  if(block!=Teuchos::null) {
85  RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
86  rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
87 
88  RCP<const MapType> map_i = t_block->getRangeMap();
89  RCP<const MapType> map_j = t_block->getDomainMap();
90 
91  RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
92  rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
93 
94  mat->resumeFill();
95  mat->setAllToScalar(0.0);
96  mat->fillComplete(map_j,map_i);
97  }
98  }
99  }
100  }
101 }
102 
103 template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
104 void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
105 initializeMatrix(ScalarT value)
106 {
107  using Thyra::LinearOpBase;
108  using Thyra::PhysicallyBlockedLinearOpBase;
109  using Thyra::ProductVectorSpaceBase;
110  using Teuchos::RCP;
111  using Teuchos::rcp_dynamic_cast;
112 
113  if(get_A()!=Teuchos::null) {
114  RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
115  = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
116  RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
117  RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
118 
119  // loop over block entries
120  for(int i=0;i<range->numBlocks();i++) {
121  for(int j=0;j<domain->numBlocks();j++) {
122  RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
123  if(block!=Teuchos::null) {
124  RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
125  rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
126 
127  // why do I have to do this?
128  RCP<const MapType> map_i = t_block->getRangeMap();
129  RCP<const MapType> map_j = t_block->getDomainMap();
130 
131  RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
132  rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
133 
134  mat->resumeFill();
135  mat->setAllToScalar(value);
136  mat->fillComplete(map_j,map_i);
137  }
138  }
139  }
140  }
141 }
142 
143 template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
144 void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
145 clear()
146 {
147  set_x(Teuchos::null);
148  set_dxdt(Teuchos::null);
149  set_f(Teuchos::null);
150  set_A(Teuchos::null);
151 }
152 
153 template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
154 void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
155 beginFill()
156 {
157  using Thyra::LinearOpBase;
158  using Thyra::PhysicallyBlockedLinearOpBase;
159  using Thyra::ProductVectorSpaceBase;
160  using Teuchos::RCP;
161  using Teuchos::rcp_dynamic_cast;
162 
163  if(get_A()!=Teuchos::null) {
164  RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
165  = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
166  RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
167  RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
168 
169  // loop over block entries
170  for(int i=0;i<range->numBlocks();i++) {
171  for(int j=0;j<domain->numBlocks();j++) {
172  RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
173  if(block!=Teuchos::null) {
174  RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
175  rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
176 
177  RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
178  rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
179 
180  mat->resumeFill();
181  }
182  }
183  }
184  }
185 }
186 
187 template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
188 void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
189 endFill()
190 {
191  using Thyra::LinearOpBase;
192  using Thyra::PhysicallyBlockedLinearOpBase;
193  using Thyra::ProductVectorSpaceBase;
194  using Teuchos::RCP;
195  using Teuchos::rcp_dynamic_cast;
196 
197  if(get_A()!=Teuchos::null) {
198  RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
199  = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
200  RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
201  RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
202 
203  // loop over block entries
204  for(int i=0;i<range->numBlocks();i++) {
205  for(int j=0;j<domain->numBlocks();j++) {
206  RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
207  if(block!=Teuchos::null) {
208  RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
209  rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
210 
211  RCP<const MapType> map_i = t_block->getRangeMap();
212  RCP<const MapType> map_j = t_block->getDomainMap();
213 
214  RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
215  rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
216 
217  mat->fillComplete(map_j,map_i);
218  }
219  }
220  }
221  }
222 }
223 
224 }
panzer::BlockedTpetraLinearObjContainer::checkCompatibility
bool checkCompatibility() const
Make sure row and column spaces match up.
Definition: Panzer_BlockedTpetraLinearObjContainer_impl.hpp:22
panzer::BlockedTpetraLinearObjContainer::initializeMatrix
void initializeMatrix(ScalarT value)
Put a particular scalar in the matrix.
Definition: Panzer_BlockedTpetraLinearObjContainer_impl.hpp:105
Generated by   doxygen 1.8.5