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
17namespace panzer {
18
20template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
21bool BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
22checkCompatibility() 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
61template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
62void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
63initialize()
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
103template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
104void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
105initializeMatrix(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
143template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
144void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
145clear()
146{
147 set_x(Teuchos::null);
148 set_dxdt(Teuchos::null);
149 set_f(Teuchos::null);
150 set_A(Teuchos::null);
151}
152
153template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
154void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
155beginFill()
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
187template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
188void BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
189endFill()
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.9.8