Panzer Version of the Day
Loading...
Searching...
No Matches
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"
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
69namespace panzer {
70
71template<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::GlobalIndexer> & 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);
116
117template<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
129template<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
260template<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
286template <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
296 Base::evaluateFields(d);
297
298 // grab local data for inputing
299 Teuchos::ArrayRCP<double> local_dgdx;
300 RCP<SpmdVectorBase<double> > 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
static std::string buildLookupName(const std::string &responseName)
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::GlobalIndexer > &indexer, const Teuchos::RCP< ProbeScatterBase > &probeScatter)
A constructor with concrete arguments instead of a parameter list.
int num_cells
DEPRECATED - use: numCells()
Teuchos::RCP< GlobalEvaluationDataContainer > gedc