Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
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+"\"!");
144  Teuchos::RCP<const panzer::PureBasis> basis = found->second;
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
160  if(basis->getElementSpace()==panzer::PureBasis::HGRAD ||
161  basis->getElementSpace()==panzer::PureBasis::CONST) {
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  }
176  else if(basis->getElementSpace()==panzer::PureBasis::HCURL) {
177  TEUCHOS_ASSERT(centroidRule!=Teuchos::null);
178 
179  // register basis values evaluator
180  {
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  }
211  else if(basis->getElementSpace()==panzer::PureBasis::HDIV) {
212  TEUCHOS_ASSERT(centroidRule!=Teuchos::null);
213 
214  // register basis values evaluator
215  {
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  RCP<Intrepid2::Basis<PHX::exec_space,double,double>> intrepidBasis = itr->second->getIntrepid2Basis();
289 
290  // we've got coordinates, lets commpute the "centroid"
291  Kokkos::DynRankView<double,PHX::Device> coords("coords",intrepidBasis->getCardinality(),
292  intrepidBasis->getBaseCellTopology().getDimension());
293  intrepidBasis->getDofCoords(coords);
294  TEUCHOS_ASSERT(coords.rank()==2);
295  TEUCHOS_ASSERT(coords.extent_int(1)==baseDimension);
296 
297  for(int i=0;i<coords.extent_int(0);i++)
298  for(int d=0;d<coords.extent_int(1);d++)
299  centroid(0,d) += coords(i,d);
300 
301  // take the average
302  for(int d=0;d<coords.extent_int(1);d++)
303  centroid(0,d) /= coords.extent(0);
304 
305  return;
306  }
307 
308  // no centroid was found...die
309  TEUCHOS_ASSERT(false);
310 }
311 
312 template <typename EvalT>
314 scaleField(const std::string & fieldName,double fieldScalar)
315 {
316  fieldToScalar_[fieldName] = fieldScalar;
317 }
318 
319 template <typename EvalT>
322 {
323  if(PHX::print<EvalT>()==PHX::print<panzer::Traits::Residual>())
324  return true;
325 
326  return false;
327 }
328 
329 template <typename EvalT>
331 addAdditionalField(const std::string & fieldName,const Teuchos::RCP<const panzer::PureBasis> & basis)
332 {
333  additionalFields_.push_back(std::make_pair(fieldName,basis));
334 }
335 
336 template <typename EvalT>
338 deleteRemovedFields(const std::vector<std::string> & removedFields,
339  std::vector<std::pair<std::string,Teuchos::RCP<const panzer::PureBasis> > > & fields) const
340 {
342  functor.removedFields_ = removedFields;
343 
344  // This is the Erase-Remove Idiom: see http://en.wikipedia.org/wiki/Erase-remove_idiom
345  fields.erase(std::remove_if(fields.begin(),fields.end(),functor),fields.end());
346 }
347 
348 }
349 
350 #endif
std::string name() const
A unique key that is the combination of the basis type and basis order.
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.
RCP< const T > getConst() const
int baseCellDimension() const
Dimension of the base cell. NOT the dimension of the local side, even if the side() method returns tr...
Object that contains information on the physics and discretization of a block of elements with the SA...
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...
const panzer::CellData & cellData() const
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)
virtual Teuchos::RCP< panzer::ResponseBase > buildResponseObject(const std::string &responseName) const
Interpolates basis DOF values to IP DOF values.
void addAdditionalField(const std::string &fieldName, const Teuchos::RCP< const panzer::PureBasis > &basis)
Interpolates basis DOF values to IP DOF Curl values.
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.
const std::vector< std::vector< std::string > > & getCoordinateDOFs() const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
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 > & getTangentFields() const
Returns list of tangent fields from DOFs and tangent param names.
Teuchos::RCP< const FieldLibraryBase > getFieldLibraryBase() const
virtual void buildAndRegisterEvaluators(const std::string &responseName, PHX::FieldManager< panzer::Traits > &fm, const panzer::PhysicsBlock &physicsBlock, const Teuchos::ParameterList &user_data) const
Description and data layouts associated with a particular basis.
#define TEUCHOS_ASSERT(assertion_test)
Interpolates basis DOF values to IP DOF values.
const std::vector< StrPureBasisPair > & getProvidedDOFs() const
void computeReferenceCentroid(const std::map< std::string, Teuchos::RCP< const panzer::PureBasis > > &bases, int baseDimension, Kokkos::DynRankView< double, PHX::Device > &centroid) const