Panzer Version of the Day
Loading...
Searching...
No Matches
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
14namespace panzer_stk {
15
16namespace {
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
30template <typename EvalT>
31Teuchos::RCP<panzer::ResponseBase> ResponseEvaluatorFactory_SolutionWriter<EvalT>::
32buildResponseObject(const std::string & responseName) const
33{
34 return Teuchos::rcp(new Response_STKDummy(responseName));
35}
36
37template <typename EvalT>
39buildAndRegisterEvaluators(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_) {
79 Teuchos::RCP<const panzer::FieldLibraryBase> fieldLib = physicsBlock.getFieldLibraryBase();
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+"\"!");
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
168 Teuchos::RCP<PHX::Evaluator<panzer::Traits> > eval =
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 {
181 Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator
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++) {
189 Teuchos::ParameterList p;
190 p.set("Name",fields[f]);
191 p.set("Basis",basis);
192 p.set("Point Rule",centroidRule.getConst());
193 Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator
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 {
203 Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator
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 {
216 Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator
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++) {
224 Teuchos::ParameterList p;
225 p.set("Name",fields[f]);
226 p.set("Basis",basis);
227 p.set("Point Rule",centroidRule.getConst());
228 Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator
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 {
238 Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator
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
259template <typename EvalT>
261bucketByBasisType(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
273template <typename EvalT>
275computeReferenceCentroid(const std::map<std::string,Teuchos::RCP<const panzer::PureBasis> > & bases,
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 auto l_centroid = centroid;
284
285 // loop over each possible basis
286 for(std::map<std::string,RCP<const panzer::PureBasis> >::const_iterator itr=bases.begin();
287 itr!=bases.end();++itr) {
288
289 RCP<Intrepid2::Basis<PHX::exec_space,double,double>> intrepidBasis = itr->second->getIntrepid2Basis();
290
291 // we've got coordinates, lets commpute the "centroid"
292 Kokkos::DynRankView<double,PHX::Device> coords("coords",intrepidBasis->getCardinality(),
293 intrepidBasis->getBaseCellTopology().getDimension());
294 intrepidBasis->getDofCoords(coords);
295 TEUCHOS_ASSERT(coords.rank()==2);
296 TEUCHOS_ASSERT(coords.extent_int(1)==baseDimension);
297
298 Kokkos::parallel_for(coords.extent_int(0), KOKKOS_LAMBDA (int i) {
299 for(int d=0;d<coords.extent_int(1);d++)
300 l_centroid(0,d) += coords(i,d)/coords.extent(0);
301 });
302 Kokkos::fence();
303
304 return;
305 }
306
307 // no centroid was found...die
308 TEUCHOS_ASSERT(false);
309}
310
311template <typename EvalT>
313scaleField(const std::string & fieldName,double fieldScalar)
314{
315 fieldToScalar_[fieldName] = fieldScalar;
316 scaledFieldsHash_.insert(fieldName);
317}
318
319template <typename EvalT>
321typeSupported() const
322{
323 if(PHX::print<EvalT>()==PHX::print<panzer::Traits::Residual>())
324 return true;
325
326 return false;
327}
328
329template <typename EvalT>
331addAdditionalField(const std::string & fieldName,const Teuchos::RCP<const panzer::PureBasis> & basis)
332{
333 additionalFields_.push_back(std::make_pair(fieldName,basis));
334}
335
336template <typename EvalT>
338deleteRemovedFields(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
Interpolates basis DOF values to IP DOF values.
int baseCellDimension() const
Dimension of the base cell. NOT the dimension of the local side, even if the side() method returns tr...
Interpolates basis DOF values to IP DOF Curl values.
Object that contains information on the physics and discretization of a block of elements with the SA...
const std::vector< StrPureBasisPair > & getTangentFields() const
Returns list of tangent fields from DOFs and tangent param names.
Teuchos::RCP< const FieldLibraryBase > getFieldLibraryBase() 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.
const panzer::CellData & cellData() const
const std::vector< StrPureBasisPair > & getProvidedDOFs() const
const std::vector< std::vector< std::string > > & getCoordinateDOFs() const
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.
void addAdditionalField(const std::string &fieldName, const Teuchos::RCP< const panzer::PureBasis > &basis)
virtual void buildAndRegisterEvaluators(const std::string &responseName, PHX::FieldManager< panzer::Traits > &fm, const panzer::PhysicsBlock &physicsBlock, const Teuchos::ParameterList &user_data) const
void computeReferenceCentroid(const std::map< std::string, Teuchos::RCP< const panzer::PureBasis > > &bases, int baseDimension, Kokkos::DynRankView< double, PHX::Device > &centroid) 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)
virtual Teuchos::RCP< panzer::ResponseBase > buildResponseObject(const std::string &responseName) const