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