Panzer  Version of the Day
Panzer_ResponseScatterEvaluator_Probe_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 #ifndef PANZER_RESPONSE_SCATTER_EVALUATOR_EXTREMEVALUE_IMPL_HPP
44 #define PANZER_RESPONSE_SCATTER_EVALUATOR_EXTREMEVALUE_IMPL_HPP
45 
46 #include <iostream>
47 #include <string>
48 
49 #include "PanzerDiscFE_config.hpp"
50 
51 #include "Phalanx_Evaluator_Macros.hpp"
52 #include "Phalanx_MDField.hpp"
53 #include "Phalanx_DataLayout_MDALayout.hpp"
54 
55 #include "Panzer_CellData.hpp"
56 #include "Panzer_PointRule.hpp"
58 #include "Panzer_ResponseBase.hpp"
59 #include "Panzer_Dimension.hpp"
60 
61 #include "Intrepid2_FunctionSpaceTools.hpp"
62 
63 #include "Thyra_SpmdVectorBase.hpp"
64 #include "Teuchos_ArrayRCP.hpp"
65 
66 #include "Kokkos_ViewFactory.hpp"
67 
68 namespace panzer {
69 
70 template<typename EvalT, typename Traits, typename LO, typename GO>
73  const std::string & responseName,
74  const std::string & fieldName,
75  const int fieldComponent,
76  const Teuchos::Array<double>& point,
77  const IntegrationRule & ir,
78  const Teuchos::RCP<const PureBasis>& basis,
79  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & indexer,
80  const Teuchos::RCP<ProbeScatterBase> & probeScatter)
81  : responseName_(responseName)
82  , fieldName_(fieldName)
83  , fieldComponent_(fieldComponent)
84  , point_(point)
85  , basis_(basis)
86  , topology_(ir.topology)
87  , globalIndexer_(indexer)
88  , scatterObj_(probeScatter)
89  , cellIndex_(0)
90 {
91  using Teuchos::RCP;
92  using Teuchos::rcp;
93 
94  // the field manager will allocate all of these fields
95  field_ = PHX::MDField<const ScalarT,Cell,BASIS>(fieldName,basis_->functional);
96  this->addDependentField(field_);
97 
98  num_basis = basis->cardinality();
99  num_dim = basis->dimension();
100  TEUCHOS_ASSERT(num_dim == static_cast<size_t>(point_.size()));
101 
102  basis_values_ = Kokkos::DynRankView<double,PHX::Device>(
103  "basis_values", 1, num_basis, 1); // Cell, Basis, Point
104 
105  // build dummy target tag
106  std::string dummyName =
107  ResponseBase::buildLookupName(responseName) + " dummy target";
108  RCP<PHX::DataLayout> dl_dummy = rcp(new PHX::MDALayout<panzer::Dummy>(0));
109  scatterHolder_ = rcp(new PHX::Tag<ScalarT>(dummyName,dl_dummy));
110  this->addEvaluatedField(*scatterHolder_);
111 
112  std::string n = "Probe Response Scatter: " + responseName;
113  this->setName(n);
114 }
115 
116 template<typename EvalT, typename Traits, typename LO, typename GO>
119 {
120  // extract linear object container
121  responseObj_ =
122  Teuchos::rcp_dynamic_cast<Response_Probe<EvalT> >(
123  d.gedc.getDataObject(ResponseBase::buildLookupName(responseName_)),
124  true);
125 }
126 
127 
128 template<typename EvalT, typename Traits, typename LO, typename GO>
132 {
133  this->utils.setFieldData(field_,fm);
134 }
135 
136 template<typename EvalT, typename Traits, typename LO, typename GO>
139 {
140  typedef Intrepid2::CellTools<double> CTD;
141  typedef Intrepid2::FunctionSpaceTools FST;
142 
143  Kokkos::DynRankView<int,PHX::Device> inCell("inCell", 1);
144  Kokkos::DynRankView<double,PHX::Device> physical_points_cell(
145  "physical_points_cell", 1, num_dim);
146  for (size_t i=0; i<num_dim; ++i)
147  physical_points_cell(0,i) = point_[i];
148 
149  // Find which cell contains our point
150  cellIndex_ = -1;
151  bool haveProbe = false;
152  for (index_t cell=0; cell<static_cast<int>(d.num_cells); ++cell) {
153  CTD::checkPointwiseInclusion(inCell,
154  physical_points_cell,
155  this->wda(d).cell_vertex_coordinates,
156  *topology_,
157  cell);
158  if (inCell(0) == 1) {
159  cellIndex_ = cell;
160  haveProbe = true;
161  break;
162  }
163  }
164 
165  // If no cell does, we're done
166  if (!haveProbe) {
167  return false;
168  }
169 
170  // Map point to reference frame
171  const size_t num_vertex = this->wda(d).cell_vertex_coordinates.dimension_1();
172  Kokkos::DynRankView<double,PHX::Device> cell_coords(
173  "cell_coords", 1, num_vertex, num_dim); // Cell, Basis, Dim
174  for (size_t i=0; i<num_vertex; ++i) {
175  for (size_t j=0; j<num_dim; ++j) {
176  cell_coords(0,i,j) = this->wda(d).cell_vertex_coordinates(cellIndex_,i,j);
177  }
178  }
179  Kokkos::DynRankView<double,PHX::Device> physical_points(
180  "physical_points", 1, 1, num_dim); // Cell, Point, Dim
181  for (size_t i=0; i<num_dim; ++i)
182  physical_points(0,0,i) = physical_points_cell(0,i);
183  Kokkos::DynRankView<double,PHX::Device> reference_points(
184  "reference_points", 1, 1, num_dim); // Cell, Point, Dim
185  CTD::mapToReferenceFrame(reference_points, physical_points, cell_coords,
186  *topology_);
187  Kokkos::DynRankView<double,PHX::Device> reference_points_cell(
188  "reference_points_cell", 1, num_dim); // Point, Dim
189  for (size_t i=0; i<num_dim; ++i)
190  reference_points_cell(0,i) = reference_points(0,0,i);
191 
192  // Compute basis functions at point
193  if (basis_->getElementSpace() == PureBasis::CONST ||
194  basis_->getElementSpace() == PureBasis::HGRAD) {
195 
196  // Evaluate basis at reference values
197  Kokkos::DynRankView<double,PHX::Device>
198  ref_basis_values("ref_basis_values", num_basis, 1); // Basis, Point
199  basis_->getIntrepid2Basis()->getValues(ref_basis_values,
200  reference_points_cell,
201  Intrepid2::OPERATOR_VALUE);
202 
203  // Apply transformation to physical frame
204  FST::HGRADtransformVALUE<double>(basis_values_, ref_basis_values);
205 
206  }
207  else if (basis_->getElementSpace() == PureBasis::HCURL ||
208  basis_->getElementSpace() == PureBasis::HDIV) {
209 
210  // Evaluate basis at reference values
211  Kokkos::DynRankView<double,PHX::Device> ref_basis_values(
212  "ref_basis_values", num_basis, 1, num_dim); // Basis, Point, Dim
213  basis_->getIntrepid2Basis()->getValues(ref_basis_values,
214  reference_points_cell,
215  Intrepid2::OPERATOR_VALUE);
216 
217  // Apply transformation to physical frame
218  Kokkos::DynRankView<double,PHX::Device> jac
219  ("jac", 1, 1, num_dim, num_dim); // Cell, Point, Dim, Dim
220  CTD::setJacobian(jac, reference_points, cell_coords, *topology_);
221  Kokkos::DynRankView<double,PHX::Device> basis_values_vec(
222  "basis_values_vec", 1, num_basis, 1, num_dim); // Cell, Basis, Point, Dim
223  if (basis_->getElementSpace() == PureBasis::HCURL) {
224  Kokkos::DynRankView<double,PHX::Device> jac_inv(
225  "jac_inv", 1, 1, num_dim, num_dim); // Cell, Point, Dim, Dim
226  CTD::setJacobianInv(jac_inv, jac);
227  FST::HCURLtransformVALUE<double>(basis_values_vec, jac_inv,
228  ref_basis_values);
229  }
230  else {
231  Kokkos::DynRankView<double,PHX::Device> jac_det(
232  "jac_det", 1, 1); // Cell Point
233  CTD::setJacobianDet(jac_det, jac);
234  FST::HDIVtransformVALUE<double>(basis_values_vec, jac, jac_det,
235  ref_basis_values);
236  }
237 
238  // Compute element orientations
239  std::vector<double> orientation;
240  globalIndexer_->getElementOrientation(cellIndex_, orientation);
241  std::string blockId = this->wda(d).block_id;
242  int fieldNum = globalIndexer_->getFieldNum(fieldName_);
243  const std::vector<int> & elmtOffset =
244  globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
245 
246  // Extract component of basis
247  for (size_t i=0; i<num_basis; ++i) {
248  int offset = elmtOffset[i];
249  basis_values_(0,i,0) =
250  orientation[offset] * basis_values_vec(0,i,0,fieldComponent_);
251  }
252 
253  }
254 
255  return true;
256 }
257 
258 template<typename EvalT, typename Traits, typename LO, typename GO>
261 {
262  // Compute basis values at point
263  const bool haveProbe = computeBasisValues(d);
264 
265  if (!haveProbe)
266  return;
267 
268  // Get field coefficients for cell
269  Kokkos::DynRankView<ScalarT,PHX::Device> field_coeffs =
270  Kokkos::createDynRankView(field_.get_static_view(), "field_val",
271  1, num_basis); // Cell, Basis
272  for (size_t i=0; i<num_basis; ++i)
273  field_coeffs(0,i) = field_(cellIndex_,i);
274 
275  // Evaluate FE interpolant at point
276  Kokkos::DynRankView<ScalarT,PHX::Device> field_val =
277  Kokkos::createDynRankView(field_coeffs, "field_val", 1, 1); // Cell, Point
278  Intrepid2::FunctionSpaceTools::evaluate<ScalarT>(
279  field_val, field_coeffs, basis_values_);
280  responseObj_->value = field_val(0,0);
281  responseObj_->have_probe = true;
282 }
283 
284 template <typename LO, typename GO>
287 {
288  using Teuchos::RCP;
289  using Teuchos::rcp_dynamic_cast;
290  using Thyra::SpmdVectorBase;
291 
292  TEUCHOS_ASSERT(this->scatterObj_!=Teuchos::null);
293 
294  Base::evaluateFields(d);
295 
296  // grab local data for inputing
297  Teuchos::ArrayRCP<double> local_dgdx;
298  RCP<SpmdVectorBase<double> > dgdx =
299  rcp_dynamic_cast<SpmdVectorBase<double> >(this->responseObj_->getGhostedVector());
300  dgdx->getNonconstLocalData(ptrFromRef(local_dgdx));
301  TEUCHOS_ASSERT(!local_dgdx.is_null());
302 
303  this->scatterObj_->scatterDerivative(this->responseObj_->value,
304  this->cellIndex_,
305  this->responseObj_->have_probe,
306  d,this->wda,
307  local_dgdx);
308 }
309 
310 }
311 
312 #endif
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > orientation
Teuchos::RCP< GlobalEvaluationData > getDataObject(const std::string &key) const
Kokkos::DynRankView< typename InputArray::value_type, PHX::Device > createDynRankView(const InputArray &a, const std::string &name, const DimensionPack... dims)
Wrapper to simplify Panzer use of Sacado ViewFactory.
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
GlobalEvaluationDataContainer gedc
ResponseScatterEvaluator_ProbeBase(const std::string &responseName, const std::string &fieldName, const int fieldComponent, const Teuchos::Array< double > &point, const IntegrationRule &ir, const Teuchos::RCP< const PureBasis > &basis, const Teuchos::RCP< const panzer::UniqueGlobalIndexer< LO, GO > > &indexer, const Teuchos::RCP< ProbeScatterBase > &probeScatter)
A constructor with concrete arguments instead of a parameter list.