Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_ParameterListCallbackBlocked.cpp
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 "PanzerAdaptersSTK_config.hpp"
44#ifdef PANZER_HAVE_TEKO
45
47
48namespace panzer_stk {
49
50using Teuchos::RCP;
51using Teuchos::rcp;
52
53ParameterListCallbackBlocked::ParameterListCallbackBlocked(
54 const Teuchos::RCP<const panzer_stk::STKConnManager> & connManager,
55 const Teuchos::RCP<const panzer::BlockedDOFManager> & blocked_ugi,
56 const Teuchos::RCP<const panzer::BlockedDOFManager> & aux_blocked_ugi)
57 : connManager_(connManager), blocked_ugi_(blocked_ugi), aux_blocked_ugi_(aux_blocked_ugi)
58{}
59
60Teuchos::RCP<Teuchos::ParameterList>
61ParameterListCallbackBlocked::request(const Teko::RequestMesg & rm)
62{
63 TEUCHOS_ASSERT(handlesRequest(rm)); // design by contract
64
65 // loop over parameter list and set the field by a particular key
66 Teuchos::RCP<Teuchos::ParameterList> outputPL = rcp(new Teuchos::ParameterList);
67 Teuchos::RCP<const Teuchos::ParameterList> inputPL = rm.getParameterList();
68 Teuchos::ParameterList::ConstIterator itr;
69 for(itr=inputPL->begin();itr!=inputPL->end();++itr) {
70 std::string * str_ptr = 0; // just used as a template specifier
71 std::string field = inputPL->entry(itr).getValue(str_ptr);
72 setFieldByKey(itr->first,field,*outputPL);
73 }
74
75 return outputPL;
76}
77
78bool ParameterListCallbackBlocked::handlesRequest(const Teko::RequestMesg & rm)
79{
80 // check if is a parameter list message, and that the parameter
81 // list contains the right fields
82 if(rm.getName()=="Parameter List") {
83 bool isHandled = true;
84 Teuchos::RCP<const Teuchos::ParameterList> pl = rm.getParameterList();
85 std::string field;
86 if(pl->isType<std::string>("x-coordinates")) {
87 field = pl->get<std::string>("x-coordinates");
88 if(!isField(field)) {
89 return false;
90 }
91 }
92 if(pl->isType<std::string>("y-coordinates")) {
93 // we assume that the fields must be the same
94 if(field != pl->get<std::string>("y-coordinates")) {
95 return false;
96 }
97 }
98 if(pl->isType<std::string>("z-coordinates")) {
99 // we assume that the fields must be the same
100 if(field != pl->get<std::string>("z-coordinates")) {
101 return false;
102 }
103 }
104 if(pl->isType<std::string>("Coordinates")){
105 field = pl->get<std::string>("Coordinates");
106 }
107#ifdef PANZER_HAVE_EPETRA_STACK
108 if(pl->isType<std::string>("Coordinates-Epetra")){
109 field = pl->get<std::string>("Coordinates-Epetra");
110 }
111#endif
112
113 return isHandled;
114 }
115 else return false;
116}
117
118void ParameterListCallbackBlocked::preRequest(const Teko::RequestMesg & rm)
119{
120 TEUCHOS_ASSERT(handlesRequest(rm)); // design by contract
121
122 const std::string& field(getHandledField(*rm.getParameterList()));
123
124 // Check if the field is in the main UGI. If it's not, assume it's in the
125 // auxiliary UGI.
126 bool useAux(true);
127 std::vector<Teuchos::RCP<panzer::GlobalIndexer>>
128 fieldDOFMngrs = blocked_ugi_->getFieldDOFManagers();
129 for (int b(0); b < static_cast<int>(fieldDOFMngrs.size()); ++b)
130 {
131 for (int f(0); f < fieldDOFMngrs[b]->getNumFields(); ++f)
132 {
133 if (fieldDOFMngrs[b]->getFieldString(f) == field)
134 useAux = false;
135 }
136 }
137
138 int block(-1);
139 if (useAux)
140 block =
141 aux_blocked_ugi_->getFieldBlock(aux_blocked_ugi_->getFieldNum(field));
142 else
143 block = blocked_ugi_->getFieldBlock(blocked_ugi_->getFieldNum(field));
144
145 // Empty... Nothing to do.
146#ifdef PANZER_HAVE_EPETRA_STACK
147 if (rm.getParameterList()->isType<std::string>("Coordinates-Epetra")) {
148 buildArrayToVectorEpetra(block, field, useAux);
149 buildCoordinatesEpetra(field, useAux);
150 } else
151#endif
152 {
153 buildArrayToVectorTpetra(block, field, useAux);
154 buildCoordinatesTpetra(field, useAux);
155 }
156}
157
158void ParameterListCallbackBlocked::setFieldByKey(const std::string & key,const std::string & field,Teuchos::ParameterList & pl) const
159{
160 // x-, y-, z-coordinates needed for ML, Coordinates needed for MueLu
161 if(key=="x-coordinates") {
162 double * x = const_cast<double *>(&getCoordinateByField(0,field)[0]);
163 pl.set<double*>(key,x);
164 }
165 else if(key=="y-coordinates") {
166 double * y = const_cast<double *>(&getCoordinateByField(1,field)[0]);
167 pl.set<double*>(key,y);
168 }
169 else if(key=="z-coordinates") {
170 double * z = const_cast<double *>(&getCoordinateByField(2,field)[0]);
171 pl.set<double*>(key,z);
172 } else if(key == "Coordinates") {
173 pl.set<Teuchos::RCP<Tpetra::MultiVector<double,int,panzer::GlobalOrdinal,panzer::TpetraNodeType> > >(key,coordsVecTp_);
174#ifdef PANZER_HAVE_EPETRA_STACK
175 } else if(key == "Coordinates-Epetra") {
176 pl.set<Teuchos::RCP<Epetra_MultiVector> >("Coordinates",coordsVecEp_);
177 // pl.remove("Coordinates-Epetra");
178#endif
179 }
180 else
181 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
182 "ParameterListCallback cannot handle key=\"" << key << "\"");
183}
184
185void ParameterListCallbackBlocked::buildArrayToVectorTpetra(int block,const std::string & field, const bool useAux)
186{
187 if(arrayToVectorTpetra_[field]==Teuchos::null) {
188 Teuchos::RCP<const panzer::GlobalIndexer> ugi;
189 if(useAux)
190 ugi = aux_blocked_ugi_->getFieldDOFManagers()[block];
191 else
192 ugi = blocked_ugi_->getFieldDOFManagers()[block];
193 arrayToVectorTpetra_[field] = Teuchos::rcp(new panzer::ArrayToFieldVector(ugi));
194 }
195}
196
197#ifdef PANZER_HAVE_EPETRA_STACK
198void ParameterListCallbackBlocked::buildArrayToVectorEpetra(int block,const std::string & field, const bool useAux)
199{
200 if(arrayToVectorEpetra_[field]==Teuchos::null) {
201 Teuchos::RCP<const panzer::GlobalIndexer> ugi;
202 if(useAux)
203 ugi = aux_blocked_ugi_->getFieldDOFManagers()[block];
204 else
205 ugi = blocked_ugi_->getFieldDOFManagers()[block];
206 arrayToVectorEpetra_[field] = Teuchos::rcp(new panzer::ArrayToFieldVectorEpetra(ugi));
207 }
208}
209#endif
210
211void ParameterListCallbackBlocked::buildCoordinatesTpetra(const std::string & field, const bool useAux)
212{
213 std::map<std::string,Kokkos::DynRankView<double,PHX::Device> > data;
214
215 Teuchos::RCP<const panzer::Intrepid2FieldPattern> fieldPattern = getFieldPattern(field,useAux);
216
217 std::vector<std::string> elementBlocks;
218 if(useAux)
219 aux_blocked_ugi_->getElementBlockIds(elementBlocks);
220 else
221 blocked_ugi_->getElementBlockIds(elementBlocks);
222 for(std::size_t i=0;i<elementBlocks.size();++i) {
223 std::string blockId = elementBlocks[i];
224 std::vector<std::size_t> localCellIds;
225
226 // allocate block of data to store coordinates
227 Kokkos::DynRankView<double,PHX::Device> & fieldData = data[blockId];
228 fieldData = Kokkos::DynRankView<double,PHX::Device>("fieldData",connManager_->getElementBlock(blockId).size(),fieldPattern->numberIds());
229
230 if(fieldPattern->supportsInterpolatoryCoordinates()) {
231 // get degree of freedom coordiantes
232 connManager_->getDofCoords(blockId,*fieldPattern,localCellIds,fieldData);
233 }
234 else {
235 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
236 out.setOutputToRootOnly(-1);
237 out << "WARNING: In ParameterListCallback::buildCoordinates(), the Intrepid2::FieldPattern in "
238 << "block \"" << blockId << "\" does not support interpolatory coordinates. "
239 << "This may be fine if coordinates are not actually needed. However if they are then bad things "
240 << "will happen. Enjoy!" << std::endl;
241
242 return;
243 }
244 }
245
246 coordsVecTp_ = arrayToVectorTpetra_[field]->template getDataVector<double>(field,data);
247
248 switch(coordsVecTp_->getNumVectors()) {
249 case 3:
250 zcoords_[field].resize(coordsVecTp_->getLocalLength());
251 coordsVecTp_->getVector(2)->get1dCopy(Teuchos::arrayViewFromVector(zcoords_[field]));
252 // Intentional fall-through.
253 case 2:
254 ycoords_[field].resize(coordsVecTp_->getLocalLength());
255 coordsVecTp_->getVector(1)->get1dCopy(Teuchos::arrayViewFromVector(ycoords_[field]));
256 // Intentional fall-through.
257 case 1:
258 xcoords_[field].resize(coordsVecTp_->getLocalLength());
259 coordsVecTp_->getVector(0)->get1dCopy(Teuchos::arrayViewFromVector(xcoords_[field]));
260 break;
261 default:
262 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
263 "ParameterListCallback::buildCoordinates: Constructed multivector has nonphysical dimensions.");
264 break;
265 }
266}
267
268#ifdef PANZER_HAVE_EPETRA_STACK
269void ParameterListCallbackBlocked::buildCoordinatesEpetra(const std::string & field, const bool useAux)
270{
271 std::map<std::string,Kokkos::DynRankView<double,PHX::Device> > data;
272
273 Teuchos::RCP<const panzer::Intrepid2FieldPattern> fieldPattern = getFieldPattern(field,useAux);
274
275 std::vector<std::string> elementBlocks;
276 if(useAux)
277 aux_blocked_ugi_->getElementBlockIds(elementBlocks);
278 else
279 blocked_ugi_->getElementBlockIds(elementBlocks);
280 for(std::size_t i=0;i<elementBlocks.size();++i) {
281 std::string blockId = elementBlocks[i];
282 std::vector<std::size_t> localCellIds;
283
284 // allocate block of data to store coordinates
285 Kokkos::DynRankView<double,PHX::Device> & fieldData = data[blockId];
286 fieldData = Kokkos::DynRankView<double,PHX::Device>("fieldData",connManager_->getElementBlock(blockId).size(),fieldPattern->numberIds());
287
288 if(fieldPattern->supportsInterpolatoryCoordinates()) {
289 // get degree of freedom coordiantes
290 connManager_->getDofCoords(blockId,*fieldPattern,localCellIds,fieldData);
291 }
292 else {
293 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
294 out.setOutputToRootOnly(-1);
295 out << "WARNING: In ParameterListCallback::buildCoordinates(), the Intrepid2::FieldPattern in "
296 << "block \"" << blockId << "\" does not support interpolatory coordinates. "
297 << "This may be fine if coordinates are not actually needed. However if they are then bad things "
298 << "will happen. Enjoy!" << std::endl;
299
300 return;
301 }
302 }
303
304 coordsVecEp_ = arrayToVectorEpetra_[field]->template getDataVector<double>(field,data);
305}
306#endif
307
308std::string ParameterListCallbackBlocked::
309getHandledField(const Teuchos::ParameterList & pl) const
310{
311 // because this method assumes handlesRequest is true, this call will succeed
312 if(pl.isType<std::string>("x-coordinates"))
313 return pl.get<std::string>("x-coordinates");
314 else if(pl.isType<std::string>("Coordinates"))
315 return pl.get<std::string>("Coordinates");
316#ifdef PANZER_HAVE_EPETRA_STACK
317 else if(pl.isType<std::string>("Coordinates-Epetra"))
318 return pl.get<std::string>("Coordinates-Epetra");
319#endif
320 else
321#ifdef PANZER_HAVE_EPETRA_STACK
322 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Neither x-coordinates nor Coordinates or Coordinates-Epetra field provided.");
323#else
324 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Neither x-coordinates or Coordinates field provided.");
325#endif
326}
327
328const std::vector<double> & ParameterListCallbackBlocked::
329getCoordinateByField(int dim,const std::string & field) const
330{
331 TEUCHOS_ASSERT(dim>=0);
332 TEUCHOS_ASSERT(dim<=2);
333
334 // get the coordinate vector you want
335 const std::map<std::string,std::vector<double> > * coord = nullptr;
336 if(dim==0) coord = &xcoords_;
337 else if(dim==1) coord = &ycoords_;
338 else if(dim==2) coord = &zcoords_;
339 else {
340 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
341 "ParameterListCallbackBlocked::getCoordinateByField: dim is not in range of [0,2] for field \""
342 << field << "\".");
343 }
344
345 std::map<std::string,std::vector<double> >::const_iterator itr;
346 itr = coord->find(field);
347
348 TEUCHOS_TEST_FOR_EXCEPTION(itr==coord->end(),std::runtime_error,
349 "ParameterListCallbackBlocked::getCoordinateByField: Coordinates for field \"" + field +
350 "\" dimension " << dim << " have not been built!");
351
352 return itr->second;
353}
354
355Teuchos::RCP<const panzer::Intrepid2FieldPattern> ParameterListCallbackBlocked
356::getFieldPattern(const std::string & fieldName, const bool useAux) const
357{
358 std::vector<std::string> elementBlocks;
359 if(useAux)
360 aux_blocked_ugi_->getElementBlockIds(elementBlocks);
361 else
362 blocked_ugi_->getElementBlockIds(elementBlocks);
363
364 for(std::size_t e=0;e<elementBlocks.size();e++) {
365 std::string blockId = elementBlocks[e];
366
367 if(blocked_ugi_->fieldInBlock(fieldName,blockId))
368 return Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(blocked_ugi_->getFieldPattern(blockId,fieldName),true);
369
370 if(aux_blocked_ugi_->fieldInBlock(fieldName,blockId))
371 return Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(aux_blocked_ugi_->getFieldPattern(blockId,fieldName),true);
372 }
373
374 return Teuchos::null;
375}
376
377
378}
379
380#endif
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.