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