Panzer Version of the Day
Loading...
Searching...
No Matches
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#include "Tpetra_FECrsMatrix.hpp"
17
18namespace panzer {
19
21template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
22bool BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
23checkCompatibility() const
24{
25 using Thyra::VectorSpaceBase;
26 using Teuchos::RCP;
27 using Teuchos::null;
28
29 bool x_matches=false, f_matches=false, dxdt_matches=false;
30
31 if(get_A()!=null) {
32 RCP<const VectorSpaceBase<ScalarT> > range = get_A()->range();
33 RCP<const VectorSpaceBase<ScalarT> > domain = get_A()->domain();
34
35 if(get_x()!=null)
36 x_matches = range->isCompatible(*get_x()->space());
37 else
38 x_matches = true; // nothing to compare
39
40 if(get_dxdt()!=null)
41 dxdt_matches = range->isCompatible(*get_dxdt()->space());
42 else
43 dxdt_matches = true; // nothing to compare
44
45 if(get_f()!=null)
46 f_matches = range->isCompatible(*get_f()->space());
47 else
48 f_matches = true; // nothing to compare
49 }
50 else if(get_x()!=null && get_dxdt()!=null) {
51 f_matches = true; // nothing to compare f to
52 x_matches = get_x()->space()->isCompatible(*get_dxdt()->space()); // dxdt and x are in the same space
53 dxdt_matches = x_matches;
54 }
55 else {
56 f_matches = x_matches = dxdt_matches = true; // nothing to compare to
57 }
58
59 return x_matches && dxdt_matches && f_matches;
60}
61
62template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
63void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
64initialize()
65{
66 using Thyra::LinearOpBase;
67 using Thyra::PhysicallyBlockedLinearOpBase;
68 using Thyra::ProductVectorSpaceBase;
69 using Teuchos::RCP;
70 using Teuchos::rcp_dynamic_cast;
71
72 if(get_x()!=Teuchos::null) Thyra::assign<ScalarT>(x.ptr(),0.0);
73 if(get_dxdt()!=Teuchos::null) Thyra::assign<ScalarT>(get_dxdt().ptr(),0.0);
74 if(get_f()!=Teuchos::null) Thyra::assign<ScalarT>(get_f().ptr(),0.0);
75 if(get_A()!=Teuchos::null) {
76 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
77 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
78 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
79 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
80
81 // loop over block entries
82 for(int i=0;i<range->numBlocks();i++) {
83 for(int j=0;j<domain->numBlocks();j++) {
84 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
85 if(block!=Teuchos::null) {
86 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
87 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
88
89 RCP<const MapType> map_i = t_block->getRangeMap();
90 RCP<const MapType> map_j = t_block->getDomainMap();
91
92 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
93 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
94
95 // An FECrsMatrix block must not be wrapped in resumeFill/fillComplete: those
96 // are the plain inherited CrsMatrix calls and would corrupt its owned /
97 // owned+shared fill state. setAllToScalar alone is safe and does exactly what
98 // is wanted here (see also the ghost-row clearing in the factory's
99 // beginFill(), which covers the rows this cannot reach).
100 const bool isFE = Teuchos::nonnull(
101 rcp_dynamic_cast<Tpetra::FECrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(mat));
102
103 if(!isFE) mat->resumeFill();
104 mat->setAllToScalar(0.0);
105 if(!isFE) mat->fillComplete(map_j,map_i);
106 }
107 }
108 }
109 }
110}
111
112template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
113void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
114initializeMatrix(ScalarT value)
115{
116 using Thyra::LinearOpBase;
117 using Thyra::PhysicallyBlockedLinearOpBase;
118 using Thyra::ProductVectorSpaceBase;
119 using Teuchos::RCP;
120 using Teuchos::rcp_dynamic_cast;
121
122 if(get_A()!=Teuchos::null) {
123 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
124 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
125 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
126 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
127
128 // loop over block entries
129 for(int i=0;i<range->numBlocks();i++) {
130 for(int j=0;j<domain->numBlocks();j++) {
131 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
132 if(block!=Teuchos::null) {
133 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
134 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
135
136 // why do I have to do this?
137 RCP<const MapType> map_i = t_block->getRangeMap();
138 RCP<const MapType> map_j = t_block->getDomainMap();
139
140 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
141 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
142
143 // See initialize(): never resumeFill/fillComplete an FECrsMatrix block.
144 const bool isFE = Teuchos::nonnull(
145 rcp_dynamic_cast<Tpetra::FECrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(mat));
146
147 if(!isFE) mat->resumeFill();
148 mat->setAllToScalar(value);
149 if(!isFE) mat->fillComplete(map_j,map_i);
150 }
151 }
152 }
153 }
154}
155
156template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
157void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
158clear()
159{
160 set_x(Teuchos::null);
161 set_dxdt(Teuchos::null);
162 set_f(Teuchos::null);
163 set_A(Teuchos::null);
164}
165
166template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
167void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
168beginFill()
169{
170 using Thyra::LinearOpBase;
171 using Thyra::PhysicallyBlockedLinearOpBase;
172 using Thyra::ProductVectorSpaceBase;
173 using Teuchos::RCP;
174 using Teuchos::rcp_dynamic_cast;
175
176 if(get_A()!=Teuchos::null) {
177 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
178 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
179 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
180 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
181
182 // loop over block entries
183 for(int i=0;i<range->numBlocks();i++) {
184 for(int j=0;j<domain->numBlocks();j++) {
185 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
186 if(block!=Teuchos::null) {
187 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
188 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
189
190 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
191 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
192
193 // FE blocks are driven by BlockedTpetraLinearObjFactory::beginFill(), which
194 // calls beginAssembly(). resumeFill() here would bypass that state machine.
195 if(Teuchos::nonnull(rcp_dynamic_cast<Tpetra::FECrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(mat)))
196 continue;
197
198 mat->resumeFill();
199 }
200 }
201 }
202 }
203}
204
205template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
206void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
207endFill()
208{
209 using Thyra::LinearOpBase;
210 using Thyra::PhysicallyBlockedLinearOpBase;
211 using Thyra::ProductVectorSpaceBase;
212 using Teuchos::RCP;
213 using Teuchos::rcp_dynamic_cast;
214
215 if(get_A()!=Teuchos::null) {
216 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
217 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
218 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
219 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
220
221 // loop over block entries
222 for(int i=0;i<range->numBlocks();i++) {
223 for(int j=0;j<domain->numBlocks();j++) {
224 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
225 if(block!=Teuchos::null) {
226 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
227 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
228
229 RCP<const MapType> map_i = t_block->getRangeMap();
230 RCP<const MapType> map_j = t_block->getDomainMap();
231
232 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
233 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
234
235 // See beginFill(): FE blocks are closed by the factory via endAssembly().
236 if(Teuchos::nonnull(rcp_dynamic_cast<Tpetra::FECrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(mat)))
237 continue;
238
239 mat->fillComplete(map_j,map_i);
240 }
241 }
242 }
243 }
244}
245
246}
panzer::BlockedTpetraLinearObjContainer::initialize
virtual void initialize()
Zeroes out all linear objects currently held by this container.
Definition Panzer_BlockedTpetraLinearObjContainer_impl.hpp:64
panzer::BlockedTpetraLinearObjContainer::checkCompatibility
bool checkCompatibility() const
Make sure row and column spaces match up.
Definition Panzer_BlockedTpetraLinearObjContainer_impl.hpp:23
panzer::BlockedTpetraLinearObjContainer::initializeMatrix
void initializeMatrix(ScalarT value)
Put a particular scalar in the matrix.
Definition Panzer_BlockedTpetraLinearObjContainer_impl.hpp:114
Generated by doxygen 1.9.8