Panzer  Version of the Day
Panzer_BasisValues2_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 #include "PanzerDiscFE_config.hpp"
44 #include "Panzer_Traits.hpp"
45 
47 #include "Kokkos_ViewFactory.hpp"
49 
50 #include "Intrepid2_Utils.hpp"
51 #include "Intrepid2_FunctionSpaceTools.hpp"
52 #include "Intrepid2_Orientation.hpp"
53 #include "Intrepid2_OrientationTools.hpp"
54 
55 namespace panzer {
56 namespace {
57 
58 template<typename Scalar>
59 void
60 applyOrientationsImpl(const int num_cells,
61  Kokkos::DynRankView<Scalar, PHX::Device> view,
62  Kokkos::DynRankView<Intrepid2::Orientation, PHX::Device> orientations,
63  const typename BasisValues2<Scalar>::IntrepidBasis & basis)
64 {
65  using ots=Intrepid2::OrientationTools<PHX::Device>;
66 
67  auto sub_orientations = Kokkos::subview(orientations, std::make_pair(0,num_cells));
68 
69  // There are two options based on rank
70  if(view.rank() == 3){
71  // Grab subview of object to re-orient and create a copy of it
72  auto sub_view = Kokkos::subview(view, std::make_pair(0,num_cells), Kokkos::ALL(), Kokkos::ALL());
73  auto sub_view_clone = Kokkos::createDynRankView(view, "sub_view_clone", sub_view.extent(0), sub_view.extent(1), sub_view.extent(2));
74  Kokkos::deep_copy(sub_view_clone, sub_view);
75 
76  // Apply the orientations to the subview
77  ots::modifyBasisByOrientation(sub_view, sub_view_clone, sub_orientations, &basis);
78  } else if (view.rank() == 4){
79  // Grab subview of object to re-orient and create a copy of it
80  auto sub_view = Kokkos::subview(view, std::make_pair(0,num_cells), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
81  auto sub_view_clone = Kokkos::createDynRankView(view, "sub_view_clone", sub_view.extent(0), sub_view.extent(1), sub_view.extent(2), sub_view.extent(3));
82  Kokkos::deep_copy(sub_view_clone, sub_view);
83 
84  // Apply the orientations to the subview
85  ots::modifyBasisByOrientation(sub_view, sub_view_clone, sub_orientations, &basis);
86  } else
87  throw std::logic_error("applyOrientationsImpl : Unknown view of rank " + std::to_string(view.rank()));
88 }
89 
90 template<typename Scalar>
91 void
92 applyOrientationsImpl(const int num_cells,
93  Kokkos::DynRankView<Scalar, PHX::Device> view,
94  const std::vector<Intrepid2::Orientation> & orientations,
95  const typename BasisValues2<Scalar>::IntrepidBasis & basis)
96 {
97 
98  // Move orientations vector to device
99  Kokkos::DynRankView<Intrepid2::Orientation,PHX::Device> device_orientations("drv_orts", num_cells);
100  auto host_orientations = Kokkos::create_mirror_view(device_orientations);
101  for(int i=0; i < num_cells; ++i)
102  host_orientations(i) = orientations[i];
103  Kokkos::deep_copy(device_orientations,host_orientations);
104 
105  // Call the device orientations applicator
106  applyOrientationsImpl(num_cells, view, device_orientations, basis);
107 }
108 
109 }
110 
111 
112 template <typename Scalar>
114 BasisValues2(const std::string & pre,
115  const bool allocArrays,
116  const bool buildWeighted)
117  : compute_derivatives(true)
118  , build_weighted(buildWeighted)
119  , alloc_arrays(allocArrays)
120  , prefix(pre)
121  , ddims_(1,0)
122  , references_evaluated(false)
123  , orientations_applied_(false)
124  , num_cells_(0)
125  , num_evaluate_cells_(0)
126  , is_uniform_(false)
127 
128 {
129  // Default all lazy evaluated components to not-evaluated
131  basis_scalar_evaluated_ = false;
133  basis_vector_evaluated_ = false;
135  grad_basis_evaluated_ = false;
140  div_basis_ref_evaluated_ = false;
141  div_basis_evaluated_ = false;
150 }
151 
152 template <typename Scalar>
154 evaluateValues(const PHX::MDField<Scalar,IP,Dim> & cub_points,
155  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
156  const PHX::MDField<Scalar,Cell,IP> & jac_det,
157  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
158  const int in_num_cells)
159 {
160 
161  build_weighted = false;
162 
163  setupUniform(basis_layout, cub_points, jac, jac_det, jac_inv, in_num_cells);
164 
165  const auto elmtspace = getElementSpace();
166  const int num_dims = jac.extent(2);
167 
168  // Evaluate Reference values
169  if(elmtspace == PureBasis::HGRAD or elmtspace == PureBasis::CONST or elmtspace == PureBasis::HVOL)
170  getBasisValuesRef(true,true);
171 
172  if(elmtspace == PureBasis::HDIV or elmtspace == PureBasis::HCURL)
173  getVectorBasisValuesRef(true,true);
174 
175  if(elmtspace == PureBasis::HGRAD and compute_derivatives)
176  getGradBasisValuesRef(true,true);
177 
178  if(elmtspace == PureBasis::HCURL and compute_derivatives){
179  if(num_dims == 2)
180  getCurl2DVectorBasisRef(true,true);
181  else if(num_dims == 3)
182  getCurlVectorBasisRef(true,true);
183  }
184 
185  if(elmtspace == PureBasis::HDIV and compute_derivatives)
186  getDivVectorBasisRef(true, true);
187 
188  references_evaluated = true;
189 
190  // Evaluate real space values
191  if(elmtspace == PureBasis::HGRAD or elmtspace == PureBasis::CONST or elmtspace == PureBasis::HVOL)
192  getBasisValues(false,true,true);
193 
194  if(elmtspace == PureBasis::HDIV or elmtspace == PureBasis::HCURL)
195  getVectorBasisValues(false,true,true);
196 
197  if(elmtspace == PureBasis::HGRAD and compute_derivatives)
198  getGradBasisValues(false,true,true);
199 
200  if(elmtspace == PureBasis::HCURL and compute_derivatives){
201  if(num_dims == 2)
202  getCurl2DVectorBasis(false,true,true);
203  else if(num_dims == 3)
204  getCurlVectorBasis(false,true,true);
205  }
206 
207  if(elmtspace == PureBasis::HDIV and compute_derivatives)
208  getDivVectorBasis(false,true,true);
209 
210  // Orientations have been applied only if this pointer is valid
211  orientations_applied_ = orientations_ != Teuchos::null;
212 }
213 
214 template <typename Scalar>
216 evaluateValues(const PHX::MDField<Scalar,IP,Dim> & cub_points,
217  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
218  const PHX::MDField<Scalar,Cell,IP> & jac_det,
219  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
220  const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
221  const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
222  bool use_vertex_coordinates,
223  const int in_num_cells)
224 {
225 
226  // This is the base call that will add all the non-weighted versions
227  evaluateValues(cub_points, jac, jac_det, jac_inv, in_num_cells);
228  if(weighted_measure.size() > 0)
229  setWeightedMeasure(weighted_measure);
230 
231  cell_vertex_coordinates_ = vertex_coordinates;
232 
233  // Add the weighted versions of all basis components - this will add to the non-weighted versions generated by evaluateValues above
234 
235  const auto elmtspace = getElementSpace();
236  const int num_dims = jac.extent(2);
237 
238  if(elmtspace == PureBasis::HGRAD or elmtspace == PureBasis::CONST or elmtspace == PureBasis::HVOL)
239  getBasisValues(true,true,true);
240 
241  if(elmtspace == PureBasis::HDIV or elmtspace == PureBasis::HCURL)
242  getVectorBasisValues(true,true,true);
243 
244  if(elmtspace == PureBasis::HGRAD and compute_derivatives)
245  getGradBasisValues(true,true,true);
246 
247  if(elmtspace == PureBasis::HCURL and compute_derivatives){
248  if(num_dims == 2)
249  getCurl2DVectorBasis(true,true,true);
250  else if(num_dims == 3)
251  getCurlVectorBasis(true,true,true);
252  }
253 
254  if(elmtspace == PureBasis::HDIV and compute_derivatives)
255  getDivVectorBasis(true,true,true);
256 
257  // Add the vertex components
258  if(use_vertex_coordinates){
259  getBasisCoordinatesRef(true,true);
260  getBasisCoordinates(true,true);
261  }
262 
263  // Orientations have been applied only if this pointer is valid
264  orientations_applied_ = orientations_ != Teuchos::null;
265 }
266 
267 template <typename Scalar>
269 evaluateValues(const PHX::MDField<Scalar,Cell,IP,Dim> & cub_points,
270  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
271  const PHX::MDField<Scalar,Cell,IP> & jac_det,
272  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
273  const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
274  const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
275  bool use_vertex_coordinates,
276  const int in_num_cells)
277 {
278 
279  cell_vertex_coordinates_ = vertex_coordinates;
280 
281  setup(basis_layout, cub_points, jac, jac_det, jac_inv, in_num_cells);
282  if(weighted_measure.size() > 0)
283  setWeightedMeasure(weighted_measure);
284 
285  const auto elmtspace = getElementSpace();
286  const int num_dims = jac.extent(2);
287 
288  if(elmtspace == PureBasis::HGRAD or elmtspace == PureBasis::CONST or elmtspace == PureBasis::HVOL){
289  getBasisValues(false,true,true);
290  if(build_weighted) getBasisValues(true,true,true);
291  }
292 
293  if(elmtspace == PureBasis::HDIV or elmtspace == PureBasis::HCURL){
294  getVectorBasisValues(false,true,true);
295  if(build_weighted) getVectorBasisValues(true,true,true);
296  }
297 
298  if(elmtspace == PureBasis::HGRAD and compute_derivatives){
299  getGradBasisValues(false,true,true);
300  if(build_weighted) getGradBasisValues(true,true,true);
301  }
302 
303  if(elmtspace == PureBasis::HCURL and compute_derivatives){
304  if(num_dims == 2){
305  getCurl2DVectorBasis(false,true,true);
306  if(build_weighted) getCurl2DVectorBasis(true,true,true);
307  } else if(num_dims == 3) {
308  getCurlVectorBasis(false,true,true);
309  if(build_weighted) getCurlVectorBasis(true,true,true);
310  }
311  }
312 
313  if(elmtspace == PureBasis::HDIV and compute_derivatives){
314  getDivVectorBasis(false,true,true);
315  if(build_weighted) getDivVectorBasis(true,true,true);
316  }
317 
318  // Add the vertex components
319  if(use_vertex_coordinates){
320  getBasisCoordinatesRef(true,true);
321  getBasisCoordinates(true,true);
322  }
323 
324  // Orientations have been applied only if this pointer is valid
325  orientations_applied_ = orientations_ != Teuchos::null;
326 
327 }
328 
329 template <typename Scalar>
331 evaluateValuesCV(const PHX::MDField<Scalar,Cell,IP,Dim> & cub_points,
332  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
333  const PHX::MDField<Scalar,Cell,IP> & jac_det,
334  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv)
335 {
336 
337  PHX::MDField<Scalar,Cell,IP> weighted_measure;
338  PHX::MDField<Scalar,Cell,NODE,Dim> vertex_coordinates;
339  evaluateValues(cub_points,jac,jac_det,jac_inv,weighted_measure,vertex_coordinates,false,jac.extent(0));
340 
341 }
342 
343 template <typename Scalar>
345 evaluateValuesCV(const PHX::MDField<Scalar,Cell,IP,Dim> & cell_cub_points,
346  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
347  const PHX::MDField<Scalar,Cell,IP> & jac_det,
348  const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
349  const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
350  bool use_vertex_coordinates,
351  const int in_num_cells)
352 {
353  PHX::MDField<Scalar,Cell,IP> weighted_measure;
354  evaluateValues(cell_cub_points,jac,jac_det,jac_inv,weighted_measure,vertex_coordinates,use_vertex_coordinates,in_num_cells);
355 }
356 
357 template <typename Scalar>
359 evaluateBasisCoordinates(const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
360  const int in_num_cells)
361 {
362  num_evaluate_cells_ = in_num_cells < 0 ? vertex_coordinates.extent(0) : in_num_cells;
363  cell_vertex_coordinates_ = vertex_coordinates;
364 
365  getBasisCoordinates(true,true);
366 }
367 
368 // method for applying orientations
369 template <typename Scalar>
371 applyOrientations(const std::vector<Intrepid2::Orientation> & orientations,
372  const int in_num_cells)
373 {
374  if (!intrepid_basis->requireOrientation())
375  return;
376 
377  // We only allow the orientations to be applied once
378  TEUCHOS_ASSERT(not orientations_applied_);
379 
380  const PureBasis::EElementSpace elmtspace = getElementSpace();
381 
382  const int num_cell_basis_layout = in_num_cells < 0 ? basis_layout->numCells() : in_num_cells;
383  const int num_cell_orientation = orientations.size();
384  const int num_cells = num_cell_basis_layout < num_cell_orientation ? num_cell_basis_layout : num_cell_orientation;
385  const int num_dim = basis_layout->dimension();
386 
387  // Copy orientations to device
388  Kokkos::DynRankView<Intrepid2::Orientation,PHX::Device> device_orientations("device_orientations", num_cells);
389  auto host_orientations = Kokkos::create_mirror_view(device_orientations);
390  for(int i=0; i < num_cells; ++i)
391  host_orientations(i) = orientations[i];
392  Kokkos::deep_copy(device_orientations,host_orientations);
393 
394  if (elmtspace == PureBasis::HGRAD){
395  applyOrientationsImpl<Scalar>(num_cells, basis_scalar.get_view(), device_orientations, *intrepid_basis);
396  if(build_weighted) applyOrientationsImpl<Scalar>(num_cells, weighted_basis_scalar.get_view(), device_orientations, *intrepid_basis);
397  if(compute_derivatives) applyOrientationsImpl<Scalar>(num_cells, grad_basis.get_view(), device_orientations, *intrepid_basis);
398  if(compute_derivatives and build_weighted) applyOrientationsImpl<Scalar>(num_cells, weighted_grad_basis.get_view(), device_orientations, *intrepid_basis);
399  }
400 
401  if(elmtspace == PureBasis::HCURL){
402  applyOrientationsImpl<Scalar>(num_cells, basis_vector.get_view(), device_orientations, *intrepid_basis);
403  if(build_weighted) applyOrientationsImpl<Scalar>(num_cells, weighted_basis_vector.get_view(), device_orientations, *intrepid_basis);
404  if(num_dim == 2){
405  if(compute_derivatives) applyOrientationsImpl<Scalar>(num_cells, curl_basis_scalar.get_view(), device_orientations, *intrepid_basis);
406  if(compute_derivatives and build_weighted) applyOrientationsImpl<Scalar>(num_cells, weighted_curl_basis_scalar.get_view(), device_orientations, *intrepid_basis);
407  }
408  if(num_dim == 3){
409  if(compute_derivatives) applyOrientationsImpl<Scalar>(num_cells, curl_basis_vector.get_view(), device_orientations, *intrepid_basis);
410  if(compute_derivatives and build_weighted) applyOrientationsImpl<Scalar>(num_cells, weighted_curl_basis_vector.get_view(), device_orientations, *intrepid_basis);
411  }
412  }
413 
414  if(elmtspace == PureBasis::HDIV){
415  applyOrientationsImpl<Scalar>(num_cells, basis_vector.get_view(), device_orientations, *intrepid_basis);
416  if(build_weighted) applyOrientationsImpl<Scalar>(num_cells, weighted_basis_vector.get_view(), device_orientations, *intrepid_basis);
417  if(compute_derivatives) applyOrientationsImpl<Scalar>(num_cells, div_basis.get_view(), device_orientations, *intrepid_basis);
418  if(compute_derivatives and build_weighted) applyOrientationsImpl<Scalar>(num_cells, weighted_div_basis.get_view(), device_orientations, *intrepid_basis);
419  }
420 
421  orientations_applied_ = true;
422 }
423 
424 // method for applying orientations
425 template <typename Scalar>
427 applyOrientations(const PHX::MDField<const Scalar,Cell,BASIS> & orientations)
428 {
429  TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::applyOrientations : this should not be called.");
430 }
431 
432 template <typename Scalar>
434 { return basis_layout->getBasis()->getElementSpace(); }
435 
436 template <typename Scalar>
438 setupArrays(const Teuchos::RCP<const panzer::BasisIRLayout>& layout,
439  bool computeDerivatives)
440 {
441  MDFieldArrayFactory af(prefix,alloc_arrays);
442 
443  compute_derivatives = computeDerivatives;
444  basis_layout = layout;
445  num_cells_ = basis_layout->numCells();
446  Teuchos::RCP<const panzer::PureBasis> basisDesc = layout->getBasis();
447 
448  // for convience pull out basis and quadrature information
449  int num_quad = layout->numPoints();
450  int dim = basisDesc->dimension();
451  int card = basisDesc->cardinality();
452  int numcells = basisDesc->numCells();
453  panzer::PureBasis::EElementSpace elmtspace = basisDesc->getElementSpace();
454  Teuchos::RCP<const shards::CellTopology> cellTopo = basisDesc->getCellTopology();
455 
456  intrepid_basis = basisDesc->getIntrepid2Basis<PHX::Device::execution_space,Scalar,Scalar>();
457 
458  // allocate field containers
459  // field sizes defined by http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/doc/html/basis_page.html#basis_md_array_sec
460 
461  // compute basis fields
462  if(elmtspace==panzer::PureBasis::HGRAD) {
463  // HGRAD is a nodal field
464 
465  // build values
467  basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("basis_ref",card,num_quad); // F, P
468  basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("basis",numcells,card,num_quad);
469 
470  if(build_weighted)
471  weighted_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_basis",numcells,card,num_quad);
472 
473  // build gradients
475 
476  if(compute_derivatives) {
477  grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("grad_basis_ref",card,num_quad,dim); // F, P, D
478  grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("grad_basis",numcells,card,num_quad,dim);
479 
480  if(build_weighted)
481  weighted_grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_grad_basis",numcells,card,num_quad,dim);
482  }
483 
484  // build curl
486 
487  // nothing - HGRAD does not support CURL operation
488  }
489  else if(elmtspace==panzer::PureBasis::HCURL) {
490  // HCURL is a vector field
491 
492  // build values
494 
495  basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("basis_ref",card,num_quad,dim); // F, P, D
496  basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis",numcells,card,num_quad,dim);
497 
498  if(build_weighted)
499  weighted_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_basis",numcells,card,num_quad,dim);
500 
501  // build gradients
503 
504  // nothing - HCURL does not support GRAD operation
505 
506  // build curl
508 
509  if(compute_derivatives) {
510  if(dim==2) {
511  // curl of HCURL basis is not dimension dependent
512  curl_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("curl_basis_ref",card,num_quad); // F, P
513  curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("curl_basis",numcells,card,num_quad);
514 
515  if(build_weighted)
516  weighted_curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_curl_basis",numcells,card,num_quad);
517  }
518  else if(dim==3){
519  curl_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("curl_basis_ref",card,num_quad,dim); // F, P, D
520  curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("curl_basis",numcells,card,num_quad,dim);
521 
522  if(build_weighted)
523  weighted_curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_curl_basis",numcells,card,num_quad,dim);
524  }
525  else { TEUCHOS_ASSERT(false); } // what do I do with 1D?
526  }
527  }
528  else if(elmtspace==panzer::PureBasis::HDIV) {
529  // HDIV is a vector field
530 
531  // build values
533 
534  basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("basis_ref",card,num_quad,dim); // F, P, D
535  basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis",numcells,card,num_quad,dim);
536 
537  if(build_weighted)
538  weighted_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_basis",numcells,card,num_quad,dim);
539 
540  // build gradients
542 
543  // nothing - HCURL does not support GRAD operation
544 
545  // build curl
547 
548  // nothing - HDIV does not support CURL operation
549 
550  // build div
552 
553  if(compute_derivatives) {
554  div_basis_ref = af.buildStaticArray<Scalar,BASIS,IP>("div_basis_ref",card,num_quad); // F, P
555  div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("div_basis",numcells,card,num_quad);
556 
557  if(build_weighted)
558  weighted_div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_div_basis",numcells,card,num_quad);
559  }
560  }
561  else if(elmtspace==panzer::PureBasis::CONST || elmtspace==panzer::PureBasis::HVOL) {
562  // CONST is a nodal field
563 
564  // build values
566  basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("basis_ref",card,num_quad); // F, P
567  basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("basis",numcells,card,num_quad);
568 
569  if(build_weighted)
570  weighted_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_basis",numcells,card,num_quad);
571 
572  // build gradients
574 
575  // nothing - CONST does not support GRAD operation
576 
577  // build curl
579 
580  // nothing - CONST does not support CURL operation
581 
582  // build div
584 
585  // nothing - CONST does not support DIV operation
586  }
587  else { TEUCHOS_ASSERT(false); }
588 
589  basis_coordinates_ref = af.buildStaticArray<Scalar,BASIS,Dim>("basis_coordinates_ref",card,dim);
590  basis_coordinates = af.buildStaticArray<Scalar,Cell,BASIS,Dim>("basis_coordinates",numcells,card,dim);
591 }
592 
593 
594 //=======================================================================================================
595 // New Interface
596 
597 
598 template <typename Scalar>
599 void
601 setup(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
602  PHX::MDField<Scalar, Cell, IP, Dim> reference_points,
603  PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian,
604  PHX::MDField<Scalar, Cell, IP> point_jacobian_determinant,
605  PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
606  const int num_evaluated_cells)
607 {
608  basis_layout = basis;
609  intrepid_basis = basis->getBasis()->getIntrepid2Basis<PHX::Device::execution_space,Scalar,Scalar>();
610  num_cells_ = basis_layout->numCells();
611  num_evaluate_cells_ = num_evaluated_cells >= 0 ? num_evaluated_cells : num_cells_;
612  build_weighted = false;
613  is_uniform_ = false;
614 
615  cubature_points_ref_ = reference_points;
616  cubature_jacobian_ = point_jacobian;
617  cubature_jacobian_determinant_ = point_jacobian_determinant;
618  cubature_jacobian_inverse_ = point_jacobian_inverse;
619 
620  // Reset internal data
621  resetArrays();
622 
623 }
624 
625 template <typename Scalar>
626 void
628 setupUniform(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
629  PHX::MDField<Scalar, IP, Dim> reference_points,
630  PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian,
631  PHX::MDField<Scalar, Cell, IP> point_jacobian_determinant,
632  PHX::MDField<Scalar, Cell, IP, Dim, Dim> point_jacobian_inverse,
633  const int num_evaluated_cells)
634 {
635  basis_layout = basis;
636  intrepid_basis = basis->getBasis()->getIntrepid2Basis<PHX::Device::execution_space,Scalar,Scalar>();
637  num_cells_ = basis_layout->numCells();
638  num_evaluate_cells_ = num_evaluated_cells >= 0 ? num_evaluated_cells : num_cells_;
639  cubature_points_uniform_ref_ = reference_points;
640  build_weighted = false;
641  is_uniform_ = true;
642 
643  cubature_jacobian_ = point_jacobian;
644  cubature_jacobian_determinant_ = point_jacobian_determinant;
645  cubature_jacobian_inverse_ = point_jacobian_inverse;
646 
647  // Reset internal data
648  resetArrays();
649 
650 }
651 
652 template <typename Scalar>
653 void
655 setOrientations(const Teuchos::RCP<const OrientationsInterface> & orientations,
656  const int num_orientations_cells)
657 {
658  if(num_orientations_cells < 0)
659  num_orientations_cells_ = num_evaluate_cells_;
660  else
661  num_orientations_cells_ = num_orientations_cells;
662  if(orientations == Teuchos::null){
663  orientations_applied_ = false;
664  orientations_ = Teuchos::null;
665  // Normally we would reset arrays here, but it seems to causes a lot of problems
666  } else {
667  orientations_ = orientations;
668  orientations_applied_ = true;
669  // Setting orientations means that we need to reset our arrays
670  resetArrays();
671  }
672 }
673 
674 template <typename Scalar>
675 void
677 setCellVertexCoordinates(PHX::MDField<Scalar,Cell,NODE,Dim> vertex_coordinates)
678 {
679  cell_vertex_coordinates_ = vertex_coordinates;
680 }
681 
682 template <typename Scalar>
683 void
686 {
687  // Turn off all evaluated fields (forces re-evaluation)
688  basis_ref_scalar_evaluated_ = false;
689  basis_scalar_evaluated_ = false;
690  basis_ref_vector_evaluated_ = false;
691  basis_vector_evaluated_ = false;
692  grad_basis_ref_evaluated_ = false;
693  grad_basis_evaluated_ = false;
694  curl_basis_ref_scalar_evaluated_ = false;
695  curl_basis_scalar_evaluated_ = false;
696  curl_basis_ref_vector_evaluated_ = false;
697  curl_basis_vector_evaluated_ = false;
698  div_basis_ref_evaluated_ = false;
699  div_basis_evaluated_ = false;
700  weighted_basis_scalar_evaluated_ = false;
701  weighted_basis_vector_evaluated_ = false;
702  weighted_grad_basis_evaluated_ = false;
703  weighted_curl_basis_scalar_evaluated_ = false;
704  weighted_curl_basis_vector_evaluated_ = false;
705  weighted_div_basis_evaluated_ = false;
706  basis_coordinates_ref_evaluated_ = false;
707  basis_coordinates_evaluated_ = false;
708 
709  // TODO: Enable this feature - requires the old interface to go away
710  // De-allocate arrays if necessary
711 // if(not alloc_arrays){
712 // basis_ref_scalar = Array_BasisIP();
713 // basis_scalar = Array_CellBasisIP();
714 // basis_ref_vector = Array_BasisIPDim();
715 // basis_vector = Array_CellBasisIPDim();
716 // grad_basis_ref = Array_BasisIPDim();
717 // grad_basis = Array_CellBasisIPDim();
718 // curl_basis_ref_scalar = Array_BasisIP();
719 // curl_basis_scalar = Array_CellBasisIP();
720 // curl_basis_ref_vector = Array_BasisIPDim();
721 // curl_basis_vector = Array_CellBasisIPDim();
722 // div_basis_ref = Array_BasisIP();
723 // div_basis = Array_CellBasisIP();
724 // weighted_basis_scalar = Array_CellBasisIP();
725 // weighted_basis_vector = Array_CellBasisIPDim();
726 // weighted_grad_basis = Array_CellBasisIPDim();
727 // weighted_curl_basis_scalar = Array_CellBasisIP();
728 // weighted_curl_basis_vector = Array_CellBasisIPDim();
729 // weighted_div_basis = Array_CellBasisIP();
730 // basis_coordinates_ref = Array_BasisDim();
731 // basis_coordinates = Array_CellBasisDim();
732 // }
733 }
734 
735 template <typename Scalar>
736 void
738 setWeightedMeasure(PHX::MDField<Scalar, Cell, IP> weighted_measure)
739 {
740  TEUCHOS_TEST_FOR_EXCEPT_MSG(build_weighted,
741  "BasisValues2::setWeightedMeasure : Weighted measure already set. Can only set weighted measure once after setup or setupUniform have beens called.");
742  cubature_weights_ = weighted_measure;
743  build_weighted = true;
744 }
745 
746 // If the array is allocated we can not reassign it - this means we
747 // have to deep copy into it. The use of deep_copy is need when an
748 // applicaiton or the library has cached a BasisValues2 member view
749 // outside of the BasisValues object. This could happen when the basis
750 // values object is embedded in an evalautor for mesh motion.
751 #define PANZER_CACHE_DATA(name) \
752  if(cache) { \
753  if(name.size()==tmp_##name.size()){ \
754  Kokkos::deep_copy(name.get_view(), tmp_##name.get_view()); \
755  } else { \
756  name = tmp_##name; \
757  } \
758  name##_evaluated_ = true; \
759  }
760 
761 template <typename Scalar>
764 getBasisCoordinatesRef(const bool cache,
765  const bool force) const
766 {
767  // Check if array already exists
768  if(basis_coordinates_ref_evaluated_ and not force)
769  return basis_coordinates_ref;
770 
771  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
772 
773  const int num_card = basis_layout->cardinality();
774  const int num_dim = basis_layout->dimension();
775 
776  using coordsScalarType = typename Intrepid2::Basis<PHX::Device::execution_space,Scalar,Scalar>::scalarType;
777  auto tmp_basis_coordinates_ref = af.buildStaticArray<coordsScalarType,BASIS,Dim>("basis_coordinates_ref", num_card, num_dim);
778  intrepid_basis->getDofCoords(tmp_basis_coordinates_ref.get_view());
779  PHX::Device().fence();
780 
781  // Store for later if cache is enabled
782  PANZER_CACHE_DATA(basis_coordinates_ref)
783 
784  return tmp_basis_coordinates_ref;
785 }
786 
787 template <typename Scalar>
790 getBasisValuesRef(const bool cache,
791  const bool force) const
792 {
793  // Check if array already exists
794  if(basis_ref_scalar_evaluated_ and not force)
795  return basis_ref_scalar;
796 
797  // Reference quantities only exist for a uniform reference space
798  TEUCHOS_ASSERT(hasUniformReferenceSpace());
799 
800  // Make sure basis is valid
801  PureBasis::EElementSpace elmtspace = getElementSpace();
802  TEUCHOS_ASSERT(elmtspace==PureBasis::HGRAD || elmtspace==PureBasis::CONST || elmtspace==PureBasis::HVOL);
803 
804  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
805 
806  const int num_quad = basis_layout->numPoints();
807  const int num_card = basis_layout->cardinality();
808 
809  auto tmp_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_quad);
810  intrepid_basis->getValues(tmp_basis_ref_scalar.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_VALUE);
811  PHX::Device().fence();
812 
813  // Store for later if cache is enabled
814  PANZER_CACHE_DATA(basis_ref_scalar);
815 
816  return tmp_basis_ref_scalar;
817 }
818 
819 template <typename Scalar>
822 getVectorBasisValuesRef(const bool cache,
823  const bool force) const
824 {
825  // Check if array already exists
826  if(basis_ref_vector_evaluated_ and not force)
827  return basis_ref_vector;
828 
829  // Reference quantities only exist for a uniform reference space
830  TEUCHOS_ASSERT(hasUniformReferenceSpace());
831 
832  // Make sure basis is valid
833  PureBasis::EElementSpace elmtspace = getElementSpace();
834  TEUCHOS_ASSERT(elmtspace==PureBasis::HDIV || elmtspace==PureBasis::HCURL);
835 
836  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
837 
838  const int num_quad = basis_layout->numPoints();
839  const int num_card = basis_layout->cardinality();
840  const int num_dim = basis_layout->dimension();
841 
842  auto tmp_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
843  intrepid_basis->getValues(tmp_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE);
844  PHX::Device().fence();
845 
846  // Store for later if cache is enabled
847  PANZER_CACHE_DATA(basis_ref_vector);
848 
849  return tmp_basis_ref_vector;
850 }
851 
852 template <typename Scalar>
855 getGradBasisValuesRef(const bool cache,
856  const bool force) const
857 {
858  // Check if array already exists
859  if(grad_basis_ref_evaluated_ and not force)
860  return grad_basis_ref;
861 
862  // Reference quantities only exist for a uniform reference space
863  TEUCHOS_ASSERT(hasUniformReferenceSpace());
864 
865  // Make sure basis is valid
866  PureBasis::EElementSpace elmtspace = getElementSpace();
867  TEUCHOS_ASSERT(elmtspace==PureBasis::HGRAD);
868 
869  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
870 
871  const int num_quad = basis_layout->numPoints();
872  const int num_card = basis_layout->cardinality();
873  const int num_dim = basis_layout->dimension();
874 
875  auto tmp_grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
876  intrepid_basis->getValues(tmp_grad_basis_ref.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_GRAD);
877  PHX::Device().fence();
878 
879  // Store for later if cache is enabled
880  PANZER_CACHE_DATA(grad_basis_ref);
881 
882  return tmp_grad_basis_ref;
883 }
884 
885 template <typename Scalar>
888 getCurl2DVectorBasisRef(const bool cache,
889  const bool force) const
890 {
891  // Check if array already exists
892  if(curl_basis_ref_scalar_evaluated_ and not force)
893  return curl_basis_ref_scalar;
894 
895  // Reference quantities only exist for a uniform reference space
896  TEUCHOS_ASSERT(hasUniformReferenceSpace());
897  TEUCHOS_ASSERT(basis_layout->dimension() == 2);
898 
899  // Make sure basis is valid
900  PureBasis::EElementSpace elmtspace = getElementSpace();
901  TEUCHOS_ASSERT(elmtspace==PureBasis::HCURL);
902 
903  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
904 
905  const int num_quad = basis_layout->numPoints();
906  const int num_card = basis_layout->cardinality();
907 
908  auto tmp_curl_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("dyn_curl_basis_ref_scalar",num_card,num_quad);
909  intrepid_basis->getValues(tmp_curl_basis_ref_scalar.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_CURL);
910  PHX::Device().fence();
911 
912  // Store for later if cache is enabled
913  PANZER_CACHE_DATA(curl_basis_ref_scalar);
914 
915  return tmp_curl_basis_ref_scalar;
916 }
917 
918 template <typename Scalar>
921 getCurlVectorBasisRef(const bool cache,
922  const bool force) const
923 {
924  // Check if array already exists
925  if(curl_basis_ref_vector_evaluated_ and not force)
926  return curl_basis_ref_vector;
927 
928  // Reference quantities only exist for a uniform reference space
929  TEUCHOS_ASSERT(hasUniformReferenceSpace());
930  TEUCHOS_ASSERT(basis_layout->dimension() == 3);
931 
932  // Make sure basis is valid
933  PureBasis::EElementSpace elmtspace = getElementSpace();
934  TEUCHOS_ASSERT(elmtspace==PureBasis::HCURL);
935 
936  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
937 
938  const int num_quad = basis_layout->numPoints();
939  const int num_card = basis_layout->cardinality();
940  const int num_dim = basis_layout->dimension();
941 
942  auto tmp_curl_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("dyn_curl_basis_ref_vector",num_card,num_quad,num_dim);
943  intrepid_basis->getValues(tmp_curl_basis_ref_vector.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_CURL);
944  PHX::Device().fence();
945 
946  // Store for later if cache is enabled
947  PANZER_CACHE_DATA(curl_basis_ref_vector);
948 
949  return tmp_curl_basis_ref_vector;
950 }
951 
952 template <typename Scalar>
955 getDivVectorBasisRef(const bool cache,
956  const bool force) const
957 {
958  // Check if array already exists
959  if(div_basis_ref_evaluated_ and not force)
960  return div_basis_ref;
961 
962  // Reference quantities only exist for a uniform reference space
963  TEUCHOS_ASSERT(hasUniformReferenceSpace());
964 
965  // Make sure basis is valid
966  PureBasis::EElementSpace elmtspace = getElementSpace();
967  TEUCHOS_ASSERT(elmtspace==PureBasis::HDIV);
968 
969  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
970 
971  const int num_quad = basis_layout->numPoints();
972  const int num_card = basis_layout->cardinality();
973 
974  auto tmp_div_basis_ref = af.buildStaticArray<Scalar,BASIS,IP>("dyn_div_basis_ref_scalar",num_card,num_quad);
975  intrepid_basis->getValues(tmp_div_basis_ref.get_view(), cubature_points_uniform_ref_.get_view(), Intrepid2::OPERATOR_DIV);
976  PHX::Device().fence();
977 
978  // Store for later if cache is enabled
979  PANZER_CACHE_DATA(div_basis_ref);
980 
981  return tmp_div_basis_ref;
982 }
983 
984 template <typename Scalar>
987 getBasisCoordinates(const bool cache,
988  const bool force) const
989 {
990  // Check if array already exists
991  if(basis_coordinates_evaluated_ and not force)
992  return basis_coordinates;
993 
994  TEUCHOS_ASSERT(cell_vertex_coordinates_.size() > 0);
995 
996  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
997  const auto s_vertex_coordinates = Kokkos::subview(cell_vertex_coordinates_.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
998 
999  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
1000 
1001  const int num_card = basis_layout->cardinality();
1002  const int num_dim = basis_layout->dimension();
1003 
1004  auto tmp_basis_coordinates = af.buildStaticArray<Scalar, Cell, BASIS, IP>("basis_coordinates", num_cells_, num_card, num_dim);
1005  auto s_aux = Kokkos::subview(tmp_basis_coordinates.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1006 
1007  // Some kind of type bug in Intrepid's mapToPhysicalFrame call? Can't let bcr be const
1008  auto bcr = getBasisCoordinatesRef(false);
1009  Kokkos::DynRankView<double> bcr_copy("nonconst_bcr",bcr.extent(0),bcr.extent(1));
1010  Kokkos::deep_copy(bcr_copy,bcr.get_view());
1011 
1012  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1013  cell_tools.mapToPhysicalFrame(s_aux, bcr_copy, s_vertex_coordinates, intrepid_basis->getBaseCellTopology());
1014  PHX::Device().fence();
1015 
1016  // Store for later if cache is enabled
1017  PANZER_CACHE_DATA(basis_coordinates);
1018 
1019  return tmp_basis_coordinates;
1020 }
1021 
1022 template <typename Scalar>
1025 getBasisValues(const bool weighted,
1026  const bool cache,
1027  const bool force) const
1028 {
1029  if(weighted){
1030  if(weighted_basis_scalar_evaluated_ and not force)
1031  return weighted_basis_scalar;
1032  } else
1033  if(basis_scalar_evaluated_ and not force)
1034  return basis_scalar;
1035 
1036  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
1037 
1038  const int num_cells = num_cells_;
1039  const int num_points = basis_layout->numPoints();
1040  const int num_card = basis_layout->cardinality();
1041  const int num_dim = basis_layout->dimension();
1042 
1043  if(weighted){
1044 
1045  TEUCHOS_ASSERT(cubature_weights_.size() > 0);
1046 
1047  // Get the basis_scalar values - do not cache
1048  const auto bv = getBasisValues(false, force);
1049 
1050  // Apply the weighted measure (cubature weights)
1051  auto tmp_weighted_basis_scalar = af.buildStaticArray<Scalar, Cell, BASIS, IP>("weighted_basis_scalar", num_cells, num_card, num_points);
1052 
1053  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1054  auto s_aux = Kokkos::subview(tmp_weighted_basis_scalar.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1055  auto s_cw = Kokkos::subview(cubature_weights_.get_view(),cell_range,Kokkos::ALL());
1056  auto s_bv = Kokkos::subview(bv.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1057 
1058  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1059  fst::multiplyMeasure(s_aux,s_cw,s_bv);
1060 
1061  // NOTE: Weighted path has orientations already applied so doesn't
1062  // need the applyOrientations call at the bottom of this function.
1063 
1064  // Store for later if cache is enabled
1065  PANZER_CACHE_DATA(weighted_basis_scalar);
1066 
1067  return tmp_weighted_basis_scalar;
1068 
1069  } else {
1070 
1071  const auto element_space = getElementSpace();
1072  TEUCHOS_ASSERT(element_space == PureBasis::HVOL || element_space == PureBasis::HGRAD || element_space == PureBasis::CONST);
1073 
1074  // HVol requires the jacobian determinant
1075  if(element_space == PureBasis::HVOL){
1076  TEUCHOS_ASSERT(cubature_jacobian_determinant_.size() > 0);
1077  }
1078 
1079  auto cell_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("cell_basis_ref_scalar",num_card,num_points);
1080  auto tmp_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("basis_scalar",num_cells,num_card,num_points);
1081 
1082  if(hasUniformReferenceSpace()){
1083 
1084  // Apply a single reference representation to all cells
1085  intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE);
1086 
1087  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1088  auto s_aux = Kokkos::subview(tmp_basis_scalar.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1089 
1090  // Apply transformation (HGRAD version is just a copy operation)
1091  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1092  if(element_space == PureBasis::HVOL){
1093  auto s_cjd = Kokkos::subview(cubature_jacobian_determinant_.get_view(), cell_range, Kokkos::ALL());
1094  fst::HVOLtransformVALUE(s_aux,s_cjd,cell_basis_ref_scalar.get_view());
1095  }else if(element_space == PureBasis::HGRAD || element_space == PureBasis::CONST)
1096  fst::HGRADtransformVALUE(s_aux,cell_basis_ref_scalar.get_view());
1097 
1098  PHX::Device().fence();
1099 
1100  } else {
1101 
1102  // This is ugly. The algorithm is restricted to host/serial due
1103  // to a call to intrepid tools that require uniform reference
1104  // representation. For DG, CVFEM and sidesets this reference is
1105  // nonuniform.
1106 
1107  // Local allocation used for each cell
1108  auto cell_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("cell_basis_scalar",1,num_card,num_points);
1109  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
1110  auto cell_jac_det = af.buildStaticArray<Scalar,Cell,IP>("cell_jac_det",1,num_points);
1111 
1112  // The array factory is difficult to extend to host space
1113  // without extra template magic and changing a ton of code in a
1114  // non-backwards compatible way, so we use some of the arrays
1115  // above only to get derivative array sized correctly and then
1116  // create the mirror on host.
1117  auto cell_basis_scalar_host = Kokkos::create_mirror_view(cell_basis_scalar.get_view());
1118  auto cell_cub_points_host = Kokkos::create_mirror_view(cell_cub_points.get_view());
1119  auto cell_jac_det_host = Kokkos::create_mirror_view(cell_jac_det.get_view());
1120  auto cell_basis_ref_scalar_host = Kokkos::create_mirror_view(cell_basis_ref_scalar.get_view());
1121  auto cubature_points_ref_host = Kokkos::create_mirror_view(cubature_points_ref_.get_view());
1122  Kokkos::deep_copy(cubature_points_ref_host,cubature_points_ref_.get_view());
1123  auto cubature_jacobian_determinant_host = Kokkos::create_mirror_view(cubature_jacobian_determinant_.get_view());
1124  Kokkos::deep_copy(cubature_jacobian_determinant_host,cubature_jacobian_determinant_.get_view());
1125  auto tmp_basis_scalar_host = Kokkos::create_mirror_view(tmp_basis_scalar.get_view());
1126 
1127  // We have to iterate through cells and apply a separate reference representation for each
1128  for(int cell=0; cell<num_evaluate_cells_; ++cell){
1129 
1130  // Load external into cell-local arrays
1131  for(int p=0;p<num_points;++p)
1132  for(int d=0;d<num_dim;++d)
1133  cell_cub_points_host(p,d)=cubature_points_ref_host(cell,p,d);
1134  Kokkos::deep_copy(cell_cub_points.get_view(),cell_cub_points_host);
1135 
1136  // Convert to mesh space. The intrepid_basis is virtual and
1137  // built in another class, short of adding a "clone on host"
1138  // method, we will have to make this call on device. This is a
1139  // major performance hit.
1140  intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE);
1141  Kokkos::deep_copy(cell_basis_ref_scalar_host,cell_basis_ref_scalar.get_view());
1142 
1143  using fst=Intrepid2::FunctionSpaceTools<Kokkos::HostSpace>;
1144 
1145  if(element_space == PureBasis::HVOL){
1146  // Need the jacobian determinant for HVOL
1147  for(int p=0;p<num_points;++p)
1148  cell_jac_det_host(0,p)=cubature_jacobian_determinant_host(cell,p);
1149 
1150  fst::HVOLtransformVALUE(cell_basis_scalar_host,cell_jac_det_host,cell_basis_ref_scalar_host);
1151  } else if(element_space == PureBasis::HGRAD || element_space == PureBasis::CONST)
1152  fst::HGRADtransformVALUE(cell_basis_scalar_host,cell_basis_ref_scalar_host);
1153 
1154  // Copy cell quantity back into main array
1155  for(int b=0; b<num_card; ++b)
1156  for(int p=0; p<num_points; ++p)
1157  tmp_basis_scalar_host(cell,b,p) = cell_basis_scalar_host(0,b,p);
1158 
1159  Kokkos::deep_copy(tmp_basis_scalar.get_view(),tmp_basis_scalar_host);
1160  }
1161  }
1162 
1163  // NOTE: weighted already has orientations applied, so this code
1164  // should only be reached if non-weighted. This is a by-product of
1165  // the two construction paths in panzer. Unification should help
1166  // fix the logic.
1167 
1168  if(orientations_ != Teuchos::null)
1169  applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_basis_scalar.get_view(), *orientations_->getOrientations(), *intrepid_basis);
1170 
1171  // Store for later if cache is enabled
1172  PANZER_CACHE_DATA(basis_scalar);
1173 
1174  return tmp_basis_scalar;
1175 
1176  }
1177 
1178 }
1179 
1180 template <typename Scalar>
1183 getVectorBasisValues(const bool weighted,
1184  const bool cache,
1185  const bool force) const
1186 {
1187  if(weighted){
1188  if(weighted_basis_vector_evaluated_ and not force)
1189  return weighted_basis_vector;
1190  } else
1191  if(basis_vector_evaluated_ and not force)
1192  return basis_vector;
1193 
1194  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
1195 
1196  const int num_cells = num_cells_;
1197  const int num_points = basis_layout->numPoints();
1198  const int num_card = basis_layout->cardinality();
1199  const int num_dim = basis_layout->dimension();
1200 
1201  if(weighted){
1202 
1203  TEUCHOS_ASSERT(cubature_weights_.size() > 0);
1204 
1205  // Get the basis_scalar values - do not cache
1206  const auto bv = getVectorBasisValues(false, force);
1207 
1208  // Apply the weighted measure (cubature weights)
1209  auto tmp_weighted_basis_vector = af.buildStaticArray<Scalar, Cell, BASIS, IP, Dim>("weighted_basis_vector", num_cells, num_card, num_points, num_dim);
1210 
1211  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1212  auto s_aux = Kokkos::subview(tmp_weighted_basis_vector.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1213  auto s_cw = Kokkos::subview(cubature_weights_.get_view(),cell_range,Kokkos::ALL());
1214  auto s_bv = Kokkos::subview(bv.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1215 
1216  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1217  fst::multiplyMeasure(s_aux, s_cw, s_bv);
1218 
1219  // Store for later if cache is enabled
1220  PANZER_CACHE_DATA(weighted_basis_vector);
1221 
1222  return tmp_weighted_basis_vector;
1223 
1224  } else { // non-weighted
1225 
1226  const auto element_space = getElementSpace();
1227  TEUCHOS_ASSERT(element_space == PureBasis::HCURL || element_space == PureBasis::HDIV);
1228  TEUCHOS_ASSERT(num_dim != 1);
1229 
1230  // HDIV and HCURL have unique jacobian requirements
1231  if(element_space == PureBasis::HCURL){
1232  TEUCHOS_ASSERT(cubature_jacobian_inverse_.size() > 0);
1233  } else if(element_space == PureBasis::HDIV){
1234  TEUCHOS_ASSERT(cubature_jacobian_.size() > 0 && cubature_jacobian_determinant_.size() > 0);
1235  }
1236 
1237  auto cell_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("cell_basis_ref_scalar",num_card,num_points,num_dim);
1238  auto tmp_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis_vector",num_cells,num_card,num_points,num_dim);
1239 
1240  if(hasUniformReferenceSpace()){
1241 
1242  // Apply a single reference representation to all cells
1243  intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_VALUE);
1244 
1245  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1246  auto s_aux = Kokkos::subview(tmp_basis_vector.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1247 
1248  // Apply transformation (HGRAD version is just a copy operation)
1249  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1250  if(element_space == PureBasis::HCURL){
1251 
1252  auto s_jac_inv = Kokkos::subview(cubature_jacobian_inverse_.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1253 
1254  fst::HCURLtransformVALUE(s_aux,s_jac_inv,cell_basis_ref_vector.get_view());
1255  } else if(element_space == PureBasis::HDIV){
1256 
1257  auto s_jac = Kokkos::subview(cubature_jacobian_.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1258  auto s_jac_det = Kokkos::subview(cubature_jacobian_determinant_.get_view(), cell_range, Kokkos::ALL());
1259 
1260  fst::HDIVtransformVALUE(s_aux,s_jac, s_jac_det, cell_basis_ref_vector.get_view());
1261  }
1262 
1263  PHX::Device().fence();
1264 
1265  } else {
1266 
1267  // This is ugly. The algorithm is restricted to host/serial due
1268  // to intrepid tools that requiring uniform reference
1269  // representation. For DG, CVFEM and sidesets this reference is
1270  // nonuniform.
1271 
1272  // Local allocation used for each cell
1273  auto cell_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_basis_vector",1,num_card,num_points,num_dim);
1274  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
1275  auto cell_jac_det = af.buildStaticArray<Scalar,Cell,IP>("cell_jac_det",1,num_points);
1276  auto cell_jac = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac",1,num_points,num_dim,num_dim);
1277  auto cell_jac_inv = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac_inv",1,num_points,num_dim,num_dim);
1278 
1279  // The array factory is difficult to extend to host space
1280  // without extra template magic and changing a ton of code in a
1281  // non-backwards compatible way, so we use some of the arrays
1282  // above only to get derivative array sized correctly and then
1283  // create the mirror on host.
1284  auto cell_basis_vector_host = Kokkos::create_mirror_view(cell_basis_vector.get_view());
1285  auto cell_cub_points_host = Kokkos::create_mirror_view(cell_cub_points.get_view());
1286  auto cell_jac_det_host = Kokkos::create_mirror_view(cell_jac_det.get_view());
1287  auto cell_jac_host = Kokkos::create_mirror_view(cell_jac.get_view());
1288  auto cell_jac_inv_host = Kokkos::create_mirror_view(cell_jac_inv.get_view());
1289  auto cell_basis_ref_vector_host = Kokkos::create_mirror_view(cell_basis_ref_vector.get_view());
1290  auto cubature_points_ref_host = Kokkos::create_mirror_view(cubature_points_ref_.get_view());
1291  Kokkos::deep_copy(cubature_points_ref_host,cubature_points_ref_.get_view());
1292  auto cubature_jacobian_host = Kokkos::create_mirror_view(cubature_jacobian_.get_view());
1293  Kokkos::deep_copy(cubature_jacobian_host,cubature_jacobian_.get_view());
1294  auto cubature_jacobian_inverse_host = Kokkos::create_mirror_view(cubature_jacobian_inverse_.get_view());
1295  Kokkos::deep_copy(cubature_jacobian_inverse_host,cubature_jacobian_inverse_.get_view());
1296  auto cubature_jacobian_determinant_host = Kokkos::create_mirror_view(cubature_jacobian_determinant_.get_view());
1297  Kokkos::deep_copy(cubature_jacobian_determinant_host,cubature_jacobian_determinant_.get_view());
1298  auto tmp_basis_vector_host = Kokkos::create_mirror_view(tmp_basis_vector.get_view());
1299 
1300  // We have to iterate through cells and apply a separate reference representation for each
1301  for(int cell=0; cell<num_evaluate_cells_; ++cell){
1302 
1303  // Load external into cell-local arrays
1304  for(int p=0;p<num_points;++p)
1305  for(int d=0;d<num_dim;++d)
1306  cell_cub_points_host(p,d)=cubature_points_ref_host(cell,p,d);
1307  Kokkos::deep_copy(cell_cub_points.get_view(),cell_cub_points_host);
1308 
1309  // Convert to mesh space. The intrepid_basis is virtual and
1310  // built in another class, short of adding a "clone on host"
1311  // method, we will have to make this call on device. This is a
1312  // major performance hit.
1313  intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE);
1314  Kokkos::deep_copy(cell_basis_ref_vector_host,cell_basis_ref_vector.get_view());
1315 
1316  using fst=Intrepid2::FunctionSpaceTools<Kokkos::HostSpace>;
1317 
1318  if(element_space == PureBasis::HCURL){
1319  // Need the jacobian inverse for HCurl
1320  for(int p=0;p<num_points;++p)
1321  for(int i=0; i<num_dim; ++i)
1322  for(int j=0; j<num_dim; ++j)
1323  cell_jac_inv_host(0,p,i,j)=cubature_jacobian_inverse_host(cell,p,i,j);
1324 
1325  fst::HCURLtransformVALUE(cell_basis_vector_host,cell_jac_inv_host,cell_basis_ref_vector_host);
1326  } else {
1327  // Need the jacobian for HDiv
1328  for(int p=0;p<num_points;++p)
1329  for(int i=0; i<num_dim; ++i)
1330  for(int j=0; j<num_dim; ++j)
1331  cell_jac_host(0,p,i,j)=cubature_jacobian_host(cell,p,i,j);
1332 
1333  for(int p=0;p<num_points;++p)
1334  cell_jac_det_host(0,p)=cubature_jacobian_determinant_host(cell,p);
1335 
1336  fst::HDIVtransformVALUE(cell_basis_vector_host,cell_jac_host,cell_jac_det_host,cell_basis_ref_vector_host);
1337  }
1338 
1339  // Copy cell quantity back into main array
1340  for(int b=0; b<num_card; ++b)
1341  for(int p=0; p<num_points; ++p)
1342  for(int d=0; d<num_dim; ++d)
1343  tmp_basis_vector_host(cell,b,p,d) = cell_basis_vector_host(0,b,p,d);
1344 
1345  Kokkos::deep_copy(tmp_basis_vector.get_view(),tmp_basis_vector_host);
1346  }
1347  }
1348 
1349  if(orientations_ != Teuchos::null)
1350  applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_basis_vector.get_view(), *orientations_->getOrientations(), *intrepid_basis);
1351 
1352  // Store for later if cache is enabled
1353  PANZER_CACHE_DATA(basis_vector);
1354 
1355  return tmp_basis_vector;
1356 
1357  }
1358 
1359 }
1360 
1361 template <typename Scalar>
1364 getGradBasisValues(const bool weighted,
1365  const bool cache,
1366  const bool force) const
1367 {
1368  if(weighted){
1369  if(weighted_grad_basis_evaluated_ and not force)
1370  return weighted_grad_basis;
1371  } else
1372  if(grad_basis_evaluated_ and not force)
1373  return grad_basis;
1374 
1375  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
1376 
1377  const int num_cells = num_cells_;
1378  const int num_points = basis_layout->numPoints();
1379  const int num_card = basis_layout->cardinality();
1380  const int num_dim = basis_layout->dimension();
1381 
1382  if(weighted){
1383 
1384  TEUCHOS_ASSERT(cubature_weights_.size() > 0);
1385 
1386  // Get the basis_scalar values - do not cache
1387  const auto bv = getGradBasisValues(false, force);
1388 
1389  // Apply the weighted measure (cubature weights)
1390  auto tmp_weighted_grad_basis = af.buildStaticArray<Scalar, Cell, BASIS, IP, Dim>("weighted_grad_basis", num_cells, num_card, num_points, num_dim);
1391 
1392  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1393  auto s_aux = Kokkos::subview(tmp_weighted_grad_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1394  auto s_cw = Kokkos::subview(cubature_weights_.get_view(), cell_range, Kokkos::ALL());
1395  auto s_bv = Kokkos::subview(bv.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1396 
1397  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1398  fst::multiplyMeasure(s_aux,s_cw,s_bv);
1399 
1400  // Store for later if cache is enabled
1401  PANZER_CACHE_DATA(weighted_grad_basis);
1402 
1403  return tmp_weighted_grad_basis;
1404 
1405  } else {
1406 
1407  TEUCHOS_ASSERT(cubature_jacobian_inverse_.size() > 0);
1408 
1409  const auto element_space = getElementSpace();
1410  TEUCHOS_ASSERT(element_space == PureBasis::CONST || element_space == PureBasis::HGRAD);
1411 
1412  auto cell_grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("cell_grad_basis_ref",num_card,num_points,num_dim);
1413  auto tmp_grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis_scalar",num_cells,num_card,num_points,num_dim);
1414 
1415  if(hasUniformReferenceSpace()){
1416 
1417  // Apply a single reference representation to all cells
1418  intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_GRAD);
1419 
1420  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1421  auto s_aux = Kokkos::subview(tmp_grad_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1422  auto s_jac_inv = Kokkos::subview(cubature_jacobian_inverse_.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1423 
1424  // Apply transformation
1425  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1426  fst::HGRADtransformGRAD(s_aux, s_jac_inv,cell_grad_basis_ref.get_view());
1427 
1428  PHX::Device().fence();
1429 
1430  } else {
1431 
1432  // This is ugly. The algorithm is restricted to host/serial due
1433  // to intrepid tools that requiring uniform reference
1434  // representation. For DG, CVFEM and sidesets this reference is
1435  // nonuniform.
1436 
1437  // Local allocation used for each cell
1438  auto cell_grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_grad_basis",1,num_card,num_points,num_dim);
1439  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
1440  auto cell_jac_inv = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac_inv",1,num_points,num_dim,num_dim);
1441 
1442  // The array factory is difficult to extend to host space
1443  // without extra template magic and changing a ton of code in a
1444  // non-backwards compatible way, so we use some of the arrays
1445  // above only to get derivative array sized correctly and then
1446  // create the mirror on host.
1447 
1448  // auto cell_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("cell_basis_ref_scalar",num_card,num_points,num_dim);
1449 
1450  auto cell_cub_points_host = Kokkos::create_mirror_view(cell_cub_points.get_view());
1451  auto cubature_points_ref_host = Kokkos::create_mirror_view(cubature_points_ref_.get_view());
1452  Kokkos::deep_copy(cubature_points_ref_host,cubature_points_ref_.get_view());
1453  auto cell_jac_inv_host = Kokkos::create_mirror_view(cell_jac_inv.get_view());
1454  auto cubature_jacobian_inverse_host = Kokkos::create_mirror_view(cubature_jacobian_inverse_.get_view());
1455  Kokkos::deep_copy(cubature_jacobian_inverse_host,cubature_jacobian_inverse_.get_view());
1456  auto cell_grad_basis_ref_host = Kokkos::create_mirror_view(cell_grad_basis_ref.get_view());
1457  auto cell_grad_basis_host = Kokkos::create_mirror_view(cell_grad_basis.get_view());
1458  auto tmp_grad_basis_host = Kokkos::create_mirror_view(tmp_grad_basis.get_view());
1459 
1460  // We have to iterate through cells and apply a separate reference representation for each
1461  for(int cell=0; cell<num_evaluate_cells_; ++cell){
1462 
1463  // Load external into cell-local arrays
1464  for(int p=0;p<num_points;++p)
1465  for(int d=0;d<num_dim;++d)
1466  cell_cub_points_host(p,d)=cubature_points_ref_host(cell,p,d);
1467 
1468  // Convert to mesh space. The intrepid_basis is virtual and
1469  // built in another class, short of adding a "clone on host"
1470  // method, we will have to make this call on device. This is a
1471  // major performance hit.
1472  Kokkos::deep_copy(cell_cub_points.get_view(),cell_cub_points_host);
1473  intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_GRAD);
1474  Kokkos::deep_copy(cell_grad_basis_ref_host,cell_grad_basis_ref.get_view());
1475 
1476  for(int p=0;p<num_points;++p)
1477  for(int d=0;d<num_dim;++d)
1478  for(int d2=0;d2<num_dim;++d2)
1479  cell_jac_inv_host(0,p,d,d2)=cubature_jacobian_inverse_host(cell,p,d,d2);
1480 
1481  using fst=Intrepid2::FunctionSpaceTools<Kokkos::HostSpace>;
1482  fst::HGRADtransformGRAD(cell_grad_basis_host,cell_jac_inv_host,cell_grad_basis_ref_host);
1483  // PHX::Device().fence();
1484 
1485  // Copy cell quantity back into main array
1486  for(int b=0; b<num_card; ++b)
1487  for(int p=0; p<num_points; ++p)
1488  for(int d=0; d<num_dim; ++d)
1489  tmp_grad_basis_host(cell,b,p,d) = cell_grad_basis_host(0,b,p,d);
1490  Kokkos::deep_copy(tmp_grad_basis.get_view(),tmp_grad_basis_host);
1491  }
1492  }
1493 
1494  if(orientations_ != Teuchos::null)
1495  applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_grad_basis.get_view(), *orientations_->getOrientations(), *intrepid_basis);
1496 
1497  // Store for later if cache is enabled
1498  PANZER_CACHE_DATA(grad_basis);
1499 
1500  return tmp_grad_basis;
1501 
1502  }
1503 
1504 }
1505 
1506 template <typename Scalar>
1509 getCurl2DVectorBasis(const bool weighted,
1510  const bool cache,
1511  const bool force) const
1512 {
1513  if(weighted){
1514  if(weighted_curl_basis_scalar_evaluated_ and not force)
1515  return weighted_curl_basis_scalar;
1516  } else
1517  if(curl_basis_scalar_evaluated_ and not force)
1518  return curl_basis_scalar;
1519 
1520  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
1521 
1522  const int num_cells = num_cells_;
1523  const int num_points = basis_layout->numPoints();
1524  const int num_card = basis_layout->cardinality();
1525  const int num_dim = basis_layout->dimension();
1526 
1527  if(weighted){
1528 
1529  TEUCHOS_ASSERT(cubature_weights_.size() > 0);
1530 
1531  // Get the basis_scalar values - do not cache
1532  const auto bv = getCurl2DVectorBasis(false, force);
1533 
1534  // Apply the weighted measure (cubature weights)
1535  auto tmp_weighted_curl_basis_scalar = af.buildStaticArray<Scalar, Cell, BASIS, IP>("weighted_curl_basis_scalar", num_cells, num_card, num_points);
1536 
1537  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1538  auto s_aux = Kokkos::subview(tmp_weighted_curl_basis_scalar.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1539  auto s_cw = Kokkos::subview(cubature_weights_.get_view(), cell_range, Kokkos::ALL());
1540  auto s_bv = Kokkos::subview(bv.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1541 
1542  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1543  fst::multiplyMeasure(s_aux,s_cw,s_bv);
1544 
1545  // Store for later if cache is enabled
1546  PANZER_CACHE_DATA(weighted_curl_basis_scalar);
1547 
1548  return tmp_weighted_curl_basis_scalar;
1549 
1550  } else {
1551 
1552  TEUCHOS_ASSERT(cubature_jacobian_determinant_.size() > 0);
1553  TEUCHOS_ASSERT(num_dim == 2);
1554 
1555  const auto element_space = getElementSpace();
1556  TEUCHOS_ASSERT(element_space == PureBasis::HCURL);
1557 
1558  auto cell_curl_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("cell_curl_basis_ref_scalar",num_card,num_points);
1559  auto tmp_curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("curl_basis_scalar",num_cells,num_card,num_points);
1560 
1561  if(hasUniformReferenceSpace()){
1562 
1563  intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_CURL);
1564 
1565  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1566  auto s_aux = Kokkos::subview(tmp_curl_basis_scalar.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1567  auto s_jac_det = Kokkos::subview(cubature_jacobian_determinant_.get_view(), cell_range, Kokkos::ALL());
1568 
1569  // note only volume deformation is needed!
1570  // this relates directly to this being in
1571  // the divergence space in 2D!
1572  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1573  fst::HDIVtransformDIV(s_aux,s_jac_det,cell_curl_basis_ref_scalar.get_view());
1574  PHX::Device().fence();
1575 
1576  } else {
1577 
1578  // This is ugly. The algorithm is restricted to host/serial due
1579  // to intrepid tools that requiring uniform reference
1580  // representation. For DG, CVFEM and sidesets this reference is
1581  // nonuniform.
1582 
1583  // Local allocation used for each cell
1584  auto cell_curl = af.buildStaticArray<Scalar,Cell,BASIS,IP>("cell_curl",1,num_card,num_points);
1585  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
1586  auto cell_jac_det = af.buildStaticArray<Scalar,Cell,IP>("cell_jac_det",1,num_points);
1587 
1588  // The array factory is difficult to extend to host space
1589  // without extra template magic and changing a ton of code in a
1590  // non-backwards compatible way, so we use some of the arrays
1591  // above only to get derivative array sized correctly and then
1592  // create the mirror on host.
1593  auto cell_cub_points_host = Kokkos::create_mirror_view(cell_cub_points.get_view());
1594  auto cubature_points_ref_host = Kokkos::create_mirror_view(cubature_points_ref_.get_view());
1595  Kokkos::deep_copy(cubature_points_ref_host,cubature_points_ref_.get_view());
1596  auto cell_jac_det_host = Kokkos::create_mirror_view(cell_jac_det.get_view());
1597  auto cubature_jacobian_determinant_host = Kokkos::create_mirror_view(cubature_jacobian_determinant_.get_view());
1598  Kokkos::deep_copy(cubature_jacobian_determinant_host,cubature_jacobian_determinant_.get_view());
1599  auto cell_curl_basis_ref_scalar_host = Kokkos::create_mirror_view(cell_curl_basis_ref_scalar.get_view());
1600  auto cell_curl_host = Kokkos::create_mirror_view(cell_curl.get_view());
1601  auto tmp_curl_basis_scalar_host = Kokkos::create_mirror_view(tmp_curl_basis_scalar.get_view());
1602 
1603  // We have to iterate through cells and apply a separate reference representation for each
1604  for(int cell=0; cell<num_evaluate_cells_; ++cell){
1605 
1606  // Load external into cell-local arrays
1607  for(int p=0;p<num_points;++p)
1608  for(int d=0;d<num_dim;++d)
1609  cell_cub_points_host(p,d)=cubature_points_ref_host(cell,p,d);
1610  for(int p=0;p<num_points;++p)
1611  cell_jac_det_host(0,p)=cubature_jacobian_determinant_host(cell,p);
1612 
1613 
1614  // Convert to mesh space. The intrepid_basis is virtual and
1615  // built in another class, short of adding a "clone on host"
1616  // method, we will have to make this call on device. This is a
1617  // major performance hit.
1618  Kokkos::deep_copy(cell_cub_points.get_view(),cell_cub_points_host);
1619  intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL);
1620  Kokkos::deep_copy(cell_curl_basis_ref_scalar_host,cell_curl_basis_ref_scalar.get_view());
1621 
1622  // note only volume deformation is needed!
1623  // this relates directly to this being in
1624  // the divergence space in 2D!
1625  using fst=Intrepid2::FunctionSpaceTools<Kokkos::HostSpace>;
1626  fst::HDIVtransformDIV(cell_curl_host,cell_jac_det_host,cell_curl_basis_ref_scalar_host);
1627  PHX::Device().fence();
1628 
1629  // Copy cell quantity back into main array
1630  for(int b=0; b<num_card; ++b)
1631  for(int p=0; p<num_points; ++p)
1632  tmp_curl_basis_scalar_host(cell,b,p) = cell_curl_host(0,b,p);
1633  Kokkos::deep_copy(tmp_curl_basis_scalar.get_view(),tmp_curl_basis_scalar_host);
1634  }
1635  }
1636 
1637  if(orientations_ != Teuchos::null)
1638  applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_curl_basis_scalar.get_view(), *orientations_->getOrientations(), *intrepid_basis);
1639 
1640  // Store for later if cache is enabled
1641  PANZER_CACHE_DATA(curl_basis_scalar);
1642 
1643  return tmp_curl_basis_scalar;
1644 
1645  }
1646 
1647 }
1648 
1649 template <typename Scalar>
1652 getCurlVectorBasis(const bool weighted,
1653  const bool cache,
1654  const bool force) const
1655 {
1656  if(weighted){
1657  if(weighted_curl_basis_vector_evaluated_ and not force)
1658  return weighted_curl_basis_vector;
1659  } else
1660  if(curl_basis_vector_evaluated_ and not force)
1661  return curl_basis_vector;
1662 
1663  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
1664 
1665  const int num_cells = num_cells_;
1666  const int num_points = basis_layout->numPoints();
1667  const int num_card = basis_layout->cardinality();
1668  const int num_dim = basis_layout->dimension();
1669 
1670  if(weighted){
1671 
1672  TEUCHOS_ASSERT(cubature_weights_.size() > 0);
1673 
1674  // Get the basis_scalar values - do not cache
1675  const auto bv = getCurlVectorBasis(false, force);
1676 
1677  // Apply the weighted measure (cubature weights)
1678  auto tmp_weighted_curl_basis_vector = af.buildStaticArray<Scalar, Cell, BASIS, IP, Dim>("weighted_curl_basis_vector", num_cells, num_card, num_points, num_dim);
1679 
1680  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1681  auto s_aux = Kokkos::subview(tmp_weighted_curl_basis_vector.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1682  auto s_cw = Kokkos::subview(cubature_weights_.get_view(), cell_range, Kokkos::ALL());
1683  auto s_bv = Kokkos::subview(bv.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1684 
1685  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1686  fst::multiplyMeasure(s_aux, s_cw, s_bv);
1687 
1688  // Store for later if cache is enabled
1689  PANZER_CACHE_DATA(weighted_curl_basis_vector);
1690 
1691  return tmp_weighted_curl_basis_vector;
1692 
1693  } else {
1694 
1695  TEUCHOS_ASSERT(cubature_jacobian_determinant_.size() > 0);
1696  TEUCHOS_ASSERT(num_dim == 3);
1697 
1698  const auto element_space = getElementSpace();
1699  TEUCHOS_ASSERT(element_space == PureBasis::HCURL);
1700 
1701  auto cell_curl_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("cell_curl_basis_ref_vector",num_card,num_points,num_dim);
1702  auto tmp_curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("curl_basis_vector",num_cells,num_card,num_points,num_dim);
1703 
1704  if(hasUniformReferenceSpace()){
1705 
1706  intrepid_basis->getValues(cell_curl_basis_ref_vector.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_CURL);
1707 
1708  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1709  auto s_aux = Kokkos::subview(tmp_curl_basis_vector.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1710  auto s_jac = Kokkos::subview(cubature_jacobian_.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1711  auto s_jac_det = Kokkos::subview(cubature_jacobian_determinant_.get_view(), cell_range, Kokkos::ALL());
1712 
1713  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1714  fst::HCURLtransformCURL(s_aux, s_jac, s_jac_det,cell_curl_basis_ref_vector.get_view());
1715  PHX::Device().fence();
1716 
1717  } else {
1718 
1719  // This is ugly. The algorithm is restricted to host/serial due
1720  // to intrepid tools that requiring uniform reference
1721  // representation. For DG, CVFEM and sidesets this reference is
1722  // nonuniform.
1723 
1724  // Local allocation used for each cell
1725  auto cell_curl = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_curl",1,num_card,num_points,num_dim);
1726  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
1727  auto cell_jac = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac",1,num_points,num_dim,num_dim);
1728  auto cell_jac_det = af.buildStaticArray<Scalar,Cell,IP>("cell_jac_det",1,num_points);
1729 
1730  // The array factory is difficult to extend to host space
1731  // without extra template magic and changing a ton of code in a
1732  // non-backwards compatible way, so we use some of the arrays
1733  // above only to get derivative array sized correctly and then
1734  // create the mirror on host.
1735  auto cell_cub_points_host = Kokkos::create_mirror_view(cell_cub_points.get_view());
1736  auto cubature_points_ref_host = Kokkos::create_mirror_view(cubature_points_ref_.get_view());
1737  Kokkos::deep_copy(cubature_points_ref_host,cubature_points_ref_.get_view());
1738  auto cell_jac_host = Kokkos::create_mirror_view(cell_jac.get_view());
1739  auto cubature_jacobian_host = Kokkos::create_mirror_view(cubature_jacobian_.get_view());
1740  Kokkos::deep_copy(cubature_jacobian_host,cubature_jacobian_.get_view());
1741  auto cell_jac_det_host = Kokkos::create_mirror_view(cell_jac_det.get_view());
1742  auto cubature_jacobian_determinant_host = Kokkos::create_mirror_view(cubature_jacobian_determinant_.get_view());
1743  Kokkos::deep_copy(cubature_jacobian_determinant_host,cubature_jacobian_determinant_.get_view());
1744  auto cell_curl_basis_ref_vector_host = Kokkos::create_mirror_view(cell_curl_basis_ref_vector.get_view());
1745  auto cell_curl_host = Kokkos::create_mirror_view(cell_curl.get_view());
1746  auto tmp_curl_basis_vector_host = Kokkos::create_mirror_view(tmp_curl_basis_vector.get_view());
1747 
1748  // We have to iterate through cells and apply a separate reference representation for each
1749  for(int cell=0; cell<num_evaluate_cells_; ++cell){
1750 
1751  // Load external into cell-local arrays
1752  for(int p=0;p<num_points;++p)
1753  for(int d=0;d<num_dim;++d)
1754  cell_cub_points_host(p,d)=cubature_points_ref_host(cell,p,d);
1755  for(int p=0;p<num_points;++p)
1756  for(int d=0;d<num_dim;++d)
1757  for(int d2=0;d2<num_dim;++d2)
1758  cell_jac_host(0,p,d,d2)=cubature_jacobian_host(cell,p,d,d2);
1759  for(int p=0;p<num_points;++p)
1760  cell_jac_det_host(0,p)=cubature_jacobian_determinant_host(cell,p);
1761 
1762  // Convert to mesh space. The intrepid_basis is virtual and
1763  // built in another class, short of adding a "clone on host"
1764  // method, we will have to make this call on device. This is a
1765  // major performance hit.
1766  Kokkos::deep_copy(cell_cub_points.get_view(),cell_cub_points_host);
1767  intrepid_basis->getValues(cell_curl_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL);
1768  Kokkos::deep_copy(cell_curl_basis_ref_vector_host,cell_curl_basis_ref_vector.get_view());
1769 
1770  using fst=Intrepid2::FunctionSpaceTools<Kokkos::HostSpace>;
1771  fst::HCURLtransformCURL(cell_curl_host,cell_jac_host,cell_jac_det_host,cell_curl_basis_ref_vector_host);
1772  // PHX::Device().fence();
1773 
1774  // Copy cell quantity back into main array
1775  for(int b=0; b<num_card; ++b)
1776  for(int p=0; p<num_points; ++p)
1777  for(int d=0;d<num_dim;++d)
1778  tmp_curl_basis_vector_host(cell,b,p,d) = cell_curl_host(0,b,p,d);
1779  Kokkos::deep_copy(tmp_curl_basis_vector.get_view(),tmp_curl_basis_vector_host);
1780  }
1781  }
1782 
1783  if(orientations_ != Teuchos::null)
1784  applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_curl_basis_vector.get_view(), *orientations_->getOrientations(), *intrepid_basis);
1785 
1786  // Store for later if cache is enabled
1787  PANZER_CACHE_DATA(curl_basis_vector);
1788 
1789  return tmp_curl_basis_vector;
1790 
1791  }
1792 
1793 }
1794 
1795 template <typename Scalar>
1798 getDivVectorBasis(const bool weighted,
1799  const bool cache,
1800  const bool force) const
1801 {
1802  if(weighted){
1803  if(weighted_div_basis_evaluated_ and not force)
1804  return weighted_div_basis;
1805  } else
1806  if(div_basis_evaluated_ and not force)
1807  return div_basis;
1808 
1809  MDFieldArrayFactory af(prefix,getExtendedDimensions(),true);
1810 
1811  const int num_cells = num_cells_;
1812  const int num_points = basis_layout->numPoints();
1813  const int num_card = basis_layout->cardinality();
1814  const int num_dim = basis_layout->dimension();
1815 
1816  if(weighted){
1817 
1818  TEUCHOS_ASSERT(cubature_weights_.size() > 0);
1819 
1820  // Get the basis_scalar values - do not cache
1821  const auto bv = getDivVectorBasis(false, force);
1822 
1823  // Apply the weighted measure (cubature weights)
1824  auto tmp_weighted_div_basis = af.buildStaticArray<Scalar, Cell, BASIS, IP>("weighted_div_basis", num_cells, num_card, num_points);
1825 
1826  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1827  auto s_aux = Kokkos::subview(tmp_weighted_div_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1828  auto s_cw = Kokkos::subview(cubature_weights_.get_view(), cell_range, Kokkos::ALL());
1829  auto s_bv = Kokkos::subview(bv.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1830 
1831  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1832  fst::multiplyMeasure(s_aux, s_cw, s_bv);
1833 
1834  // Store for later if cache is enabled
1835  PANZER_CACHE_DATA(weighted_div_basis);
1836 
1837  return tmp_weighted_div_basis;
1838 
1839  } else {
1840 
1841  TEUCHOS_ASSERT(cubature_jacobian_determinant_.size() > 0);
1842 
1843  const auto element_space = getElementSpace();
1844  TEUCHOS_ASSERT(element_space == PureBasis::HDIV);
1845 
1846  auto cell_div_basis_ref = af.buildStaticArray<Scalar,BASIS,IP>("cell_div_basis_ref",num_card,num_points);
1847  auto tmp_div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("div_basis",num_cells,num_card,num_points);
1848 
1849  if(hasUniformReferenceSpace()){
1850 
1851  intrepid_basis->getValues(cell_div_basis_ref.get_view(),cubature_points_uniform_ref_.get_view(),Intrepid2::OPERATOR_DIV);
1852 
1853  const std::pair<int,int> cell_range(0,num_evaluate_cells_);
1854  auto s_aux = Kokkos::subview(tmp_div_basis.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1855  auto s_jac_det = Kokkos::subview(cubature_jacobian_determinant_.get_view(), cell_range, Kokkos::ALL());
1856 
1857  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1858  fst::HDIVtransformDIV(s_aux,s_jac_det,cell_div_basis_ref.get_view());
1859  PHX::Device().fence();
1860 
1861  } else {
1862 
1863  // This is ugly. The algorithm is restricted to host/serial due
1864  // to intrepid tools that requiring uniform reference
1865  // representation. For DG, CVFEM and sidesets this reference is
1866  // nonuniform.
1867 
1868  // Local allocation used for each cell
1869  auto cell_div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("cell_div_basis",1,num_card,num_points);
1870  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
1871  auto cell_jac_det = af.buildStaticArray<Scalar,Cell,IP>("cell_jac_det",1,num_points);
1872 
1873  // The array factory is difficult to extend to host space
1874  // without extra template magic and changing a ton of code in a
1875  // non-backwards compatible way, so we use some of the arrays
1876  // above only to get derivative array sized correctly and then
1877  // create the mirror on host.
1878  auto cell_cub_points_host = Kokkos::create_mirror_view(cell_cub_points.get_view());
1879  auto cubature_points_ref_host = Kokkos::create_mirror_view(cubature_points_ref_.get_view());
1880  Kokkos::deep_copy(cubature_points_ref_host,cubature_points_ref_.get_view());
1881  auto cell_jac_det_host = Kokkos::create_mirror_view(cell_jac_det.get_view());
1882  auto cubature_jacobian_determinant_host = Kokkos::create_mirror_view(cubature_jacobian_determinant_.get_view());
1883  Kokkos::deep_copy(cubature_jacobian_determinant_host,cubature_jacobian_determinant_.get_view());
1884  auto cell_div_basis_ref_host = Kokkos::create_mirror_view(cell_div_basis_ref.get_view());
1885  auto cell_div_basis_host = Kokkos::create_mirror_view(cell_div_basis.get_view());
1886  auto tmp_div_basis_host = Kokkos::create_mirror_view(tmp_div_basis.get_view());
1887 
1888  // We have to iterate through cells and apply a separate reference representation for each
1889  for(int cell=0; cell<num_evaluate_cells_; ++cell){
1890 
1891  // Load external into cell-local arrays
1892  for(int p=0;p<num_points;++p)
1893  for(int d=0;d<num_dim;++d)
1894  cell_cub_points_host(p,d)=cubature_points_ref_host(cell,p,d);
1895  for(int p=0;p<num_points;++p)
1896  cell_jac_det_host(0,p)=cubature_jacobian_determinant_host(cell,p);
1897 
1898  // Convert to mesh space. The intrepid_basis is virtual and
1899  // built in another class, short of adding a "clone on host"
1900  // method, we will have to make this call on device. This is a
1901  // major performance hit.
1902  Kokkos::deep_copy(cell_cub_points.get_view(),cell_cub_points_host);
1903  intrepid_basis->getValues(cell_div_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_DIV);
1904  Kokkos::deep_copy(cell_div_basis_ref_host,cell_div_basis_ref.get_view());
1905 
1906  using fst=Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>;
1907  fst::HDIVtransformDIV(cell_div_basis.get_view(),cell_jac_det.get_view(),cell_div_basis_ref.get_view());
1908  Kokkos::deep_copy(cell_div_basis_host, cell_div_basis.get_static_view());
1909 
1910 
1911  // Copy cell quantity back into main array
1912  for(int b=0; b<num_card; ++b)
1913  for(int p=0; p<num_points; ++p)
1914  tmp_div_basis_host(cell,b,p) = cell_div_basis_host(0,b,p);
1915  Kokkos::deep_copy(tmp_div_basis.get_view(),tmp_div_basis_host);
1916  }
1917  }
1918 
1919  if(orientations_ != Teuchos::null)
1920  applyOrientationsImpl<Scalar>(num_orientations_cells_, tmp_div_basis.get_view(), *orientations_->getOrientations(), *intrepid_basis);
1921 
1922  // Store for later if cache is enabled
1923  PANZER_CACHE_DATA(div_basis);
1924 
1925  return tmp_div_basis;
1926 
1927  }
1928 
1929 }
1930 
1931 //=======================================================================================================
1932 
1933 } // namespace panzer
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
void setCellVertexCoordinates(PHX::MDField< Scalar, Cell, NODE, Dim > vertex_coordinates)
Set the cell vertex coordinates (required for getBasisCoordinates())
ConstArray_BasisIPDim getVectorBasisValuesRef(const bool cache=true, const bool force=false) const
Get the vector basis values evaluated at reference points.
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, IP, Dim > &cell_cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim > &jac, const PHX::MDField< Scalar, Cell, IP > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim > &jac_inv)
PHX::MDField< const Scalar, Cell, BASIS, IP > ConstArray_CellBasisIP
bool basis_ref_scalar_evaluated_
Used to check if arrays have been cached.
void setWeightedMeasure(PHX::MDField< Scalar, Cell, IP > weighted_measure)
Set the cubature weights (weighted measure) for the basis values object - required to get weighted ba...
PHX::MDField< const Scalar, Cell, BASIS, Dim > ConstArray_CellBasisDim
void applyOrientations(const PHX::MDField< const Scalar, Cell, BASIS > &orientations)
Method to apply orientations to a basis values container.
PHX::MDField< const Scalar, BASIS, IP > ConstArray_BasisIP
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.
void setOrientations(const Teuchos::RCP< const OrientationsInterface > &orientations, const int num_orientations_cells=-1)
Set the orientations object for applying orientations using the lazy evaluation path - required for c...
void evaluateBasisCoordinates(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int in_num_cells=-1)
ConstArray_CellBasisIP getDivVectorBasis(const bool weighted, const bool cache=true, const bool force=false) const
Get the divergence of a vector basis evaluated at mesh points.
ConstArray_CellBasisIPDim getVectorBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the vector basis values evaluated at mesh points.
ConstArray_CellBasisIPDim getGradBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the gradient of the basis evaluated at mesh points.
ConstArray_CellBasisIPDim getCurlVectorBasis(const bool weighted, const bool cache=true, const bool force=false) const
Get the curl of a 3D vector basis evaluated at mesh points.
void setupUniform(const Teuchos::RCP< const panzer::BasisIRLayout > &basis, PHX::MDField< Scalar, IP, Dim > reference_points, PHX::MDField< Scalar, Cell, IP, Dim, Dim > point_jacobian, PHX::MDField< Scalar, Cell, IP > point_jacobian_determinant, PHX::MDField< Scalar, Cell, IP, Dim, Dim > point_jacobian_inverse, const int num_evaluated_cells=-1)
Setup for lazy evaluation for uniform point layout.
ConstArray_BasisIP getDivVectorBasisRef(const bool cache=true, const bool force=false) const
Get the divergence of a vector basis evaluated at reference points.
#define PANZER_CACHE_DATA(name)
void setupArrays(const Teuchos::RCP< const panzer::BasisIRLayout > &basis, bool computeDerivatives=true)
Sizes/allocates memory for arrays.
ConstArray_BasisDim getBasisCoordinatesRef(const bool cache=true, const bool force=false) const
Get the reference coordinates for basis.
PHX::MDField< const Scalar, Cell, BASIS, IP, Dim > ConstArray_CellBasisIPDim
Intrepid2::Basis< PHX::Device::execution_space, Scalar, Scalar > IntrepidBasis
PHX::MDField< const Scalar, BASIS, IP, Dim > ConstArray_BasisIPDim
BasisValues2(const std::string &prefix="", const bool allocArrays=false, const bool buildWeighted=false)
Main constructor.
ConstArray_BasisIPDim getGradBasisValuesRef(const bool cache=true, const bool force=false) const
Get the gradient of the basis evaluated at reference points.
ConstArray_BasisIP getCurl2DVectorBasisRef(const bool cache=true, const bool force=false) const
Get the curl of a vector basis evaluated at reference points.
ConstArray_CellBasisIP getCurl2DVectorBasis(const bool weighted, const bool cache=true, const bool force=false) const
Get the curl of a 2D vector basis evaluated at mesh points.
ConstArray_CellBasisIP getBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the basis values evaluated at mesh points.
ConstArray_CellBasisDim getBasisCoordinates(const bool cache=true, const bool force=false) const
Carterisan coordinates for basis coefficients in mesh space.
ConstArray_BasisIPDim getCurlVectorBasisRef(const bool cache=true, const bool force=false) const
Get the curl of a vector basis evaluated at reference points.
ConstArray_BasisIP getBasisValuesRef(const bool cache=true, const bool force=false) const
Get the basis values evaluated at reference points.
PureBasis::EElementSpace getElementSpace() const
PHX::MDField< const Scalar, BASIS, Dim > ConstArray_BasisDim
void setup(const Teuchos::RCP< const panzer::BasisIRLayout > &basis, PHX::MDField< Scalar, Cell, IP, Dim > reference_points, PHX::MDField< Scalar, Cell, IP, Dim, Dim > point_jacobian, PHX::MDField< Scalar, Cell, IP > point_jacobian_determinant, PHX::MDField< Scalar, Cell, IP, Dim, Dim > point_jacobian_inverse, const int num_evaluated_cells=-1)
Setup for lazy evaluation for non-uniform point layout.
void evaluateValues(const PHX::MDField< Scalar, IP, Dim > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim > &jac, const PHX::MDField< Scalar, Cell, IP > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim > &jac_inv, const int in_num_cells=-1)