Panzer  Version of the Day
Panzer_STK_ParameterListCallbackBlocked_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 #ifdef PANZER_HAVE_TEKO
44 
45 using Teuchos::RCP;
46 using Teuchos::rcp;
47 
48 namespace panzer_stk {
49 
50 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
51 ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::ParameterListCallbackBlocked(
52  const Teuchos::RCP<const panzer_stk::STKConnManager<GlobalOrdinalT> > & connManager,
53  const Teuchos::RCP<const panzer::BlockedDOFManager<int,GlobalOrdinalT> > & blocked_ugi)
54  : connManager_(connManager), blocked_ugi_(blocked_ugi)
55 {
56 }
57 
58 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
59 Teuchos::RCP<Teuchos::ParameterList>
60 ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::request(const Teko::RequestMesg & rm)
61 {
62  TEUCHOS_ASSERT(handlesRequest(rm)); // design by contract
63 
64  // loop over parameter list and set the field by a particular key
65  Teuchos::RCP<Teuchos::ParameterList> outputPL = rcp(new Teuchos::ParameterList);
66  Teuchos::RCP<const Teuchos::ParameterList> inputPL = rm.getParameterList();
67  Teuchos::ParameterList::ConstIterator itr;
68  for(itr=inputPL->begin();itr!=inputPL->end();++itr) {
69  std::string * str_ptr = 0; // just used as a template specifier
70  std::string field = inputPL->entry(itr).getValue(str_ptr);
71  setFieldByKey(itr->first,field,*outputPL);
72  }
73 
74  return outputPL;
75 }
76 
77 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
78 bool ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::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 
108  return isHandled;
109  }
110  else return false;
111 }
112 
113 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
114 void ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::preRequest(const Teko::RequestMesg & rm)
115 {
116  TEUCHOS_ASSERT(handlesRequest(rm)); // design by contract
117 
118  const std::string & field = getHandledField(*rm.getParameterList());
119  int block = blocked_ugi_->getFieldBlock(blocked_ugi_->getFieldNum(field));
120 
121  // empty...nothing to do
122  buildArrayToVector(block,field);
123  buildCoordinates(field);
124 }
125 
126 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
127 void ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::setFieldByKey(const std::string & key,const std::string & field,Teuchos::ParameterList & pl) const
128 {
129  // x-, y-, z-coordinates needed for ML, Coordinates needed for MueLu
130  if(key=="x-coordinates") {
131  double * x = const_cast<double *>(&getCoordinateByField(0,field)[0]);
132  pl.set<double*>(key,x);
133  }
134  else if(key=="y-coordinates") {
135  double * y = const_cast<double *>(&getCoordinateByField(1,field)[0]);
136  pl.set<double*>(key,y);
137  }
138  else if(key=="z-coordinates") {
139  double * z = const_cast<double *>(&getCoordinateByField(2,field)[0]);
140  pl.set<double*>(key,z);
141  } else if(key == "Coordinates") {
142  pl.set<Teuchos::RCP<Tpetra::MultiVector<double,int,GlobalOrdinalT,Node> > >("Coordinates",coordsVec_);
143  }
144  else
145  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
146  "ParameterListCallback cannot handle key=\"" << key << "\"");
147 }
148 
149 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
150 void ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::buildArrayToVector(int block,const std::string & field)
151 {
152  if(arrayToVector_[field]==Teuchos::null) {
153  Teuchos::RCP<const panzer::UniqueGlobalIndexer<LocalOrdinalT,GlobalOrdinalT> > ugi = blocked_ugi_->getFieldDOFManagers()[block];
154  arrayToVector_[field] = Teuchos::rcp(new panzer::ArrayToFieldVector<LocalOrdinalT,GlobalOrdinalT,Node>(ugi));
155  }
156 }
157 
158 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
159 void ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::buildCoordinates(const std::string & field)
160 {
161  std::map<std::string,Kokkos::DynRankView<double,PHX::Device> > data;
162 
163  Teuchos::RCP<const panzer::Intrepid2FieldPattern> fieldPattern = getFieldPattern(field);
164 
165  std::vector<std::string> elementBlocks;
166  blocked_ugi_->getElementBlockIds(elementBlocks);
167  for(std::size_t i=0;i<elementBlocks.size();++i) {
168  std::string blockId = elementBlocks[i];
169  std::vector<std::size_t> localCellIds;
170 
171  // allocate block of data to store coordinates
172  Kokkos::DynRankView<double,PHX::Device> & fieldData = data[blockId];
173  fieldData = Kokkos::DynRankView<double,PHX::Device>("fieldData",connManager_->getElementBlock(blockId).size(),fieldPattern->numberIds());
174 
175  if(fieldPattern->supportsInterpolatoryCoordinates()) {
176  // get degree of freedom coordiantes
177  connManager_->getDofCoords(blockId,*fieldPattern,localCellIds,fieldData);
178  }
179  else {
180  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
181  out.setOutputToRootOnly(-1);
182  out << "WARNING: In ParameterListCallback::buildCoordinates(), the Intrepid2::FieldPattern in "
183  << "block \"" << blockId << "\" does not support interpolatory coordinates. "
184  << "This may be fine if coordinates are not actually needed. However if they are then bad things "
185  << "will happen. Enjoy!" << std::endl;
186 
187  return;
188  }
189  }
190 
191  coordsVec_ = arrayToVector_[field]->template getDataVector<double>(field,data);
192 
193  switch(coordsVec_->getNumVectors()) {
194  case 3:
195  zcoords_[field].resize(coordsVec_->getLocalLength());
196  coordsVec_->getVector(2)->get1dCopy(Teuchos::arrayViewFromVector(zcoords_[field]));
197  case 2:
198  ycoords_[field].resize(coordsVec_->getLocalLength());
199  coordsVec_->getVector(1)->get1dCopy(Teuchos::arrayViewFromVector(ycoords_[field]));
200  case 1:
201  xcoords_[field].resize(coordsVec_->getLocalLength());
202  coordsVec_->getVector(0)->get1dCopy(Teuchos::arrayViewFromVector(xcoords_[field]));
203  break;
204  default:
205  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
206  "ParameterListCallback::buildCoordinates: Constructed multivector has nonphysical dimensions.");
207  break;
208  }
209 }
210 
211 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
212 std::string ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::
213 getHandledField(const Teuchos::ParameterList & pl) const
214 {
215  // because this method assumes handlesRequest is true, this call will succeed
216  if(pl.isType<std::string>("x-coordinates"))
217  return pl.get<std::string>("x-coordinates");
218  else if(pl.isType<std::string>("Coordinates"))
219  return pl.get<std::string>("Coordinates");
220  else
221  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Neither x-coordinates not Coordinates field provided.");
222 
223 }
224 
225 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
226 const std::vector<double> & ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>::
227 getCoordinateByField(int dim,const std::string & field) const
228 {
229  TEUCHOS_ASSERT(dim>=0);
230  TEUCHOS_ASSERT(dim<=2);
231 
232  // get the coordinate vector you want
233  const std::map<std::string,std::vector<double> > * coord;
234  if(dim==0) coord = &xcoords_;
235  else if(dim==1) coord = &ycoords_;
236  else if(dim==2) coord = &zcoords_;
237 
238  std::map<std::string,std::vector<double> >::const_iterator itr;
239  itr = coord->find(field);
240 
241  TEUCHOS_TEST_FOR_EXCEPTION(itr==coord->end(),std::runtime_error,
242  "ParameterListCallbackBlocked::getCoordinateByField: Coordinates for field \"" + field +
243  "\" dimension " << dim << " have not been built!");
244 
245  return itr->second;
246 }
247 
248 template <typename LocalOrdinalT,typename GlobalOrdinalT,typename Node>
249 Teuchos::RCP<const panzer::Intrepid2FieldPattern> ParameterListCallbackBlocked<LocalOrdinalT,GlobalOrdinalT,Node>
250 ::getFieldPattern(const std::string & fieldName) const
251 {
252  std::vector<std::string> elementBlocks;
253  blocked_ugi_->getElementBlockIds(elementBlocks);
254 
255  for(std::size_t e=0;e<elementBlocks.size();e++) {
256  std::string blockId = elementBlocks[e];
257 
258  if(blocked_ugi_->fieldInBlock(fieldName,blockId))
259  return Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(blocked_ugi_->getFieldPattern(blockId,fieldName),true);
260  }
261 
262  return Teuchos::null;
263 }
264 
265 
266 }
267 
268 #endif
PHX::MDField< const ScalarT, Cell, IP > field