Panzer  Version of the Day
Panzer_BasisValues2.cpp
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 #include "Panzer_BasisValues2.hpp"
46 
48 #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 
56 namespace panzer {
57 
58 template <typename Scalar>
60 evaluateValues(const PHX::MDField<Scalar,IP,Dim,void,void,void,void,void,void> & cub_points,
61  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
62  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
63  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv)
64 {
65  PHX::MDField<Scalar,Cell,IP> weighted_measure;
66  PHX::MDField<Scalar,Cell,NODE,Dim> vertex_coordinates;
67  build_weighted = false;
68  evaluateValues(cub_points,jac,jac_det,jac_inv,weighted_measure,vertex_coordinates,false);
69 }
70 
71 template <typename Scalar>
73 evaluateValues(const PHX::MDField<Scalar,IP,Dim,void,void,void,void,void,void> & cub_points,
74  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
75  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
76  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv,
77  const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
78  const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
79  bool use_vertex_coordinates)
80 {
81  MDFieldArrayFactory af("",ddims_,true);
82 
83  int num_dim = basis_layout->dimension();
84 
85  // currently this just copies on the basis objects are converted
86  // in intrepid
87  evaluateReferenceValues(cub_points,compute_derivatives,use_vertex_coordinates);
88 
89  PureBasis::EElementSpace elmtspace = getElementSpace();
90  if(elmtspace==PureBasis::CONST ||
91  elmtspace==PureBasis::HGRAD) {
92  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
93  HGRADtransformVALUE(basis_scalar.get_view(),
94  basis_ref_scalar.get_view());
95 
96  if(build_weighted) {
97  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
98  multiplyMeasure(weighted_basis_scalar.get_view(),
99  weighted_measure.get_view(),
100  basis_scalar.get_view());
101  }
102  }
103  else if(elmtspace==PureBasis::HCURL) {
104  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
105  HCURLtransformVALUE(basis_vector.get_view(),
106  jac_inv.get_view(),
107  basis_ref_vector.get_view());
108 
109  if(build_weighted) {
110  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
111  multiplyMeasure(weighted_basis_vector.get_view(),
112  weighted_measure.get_view(),
113  basis_vector.get_view());
114  }
115  }
116  else if(elmtspace==PureBasis::HDIV)
117  {
118  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
119  HDIVtransformVALUE(basis_vector.get_view(),
120  jac.get_view(),
121  jac_det.get_view(),
122  basis_ref_vector.get_view());
123 
124  if(build_weighted) {
125  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
126  multiplyMeasure(weighted_basis_vector.get_view(),
127  weighted_measure.get_view(),
128  basis_vector.get_view());
129  }
130  }
131  else { TEUCHOS_ASSERT(false); }
132 
133  if(elmtspace==PureBasis::HGRAD && compute_derivatives) {
134  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
135  HGRADtransformGRAD(grad_basis.get_view(),
136  jac_inv.get_view(),
137  grad_basis_ref.get_view());
138 
139  if(build_weighted) {
140  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
141  multiplyMeasure(weighted_grad_basis.get_view(),
142  weighted_measure.get_view(),
143  grad_basis.get_view());
144  }
145  }
146  else if(elmtspace==PureBasis::HCURL && num_dim==2 && compute_derivatives) {
147  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
148  HDIVtransformDIV(curl_basis_scalar.get_view(),
149  jac_det.get_view(), // note only volume deformation is needed!
150  // this relates directly to this being in
151  // the divergence space in 2D!
152  curl_basis_ref_scalar.get_view());
153 
154  if(build_weighted) {
155  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
156  multiplyMeasure(weighted_curl_basis_scalar.get_view(),
157  weighted_measure.get_view(),
158  curl_basis_scalar.get_view());
159  }
160  }
161  else if(elmtspace==PureBasis::HCURL && num_dim==3 && compute_derivatives) {
162  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
163  HCURLtransformCURL(curl_basis_vector.get_view(),
164  jac.get_view(),
165  jac_det.get_view(),
166  curl_basis_ref_vector.get_view());
167 
168  if(build_weighted) {
169  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
170  multiplyMeasure(weighted_curl_basis_vector.get_view(),
171  weighted_measure.get_view(),
172  curl_basis_vector.get_view());
173  }
174  }
175  else if(elmtspace==PureBasis::HDIV && compute_derivatives) {
176  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
177  HDIVtransformDIV(div_basis.get_view(),
178  jac_det.get_view(),
179  div_basis_ref.get_view());
180 
181  if(build_weighted) {
182  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
183  multiplyMeasure(weighted_div_basis.get_view(),
184  weighted_measure.get_view(),
185  div_basis.get_view());
186  }
187  }
188 
189  // If basis supports coordinate values at basis points, then
190  // compute these values
191  if(use_vertex_coordinates) {
192  // Teuchos::RCP<Intrepid2::DofCoordsInterface<ArrayDynamic> > coords
193  // = Teuchos::rcp_dynamic_cast<Intrepid2::DofCoordsInterface<ArrayDynamic> >(intrepid_basis);
194  // if (!Teuchos::is_null(coords)) {
195 /*
196  ArrayDynamic dyn_basis_coordinates_ref = af.buildArray<Scalar,BASIS,Dim>("basis_coordinates_ref",basis_coordinates_ref.extent(0),basis_coordinates_ref.extent(1));
197  coords->getDofCoords(dyn_basis_coordinates_ref);
198 
199  // fill in basis coordinates
200  for (size_type i = 0; i < basis_coordinates_ref.extent(0); ++i)
201  for (size_type j = 0; j < basis_coordinates_ref.extent(1); ++j)
202  basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j);
203 */
204 
205  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
206  cell_tools.mapToPhysicalFrame(basis_coordinates.get_view(),
207  basis_coordinates_ref.get_view(),
208  vertex_coordinates.get_view(),
209  intrepid_basis->getBaseCellTopology());
210  }
211 }
212 
213 
214 
215 
216 
217 
218 
219 template <typename Scalar>
221 evaluateBasisCoordinates(const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates)
222 {
223  MDFieldArrayFactory af("",ddims_,true);
224 
225  // intrepid_basis->getDofCoords requires DynRankView, but basis_coordinates_ref is more of a static View
226  // We use an auxiliary 'dyn' array to get around this
227  using coordsScalarType = typename Intrepid2::Basis<PHX::Device::execution_space,Scalar,Scalar>::scalarType;
228  auto dyn_basis_coordinates_ref = af.buildArray<coordsScalarType,BASIS,Dim>("basis_coordinates_ref",
229  basis_coordinates_ref.extent(0),
230  basis_coordinates_ref.extent(1));
231  intrepid_basis->getDofCoords(dyn_basis_coordinates_ref.get_view());
232 
233  // fill in basis coordinates
234  for (int i = 0; i < basis_coordinates_ref.extent_int(0); ++i)
235  for (int j = 0; j < basis_coordinates_ref.extent_int(1); ++j)
236  basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j);
237 
238  Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
239  cell_tools.mapToPhysicalFrame(basis_coordinates.get_view(),
240  basis_coordinates_ref.get_view(),
241  vertex_coordinates.get_view(),
242  intrepid_basis->getBaseCellTopology());
243 }
244 
245 
246 
247 
248 
249 template <typename Scalar>
251 evaluateValues(const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cub_points,
252  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
253  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
254  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv,
255  const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
256  const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
257  bool use_vertex_coordinates)
258 {
259 
260  PureBasis::EElementSpace elmtspace = getElementSpace();
261 
262  if(elmtspace == PureBasis::CONST){
263  evaluateValues_Const(cub_points,jac_inv,weighted_measure);
264  } else if(elmtspace == PureBasis::HGRAD){
265  evaluateValues_HGrad(cub_points,jac_inv,weighted_measure);
266  } else if(elmtspace == PureBasis::HCURL){
267  evaluateValues_HCurl(cub_points,jac,jac_det,jac_inv,weighted_measure);
268  } else if(elmtspace == PureBasis::HDIV){
269  evaluateValues_HDiv(cub_points,jac,jac_det,weighted_measure);
270  } else {
271  TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::evaluateValues : Element space not recognized.");
272  }
273 
274  if(use_vertex_coordinates) {
275  TEUCHOS_TEST_FOR_EXCEPT_MSG(elmtspace == PureBasis::CONST,"panzer::BasisValues2::evaluateValues : Const basis cannot have basis coordinates.");
276  evaluateBasisCoordinates(vertex_coordinates);
277  }
278 
279 }
280 
281 template <typename Scalar>
283 evaluateValues_Const(const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cub_points,
284  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv,
285  const PHX::MDField<Scalar,Cell,IP> & weighted_measure)
286 {
287 
288  TEUCHOS_ASSERT(getElementSpace() == PureBasis::CONST);
289 
290  typedef Intrepid2::FunctionSpaceTools<PHX::Device::execution_space> fst;
291  MDFieldArrayFactory af("",ddims_,true);
292 
293  const panzer::PureBasis & basis = *(basis_layout->getBasis());
294 
295  const int num_points = basis_layout->numPoints();
296  const int num_basis = basis.cardinality();
297  const int num_dim = basis_layout->dimension();
298  const int num_cells = basis_layout->numCells();
299 
300  auto cell_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("cell_basis_scalar",1,num_basis,num_points);
301  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
302  auto cell_grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_grad_basis",1,num_basis,num_points,num_dim);
303  auto cell_jac_inv = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac_inv",1,num_points,num_dim,num_dim);
304 
305  auto cell_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("cell_basis_ref_scalar",num_basis,num_points);
306  auto cell_grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("cell_grad_basis_ref",num_basis,num_points,num_dim);
307 
308  for(int cell=0;cell<num_cells;++cell){
309 
310  // =============================================
311  // Load external into cell-local arrays
312 
313  for(int p=0;p<num_points;++p)
314  for(int d=0;d<num_dim;++d)
315  for(int d2=0;d2<num_dim;++d2)
316  cell_jac_inv(0,p,d,d2)=jac_inv(cell,p,d,d2);
317  for(int p=0;p<num_points;++p)
318  for(int d=0;d<num_dim;++d)
319  cell_cub_points(p,d)=cub_points(cell,p,d);
320 
321  // =============================================
322  // Load Reference Values
323 
324  intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE);
325 
326  if(compute_derivatives){
327  Kokkos::deep_copy(cell_grad_basis_ref.get_view(),0.0);
328  }
329 
330  // =============================================
331  // Transform reference values to physical values
332 
333  fst::HGRADtransformVALUE(cell_basis_scalar.get_view(),cell_basis_ref_scalar.get_view());
334  for(int b=0;b<num_basis;++b)
335  for(int p=0;p<num_points;++p)
336  basis_scalar(cell,b,p)=cell_basis_scalar(0,b,p);
337 
338  if(compute_derivatives){
339  fst::HGRADtransformGRAD(cell_grad_basis.get_view(),cell_jac_inv.get_view(),cell_grad_basis_ref.get_view());
340  for(int b=0;b<num_basis;++b)
341  for(int p=0;p<num_points;++p)
342  for(int d=0;d<num_dim;++d)
343  grad_basis(cell,b,p,d)=cell_grad_basis(0,b,p,d);
344  }
345  // =============================================
346  }
347 
348 
349  if(build_weighted){
350  fst::multiplyMeasure(weighted_basis_scalar.get_view(),weighted_measure.get_view(),basis_scalar.get_view());
351  if(compute_derivatives){
352  fst::multiplyMeasure(weighted_grad_basis.get_view(),weighted_measure.get_view(),grad_basis.get_view());
353  }
354  }
355 
356 
357 }
358 
359 template <typename Scalar>
361 evaluateValues_HGrad(const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cub_points,
362  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv,
363  const PHX::MDField<Scalar,Cell,IP> & weighted_measure)
364 {
365 
366  TEUCHOS_ASSERT(getElementSpace() == PureBasis::HGRAD);
367 
368  typedef Intrepid2::FunctionSpaceTools<PHX::Device::execution_space> fst;
369  MDFieldArrayFactory af("",ddims_,true);
370 
371  const panzer::PureBasis & basis = *(basis_layout->getBasis());
372 
373  const int num_points = basis_layout->numPoints();
374  const int num_basis = basis.cardinality();
375  const int num_dim = basis_layout->dimension();
376  const int num_cells = cub_points.extent(0);
377 
378  auto cell_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("cell_basis_scalar",1,num_basis,num_points);
379  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
380  auto cell_grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_grad_basis",1,num_basis,num_points,num_dim);
381  auto cell_jac_inv = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac_inv",1,num_points,num_dim,num_dim);
382 
383  auto cell_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("cell_basis_ref_scalar",num_basis,num_points);
384  auto cell_grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("cell_grad_basis_ref",num_basis,num_points,num_dim);
385 
386  for(int cell=0;cell<num_cells;++cell){
387 
388  // =============================================
389  // Load external into cell-local arrays
390 
391  for(int p=0;p<num_points;++p)
392  for(int d=0;d<num_dim;++d)
393  for(int d2=0;d2<num_dim;++d2)
394  cell_jac_inv(0,p,d,d2)=jac_inv(cell,p,d,d2);
395  for(int p=0;p<num_points;++p)
396  for(int d=0;d<num_dim;++d)
397  cell_cub_points(p,d)=cub_points(cell,p,d);
398 
399  // =============================================
400  // Load Reference Values
401 
402  intrepid_basis->getValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE);
403 
404  if(compute_derivatives){
405  intrepid_basis->getValues(cell_grad_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_GRAD);
406  }
407 
408  // =============================================
409  // Transform reference values to physical values
410 
411  fst::HGRADtransformVALUE(cell_basis_scalar.get_view(),cell_basis_ref_scalar.get_view());
412  for(int b=0;b<num_basis;++b)
413  for(int p=0;p<num_points;++p)
414  basis_scalar(cell,b,p)=cell_basis_scalar(0,b,p);
415 
416  if(compute_derivatives){
417  fst::HGRADtransformGRAD(cell_grad_basis.get_view(),cell_jac_inv.get_view(),cell_grad_basis_ref.get_view());
418  for(int b=0;b<num_basis;++b)
419  for(int p=0;p<num_points;++p)
420  for(int d=0;d<num_dim;++d)
421  grad_basis(cell,b,p,d)=cell_grad_basis(0,b,p,d);
422  }
423  // =============================================
424  }
425 
426  if(build_weighted){
427  fst::multiplyMeasure(weighted_basis_scalar.get_view(),weighted_measure.get_view(),basis_scalar.get_view());
428  if(compute_derivatives){
429  fst::multiplyMeasure(weighted_grad_basis.get_view(),weighted_measure.get_view(),grad_basis.get_view());
430  }
431  }
432 
433 }
434 
435 
436 template <typename Scalar>
438 evaluateValues_HCurl(const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cub_points,
439  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
440  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
441  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv,
442  const PHX::MDField<Scalar,Cell,IP> & weighted_measure)
443 {
444 
445  TEUCHOS_ASSERT(getElementSpace() == PureBasis::HCURL);
446 
447 
448  typedef Intrepid2::FunctionSpaceTools<PHX::Device::execution_space> fst;
449  MDFieldArrayFactory af("",ddims_,true);
450 
451  const panzer::PureBasis & basis = *(basis_layout->getBasis());
452 
453  const int num_points = basis_layout->numPoints();
454  const int num_basis = basis.cardinality();
455  const int num_dim = basis_layout->dimension();
456  const int num_cells = basis_layout->numCells();
457 
458  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
459  auto cell_jac = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac",1,num_points,num_dim,num_dim);
460  auto cell_jac_inv = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac_inv",1,num_points,num_dim,num_dim);
461  auto cell_jac_det = af.buildStaticArray<Scalar,Cell,IP>("cell_jac_det",1,num_points);
462 
463  auto cell_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_basis_vector",1,num_basis,num_points,num_dim);
464  auto cell_curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("cell_curl_basis_scalar",1,num_basis,num_points);
465  auto cell_curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_curl_basis_vector",1,num_basis,num_points,num_dim);
466 
467  auto cell_curl_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("cell_curl_basis_ref",num_basis,num_points,num_dim);
468  auto cell_curl_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("cell_curl_basis_ref_scalar",num_basis,num_points);
469  auto cell_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("cell_basis_ref_vector",num_basis,num_points,num_dim);
470 
471  for(int cell=0;cell<num_cells;++cell){
472 
473  // =============================================
474  // Load external into cell-local arrays
475 
476  for(int p=0;p<num_points;++p)
477  for(int d=0;d<num_dim;++d)
478  for(int d2=0;d2<num_dim;++d2)
479  cell_jac(0,p,d,d2)=jac(cell,p,d,d2);
480  for(int p=0;p<num_points;++p)
481  for(int d=0;d<num_dim;++d)
482  for(int d2=0;d2<num_dim;++d2)
483  cell_jac_inv(0,p,d,d2)=jac_inv(cell,p,d,d2);
484  for(int p=0;p<num_points;++p)
485  cell_jac_det(0,p)=jac_det(cell,p);
486  for(int p=0;p<num_points;++p)
487  for(int d=0;d<num_dim;++d)
488  cell_cub_points(p,d)=cub_points(cell,p,d);
489 
490  // =============================================
491  // Load Reference Values
492 
493  intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE);
494 
495  if(compute_derivatives){
496  if(num_dim==2){
497  intrepid_basis->getValues(cell_curl_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL);
498  } else if(num_dim==3){
499  intrepid_basis->getValues(cell_curl_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_CURL);
500  }
501  }
502 
503  // =============================================
504  // Transform reference values to physical values
505 
506  fst::HCURLtransformVALUE(cell_basis_vector.get_view(),cell_jac_inv.get_view(),cell_basis_ref_vector.get_view());
507  for(int b=0;b<num_basis;++b)
508  for(int p=0;p<num_points;++p)
509  for(int d=0;d<num_dim;++d)
510  basis_vector(cell,b,p,d)=cell_basis_vector(0,b,p,d);
511 
512  if(compute_derivatives){
513  if(num_dim==2){
514  // note only volume deformation is needed!
515  // this relates directly to this being in
516  // the divergence space in 2D!
517  fst::HDIVtransformDIV(cell_curl_basis_scalar.get_view(),cell_jac_det.get_view(),cell_curl_basis_ref_scalar.get_view());
518  for(int b=0;b<num_basis;++b)
519  for(int p=0;p<num_points;++p)
520  curl_basis_scalar(cell,b,p)=cell_curl_basis_scalar(0,b,p);
521  } else if(num_dim==3) {
522  fst::HCURLtransformCURL(cell_curl_basis_vector.get_view(),cell_jac.get_view(),cell_jac_det.get_view(),cell_curl_basis_ref.get_view());
523  for(int b=0;b<num_basis;++b)
524  for(int p=0;p<num_points;++p)
525  for(int d=0;d<num_dim;++d)
526  curl_basis_vector(cell,b,p,d)=cell_curl_basis_vector(0,b,p,d);
527  } else {
528  TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::evaluateValues_HCurl : HCurl only setup for 2D and 3D.");
529  }
530  }
531  }
532 
533  if(build_weighted){
534  fst::multiplyMeasure(weighted_basis_vector.get_view(),weighted_measure.get_view(),basis_vector.get_view());
535  if(compute_derivatives){
536  if(num_dim==2){
537  fst::multiplyMeasure(weighted_curl_basis_scalar.get_view(),weighted_measure.get_view(),curl_basis_scalar.get_view());
538  } else if(num_dim==3){
539  fst::multiplyMeasure(weighted_curl_basis_vector.get_view(),weighted_measure.get_view(),curl_basis_vector.get_view());
540  }
541  }
542  }
543 
544 }
545 
546 template <typename Scalar>
548 evaluateValues_HDiv(const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cub_points,
549  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
550  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
551  const PHX::MDField<Scalar,Cell,IP> & weighted_measure)
552 {
553 
554  TEUCHOS_ASSERT(getElementSpace() == PureBasis::HDIV);
555 
556  typedef Intrepid2::FunctionSpaceTools<PHX::Device::execution_space> fst;
557  MDFieldArrayFactory af("",ddims_,true);
558 
559  const panzer::PureBasis & basis = *(basis_layout->getBasis());
560 
561  const int num_points = basis_layout->numPoints();
562  const int num_basis = basis.cardinality();
563  const int num_dim = basis_layout->dimension();
564  const int num_cells = basis_layout->numCells();
565 
566  auto cell_cub_points = af.buildStaticArray<Scalar,IP,Dim>("cell_cub_points",num_points,num_dim);
567  auto cell_jac = af.buildStaticArray<Scalar,Cell,IP,Dim,Dim>("cell_jac",1,num_points,num_dim,num_dim);
568  auto cell_jac_det = af.buildStaticArray<Scalar,Cell,IP>("cell_jac_det",1,num_points);
569 
570  auto cell_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("cell_basis_vector",1,num_basis,num_points,num_dim);
571  auto cell_div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("cell_div_basis",1,num_basis,num_points);
572 
573  auto cell_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("cell_basis_ref_vector",num_basis,num_points,num_dim);
574  auto cell_div_basis_ref = af.buildStaticArray<Scalar,BASIS,IP>("cell_div_basis_ref",num_basis,num_points);
575 
576  for(int cell=0;cell<num_cells;++cell){
577 
578  // =============================================
579  // Load external into cell-local arrays
580 
581  for(int p=0;p<num_points;++p)
582  for(int d=0;d<num_dim;++d)
583  for(int d2=0;d2<num_dim;++d2)
584  cell_jac(0,p,d,d2)=jac(cell,p,d,d2);
585  for(int p=0;p<num_points;++p)
586  cell_jac_det(0,p)=jac_det(cell,p);
587  for(int p=0;p<num_points;++p)
588  for(int d=0;d<num_dim;++d)
589  cell_cub_points(p,d)=cub_points(cell,p,d);
590  // =============================================
591  // Load Reference Values
592 
593  intrepid_basis->getValues(cell_basis_ref_vector.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE);
594 
595  if(compute_derivatives){
596  intrepid_basis->getValues(cell_div_basis_ref.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_DIV);
597  }
598 
599  // =============================================
600  // Transform reference values to physical values
601 
602  fst::HDIVtransformVALUE(cell_basis_vector.get_view(),cell_jac.get_view(),cell_jac_det.get_view(),cell_basis_ref_vector.get_view());
603  for(int b=0;b<num_basis;++b)
604  for(int p=0;p<num_points;++p)
605  for(int d=0;d<num_dim;++d)
606  basis_vector(cell,b,p,d)=cell_basis_vector(0,b,p,d);
607 
608  if(compute_derivatives){
609  fst::HDIVtransformDIV(cell_div_basis.get_view(),cell_jac_det.get_view(),cell_div_basis_ref.get_view());
610  for(int b=0;b<num_basis;++b)
611  for(int p=0;p<num_points;++p)
612  div_basis(cell,b,p)=cell_div_basis(0,b,p);
613  }
614  }
615 
616  if(build_weighted){
617  fst::multiplyMeasure(weighted_basis_vector.get_view(),weighted_measure.get_view(),basis_vector.get_view());
618  if(compute_derivatives){
619  fst::multiplyMeasure(weighted_div_basis.get_view(),weighted_measure.get_view(),div_basis.get_view());
620  }
621  }
622 
623 }
624 
625 
626 
627 
628 
629 
630 
631 
632 
633 
634 
635 
636 
637 
638 
639 
640 
641 
642 
643 
644 
645 
646 
647 
648 
649 
650 
651 
652 
653 
654 
655 
656 
657 
658 
659 
660 
661 
662 
663 
664 
665 
666 
667 template <typename Scalar>
669 evaluateValuesCV(const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cell_cub_points,
670  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
671  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
672  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv)
673 {
674  MDFieldArrayFactory af("",ddims_,true);
675 
676  int num_ip = basis_layout->numPoints();
677  int num_card = basis_layout->cardinality();
678  int num_dim = basis_layout->dimension();
679 
680  size_type num_cells = jac.extent(0);
681 
682  PureBasis::EElementSpace elmtspace = getElementSpace();
683  ArrayDynamic dyn_cub_points = af.buildArray<Scalar,IP,Dim>("dyn_cub_points", num_ip, num_dim);
684 
685  // Integration points are located on physical cells rather than reference cells,
686  // so we evaluate the basis in a loop over cells.
687  for (size_type icell = 0; icell < num_cells; ++icell)
688  {
689  for (int ip = 0; ip < num_ip; ++ip)
690  for (int d = 0; d < num_dim; ++d)
691  dyn_cub_points(ip,d) = cell_cub_points(icell,ip,d);
692 
693  if(elmtspace==PureBasis::CONST) {
694  ArrayDynamic dyn_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_ip);
695 
696  intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(),
697  dyn_cub_points.get_view(),
698  Intrepid2::OPERATOR_VALUE);
699 
700  // transform values method just transfers values to array with cell index - no need to call
701  for (int b = 0; b < num_card; ++b)
702  for (int ip = 0; ip < num_ip; ++ip)
703  basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip);
704 
705  }
706  if(elmtspace==PureBasis::HGRAD) {
707  ArrayDynamic dyn_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_ip);
708 
709  intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(),
710  dyn_cub_points.get_view(),
711  Intrepid2::OPERATOR_VALUE);
712 
713  // transform values method just transfers values to array with cell index - no need to call
714  for (int b = 0; b < num_card; ++b)
715  for (int ip = 0; ip < num_ip; ++ip)
716  basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip);
717 
718  if(compute_derivatives) {
719 
720  int one_cell = 1;
721  ArrayDynamic dyn_grad_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_grad_basis_ref",num_card,num_ip,num_dim);
722  ArrayDynamic dyn_grad_basis = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_grad_basis",one_cell,num_card,num_ip,num_dim);
723  ArrayDynamic dyn_jac_inv = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim);
724 
725  intrepid_basis->getValues(dyn_grad_basis_ref.get_view(),
726  dyn_cub_points.get_view(),
727  Intrepid2::OPERATOR_GRAD);
728 
729  int cellInd = 0;
730  for (int ip = 0; ip < num_ip; ++ip)
731  for (int d1 = 0; d1 < num_dim; ++d1)
732  for (int d2 = 0; d2 < num_dim; ++d2)
733  dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2);
734 
735  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::HGRADtransformGRAD<Scalar>(dyn_grad_basis.get_view(),
736  dyn_jac_inv.get_view(),
737  dyn_grad_basis_ref.get_view());
738 
739  for (int b = 0; b < num_card; ++b)
740  for (int ip = 0; ip < num_ip; ++ip)
741  for (int d = 0; d < num_dim; ++d)
742  grad_basis(icell,b,ip,d) = dyn_grad_basis(0,b,ip,d);
743 
744  }
745  }
746  else if(elmtspace==PureBasis::HCURL) {
747  ArrayDynamic dyn_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_ip,num_dim);
748 
749  intrepid_basis->getValues(dyn_basis_ref_vector.get_view(),
750  dyn_cub_points.get_view(),
751  Intrepid2::OPERATOR_VALUE);
752 
753  int one_cell = 1;
754  ArrayDynamic dyn_basis_vector = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_basis_vector",one_cell,num_card,num_ip,num_dim);
755  ArrayDynamic dyn_jac_inv = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim);
756 
757  int cellInd = 0;
758  for (int ip = 0; ip < num_ip; ++ip)
759  for (int d1 = 0; d1 < num_dim; ++d1)
760  for (int d2 = 0; d2 < num_dim; ++d2)
761  dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2);
762 
763  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::HCURLtransformVALUE(dyn_basis_vector.get_view(),
764  dyn_jac_inv.get_view(),
765  dyn_basis_ref_vector.get_view());
766 
767  for (int b = 0; b < num_card; ++b)
768  for (int ip = 0; ip < num_ip; ++ip)
769  for (int d = 0; d < num_dim; ++d)
770  basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d);
771 
772  if(compute_derivatives && num_dim ==2) {
773 
774  int one_cell = 1;
775  ArrayDynamic dyn_curl_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_curl_basis_ref_scalar",num_card,num_ip);
776  ArrayDynamic dyn_curl_basis_scalar = af.buildArray<Scalar,Cell,BASIS,IP>("dyn_curl_basis_scalar",one_cell,num_card,num_ip);
777  ArrayDynamic dyn_jac_det = af.buildArray<Scalar,Cell,IP>("dyn_jac_det",one_cell,num_ip);
778 
779  intrepid_basis->getValues(dyn_curl_basis_ref_scalar.get_view(),
780  dyn_cub_points.get_view(),
781  Intrepid2::OPERATOR_CURL);
782 
783  int cellInd = 0;
784  for (int ip = 0; ip < num_ip; ++ip)
785  dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
786 
787  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::HDIVtransformDIV(dyn_curl_basis_scalar.get_view(),
788  dyn_jac_det.get_view(),
789  dyn_curl_basis_ref_scalar.get_view());
790 
791  for (int b = 0; b < num_card; ++b)
792  for (int ip = 0; ip < num_ip; ++ip)
793  curl_basis_scalar(icell,b,ip) = dyn_curl_basis_scalar(0,b,ip);
794 
795  }
796  if(compute_derivatives && num_dim ==3) {
797 
798  int one_cell = 1;
799  ArrayDynamic dyn_curl_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_curl_basis_ref_vector",num_card,num_ip,num_dim);
800  ArrayDynamic dyn_curl_basis = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_curl_basis_vector",one_cell,num_card,num_ip,num_dim);
801  ArrayDynamic dyn_jac_det = af.buildArray<Scalar,Cell,IP>("dyn_jac_det",one_cell,num_ip);
802  ArrayDynamic dyn_jac = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac",one_cell,num_ip,num_dim,num_dim);
803 
804  intrepid_basis->getValues(dyn_curl_basis_ref.get_view(),
805  dyn_cub_points.get_view(),
806  Intrepid2::OPERATOR_CURL);
807 
808  int cellInd = 0;
809  for (int ip = 0; ip < num_ip; ++ip)
810  {
811  dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
812  for (int d1 = 0; d1 < num_dim; ++d1)
813  for (int d2 = 0; d2 < num_dim; ++d2)
814  dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2);
815  }
816 
817  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::HCURLtransformCURL(dyn_curl_basis.get_view(),
818  dyn_jac.get_view(),
819  dyn_jac_det.get_view(),
820  dyn_curl_basis_ref.get_view());
821 
822  for (int b = 0; b < num_card; ++b)
823  for (int ip = 0; ip < num_ip; ++ip)
824  for (int d = 0; d < num_dim; ++d)
825  curl_basis_vector(icell,b,ip,d) = dyn_curl_basis(0,b,ip,d);
826 
827  }
828 
829  }
830  else if(elmtspace==PureBasis::HDIV) {
831 
832  ArrayDynamic dyn_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_ip,num_dim);
833 
834  intrepid_basis->getValues(dyn_basis_ref_vector.get_view(),
835  dyn_cub_points.get_view(),
836  Intrepid2::OPERATOR_VALUE);
837 
838  int one_cell= 1;
839  ArrayDynamic dyn_basis_vector = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_basis_vector",one_cell,num_card,num_ip,num_dim);
840  ArrayDynamic dyn_jac = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac",one_cell,num_ip,num_dim,num_dim);
841  ArrayDynamic dyn_jac_det = af.buildArray<Scalar,Cell,IP>("dyn_jac_det",one_cell,num_ip);
842 
843  int cellInd = 0;
844  for (int ip = 0; ip < num_ip; ++ip)
845  {
846  dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
847  for (int d1 = 0; d1 < num_dim; ++d1)
848  for (int d2 = 0; d2 < num_dim; ++d2)
849  dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2);
850  }
851 
852  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::HDIVtransformVALUE(dyn_basis_vector.get_view(),
853  dyn_jac.get_view(),
854  dyn_jac_det.get_view(),
855  dyn_basis_ref_vector.get_view());
856 
857  for (int b = 0; b < num_card; ++b)
858  for (int ip = 0; ip < num_ip; ++ip)
859  for (int d = 0; d < num_dim; ++d)
860  basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d);
861 
862  if(compute_derivatives) {
863 
864  ArrayDynamic dyn_div_basis_ref = af.buildArray<Scalar,BASIS,IP>("dyn_div_basis_ref_scalar",num_card,num_ip);
865  ArrayDynamic dyn_div_basis = af.buildArray<Scalar,Cell,BASIS,IP>("dyn_div_basis_scalar",one_cell,num_card,num_ip);
866 
867  intrepid_basis->getValues(dyn_div_basis_ref.get_view(),
868  dyn_cub_points.get_view(),
869  Intrepid2::OPERATOR_DIV);
870 
871  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::HDIVtransformDIV<Scalar>(dyn_div_basis.get_view(),
872  dyn_jac_det.get_view(),
873  dyn_div_basis_ref.get_view());
874 
875  for (int b = 0; b < num_card; ++b)
876  for (int ip = 0; ip < num_ip; ++ip)
877  div_basis(icell,b,ip) = dyn_div_basis(0,b,ip);
878 
879  }
880 
881  }
882  else { TEUCHOS_ASSERT(false); }
883 
884  } // cell loop
885 
886 }
887 
888 template <typename Scalar>
890 evaluateReferenceValues(const PHX::MDField<Scalar,IP,Dim> & cub_points,bool compute_derivatives,bool use_vertex_coordinates)
891 {
892  MDFieldArrayFactory af("",ddims_,true);
893 
894  int num_quad = basis_layout->numPoints();
895  int num_dim = basis_layout->dimension();
896  int num_card = basis_layout->cardinality();
897 
898  ArrayDynamic dyn_cub_points = af.buildArray<Scalar,IP,Dim>("dyn_cub_points", num_quad,num_dim);
899 
900  for (int ip = 0; ip < num_quad; ++ip)
901  for (int d = 0; d < num_dim; ++d)
902  dyn_cub_points(ip,d) = cub_points(ip,d);
903 
904  PureBasis::EElementSpace elmtspace = getElementSpace();
905  if(elmtspace==PureBasis::HGRAD || elmtspace==PureBasis::CONST) {
906  ArrayDynamic dyn_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_quad);
907 
908  intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(),
909  dyn_cub_points.get_view(),
910  Intrepid2::OPERATOR_VALUE);
911 
912  for (int b = 0; b < num_card; ++b)
913  for (int ip = 0; ip < num_quad; ++ip)
914  basis_ref_scalar(b,ip) = dyn_basis_ref_scalar(b,ip);
915  }
916  else if(elmtspace==PureBasis::HDIV || elmtspace==PureBasis::HCURL) {
917  ArrayDynamic dyn_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
918 
919  intrepid_basis->getValues(dyn_basis_ref_vector.get_view(),
920  dyn_cub_points.get_view(),
921  Intrepid2::OPERATOR_VALUE);
922 
923  for (int b = 0; b < num_card; ++b)
924  for (int ip = 0; ip < num_quad; ++ip)
925  for (int d = 0; d < num_dim; ++d)
926  basis_ref_vector(b,ip,d) = dyn_basis_ref_vector(b,ip,d);
927  }
928  else { TEUCHOS_ASSERT(false); }
929 
930  if(elmtspace==PureBasis::HGRAD && compute_derivatives) {
931  ArrayDynamic dyn_grad_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
932 
933  intrepid_basis->getValues(dyn_grad_basis_ref.get_view(),
934  dyn_cub_points.get_view(),
935  Intrepid2::OPERATOR_GRAD);
936 
937  for (int b = 0; b < num_card; ++b)
938  for (int ip = 0; ip < num_quad; ++ip)
939  for (int d = 0; d < num_dim; ++d)
940  grad_basis_ref(b,ip,d) = dyn_grad_basis_ref(b,ip,d);
941  }
942  else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==2) {
943  ArrayDynamic dyn_curl_basis_ref = af.buildArray<Scalar,BASIS,IP>("dyn_curl_basis_ref_scalar",num_card,num_quad);
944 
945  intrepid_basis->getValues(dyn_curl_basis_ref.get_view(),
946  dyn_cub_points.get_view(),
947  Intrepid2::OPERATOR_CURL);
948 
949  for (int b = 0; b < num_card; ++b)
950  for (int ip = 0; ip < num_quad; ++ip)
951  curl_basis_ref_scalar(b,ip) = dyn_curl_basis_ref(b,ip);
952  }
953  else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==3) {
954  ArrayDynamic dyn_curl_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_curl_basis_ref_vector",num_card,num_quad,num_dim);
955 
956  intrepid_basis->getValues(dyn_curl_basis_ref.get_view(),
957  dyn_cub_points.get_view(),
958  Intrepid2::OPERATOR_CURL);
959 
960  for (int b = 0; b < num_card; ++b)
961  for (int ip = 0; ip < num_quad; ++ip)
962  for (int d = 0; d < num_dim; ++d)
963  curl_basis_ref_vector(b,ip,d) = dyn_curl_basis_ref(b,ip,d);
964  }
965  else if(elmtspace==PureBasis::HDIV && compute_derivatives) {
966  ArrayDynamic dyn_div_basis_ref = af.buildArray<Scalar,BASIS,IP>("dyn_div_basis_ref_scalar",num_card,num_quad);
967 
968  intrepid_basis->getValues(dyn_div_basis_ref.get_view(),
969  dyn_cub_points.get_view(),
970  Intrepid2::OPERATOR_DIV);
971 
972  for (int b = 0; b < num_card; ++b)
973  for (int ip = 0; ip < num_quad; ++ip)
974  div_basis_ref(b,ip) = dyn_div_basis_ref(b,ip);
975  }
976 
977 
978  if(use_vertex_coordinates) {
979  // Intrepid removes fad types from the coordinate scalar type. We
980  // pull the actual field scalar type from the basis object to be
981  // consistent.
982  if (elmtspace != PureBasis::CONST) {
983  using coordsScalarType = typename Intrepid2::Basis<PHX::Device::execution_space,Scalar,Scalar>::scalarType;
984  auto dyn_basis_coordinates_ref = af.buildArray<coordsScalarType,BASIS,Dim>("basis_coordinates_ref",
985  basis_coordinates_ref.extent(0),
986  basis_coordinates_ref.extent(1));
987  intrepid_basis->getDofCoords(dyn_basis_coordinates_ref.get_view());
988 
989  // fill in basis coordinates
990  for (int i = 0; i < basis_coordinates_ref.extent_int(0); ++i)
991  for (int j = 0; j < basis_coordinates_ref.extent_int(1); ++j)
992  basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j);
993  }
994  }
995 
996  references_evaluated = true;
997 }
998 
999 // method for applying orientations
1000 template <typename Scalar>
1002 applyOrientations(const std::vector<Intrepid2::Orientation> & orientations)
1003 {
1004  if (!intrepid_basis->requireOrientation())
1005  return;
1006 
1007  typedef Intrepid2::OrientationTools<PHX::Device> ots;
1008  const PureBasis::EElementSpace elmtspace = getElementSpace();
1009 
1010  // orientation (right now std vector) is created using push_back method.
1011  // thus, its size is the actual number of elements to be applied.
1012  // on the other hand, basis_layout num cell indicates workset size.
1013  // to get right size of cells, use minimum of them.
1014  const int num_cell_basis_layout = basis_layout->numCells(), num_cell_orientation = orientations.size();
1015  const int num_cell = num_cell_basis_layout < num_cell_orientation ? num_cell_basis_layout : num_cell_orientation;
1016  const int num_dim = basis_layout->dimension();
1017  const Kokkos::pair<int,int> range_cell(0, num_cell);
1018 
1019  Kokkos::DynRankView<Intrepid2::Orientation,PHX::Device>
1020  drv_orts((Intrepid2::Orientation*)orientations.data(), num_cell);
1021 
1025  if (elmtspace==PureBasis::HGRAD) {
1026  {
1027  {
1028  auto drv_basis_scalar = Kokkos::subview(basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL());
1029  auto drv_basis_scalar_tmp = Kokkos::createDynRankView(basis_scalar.get_view(),
1030  "drv_basis_scalar_tmp",
1031  drv_basis_scalar.extent(0), // C
1032  drv_basis_scalar.extent(1), // F
1033  drv_basis_scalar.extent(2)); // P
1034  Kokkos::deep_copy(drv_basis_scalar_tmp, drv_basis_scalar);
1035  ots::modifyBasisByOrientation(drv_basis_scalar,
1036  drv_basis_scalar_tmp,
1037  drv_orts,
1038  intrepid_basis);
1039  }
1040  if(build_weighted) {
1041  auto drv_basis_scalar = Kokkos::subview(weighted_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL());
1042  auto drv_basis_scalar_tmp = Kokkos::createDynRankView(weighted_basis_scalar.get_view(),
1043  "drv_basis_scalar_tmp",
1044  drv_basis_scalar.extent(0), // C
1045  drv_basis_scalar.extent(1), // F
1046  drv_basis_scalar.extent(2)); // P
1047  Kokkos::deep_copy(drv_basis_scalar_tmp, drv_basis_scalar);
1048  ots::modifyBasisByOrientation(drv_basis_scalar,
1049  drv_basis_scalar_tmp,
1050  drv_orts,
1051  intrepid_basis);
1052  }
1053 
1054  }
1055 
1056  if (compute_derivatives) {
1057  {
1058  auto drv_grad_basis = Kokkos::subview(grad_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1059  auto drv_grad_basis_tmp = Kokkos::createDynRankView(grad_basis.get_view(),
1060  "drv_grad_basis_tmp",
1061  drv_grad_basis.extent(0), // C
1062  drv_grad_basis.extent(1), // F
1063  drv_grad_basis.extent(2), // P
1064  drv_grad_basis.extent(3)); // D
1065  Kokkos::deep_copy(drv_grad_basis_tmp, drv_grad_basis);
1066  ots::modifyBasisByOrientation(drv_grad_basis,
1067  drv_grad_basis_tmp,
1068  drv_orts,
1069  intrepid_basis);
1070  }
1071  if(build_weighted) {
1072  auto drv_grad_basis = Kokkos::subview(weighted_grad_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1073  auto drv_grad_basis_tmp = Kokkos::createDynRankView(weighted_grad_basis.get_view(),
1074  "drv_grad_basis_tmp",
1075  drv_grad_basis.extent(0), // C
1076  drv_grad_basis.extent(1), // F
1077  drv_grad_basis.extent(2), // P
1078  drv_grad_basis.extent(3)); // D
1079  Kokkos::deep_copy(drv_grad_basis_tmp, drv_grad_basis);
1080  ots::modifyBasisByOrientation(drv_grad_basis,
1081  drv_grad_basis_tmp,
1082  drv_orts,
1083  intrepid_basis);
1084  }
1085  }
1086  }
1087 
1091  else if (elmtspace==PureBasis::HCURL && num_dim==2) {
1092  {
1093  {
1094  auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1095  auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(),
1096  "drv_basis_vector_tmp",
1097  drv_basis_vector.extent(0), // C
1098  drv_basis_vector.extent(1), // F
1099  drv_basis_vector.extent(2), // P
1100  drv_basis_vector.extent(3)); // D
1101  Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector);
1102  ots::modifyBasisByOrientation(drv_basis_vector,
1103  drv_basis_vector_tmp,
1104  drv_orts,
1105  intrepid_basis);
1106  }
1107  if(build_weighted) {
1108  auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1109  auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(),
1110  "drv_basis_vector_tmp",
1111  drv_basis_vector.extent(0), // C
1112  drv_basis_vector.extent(1), // F
1113  drv_basis_vector.extent(2), // P
1114  drv_basis_vector.extent(3)); // D
1115  Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector);
1116  ots::modifyBasisByOrientation(drv_basis_vector,
1117  drv_basis_vector_tmp,
1118  drv_orts,
1119  intrepid_basis);
1120  }
1121  }
1122 
1123  if (compute_derivatives) {
1124  {
1125  auto drv_curl_basis_scalar = Kokkos::subview(curl_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL());
1126  auto drv_curl_basis_scalar_tmp = Kokkos::createDynRankView(curl_basis_scalar.get_view(),
1127  "drv_curl_basis_scalar_tmp",
1128  drv_curl_basis_scalar.extent(0), // C
1129  drv_curl_basis_scalar.extent(1), // F
1130  drv_curl_basis_scalar.extent(2)); // P
1131  Kokkos::deep_copy(drv_curl_basis_scalar_tmp, drv_curl_basis_scalar);
1132  ots::modifyBasisByOrientation(drv_curl_basis_scalar,
1133  drv_curl_basis_scalar_tmp,
1134  drv_orts,
1135  intrepid_basis);
1136  }
1137 
1138  if(build_weighted) {
1139  auto drv_curl_basis_scalar = Kokkos::subview(weighted_curl_basis_scalar.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL());
1140  auto drv_curl_basis_scalar_tmp = Kokkos::createDynRankView(weighted_curl_basis_scalar.get_view(),
1141  "drv_curl_basis_scalar_tmp",
1142  drv_curl_basis_scalar.extent(0), // C
1143  drv_curl_basis_scalar.extent(1), // F
1144  drv_curl_basis_scalar.extent(2)); // P
1145  Kokkos::deep_copy(drv_curl_basis_scalar_tmp, drv_curl_basis_scalar);
1146  ots::modifyBasisByOrientation(drv_curl_basis_scalar,
1147  drv_curl_basis_scalar_tmp,
1148  drv_orts,
1149  intrepid_basis);
1150  }
1151  }
1152  }
1153 
1157  else if (elmtspace==PureBasis::HCURL && num_dim==3) {
1158  {
1159  {
1160  auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1161  auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(),
1162  "drv_basis_vector_tmp",
1163  drv_basis_vector.extent(0), // C
1164  drv_basis_vector.extent(1), // F
1165  drv_basis_vector.extent(2), // P
1166  drv_basis_vector.extent(3)); // D
1167  Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector);
1168  ots::modifyBasisByOrientation(drv_basis_vector,
1169  drv_basis_vector_tmp,
1170  drv_orts,
1171  intrepid_basis);
1172  }
1173  if(build_weighted) {
1174  auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1175  auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(),
1176  "drv_basis_vector_tmp",
1177  drv_basis_vector.extent(0), // C
1178  drv_basis_vector.extent(1), // F
1179  drv_basis_vector.extent(2), // P
1180  drv_basis_vector.extent(3)); // D
1181  Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector);
1182  ots::modifyBasisByOrientation(drv_basis_vector,
1183  drv_basis_vector_tmp,
1184  drv_orts,
1185  intrepid_basis);
1186  }
1187  }
1188 
1189  if (compute_derivatives) {
1190  {
1191  auto drv_curl_basis_vector = Kokkos::subview(curl_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1192  auto drv_curl_basis_vector_tmp = Kokkos::createDynRankView(curl_basis_vector.get_view(),
1193  "drv_curl_basis_vector_tmp",
1194  drv_curl_basis_vector.extent(0), // C
1195  drv_curl_basis_vector.extent(1), // F
1196  drv_curl_basis_vector.extent(2), // P
1197  drv_curl_basis_vector.extent(3)); // D
1198  Kokkos::deep_copy(drv_curl_basis_vector_tmp, drv_curl_basis_vector);
1199  ots::modifyBasisByOrientation(drv_curl_basis_vector,
1200  drv_curl_basis_vector_tmp,
1201  drv_orts,
1202  intrepid_basis);
1203  }
1204  if(build_weighted) {
1205  auto drv_curl_basis_vector = Kokkos::subview(weighted_curl_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1206  auto drv_curl_basis_vector_tmp = Kokkos::createDynRankView(weighted_curl_basis_vector.get_view(),
1207  "drv_curl_basis_vector_tmp",
1208  drv_curl_basis_vector.extent(0), // C
1209  drv_curl_basis_vector.extent(1), // F
1210  drv_curl_basis_vector.extent(2), // P
1211  drv_curl_basis_vector.extent(3)); // D
1212  Kokkos::deep_copy(drv_curl_basis_vector_tmp, drv_curl_basis_vector);
1213  ots::modifyBasisByOrientation(drv_curl_basis_vector,
1214  drv_curl_basis_vector_tmp,
1215  drv_orts,
1216  intrepid_basis);
1217  }
1218  }
1219  }
1223  else if (elmtspace==PureBasis::HDIV) {
1224  {
1225  {
1226  auto drv_basis_vector = Kokkos::subview(basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1227  auto drv_basis_vector_tmp = Kokkos::createDynRankView(basis_vector.get_view(),
1228  "drv_basis_vector_tmp",
1229  drv_basis_vector.extent(0), // C
1230  drv_basis_vector.extent(1), // F
1231  drv_basis_vector.extent(2), // P
1232  drv_basis_vector.extent(3)); // D
1233  Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector);
1234  ots::modifyBasisByOrientation(drv_basis_vector,
1235  drv_basis_vector_tmp,
1236  drv_orts,
1237  intrepid_basis);
1238  }
1239  if(build_weighted) {
1240  auto drv_basis_vector = Kokkos::subview(weighted_basis_vector.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
1241  auto drv_basis_vector_tmp = Kokkos::createDynRankView(weighted_basis_vector.get_view(),
1242  "drv_basis_vector_tmp",
1243  drv_basis_vector.extent(0), // C
1244  drv_basis_vector.extent(1), // F
1245  drv_basis_vector.extent(2), // P
1246  drv_basis_vector.extent(3)); // D
1247  Kokkos::deep_copy(drv_basis_vector_tmp, drv_basis_vector);
1248  ots::modifyBasisByOrientation(drv_basis_vector,
1249  drv_basis_vector_tmp,
1250  drv_orts,
1251  intrepid_basis);
1252  }
1253  }
1254  if (compute_derivatives) {
1255  {
1256  auto drv_div_basis = Kokkos::subview(div_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL());
1257  auto drv_div_basis_tmp = Kokkos::createDynRankView(div_basis.get_view(),
1258  "drv_div_basis_tmp",
1259  drv_div_basis.extent(0), // C
1260  drv_div_basis.extent(1), // F
1261  drv_div_basis.extent(2)); // P
1262  Kokkos::deep_copy(drv_div_basis_tmp, drv_div_basis);
1263  ots::modifyBasisByOrientation(drv_div_basis,
1264  drv_div_basis_tmp,
1265  drv_orts,
1266  intrepid_basis);
1267  }
1268  if(build_weighted) {
1269  auto drv_div_basis = Kokkos::subview(weighted_div_basis.get_view(), range_cell, Kokkos::ALL(), Kokkos::ALL());
1270  auto drv_div_basis_tmp = Kokkos::createDynRankView(weighted_div_basis.get_view(),
1271  "drv_div_basis_tmp",
1272  drv_div_basis.extent(0), // C
1273  drv_div_basis.extent(1), // F
1274  drv_div_basis.extent(2)); // P
1275  Kokkos::deep_copy(drv_div_basis_tmp, drv_div_basis);
1276  ots::modifyBasisByOrientation(drv_div_basis,
1277  drv_div_basis_tmp,
1278  drv_orts,
1279  intrepid_basis);
1280  }
1281  }
1282  }
1283 }
1284 
1285 // method for applying orientations
1286 template <typename Scalar>
1288 applyOrientations(const PHX::MDField<const Scalar,Cell,BASIS> & orientations)
1289 {
1290 
1291  TEUCHOS_TEST_FOR_EXCEPT_MSG(true,"panzer::BasisValues2::applyOrientations : this should not be called.");
1292 
1293  int num_cell = orientations.extent(0);
1294  int num_basis = orientations.extent(1);
1295  int num_dim = basis_layout->dimension();
1296  int num_ip = basis_layout->numPoints();
1297  PureBasis::EElementSpace elmtspace = getElementSpace();
1298 
1299  if(elmtspace==PureBasis::HCURL && num_dim==2) {
1300 
1301  // setup the orientations for the trial space
1302  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(basis_vector,orientations);
1303 
1304  for (int c=0; c<num_cell; c++)
1305  for (int b=0; b<num_basis; b++)
1306  for (int p=0; p<num_ip; p++)
1307  for (int d=0; d<num_dim; d++)
1308  basis_vector(c, b, p, d) *= orientations(c, b);
1309 
1310  if(compute_derivatives) {
1311  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(curl_basis_scalar,orientations);
1312  for (int c=0; c<num_cell; c++)
1313  for (int b=0; b<num_basis; b++)
1314  for (int p=0; p<num_ip; p++)
1315  curl_basis_scalar(c, b, p) *= orientations(c, b);
1316  }
1317 
1318  // setup the orientations for the test space
1319  if(build_weighted) {
1320  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view());
1321 
1322  if(compute_derivatives)
1323  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::applyFieldSigns(weighted_curl_basis_scalar.get_view(),orientations.get_view());
1324  }
1325  }
1326  else if(elmtspace==PureBasis::HCURL && num_dim==3) {
1327 
1328  // setup the orientations for the trial space
1329  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(basis_vector,orientations);
1330 
1331  for (int c=0; c<num_cell; c++)
1332  for (int b=0; b<num_basis; b++)
1333  for (int p=0; p<num_ip; p++)
1334  for (int d=0; d<num_dim; d++)
1335  basis_vector(c, b, p, d) *= orientations(c, b);
1336 
1337  if(compute_derivatives) {
1338  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(curl_basis_vector,orientations);
1339  for (int c=0; c<num_cell; c++)
1340  for (int b=0; b<num_basis; b++)
1341  for (int p=0; p<num_ip; p++)
1342  for (int d=0; d<num_dim; d++)
1343  curl_basis_vector(c, b, p,d) *= orientations(c, b);
1344  }
1345 
1346  // setup the orientations for the test space
1347  if(build_weighted) {
1348  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view());
1349 
1350  if(compute_derivatives)
1351  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::applyFieldSigns(weighted_curl_basis_vector.get_view(),orientations.get_view());
1352  }
1353  }
1354  else if(elmtspace==PureBasis::HDIV) {
1355  // setup the orientations for the trial space
1356  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(basis_vector,orientations);
1357 
1358  for (int c=0; c<num_cell; c++)
1359  for (int b=0; b<num_basis; b++)
1360  for (int p=0; p<num_ip; p++)
1361  for (int d=0; d<num_dim; d++)
1362  basis_vector(c, b, p, d) *= orientations(c, b);
1363 
1364  if(compute_derivatives) {
1365  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(div_basis,orientations);
1366 
1367  for (int c=0; c<num_cell; c++)
1368  for (int b=0; b<num_basis; b++)
1369  for (int p=0; p<num_ip; p++)
1370  div_basis(c, b, p) *= orientations(c, b);
1371  }
1372 
1373  // setup the orientations for the test space
1374  if(build_weighted) {
1375  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::applyFieldSigns(weighted_basis_vector.get_view(),orientations.get_view());
1376 
1377  if(compute_derivatives)
1378  Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::applyFieldSigns(weighted_div_basis.get_view(),orientations.get_view());
1379  }
1380  }
1381 }
1382 
1383 template <typename Scalar>
1385 { return basis_layout->getBasis()->getElementSpace(); }
1386 
1387 template <typename Scalar>
1389 setupArrays(const Teuchos::RCP<const panzer::BasisIRLayout>& layout,
1390  bool computeDerivatives)
1391 {
1392  MDFieldArrayFactory af(prefix,alloc_arrays);
1393 
1394  compute_derivatives = computeDerivatives;
1395  basis_layout = layout;
1396  Teuchos::RCP<const panzer::PureBasis> basisDesc = layout->getBasis();
1397 
1398  // for convience pull out basis and quadrature information
1399  int num_quad = layout->numPoints();
1400  int dim = basisDesc->dimension();
1401  int card = basisDesc->cardinality();
1402  int numcells = basisDesc->numCells();
1403  panzer::PureBasis::EElementSpace elmtspace = basisDesc->getElementSpace();
1404  Teuchos::RCP<const shards::CellTopology> cellTopo = basisDesc->getCellTopology();
1405 
1406  intrepid_basis = basisDesc->getIntrepid2Basis<PHX::Device::execution_space,Scalar,Scalar>();
1407 
1408  // allocate field containers
1409  // field sizes defined by http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/doc/html/basis_page.html#basis_md_array_sec
1410 
1411  // compute basis fields
1412  if(elmtspace==panzer::PureBasis::HGRAD) {
1413  // HGRAD is a nodal field
1414 
1415  // build values
1417  basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("basis_ref",card,num_quad); // F, P
1418  basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("basis",numcells,card,num_quad);
1419 
1420  if(build_weighted)
1421  weighted_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_basis",numcells,card,num_quad);
1422 
1423  // build gradients
1425 
1426  if(compute_derivatives) {
1427  grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("grad_basis_ref",card,num_quad,dim); // F, P, D
1428  grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("grad_basis",numcells,card,num_quad,dim);
1429 
1430  if(build_weighted)
1431  weighted_grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_grad_basis",numcells,card,num_quad,dim);
1432  }
1433 
1434  // build curl
1436 
1437  // nothing - HGRAD does not support CURL operation
1438  }
1439  else if(elmtspace==panzer::PureBasis::HCURL) {
1440  // HCURL is a vector field
1441 
1442  // build values
1444 
1445  basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("basis_ref",card,num_quad,dim); // F, P, D
1446  basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis",numcells,card,num_quad,dim);
1447 
1448  if(build_weighted)
1449  weighted_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_basis",numcells,card,num_quad,dim);
1450 
1451  // build gradients
1453 
1454  // nothing - HCURL does not support GRAD operation
1455 
1456  // build curl
1458 
1459  if(compute_derivatives) {
1460  if(dim==2) {
1461  // curl of HCURL basis is not dimension dependent
1462  curl_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("curl_basis_ref",card,num_quad); // F, P
1463  curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("curl_basis",numcells,card,num_quad);
1464 
1465  if(build_weighted)
1466  weighted_curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_curl_basis",numcells,card,num_quad);
1467  }
1468  else if(dim==3){
1469  curl_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("curl_basis_ref",card,num_quad,dim); // F, P, D
1470  curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("curl_basis",numcells,card,num_quad,dim);
1471 
1472  if(build_weighted)
1473  weighted_curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_curl_basis",numcells,card,num_quad,dim);
1474  }
1475  else { TEUCHOS_ASSERT(false); } // what do I do with 1D?
1476  }
1477  }
1478  else if(elmtspace==panzer::PureBasis::HDIV) {
1479  // HDIV is a vector field
1480 
1481  // build values
1483 
1484  basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("basis_ref",card,num_quad,dim); // F, P, D
1485  basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis",numcells,card,num_quad,dim);
1486 
1487  if(build_weighted)
1488  weighted_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_basis",numcells,card,num_quad,dim);
1489 
1490  // build gradients
1492 
1493  // nothing - HCURL does not support GRAD operation
1494 
1495  // build curl
1497 
1498  // nothing - HDIV does not support CURL operation
1499 
1500  // build div
1502 
1503  if(compute_derivatives) {
1504  div_basis_ref = af.buildStaticArray<Scalar,BASIS,IP>("div_basis_ref",card,num_quad); // F, P
1505  div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("div_basis",numcells,card,num_quad);
1506 
1507  if(build_weighted)
1508  weighted_div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_div_basis",numcells,card,num_quad);
1509  }
1510  }
1511  else if(elmtspace==panzer::PureBasis::CONST) {
1512  // CONST is a nodal field
1513 
1514  // build values
1516  basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("basis_ref",card,num_quad); // F, P
1517  basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("basis",numcells,card,num_quad);
1518 
1519  if(build_weighted)
1520  weighted_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_basis",numcells,card,num_quad);
1521 
1522  // build gradients
1524 
1525  // nothing - CONST does not support GRAD operation
1526 
1527  // build curl
1529 
1530  // nothing - CONST does not support CURL operation
1531 
1532  // build div
1534 
1535  // nothing - CONST does not support DIV operation
1536  }
1537  else { TEUCHOS_ASSERT(false); }
1538 
1539  basis_coordinates_ref = af.buildStaticArray<Scalar,BASIS,Dim>("basis_coordinates_ref",card,dim);
1540  basis_coordinates = af.buildStaticArray<Scalar,Cell,BASIS,Dim>("basis_coordinates",numcells,card,dim);
1541 }
1542 
1543 // do some explicit instantiation so things compile faster.
1544 
1545 #define BASIS_VALUES_INSTANTIATION(SCALAR) \
1546 template class BasisValues2<SCALAR>;
1547 
1550 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
1552 #endif
1553 
1554 } // namespace panzer
void evaluateValues_HGrad(const PHX::MDField< Scalar, Cell, IP, Dim, void, void, void, void, void > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv, const PHX::MDField< Scalar, Cell, IP > &weighted_measure)
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
void evaluateValues(const PHX::MDField< Scalar, IP, Dim, void, void, void, void, void, void > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv)
void evaluateReferenceValues(const PHX::MDField< Scalar, IP, Dim > &cub_points, bool compute_derivatives, bool use_vertex_coordinates)
ArrayTraits< Scalar, PHX::MDField< Scalar, void, void, void, void, void, void, void, void > >::size_type size_type
void evaluateBasisCoordinates(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates)
PHX::MDField< Scalar > buildArray(const std::string &str, int d0) const
void evaluateValues_Const(const PHX::MDField< Scalar, Cell, IP, Dim, void, void, void, void, void > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv, const PHX::MDField< Scalar, Cell, IP > &weighted_measure)
void applyOrientations(const PHX::MDField< const Scalar, Cell, BASIS > &orientations)
Method to apply orientations to a basis values container.
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 evaluateValues_HCurl(const PHX::MDField< Scalar, Cell, IP, Dim, void, void, void, void, void > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv, const PHX::MDField< Scalar, Cell, IP > &weighted_measure)
#define BASIS_VALUES_INSTANTIATION(SCALAR)
PANZER_FADTYPE FadType
void setupArrays(const Teuchos::RCP< const panzer::BasisIRLayout > &basis, bool computeDerivatives=true)
Sizes/allocates memory for arrays.
Sacado::Fad::DFad< Sacado::Fad::SFad< RealType, 1 > > HessianType
Description and data layouts associated with a particular basis.
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, IP, Dim, void, void, void, void, void > &cell_cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv)
PureBasis::EElementSpace getElementSpace() const
int cardinality() const
Returns the number of basis coefficients.
void evaluateValues_HDiv(const PHX::MDField< Scalar, Cell, IP, Dim, void, void, void, void, void > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP > &weighted_measure)
PHX::MDField< Scalar > ArrayDynamic