Panzer  Version of the Day
Panzer_STK_ResponseEvaluatorFactory_SolutionWriter_impl.hpp
Go to the documentation of this file.
1 #ifndef __Panzer_STK_ResponseEvaluatorFactory_SolutionWriter_impl_hpp__
2 #define __Panzer_STK_ResponseEvaluatorFactory_SolutionWriter_impl_hpp__
3 
5 #include "Panzer_STK_ScatterFields.hpp"
6 #include "Panzer_STK_ScatterVectorFields.hpp"
7 #include "Panzer_PointValues_Evaluator.hpp"
8 #include "Panzer_BasisValues_Evaluator.hpp"
9 #include "Panzer_DOF.hpp"
11 
12 #include <unordered_set>
13 
14 namespace panzer_stk {
15 
16 namespace {
18  class Response_STKDummy : public panzer::ResponseBase {
19  public:
20  Response_STKDummy(const std::string & rn)
21  : ResponseBase(rn) {}
22  virtual void scatterResponse() {}
23  virtual void initializeResponse() {}
24  private:
25  Response_STKDummy();
26  Response_STKDummy(const Response_STKDummy &);
27  };
28 }
29 
30 template <typename EvalT>
32 buildResponseObject(const std::string & responseName) const
33 {
34  return Teuchos::rcp(new Response_STKDummy(responseName));
35 }
36 
37 template <typename EvalT>
39 buildAndRegisterEvaluators(const std::string & responseName,
41  const panzer::PhysicsBlock & physicsBlock,
42  const Teuchos::ParameterList & user_data) const
43 {
44  using Teuchos::RCP;
45  using Teuchos::rcp;
46 
47  typedef std::pair<std::string,RCP<const panzer::PureBasis> > StrConstPureBasisPair;
48 
49  // this will help so we can print out any unused scaled fields as a warning
50  std::unordered_set<std::string> scaledFieldsHash = scaledFieldsHash_;
51 
52  std::map<std::string,RCP<const panzer::PureBasis> > bases;
53  std::map<std::string,std::vector<std::string> > basisBucket;
54  {
55  const std::map<std::string,RCP<panzer::PureBasis> > & nc_bases = physicsBlock.getBases();
56  bases.insert(nc_bases.begin(),nc_bases.end());
57  }
58 
59  std::vector<StrConstPureBasisPair> allFields;
60 
61  // only add in solution fields if required
62 
63  if(!addCoordinateFields_ && addSolutionFields_) {
64  // inject all the fields, including the coordinates (we will remove them shortly)
65  allFields.insert(allFields.end(),physicsBlock.getProvidedDOFs().begin(),physicsBlock.getProvidedDOFs().end());
66 
67 
68  // get a list of strings with fields to remove
69  std::vector<std::string> removedFields;
70  const std::vector<std::vector<std::string> > & coord_fields = physicsBlock.getCoordinateDOFs();
71  for(std::size_t c=0;c<coord_fields.size();c++)
72  for(std::size_t d=0;d<coord_fields[c].size();d++)
73  removedFields.push_back(coord_fields[c][d]);
74 
75  // remove all coordinate fields
76  deleteRemovedFields(removedFields,allFields);
77  }
78  else if(addCoordinateFields_ && !addSolutionFields_) {
80  const std::vector<std::vector<std::string> > & coord_fields = physicsBlock.getCoordinateDOFs();
81 
82  // get the basis and field for each coordiante
83  for(std::size_t c=0;c<coord_fields.size();c++) {
84  for(std::size_t d=0;d<coord_fields[c].size();d++) {
85  Teuchos::RCP<panzer::PureBasis> basis = // const_cast==yuck!
86  Teuchos::rcp_const_cast<panzer::PureBasis>(fieldLib->lookupBasis(coord_fields[c][d]));
87 
88  // make sure they are inserted in the allFields list
89  allFields.push_back(std::make_pair(coord_fields[c][d],basis));
90  }
91  }
92  }
93  else if(addSolutionFields_)
94  allFields.insert(allFields.end(),physicsBlock.getProvidedDOFs().begin(),physicsBlock.getProvidedDOFs().end());
95 
96  // Add in tangent fields
97  if(addSolutionFields_)
98  allFields.insert(allFields.end(),physicsBlock.getTangentFields().begin(),physicsBlock.getTangentFields().end());
99 
100  // add in bases for any addtional fields
101  for(std::size_t i=0;i<additionalFields_.size();i++)
102  bases[additionalFields_[i].second->name()] = additionalFields_[i].second;
103 
104  allFields.insert(allFields.end(),additionalFields_.begin(),additionalFields_.end());
105 
106  deleteRemovedFields(removedFields_,allFields);
107 
108  bucketByBasisType(allFields,basisBucket);
109 
110  // add this for HCURL and HDIV basis, only want to add them once: evaluate vector fields at centroid
112  RCP<panzer::PointRule> centroidRule;
113  for(std::map<std::string,Teuchos::RCP<const panzer::PureBasis> >::const_iterator itr=bases.begin();
114  itr!=bases.end();++itr) {
115 
116  if(itr->second->isVectorBasis()) {
117  centroidRule = rcp(new panzer::PointRule("Centroid",1,physicsBlock.cellData()));
118 
119  // compute centroid
120  Kokkos::DynRankView<double,PHX::Device> centroid;
121  computeReferenceCentroid(bases,physicsBlock.cellData().baseCellDimension(),centroid);
122 
123  // build pointe values evaluator
124  RCP<PHX::Evaluator<panzer::Traits> > evaluator =
125  rcp(new panzer::PointValues_Evaluator<EvalT,panzer::Traits>(centroidRule,centroid));
126  this->template registerEvaluator<EvalT>(fm, evaluator);
127 
128  break; // get out of the loop, only need one evaluator
129  }
130  }
131 
132  // add evaluators for each field
134 
135  for(std::map<std::string,std::vector<std::string> >::const_iterator itr=basisBucket.begin();
136  itr!=basisBucket.end();++itr) {
137 
138  std::string basisName = itr->first;
139  const std::vector<std::string> & fields = itr->second;
140 
141  std::map<std::string,Teuchos::RCP<const panzer::PureBasis> >::const_iterator found = bases.find(basisName);
142  TEUCHOS_TEST_FOR_EXCEPTION(found==bases.end(),std::logic_error,
143  "Could not find basis \""+basisName+"\"!");
145 
146  // determine if user has modified field scalar for each field to be written to STK
147  std::vector<double> scalars(fields.size(),1.0); // fill with 1.0
148  for(std::size_t f=0;f<fields.size();f++) {
149  std::unordered_map<std::string,double>::const_iterator f2s_itr = fieldToScalar_.find(fields[f]);
150 
151  // if scalar is found, include it in the vector and remove the field from the
152  // hash table so it won't be included in the warning message.
153  if(f2s_itr!=fieldToScalar_.end()) {
154  scalars[f] = f2s_itr->second;
155  scaledFieldsHash.erase(fields[f]);
156  }
157  }
158 
159  // write out nodal fields
162 
163  std::string fields_concat = "";
164  for(std::size_t f=0;f<fields.size();f++) {
165  fields_concat += fields[f];
166  }
167 
169  Teuchos::rcp(new ScatterFields<EvalT,panzer::Traits>("STK HGRAD Scatter Basis " +basis->name()+": "+fields_concat,
170  mesh_, basis, fields,scalars));
171 
172  // register and require evaluator fields
173  this->template registerEvaluator<EvalT>(fm, eval);
174  fm.template requireField<EvalT>(*eval->evaluatedFields()[0]);
175  }
177  TEUCHOS_ASSERT(centroidRule!=Teuchos::null);
178 
179  // register basis values evaluator
180  {
182  = Teuchos::rcp(new panzer::BasisValues_Evaluator<EvalT,panzer::Traits>(centroidRule,basis));
183  this->template registerEvaluator<EvalT>(fm, evaluator);
184  }
185 
186  // add a DOF_PointValues for each field
187  std::string fields_concat = "";
188  for(std::size_t f=0;f<fields.size();f++) {
190  p.set("Name",fields[f]);
191  p.set("Basis",basis);
192  p.set("Point Rule",centroidRule.getConst());
195 
196  this->template registerEvaluator<EvalT>(fm, evaluator);
197 
198  fields_concat += fields[f];
199  }
200 
201  // add the scatter field evaluator for this basis
202  {
204  = Teuchos::rcp(new ScatterVectorFields<EvalT,panzer::Traits>("STK HCURL Scatter Basis " +basis->name()+": "+fields_concat,
205  mesh_,centroidRule,fields,scalars));
206 
207  this->template registerEvaluator<EvalT>(fm, evaluator);
208  fm.template requireField<EvalT>(*evaluator->evaluatedFields()[0]); // require the dummy evaluator
209  }
210  }
212  TEUCHOS_ASSERT(centroidRule!=Teuchos::null);
213 
214  // register basis values evaluator
215  {
217  = Teuchos::rcp(new panzer::BasisValues_Evaluator<EvalT,panzer::Traits>(centroidRule,basis));
218  this->template registerEvaluator<EvalT>(fm, evaluator);
219  }
220 
221  // add a DOF_PointValues for each field
222  std::string fields_concat = "";
223  for(std::size_t f=0;f<fields.size();f++) {
225  p.set("Name",fields[f]);
226  p.set("Basis",basis);
227  p.set("Point Rule",centroidRule.getConst());
230 
231  this->template registerEvaluator<EvalT>(fm, evaluator);
232 
233  fields_concat += fields[f];
234  }
235 
236  // add the scatter field evaluator for this basis
237  {
239  = Teuchos::rcp(new ScatterVectorFields<EvalT,panzer::Traits>("STK HDIV Scatter Basis " +basis->name()+": "+fields_concat,
240  mesh_,centroidRule,fields,scalars));
241 
242  this->template registerEvaluator<EvalT>(fm, evaluator);
243  fm.template requireField<EvalT>(*evaluator->evaluatedFields()[0]); // require the dummy evaluator
244  }
245  }
246  }
247 
248  // print warning message for any unused scaled fields
249  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
250  out.setOutputToRootOnly(0);
251 
252  for(std::unordered_set<std::string>::const_iterator itr=scaledFieldsHash.begin();
253  itr!=scaledFieldsHash.end();itr++) {
254  out << "WARNING: STK Solution Writer did not scale the field \"" << *itr << "\" "
255  << "because it was not written." << std::endl;
256  }
257 }
258 
259 template <typename EvalT>
261 bucketByBasisType(const std::vector<std::pair<std::string,Teuchos::RCP<const panzer::PureBasis> > > & providedDofs,
262  std::map<std::string,std::vector<std::string> > & basisBucket)
263 {
264  // this should be self explanatory
265  for(std::size_t i=0;i<providedDofs.size();i++) {
266  std::string fieldName = providedDofs[i].first;
267  Teuchos::RCP<const panzer::PureBasis> basis = providedDofs[i].second;
268 
269  basisBucket[basis->name()].push_back(fieldName);
270  }
271 }
272 
273 template <typename EvalT>
276  int baseDimension,
277  Kokkos::DynRankView<double,PHX::Device> & centroid) const
278 {
279  using Teuchos::RCP;
280  using Teuchos::rcp_dynamic_cast;
281 
282  centroid = Kokkos::DynRankView<double,PHX::Device>("centroid",1,baseDimension);
283 
284  // loop over each possible basis
285  for(std::map<std::string,RCP<const panzer::PureBasis> >::const_iterator itr=bases.begin();
286  itr!=bases.end();++itr) {
287 
288  // see if this basis has coordinates
289  RCP<Intrepid2::Basis<double,Kokkos::DynRankView<double,PHX::Device> > > intrepidBasis = itr->second->getIntrepid2Basis();
290  RCP<Intrepid2::DofCoordsInterface<Kokkos::DynRankView<double,PHX::Device> > > basisCoords
291  = rcp_dynamic_cast<Intrepid2::DofCoordsInterface<Kokkos::DynRankView<double,PHX::Device> > >(intrepidBasis);
292 
293  if(basisCoords==Teuchos::null) // no coordinates...move on
294  continue;
295 
296  // we've got coordinates, lets commpute the "centroid"
297  Kokkos::DynRankView<double,PHX::Device> coords("coords",intrepidBasis->getCardinality(),
298  intrepidBasis->getBaseCellTopology().getDimension());
299  basisCoords->getDofCoords(coords);
300  TEUCHOS_ASSERT(coords.rank()==2);
301  TEUCHOS_ASSERT(coords.extent_int(1)==baseDimension);
302 
303  for(int i=0;i<coords.extent_int(0);i++)
304  for(int d=0;d<coords.extent_int(1);d++)
305  centroid(0,d) += coords(i,d);
306 
307  // take the average
308  for(int d=0;d<coords.extent_int(1);d++)
309  centroid(0,d) /= coords.dimension(0);
310 
311  return;
312  }
313 
314  // no centroid was found...die
315  TEUCHOS_ASSERT(false);
316 }
317 
318 template <typename EvalT>
320 scaleField(const std::string & fieldName,double fieldScalar)
321 {
322  fieldToScalar_[fieldName] = fieldScalar;
323 }
324 
325 template <typename EvalT>
328 {
329  if(PHX::typeAsString<EvalT>()==PHX::typeAsString<panzer::Traits::Residual>())
330  return true;
331 
332  return false;
333 }
334 
335 template <typename EvalT>
337 addAdditionalField(const std::string & fieldName,const Teuchos::RCP<const panzer::PureBasis> & basis)
338 {
339  additionalFields_.push_back(std::make_pair(fieldName,basis));
340 }
341 
342 template <typename EvalT>
344 deleteRemovedFields(const std::vector<std::string> & removedFields,
345  std::vector<std::pair<std::string,Teuchos::RCP<const panzer::PureBasis> > > & fields) const
346 {
348  functor.removedFields_ = removedFields;
349 
350  // This is the Erase-Remove Idiom: see http://en.wikipedia.org/wiki/Erase-remove_idiom
351  fields.erase(std::remove_if(fields.begin(),fields.end(),functor),fields.end());
352 }
353 
354 }
355 
356 #endif
void computeReferenceCentroid(const std::map< std::string, Teuchos::RCP< const panzer::PureBasis > > &bases, int baseDimension, Kokkos::DynRankView< double, PHX::Device > &centroid) const
std::string name() const
A unique key that is the combination of the basis type and basis order.
Object that contains information on the physics and discretization of a block of elements with the SA...
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Kokkos::View< const double *, PHX::Device > scalars
virtual Teuchos::RCP< panzer::ResponseBase > buildResponseObject(const std::string &responseName) const
virtual void buildAndRegisterEvaluators(const std::string &responseName, PHX::FieldManager< panzer::Traits > &fm, const panzer::PhysicsBlock &physicsBlock, const Teuchos::ParameterList &user_data) const
PHX::MDField< ScalarT > vector
void addAdditionalField(const std::string &fieldName, const Teuchos::RCP< const panzer::PureBasis > &basis)
Interpolates basis DOF values to IP DOF Curl values.
int baseCellDimension() const
Dimension of the base cell. NOT the dimension of the local side, even if the side() method returns tr...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Teuchos::RCP< const panzer::PureBasis > lookupBasis(const std::string &fieldName) const =0
Get the basis associated with a particular field.
Teuchos::RCP< STK_Interface > mesh_
EElementSpace getElementSpace() const
Teuchos::RCP< const FieldLibraryBase > getFieldLibraryBase() const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
const std::vector< std::vector< std::string > > & getCoordinateDOFs() const
const panzer::CellData & cellData() const
static void bucketByBasisType(const std::vector< std::pair< std::string, Teuchos::RCP< const panzer::PureBasis > > > &providedDofs, std::map< std::string, std::vector< std::string > > &basisBucket)
const std::vector< StrPureBasisPair > & getProvidedDOFs() const
const std::map< std::string, Teuchos::RCP< panzer::PureBasis > > & getBases() const
Returns the unique set of bases, key is the unique panzer::PureBasis::name() of the basis...
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
void deleteRemovedFields(const std::vector< std::string > &removedFields, std::vector< std::pair< std::string, Teuchos::RCP< const panzer::PureBasis > > > &fields) const
Delete from the argument all the fields that are in the removedFields array.
Description and data layouts associated with a particular basis.
const std::vector< StrPureBasisPair > & getTangentFields() const
Returns list of tangent fields from DOFs and tangent param names.
#define TEUCHOS_ASSERT(assertion_test)