Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_BlockedTpetraLinearObjContainer_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#include "Thyra_VectorStdOps.hpp"
44#include "Thyra_ProductVectorSpaceBase.hpp"
45#include "Thyra_TpetraLinearOp.hpp"
46
47#include "Tpetra_CrsMatrix.hpp"
48
49namespace panzer {
50
52template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
55{
57 using Teuchos::RCP;
58 using Teuchos::null;
59
60 bool x_matches=false, f_matches=false, dxdt_matches=false;
61
62 if(get_A()!=null) {
63 RCP<const VectorSpaceBase<ScalarT> > range = get_A()->range();
64 RCP<const VectorSpaceBase<ScalarT> > domain = get_A()->domain();
65
66 if(get_x()!=null)
67 x_matches = range->isCompatible(*get_x()->space());
68 else
69 x_matches = true; // nothing to compare
70
71 if(get_dxdt()!=null)
72 dxdt_matches = range->isCompatible(*get_dxdt()->space());
73 else
74 dxdt_matches = true; // nothing to compare
75
76 if(get_f()!=null)
77 f_matches = range->isCompatible(*get_f()->space());
78 else
79 f_matches = true; // nothing to compare
80 }
81 else if(get_x()!=null && get_dxdt()!=null) {
82 f_matches = true; // nothing to compare f to
83 x_matches = get_x()->space()->isCompatible(*get_dxdt()->space()); // dxdt and x are in the same space
84 dxdt_matches = x_matches;
85 }
86 else {
87 f_matches = x_matches = dxdt_matches = true; // nothing to compare to
88 }
89
90 return x_matches && dxdt_matches && f_matches;
91}
92
93template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
96{
98 using Thyra::PhysicallyBlockedLinearOpBase;
99 using Thyra::ProductVectorSpaceBase;
100 using Teuchos::RCP;
101 using Teuchos::rcp_dynamic_cast;
102
103 if(get_x()!=Teuchos::null) Thyra::assign<ScalarT>(x.ptr(),0.0);
104 if(get_dxdt()!=Teuchos::null) Thyra::assign<ScalarT>(get_dxdt().ptr(),0.0);
105 if(get_f()!=Teuchos::null) Thyra::assign<ScalarT>(get_f().ptr(),0.0);
106 if(get_A()!=Teuchos::null) {
107 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
108 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
109 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
110 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
111
112 // loop over block entries
113 for(int i=0;i<range->numBlocks();i++) {
114 for(int j=0;j<domain->numBlocks();j++) {
115 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
116 if(block!=Teuchos::null) {
117 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
118 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
119
120 RCP<const MapType> map_i = t_block->getRangeMap();
121 RCP<const MapType> map_j = t_block->getDomainMap();
122
123 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
124 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
125
126 mat->resumeFill();
127 mat->setAllToScalar(0.0);
128 mat->fillComplete(map_j,map_i);
129 }
130 }
131 }
132 }
133}
134
135template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
137initializeMatrix(ScalarT value)
138{
140 using Thyra::PhysicallyBlockedLinearOpBase;
141 using Thyra::ProductVectorSpaceBase;
142 using Teuchos::RCP;
143 using Teuchos::rcp_dynamic_cast;
144
145 if(get_A()!=Teuchos::null) {
146 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
147 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
148 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
149 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
150
151 // loop over block entries
152 for(int i=0;i<range->numBlocks();i++) {
153 for(int j=0;j<domain->numBlocks();j++) {
154 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
155 if(block!=Teuchos::null) {
156 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
157 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
158
159 // why do I have to do this?
160 RCP<const MapType> map_i = t_block->getRangeMap();
161 RCP<const MapType> map_j = t_block->getDomainMap();
162
163 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
164 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
165
166 mat->resumeFill();
167 mat->setAllToScalar(value);
168 mat->fillComplete(map_j,map_i);
169 }
170 }
171 }
172 }
173}
174
175template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
177clear()
178{
179 set_x(Teuchos::null);
180 set_dxdt(Teuchos::null);
181 set_f(Teuchos::null);
182 set_A(Teuchos::null);
183}
184
185template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
187beginFill()
188{
190 using Thyra::PhysicallyBlockedLinearOpBase;
191 using Thyra::ProductVectorSpaceBase;
192 using Teuchos::RCP;
193 using Teuchos::rcp_dynamic_cast;
194
195 if(get_A()!=Teuchos::null) {
196 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
197 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
198 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
199 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
200
201 // loop over block entries
202 for(int i=0;i<range->numBlocks();i++) {
203 for(int j=0;j<domain->numBlocks();j++) {
204 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
205 if(block!=Teuchos::null) {
206 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
207 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
208
209 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
210 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
211
212 mat->resumeFill();
213 }
214 }
215 }
216 }
217}
218
219template <typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
221endFill()
222{
224 using Thyra::PhysicallyBlockedLinearOpBase;
225 using Thyra::ProductVectorSpaceBase;
226 using Teuchos::RCP;
227 using Teuchos::rcp_dynamic_cast;
228
229 if(get_A()!=Teuchos::null) {
230 RCP<PhysicallyBlockedLinearOpBase<ScalarT> > Amat
231 = rcp_dynamic_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(get_A(),true);
232 RCP<const ProductVectorSpaceBase<ScalarT> > range = Amat->productRange();
233 RCP<const ProductVectorSpaceBase<ScalarT> > domain = Amat->productDomain();
234
235 // loop over block entries
236 for(int i=0;i<range->numBlocks();i++) {
237 for(int j=0;j<domain->numBlocks();j++) {
238 RCP<LinearOpBase<ScalarT> > block = Amat->getNonconstBlock(i,j);
239 if(block!=Teuchos::null) {
240 RCP<Tpetra::Operator<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > t_block =
241 rcp_dynamic_cast<Thyra::TpetraLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(block,true)->getTpetraOperator();
242
243 RCP<const MapType> map_i = t_block->getRangeMap();
244 RCP<const MapType> map_j = t_block->getDomainMap();
245
246 RCP<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> > mat =
247 rcp_dynamic_cast<Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >(t_block,true);
248
249 mat->fillComplete(map_j,map_i);
250 }
251 }
252 }
253 }
254}
255
256}
bool checkCompatibility() const
Make sure row and column spaces match up.
void initializeMatrix(ScalarT value)
Put a particular scalar in the matrix.