Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_IntegrationValues2.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
45
46#include "Shards_CellTopology.hpp"
47
48#include "Kokkos_DynRankView.hpp"
49#include "Intrepid2_FunctionSpaceTools.hpp"
50#include "Intrepid2_RealSpaceTools.hpp"
51#include "Intrepid2_CellTools.hpp"
52#include "Intrepid2_ArrayTools.hpp"
53#include "Intrepid2_CubatureControlVolume.hpp"
54#include "Intrepid2_CubatureControlVolumeSide.hpp"
55#include "Intrepid2_CubatureControlVolumeBoundary.hpp"
56
58#include "Panzer_Traits.hpp"
61
62// FIXME: There are some calls in Intrepid2 that require non-const arrays when they should be const - search for PHX::getNonConstDynRankViewFromConstMDField
63#include "Phalanx_GetNonConstDynRankViewFromConstMDField.hpp"
64
65namespace panzer {
66
67namespace {
68
69Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double>>
70getIntrepidCubature(const panzer::IntegrationRule & ir)
71{
73 Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double> > ic;
74
75 Intrepid2::DefaultCubatureFactory cubature_factory;
76
77 if(ir.getType() == ID::CV_SIDE){
78 ic = Teuchos::rcp(new Intrepid2::CubatureControlVolumeSide<PHX::Device::execution_space,double,double>(*ir.topology));
79 } else if(ir.getType() == ID::CV_VOLUME){
80 ic = Teuchos::rcp(new Intrepid2::CubatureControlVolume<PHX::Device::execution_space,double,double>(*ir.topology));
81 } else if(ir.getType() == ID::CV_BOUNDARY){
82 TEUCHOS_ASSERT(ir.isSide());
83 ic = Teuchos::rcp(new Intrepid2::CubatureControlVolumeBoundary<PHX::Device::execution_space,double,double>(*ir.topology,ir.getSide()));
84 } else if(ir.getType() == ID::VOLUME){
85 ic = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir.topology),ir.getOrder());
86 } else if(ir.getType() == ID::SIDE){
87 ic = cubature_factory.create<PHX::Device::execution_space,double,double>(*(ir.side_topology),ir.getOrder());
88 } else if(ir.getType() == ID::SURFACE){
89 // closed surface integrals don't exist in intrepid.
90 } else {
91 TEUCHOS_ASSERT(false);
92 }
93
94 return ic;
95}
96
97template<typename Scalar>
98void
99correctVirtualNormals(PHX::MDField<Scalar,Cell,IP,Dim> normals,
100 const int num_real_cells,
101 const int num_virtual_cells,
102 const shards::CellTopology & cell_topology,
103 const SubcellConnectivity & face_connectivity)
104{
105
106 // What we want is for the adjoining face of the virtual cell to have normals that are the negated real cell's normals.
107 // we correct the normals here:
108 const int space_dim = cell_topology.getDimension();
109 const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1);
110 const int points = normals.extent_int(1);
111 const int points_per_face = points / faces_per_cell;
112
113 Kokkos::parallel_for(num_virtual_cells, KOKKOS_LAMBDA(const int & virtual_cell_ordinal){
114
115 const int virtual_cell = virtual_cell_ordinal+num_real_cells;
116
117 // Find the face and local face ids for the given virtual cell
118 // Note that virtual cells only connect to the owned cells through a single face
119 int face, virtual_lidx;
120 for (int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
121 // Faces that exist have positive indexes
122 face = face_connectivity.subcellForCell(virtual_cell, local_face_id);
123 if (face >= 0){
124 virtual_lidx = local_face_id;
125 break;
126 }
127 }
128
129 // Indexes for real cell
130 const int real_side = (face_connectivity.cellForSubcell(face, 0) == virtual_cell) ? 1 : 0;
131 const int real_cell = face_connectivity.cellForSubcell(face,real_side);
132 const int real_lidx = face_connectivity.localSubcellForSubcell(face,real_side);
133
134 // Make sure it is a real cell (it should actually be an owned cell)
135 KOKKOS_ASSERT(real_cell < num_real_cells);
136
137 for(int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
138
139 // Point indexes for virtual and real point on face
140 const int virtual_cell_point = points_per_face * virtual_lidx + point_ordinal;
141 const int real_cell_point = points_per_face * real_lidx + point_ordinal;
142
143 for (int d=0; d<space_dim; d++)
144 normals(virtual_cell,virtual_cell_point,d) = -normals(real_cell,real_cell_point,d);
145
146 }
147
148 // Clear other normals
149 for (int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
150
151 if (local_face_id == virtual_lidx) continue;
152
153 for (int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
154 const int point = local_face_id * points_per_face + point_ordinal;
155 for (int dim=0; dim<space_dim; dim++)
156 normals(virtual_cell,point,dim) = 0.0;
157 }
158 }
159 });
160 PHX::Device::execution_space().fence();
161}
162
163
164template<typename Scalar>
165void
166correctVirtualRotationMatrices(PHX::MDField<Scalar,Cell,IP,Dim,Dim> rotation_matrices,
167 const int num_real_cells,
168 const int num_virtual_cells,
169 const shards::CellTopology & cell_topology,
170 const SubcellConnectivity & face_connectivity)
171{
172
173 // What we want is for the adjoining face of the virtual cell to have normals that are the negated real cell's normals.
174 // we correct the normals here:
175 const int space_dim = cell_topology.getDimension();
176 const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1);
177 const int points = rotation_matrices.extent_int(1);
178 const int points_per_face = points / faces_per_cell;
179
180 Kokkos::parallel_for(num_virtual_cells, KOKKOS_LAMBDA(const int & virtual_cell_ordinal){
181
182 const int virtual_cell = virtual_cell_ordinal+num_real_cells;
183
184 // Find the face and local face ids for the given virtual cell
185 // Note that virtual cells only connect to the owned cells through a single face
186 int face, virtual_lidx;
187 for (int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
188 // Faces that exist have positive indexes
189 face = face_connectivity.subcellForCell(virtual_cell, local_face_id);
190 if (face >= 0){
191 virtual_lidx = local_face_id;
192 break;
193 }
194 }
195
196 // The normals already have the correction applied, so we just need to zero out the rotation matrices on the other faces
197
198 // Clear other rotation matrices
199 for (int local_face_id=0; local_face_id<faces_per_cell; local_face_id++){
200
201 if (local_face_id == virtual_lidx) continue;
202
203 for (int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++){
204 const int point = local_face_id * points_per_face + point_ordinal;
205 for (int dim=0; dim<3; dim++)
206 for (int dim2=0; dim2<3; dim2++)
207 rotation_matrices(virtual_cell,point,dim,dim2) = 0.0;
208 }
209 }
210 });
211 PHX::Device::execution_space().fence();
212}
213
214template<typename Scalar>
215void
216applyBasePermutation(PHX::MDField<Scalar,IP> field,
217 PHX::MDField<const int,Cell,IP> permutations)
218{
219 MDFieldArrayFactory af("",true);
220
221 const int num_ip = field.extent(0);
222
223 auto scratch = af.template buildStaticArray<Scalar,IP>("scratch", num_ip);
224 scratch.deep_copy(field);
225 Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int & ){
226 for(int ip=0; ip<num_ip; ++ip)
227 if (ip != permutations(0,ip))
228 field(ip) = scratch(permutations(0,ip));
229 });
230 PHX::Device::execution_space().fence();
231}
232
233template<typename Scalar>
234void
235applyBasePermutation(PHX::MDField<Scalar,IP,Dim> field,
236 PHX::MDField<const int,Cell,IP> permutations)
237{
238 MDFieldArrayFactory af("",true);
239
240 const int num_ip = field.extent(0);
241 const int num_dim = field.extent(1);
242
243 auto scratch = af.template buildStaticArray<Scalar,IP,Dim>("scratch", num_ip,num_dim);
244 scratch.deep_copy(field);
245 Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int & ){
246 for(int ip=0; ip<num_ip; ++ip)
247 if (ip != permutations(0,ip))
248 for(int dim=0; dim<num_dim; ++dim)
249 field(ip,dim) = scratch(permutations(0,ip),dim);
250 });
251 PHX::Device::execution_space().fence();
252}
253
254template<typename Scalar>
255void
256applyPermutation(PHX::MDField<Scalar,Cell,IP> field,
257 PHX::MDField<const int,Cell,IP> permutations)
258{
259 MDFieldArrayFactory af("",true);
260
261 const int num_cells = field.extent(0);
262 const int num_ip = field.extent(1);
263
264 auto scratch = af.template buildStaticArray<Scalar,Cell,IP>("scratch", num_cells, num_ip);
265 scratch.deep_copy(field);
266 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){
267 for(int ip=0; ip<num_ip; ++ip)
268 if (ip != permutations(cell,ip))
269 field(cell,ip) = scratch(cell,permutations(cell,ip));
270 });
271 PHX::Device::execution_space().fence();
272}
273
274template<typename Scalar>
275void
276applyPermutation(PHX::MDField<Scalar,Cell,IP,Dim> field,
277 PHX::MDField<const int,Cell,IP> permutations)
278{
279 MDFieldArrayFactory af("",true);
280
281 const int num_cells = field.extent(0);
282 const int num_ip = field.extent(1);
283 const int num_dim = field.extent(2);
284
285 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim>("scratch", num_cells, num_ip, num_dim);
286 scratch.deep_copy(field);
287 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){
288 for(int ip=0; ip<num_ip; ++ip)
289 if (ip != permutations(cell,ip))
290 for(int dim=0; dim<num_dim; ++dim)
291 field(cell,ip,dim) = scratch(cell,permutations(cell,ip),dim);
292 });
293 PHX::Device::execution_space().fence();
294}
295
296template<typename Scalar>
297void
298applyPermutation(PHX::MDField<Scalar,Cell,IP,Dim,Dim> field,
299 PHX::MDField<const int,Cell,IP> permutations)
300{
301 MDFieldArrayFactory af("",true);
302
303 const int num_cells = field.extent(0);
304 const int num_ip = field.extent(1);
305 const int num_dim = field.extent(2);
306 const int num_dim2 = field.extent(3);
307
308 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("scratch", num_cells, num_ip, num_dim, num_dim2);
309 scratch.deep_copy(field);
310 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){
311 for(int ip=0; ip<num_ip; ++ip)
312 if (ip != permutations(cell,ip))
313 for(int dim=0; dim<num_dim; ++dim)
314 for(int dim2=0; dim2<num_dim2; ++dim2)
315 field(cell,ip,dim,dim2) = scratch(cell,permutations(cell,ip),dim,dim2);
316 });
317 PHX::Device::execution_space().fence();
318}
319
320
321// Find the permutation that maps the set of points coords to other_coords. To
322// avoid possible finite precision issues, == is not used, but rather
323// min(norm(.)).
324template <typename Scalar>
325PHX::MDField<const int,Cell,IP>
326generatePermutations(const int num_cells,
327 PHX::MDField<const Scalar,Cell,IP,Dim> coords,
328 PHX::MDField<const Scalar,Cell,IP,Dim> other_coords)
329{
330 const int num_ip = coords.extent(1);
331 const int num_dim = coords.extent(2);
332
333 MDFieldArrayFactory af("",true);
334
335 // This will store the permutations to go from coords to other_coords
336 auto permutation = af.template buildStaticArray<int,Cell,IP>("permutation", num_cells, num_ip);
337 permutation.deep_copy(0);
338
339 // This is scratch space - there is likely a better way to do this
340 auto taken = af.template buildStaticArray<int,Cell,IP>("taken", num_cells, num_ip);
341 taken.deep_copy(0);
342
343 Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA(const int & cell){
344
345 for (int ip = 0; ip < num_ip; ++ip) {
346 // Find an other point to associate with ip.
347 int i_min = 0;
348 Scalar d_min = -1;
349 for (int other_ip = 0; other_ip < num_ip; ++other_ip) {
350 // For speed, skip other points that are already associated.
351 if(taken(cell,other_ip)) continue;
352 // Compute the distance between the two points.
353 Scalar d(0);
354 for (int dim = 0; dim < num_dim; ++dim) {
355 const Scalar diff = coords(cell, ip, dim) - other_coords(cell, other_ip, dim);
356 d += diff*diff;
357 }
358 if (d_min < 0 || d < d_min) {
359 d_min = d;
360 i_min = other_ip;
361 }
362 }
363 // Record the match.
364 permutation(cell,ip) = i_min;
365 // This point is now taken.
366 taken(cell,i_min) = 1;
367 }
368 });
369 PHX::Device::execution_space().fence();
370
371 return permutation;
372
373}
374
375template <typename Scalar>
376PHX::MDField<const int,Cell,IP>
377generateSurfacePermutations(const int num_cells,
378 const SubcellConnectivity face_connectivity,
379 PHX::MDField<const Scalar,Cell,IP,Dim> surface_points,
380 PHX::MDField<const Scalar,Cell,IP,Dim,Dim> surface_rotation_matrices)
381
382{
383
384 // The challenge for this call is handling wedge-based periodic boundaries
385 // We need to make sure that we can align points along faces that are rotated with respect to one another.
386 // Below we will see an algorithm that rotates two boundaries into a shared frame, then figures out
387 // how to permute the points on one face to line up with the other.
388
389 const int num_points = surface_points.extent_int(1);
390 const int num_faces_per_cell = face_connectivity.numSubcellsOnCellHost(0);
391 const int num_points_per_face = num_points / num_faces_per_cell;
392 const int cell_dim = surface_points.extent(2);
393
394 MDFieldArrayFactory af("",true);
395
396 // This will store the permutations
397 auto permutation = af.template buildStaticArray<int,Cell,IP>("permutation", num_cells, num_points);
398 permutation.deep_copy(0);
399
400 // Fill permutations with trivial values (i.e. no permutation - this will get overwritten for some cells)
401 Kokkos::parallel_for(num_cells,KOKKOS_LAMBDA (const int & cell) {
402 for(int point=0; point<num_points; ++point)
403 permutation(cell,point) = point;
404 });
405
406 // Now we need to align the cubature points for each face
407 // If there is only one point there is no need to align things
408 if(num_points_per_face != 1) {
409 // To enforce that quadrature points on faces are aligned properly we will iterate through faces,
410 // map points to a plane shared by the faces, then re-order quadrature points on the "1" face to
411 // line up with the "0" face
412
413 // Utility calls
414#define PANZER_DOT(a,b) (a[0]*b[0] + a[1]*b[1] + a[2]*b[2])
415#define PANZER_CROSS(a,b,c) {c[0] = a[1]*b[2] - a[2]*b[1]; c[1] = a[2]*b[0] - a[0]*b[2]; c[2] = a[0]*b[1] - a[1]*b[0];}
416
417 // Iterate through faces
418 Kokkos::parallel_for("face iteration",face_connectivity.numSubcells(),KOKKOS_LAMBDA (const int face) {
419 // Cells for sides 0 and 1
420 const int cell_0 = face_connectivity.cellForSubcell(face,0);
421 const int cell_1 = face_connectivity.cellForSubcell(face,1);
422
423 // If this face doesn't connect to anything we don't need to worry about alignment
424 if(cell_1 < 0)
425 return;
426
427 // Local face index for sides 0 and 1
428 const int lidx_0 = face_connectivity.localSubcellForSubcell(face,0);
429 const int lidx_1 = face_connectivity.localSubcellForSubcell(face,1);
430
431 // If the cell exists, then the local face index needs to exist
432 KOKKOS_ASSERT(lidx_1 >= 0);
433
434 // To compare points on planes, it is good to center the planes around the origin
435 // Find the face center for the left and right cell (may not be the same point - e.g. periodic conditions)
436 // We also define a length scale r2 which gives an order of magnitude approximation to the cell size squared
437 Scalar xc0[3] = {0}, xc1[3] = {0}, r2 = 0;
438 for(int face_point=0; face_point<num_points_per_face; ++face_point){
439 Scalar dx2 = 0.;
440 for(int dim=0; dim<cell_dim; ++dim){
441 xc0[dim] += surface_points(cell_0,lidx_0*num_points_per_face + face_point,dim);
442 xc1[dim] += surface_points(cell_1,lidx_1*num_points_per_face + face_point,dim);
443 const Scalar dx = surface_points(cell_0,lidx_0*num_points_per_face + face_point,dim) - surface_points(cell_0,lidx_0*num_points_per_face,dim);
444 dx2 += dx*dx;
445 }
446
447 // Check if the distance squared between these two points is larger than the others (doesn't need to be perfect)
448 r2 = (r2 < dx2) ? dx2 : r2;
449 }
450 for(int dim=0; dim<cell_dim; ++dim){
451 xc0[dim] /= double(num_points_per_face);
452 xc1[dim] /= double(num_points_per_face);
453 }
454
455 // TODO: This needs to be adaptable to having curved faces - for now this will work
456
457 // We need to define a plane with two vectors (transverse and binormal)
458 // These will be used with the face centers to construct a local reference frame for aligning points
459
460 // We use the first point on the face to define the normal (may break for curved faces)
461 const int example_point_0 = lidx_0*num_points_per_face;
462 const int example_point_1 = lidx_1*num_points_per_face;
463
464 // Load the transverse and binormal for the 0 cell (default)
465 Scalar t[3] = {surface_rotation_matrices(cell_0,example_point_0,1,0), surface_rotation_matrices(cell_0,example_point_0,1,1), surface_rotation_matrices(cell_0,example_point_0,1,2)};
466 Scalar b[3] = {surface_rotation_matrices(cell_0,example_point_0,2,0), surface_rotation_matrices(cell_0,example_point_0,2,1), surface_rotation_matrices(cell_0,example_point_0,2,2)};
467
468 // In case the faces are not antiparallel (e.g. periodic wedge), we may need to change the transverse and binormal
469 {
470
471 // Get the normals for each face for constructing one of the plane vectors (transverse)
472 const Scalar n0[3] = {surface_rotation_matrices(cell_0,example_point_0,0,0), surface_rotation_matrices(cell_0,example_point_0,0,1), surface_rotation_matrices(cell_0,example_point_0,0,2)};
473 const Scalar n1[3] = {surface_rotation_matrices(cell_1,example_point_1,0,0), surface_rotation_matrices(cell_1,example_point_1,0,1), surface_rotation_matrices(cell_1,example_point_1,0,2)};
474
475 // n_0*n_1 == -1 (antiparallel), n_0*n_1 == 1 (parallel - bad), |n_0*n_1| < 1 (other)
476 const Scalar n0_dot_n1 = PANZER_DOT(n0,n1);
477
478 // FIXME: Currently virtual cells will set their surface normal along the same direction as the cell they "reflect"
479 // This causes a host of issues (e.g. identifying 180 degree periodic wedges), but we have to support virtual cells as a priority
480 // Therefore, we will just assume that the ordering is fine (not valid for 180 degree periodic wedges)
481 if(Kokkos::fabs(n0_dot_n1 - 1.) < 1.e-8)
482 return;
483
484 // Rotated faces case (e.g. periodic wedge) we need to check for non-antiparallel face normals
485 if(Kokkos::fabs(n0_dot_n1 + 1.) > 1.e-2){
486
487 // Now we need to define an arbitrary transverse and binormal in the plane across which the faces are anti-symmetric
488 // We can do this by setting t = n_0 \times n_1
489 PANZER_CROSS(n0,n1,t);
490
491 // Normalize the transverse vector
492 const Scalar mag_t = Kokkos::sqrt(PANZER_DOT(t,t));
493 t[0] /= mag_t;
494 t[1] /= mag_t;
495 t[2] /= mag_t;
496
497 // The binormal will be the sum of the normals (does not need to be a right handed system)
498 b[0] = n0[0] + n1[0];
499 b[1] = n0[1] + n1[1];
500 b[2] = n0[2] + n1[2];
501
502 // Normalize the binormal vector
503 const Scalar mag_b = Kokkos::sqrt(PANZER_DOT(b,b));
504 b[0] /= mag_b;
505 b[1] /= mag_b;
506 b[2] /= mag_b;
507
508 }
509 }
510
511 // Now that we have a reference coordinate plane in which to align our points
512 // Point on the transverse/binormal plane
513 Scalar p0[2] = {0};
514 Scalar p1[2] = {0};
515
516 // Differential position in mesh
517 Scalar x0[3] = {0};
518 Scalar x1[3] = {0};
519
520 // Iterate through points in cell 1 and find which point they align with in cell 0
521 for(int face_point_1=0; face_point_1<num_points_per_face; ++face_point_1){
522
523 // Get the point index in the 1 cell
524 const int point_1 = lidx_1*num_points_per_face + face_point_1;
525
526 // Load position shifted by face center
527 for(int dim=0; dim<cell_dim; ++dim)
528 x1[dim] = surface_points(cell_1,point_1,dim) - xc1[dim];
529
530 // Convert position to transverse/binormal plane
531 p1[0] = PANZER_DOT(x1,t);
532 p1[1] = PANZER_DOT(x1,b);
533
534 // Compare to points on the other surface
535 for(int face_point_0=0; face_point_0<num_points_per_face; ++face_point_0){
536
537 // Get the point index in the 0 cell
538 const int point_0 = lidx_0*num_points_per_face + face_point_0;
539
540 // Load position shifted by face center
541 for(int dim=0; dim<cell_dim; ++dim)
542 x0[dim] = surface_points(cell_0,point_0,dim) - xc0[dim];
543
544 // Convert position to transverse/binormal plane
545 p0[0] = PANZER_DOT(x0,t);
546 p0[1] = PANZER_DOT(x0,b);
547
548 // Find the distance squared between p0 and p1
549 const Scalar p012 = (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]);
550
551 // TODO: Should this be a constant value, or should we just find the minimum point?
552 // If the distance, compared to the size of the cell, is small, we assume these are the same points
553 if(p012 / r2 < 1.e-12){
554 permutation(cell_1,lidx_1*num_points_per_face+face_point_0) = point_1;
555 break;
556 }
557
558 }
559 }
560 });
561 PHX::Device::execution_space().fence();
562 }
563
564#undef PANZER_DOT
565#undef PANZER_CROSS
566
567 return permutation;
568
569}
570
571}
572
573//template<typename DataType>
574//using UnmanagedDynRankView = Kokkos::DynRankView<DataType,typename PHX::DevLayout<DataType>::type,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
575
576template <typename Scalar>
578IntegrationValues2(const std::string & pre,
579 const bool allocArrays):
580 num_cells_(0),
581 num_evaluate_cells_(0),
582 num_virtual_cells_(-1),
583 requires_permutation_(false),
584 alloc_arrays_(allocArrays),
585 prefix_(pre),
586 ddims_(1,0)
587{
588 resetArrays();
589}
590
591template <typename Scalar>
592void
595{
596 cub_points_evaluated_ = false;
597 side_cub_points_evaluated_ = false;
598 cub_weights_evaluated_ = false;
599 node_coordinates_evaluated_ = false;
600 jac_evaluated_ = false;
601 jac_inv_evaluated_ = false;
602 jac_det_evaluated_ = false;
603 weighted_measure_evaluated_ = false;
604 weighted_normals_evaluated_ = false;
605 surface_normals_evaluated_ = false;
606 surface_rotation_matrices_evaluated_ = false;
607 covarient_evaluated_ = false;
608 contravarient_evaluated_ = false;
609 norm_contravarient_evaluated_ = false;
610 ip_coordinates_evaluated_ = false;
611 ref_ip_coordinates_evaluated_ = false;
612
613 // TODO: We need to clear the views
614}
615
616template <typename Scalar>
618setupArrays(const Teuchos::RCP<const panzer::IntegrationRule>& ir)
619{
620 MDFieldArrayFactory af(prefix_,alloc_arrays_);
621
623
624 int_rule = ir;
625
626 const int num_nodes = ir->topology->getNodeCount();
627 const int num_cells = ir->workset_size;
628 const int num_space_dim = ir->topology->getDimension();
629
630 // Specialize content if this is quadrature at a node
631 const bool is_node_rule = (num_space_dim==1) and ir->isSide();
632 if(not is_node_rule) {
633 TEUCHOS_ASSERT(ir->getType() != ID::NONE);
634 intrepid_cubature = getIntrepidCubature(*ir);
635 }
636
637 const int num_ip = is_node_rule ? 1 : ir->num_points;
638
639 cub_points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
640
641 if (ir->isSide() && ir->cv_type == "none")
642 side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip,ir->side_topology->getDimension());
643
644 cub_weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
645
646 node_coordinates = af.template buildStaticArray<Scalar,Cell,BASIS,Dim>("node_coordinates",num_cells, num_nodes, num_space_dim);
647
648 jac = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells, num_ip, num_space_dim,num_space_dim);
649
650 jac_inv = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells, num_ip, num_space_dim,num_space_dim);
651
652 jac_det = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells, num_ip);
653
654 weighted_measure = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells, num_ip);
655
656 covarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells, num_ip, num_space_dim,num_space_dim);
657
658 contravarient = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells, num_ip, num_space_dim,num_space_dim);
659
660 norm_contravarient = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells, num_ip);
661
662 ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordiantes",num_cells, num_ip,num_space_dim);
663
664 ref_ip_coordinates = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells, num_ip,num_space_dim);
665
666 weighted_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted_normal",num_cells,num_ip,num_space_dim);
667
668 surface_normals = af.template buildStaticArray<Scalar,Cell,IP,Dim>("surface_normals",num_cells, num_ip,num_space_dim);
669
670 surface_rotation_matrices = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("surface_rotation_matrices",num_cells, num_ip,3,3);
671
672}
673
674
675// ***********************************************************
676// * Evaluation of values - NOT specialized
677// ***********************************************************
678template <typename Scalar>
680evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim> & in_node_coordinates,
681 const int in_num_cells,
682 const Teuchos::RCP<const SubcellConnectivity> & face_connectivity,
683 const int num_virtual_cells)
684{
685
686 setup(int_rule, in_node_coordinates, in_num_cells);
687
688 // Setup permutations (only required for surface integrators)
690 if((int_rule->getType() == ID::SURFACE) and (face_connectivity != Teuchos::null))
691 setupPermutations(face_connectivity, num_virtual_cells);
692
693 // Evaluate everything once permutations are generated
694 evaluateEverything();
695
696}
697
698template <typename Scalar>
699void
701evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim>& in_node_coordinates,
702 const PHX::MDField<Scalar,Cell,IP,Dim>& other_ip_coordinates,
703 const int in_num_cells)
704{
705
706 setup(int_rule, in_node_coordinates, in_num_cells);
707
708 // Setup permutations
709 setupPermutations(other_ip_coordinates);
710
711 // Evaluate everything once permutations are generated
712 evaluateEverything();
713
714}
715
716template <typename Scalar>
717void
719setupPermutations(const Teuchos::RCP<const SubcellConnectivity> & face_connectivity,
720 const int num_virtual_cells)
721{
722 TEUCHOS_ASSERT(not int_rule->isSide());
723 TEUCHOS_ASSERT(face_connectivity != Teuchos::null);
724 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() != panzer::IntegrationDescriptor::SURFACE,
725 "IntegrationValues2::setupPermutations : Face connectivity based permutations are only required for surface integration schemes");
726 TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ >= 0,
727 "IntegrationValues2::setupPermutations : Number of virtual cells for surface permuations must be >= 0");
728 resetArrays();
729 side_connectivity_ = face_connectivity;
730 num_virtual_cells_ = num_virtual_cells;
731 requires_permutation_ = false;
732 permutations_ = generateSurfacePermutations<Scalar>(num_evaluate_cells_,*face_connectivity, getCubaturePoints(false,true), getSurfaceRotationMatrices(false,true));
733 requires_permutation_ = true;
734}
735
736template <typename Scalar>
737void
739setupPermutations(const PHX::MDField<Scalar,Cell,IP,Dim> & other_ip_coordinates)
740{
741 resetArrays();
742 requires_permutation_ = false;
743 permutations_ = generatePermutations<Scalar>(num_evaluate_cells_, getCubaturePoints(false,true), other_ip_coordinates);
744 requires_permutation_ = true;
745}
746
747
748template <typename Scalar>
749void
751setup(const Teuchos::RCP<const panzer::IntegrationRule>& ir,
752 const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
753 const int num_cells)
754{
755
756 // Clear arrays just in case we are rebuilding this object
757 resetArrays();
758
759 num_cells_ = vertex_coordinates.extent(0);
760 num_evaluate_cells_ = num_cells < 0 ? vertex_coordinates.extent(0) : num_cells;
761 int_rule = ir;
762
763 TEUCHOS_ASSERT(ir->getType() != IntegrationDescriptor::NONE);
764 intrepid_cubature = getIntrepidCubature(*ir);
765
766 // Copy node coordinates into owned allocation
767 {
768 MDFieldArrayFactory af(prefix_,true);
769
770 const int num_space_dim = int_rule->topology->getDimension();
771 const int num_vertexes = int_rule->topology->getNodeCount();
772 TEUCHOS_ASSERT(static_cast<int>(vertex_coordinates.extent(1)) == num_vertexes);
773
774 auto aux = af.template buildStaticArray<Scalar,Cell,NODE,Dim>("node_coordinates",num_cells_,num_vertexes, num_space_dim);
775 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_vertexes,num_space_dim});
776 Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int & cell, const int & point, const int & dim){
777 aux(cell,point,dim) = vertex_coordinates(cell,point,dim);
778 });
779 PHX::Device::execution_space().fence();
780 node_coordinates = aux;
781 }
782
783}
784
785template <typename Scalar>
788getUniformCubaturePointsRef(const bool cache,
789 const bool force,
790 const bool apply_permutation) const
791{
792
793 if(cub_points_evaluated_ and not force)
794 return cub_points;
795
796 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
797 MDFieldArrayFactory af(prefix_,true);
798
799 int num_space_dim = int_rule->topology->getDimension();
800 int num_ip = int_rule->num_points;
801
802 auto aux = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip, num_space_dim);
803
804 if (int_rule->isSide() && num_space_dim==1) {
805 std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
806 << "non-natural integration rules.";
807 return aux;
808 }
809
810 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none",
811 "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme.");
812
813 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE,
814 "IntegrationValues2::getUniformCubaturePointsRef : Cannot build reference cubature points for surface integration scheme.");
815
816 auto weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
817
818 if (!int_rule->isSide())
819 intrepid_cubature->getCubature(aux.get_view(), weights.get_view());
820 else {
821 auto s_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip, num_space_dim-1);
822
823 intrepid_cubature->getCubature(s_cub_points.get_view(), weights.get_view());
824 cell_tools.mapToReferenceSubcell(aux.get_view(), s_cub_points.get_view(), num_space_dim-1, int_rule->getSide(), *(int_rule->topology));
825 }
826
827 PHX::Device::execution_space().fence();
828
829 if(apply_permutation and requires_permutation_)
830 applyBasePermutation(aux, permutations_);
831
832 if(cache){
833 cub_points = aux;
834 cub_points_evaluated_ = true;
835 }
836
837 return aux;
838
839}
840
841template <typename Scalar>
844getUniformSideCubaturePointsRef(const bool cache,
845 const bool force,
846 const bool apply_permutation) const
847{
848
849 if(side_cub_points_evaluated_ and not force)
850 return side_cub_points;
851
852 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
853 MDFieldArrayFactory af(prefix_,true);
854
855 int num_space_dim = int_rule->topology->getDimension();
856 int num_ip = int_rule->num_points;
857
858 auto aux = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_ip, num_space_dim-1);
859
860 if (int_rule->isSide() && num_space_dim==1) {
861 std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
862 << "non-natural integration rules.";
863 return aux;
864 }
865
866 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none",
867 "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for control volume integration scheme.");
868
869 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE,
870 "IntegrationValues2::getUniformSideCubaturePointsRef : Cannot build reference cubature points for surface integration scheme.");
871
872 TEUCHOS_TEST_FOR_EXCEPT_MSG(!int_rule->isSide(),
873 "IntegrationValues2::getUniformSideCubaturePointsRef : Requested side points, which is not supported by integration rule.");
874
875 auto weights = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
876
877 intrepid_cubature->getCubature(aux.get_view(), weights.get_view());
878
879 PHX::Device::execution_space().fence();
880
881 if(apply_permutation and requires_permutation_)
882 applyBasePermutation(aux, permutations_);
883
884 if(cache){
885 side_cub_points = aux;
886 side_cub_points_evaluated_ = true;
887 }
888
889 return aux;
890}
891
892template <typename Scalar>
895getUniformCubatureWeightsRef(const bool cache,
896 const bool force,
897 const bool apply_permutation) const
898{
899
900 if(cub_weights_evaluated_ and not force)
901 return cub_weights;
902
903 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
904 MDFieldArrayFactory af(prefix_,true);
905
906 int num_space_dim = int_rule->topology->getDimension();
907 int num_ip = int_rule->num_points;
908
909 auto aux = af.template buildStaticArray<Scalar,IP>("cub_weights",num_ip);
910
911 if (int_rule->isSide() && num_space_dim==1) {
912 std::cout << "WARNING: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
913 << "non-natural integration rules.";
914 return aux;
915 }
916
917 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none",
918 "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for control volume integration scheme.");
919
920 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE,
921 "IntegrationValues2::getUniformCubatureWeightsRef : Cannot build reference cubature weights for surface integration scheme.");
922
923 auto points = af.template buildStaticArray<Scalar,IP,Dim>("cub_points",num_ip,num_space_dim);
924
925 intrepid_cubature->getCubature(points.get_view(), aux.get_view());
926
927 PHX::Device::execution_space().fence();
928
929 if(apply_permutation and requires_permutation_)
930 applyBasePermutation(aux, permutations_);
931
932 if(cache){
933 cub_weights = aux;
934 cub_weights_evaluated_ = true;
935 }
936
937 return aux;
938
939}
940
941template <typename Scalar>
944getNodeCoordinates() const
945{
946 return node_coordinates;
947}
948
949template <typename Scalar>
952getJacobian(const bool cache,
953 const bool force) const
954{
955 if(jac_evaluated_ and not force)
956 return jac;
957
958 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
959 MDFieldArrayFactory af(prefix_,true);
960
961 int num_space_dim = int_rule->topology->getDimension();
962 int num_ip = int_rule->num_points;
963
964 // 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
965 auto const_ref_coord = getCubaturePointsRef(false,force);
966 auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord);
967 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
968 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac",num_cells_, num_ip, num_space_dim,num_space_dim);
969
970 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
971 auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
972 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
973 auto s_jac = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
974
975 cell_tools.setJacobian(s_jac, s_ref_coord, s_node_coord,*(int_rule->topology));
976
977 PHX::Device::execution_space().fence();
978
979 if(cache){
980 jac = aux;
981 jac_evaluated_ = true;
982 }
983
984 return aux;
985
986}
987
988template <typename Scalar>
991getJacobianInverse(const bool cache,
992 const bool force) const
993{
994
995 if(jac_inv_evaluated_ and not force)
996 return jac_inv;
997
998 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
999 MDFieldArrayFactory af(prefix_,true);
1000
1001 const int num_space_dim = int_rule->topology->getDimension();
1002 const int num_ip = int_rule->num_points;
1003
1004 auto jacobian = getJacobian(false,force);
1005 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("jac_inv",num_cells_, num_ip, num_space_dim,num_space_dim);
1006
1007 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1008 auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1009 auto s_jac_inv = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1010
1011 cell_tools.setJacobianInv(s_jac_inv, s_jac);
1012
1013 PHX::Device::execution_space().fence();
1014
1015 if(cache){
1016 jac_inv = aux;
1017 jac_inv_evaluated_ = true;
1018 }
1019
1020 return aux;
1021
1022}
1023
1024template <typename Scalar>
1027getJacobianDeterminant(const bool cache,
1028 const bool force) const
1029{
1030
1031 if(jac_det_evaluated_ and not force)
1032 return jac_det;
1033
1034 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1035 MDFieldArrayFactory af(prefix_,true);
1036
1037 const int num_ip = int_rule->num_points;
1038
1039 auto jacobian = getJacobian(false,force);
1040 auto aux = af.template buildStaticArray<Scalar,Cell,IP>("jac_det",num_cells_, num_ip);
1041
1042 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1043 auto s_jac = Kokkos::subview(jacobian.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1044 auto s_jac_det = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL());
1045
1046 cell_tools.setJacobianDet(s_jac_det, s_jac);
1047
1048 PHX::Device::execution_space().fence();
1049
1050 if(cache){
1051 jac_det = aux;
1052 jac_det_evaluated_ = true;
1053 }
1054
1055 return aux;
1056
1057}
1058
1059template <typename Scalar>
1062getWeightedMeasure(const bool cache,
1063 const bool force) const
1064{
1065
1066 if(weighted_measure_evaluated_ and not force)
1067 return weighted_measure;
1068
1069 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1070 MDFieldArrayFactory af(prefix_,true);
1071
1072 const int num_space_dim = int_rule->topology->getDimension();
1073 const int num_ip = int_rule->num_points;
1074
1075 auto aux = af.template buildStaticArray<Scalar,Cell,IP>("weighted_measure",num_cells_, num_ip);
1076
1077 if(int_rule->cv_type != "none"){
1078
1079 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type == "side",
1080 "IntegrationValues2::getWeightedMeasure : Weighted measure not available for side control volume integrators. Use getWeightedNormals instead.");
1081
1082 // CV integration uses a single call to map from physical space to the weighted measure - I assume this is slower than what we do with non-cv integration methods
1083
1084 auto s_cub_points = af.template buildStaticArray<Scalar, Cell, IP, Dim>("cub_points",num_evaluate_cells_,num_ip,num_space_dim);
1085
1086 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1087
1088 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1089 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1090 auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL());
1091
1092 intrepid_cubature->getCubature(s_cub_points.get_view(),s_weighted_measure,s_node_coord);
1093
1094 } else if(int_rule->getType() == IntegrationDescriptor::SURFACE){
1095
1096 const auto & cell_topology = *int_rule->topology;
1097 const int cell_dim = cell_topology.getDimension();
1098 const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount();
1099 const int cubature_order = int_rule->order();
1100 const int num_points_on_side = num_ip / num_sides;
1101
1102 Intrepid2::DefaultCubatureFactory cubature_factory;
1103 auto jacobian = getJacobian(false,force);
1104
1105 for(int side=0; side<num_sides; ++side) {
1106
1107 const int point_offset=side*num_points_on_side;
1108
1109 // Get the cubature for the side
1110 Kokkos::DynRankView<double,PHX::Device> side_cub_weights;
1111 if(cell_dim==1){
1112 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>("side_cub_weights",num_points_on_side);
1113 auto tmp_side_cub_weights_host = Kokkos::create_mirror_view(side_cub_weights);
1114 tmp_side_cub_weights_host(0)=1.;
1115 Kokkos::deep_copy(side_cub_weights,tmp_side_cub_weights_host);
1116 } else {
1117
1118 // Get the face topology from the cell topology
1119 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(cell_dim-1,side));
1120
1121 auto ic = cubature_factory.create<PHX::Device::execution_space,double,double>(face_topology,cubature_order);
1122
1123 side_cub_weights = Kokkos::DynRankView<double,PHX::Device>("side_cub_weights",num_points_on_side);
1124 auto subcell_cub_points = Kokkos::DynRankView<double,PHX::Device>("subcell_cub_points",num_points_on_side,cell_dim-1);
1125
1126 // Get the reference face points
1127 ic->getCubature(subcell_cub_points, side_cub_weights);
1128 }
1129
1130 PHX::Device::execution_space().fence();
1131
1132 // Iterating over face points
1133 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_side});
1134
1135 // Calculate measures (quadrature weights in physical space) for this side
1136 auto side_weighted_measure = Kokkos::DynRankView<Scalar,PHX::Device>("side_weighted_measure",num_evaluate_cells_,num_points_on_side);
1137 if(cell_dim == 1){
1138 auto side_cub_weights_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),side_cub_weights);
1139 Kokkos::deep_copy(side_weighted_measure, side_cub_weights_host(0));
1140 } else {
1141
1142 // Copy from complete jacobian to side jacobian
1143 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_jac",num_evaluate_cells_,num_points_on_side,cell_dim,cell_dim);
1144
1145 Kokkos::parallel_for("Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (const int cell,const int point) {
1146 for(int dim=0;dim<cell_dim;++dim)
1147 for(int dim1=0;dim1<cell_dim;++dim1)
1148 side_jacobian(cell,point,dim,dim1) = jacobian(cell,point_offset+point,dim,dim1);
1149 });
1150
1151 auto scratch = af.template buildStaticArray<Scalar,Point>("scratch_for_compute_measure", num_evaluate_cells_*num_points_on_side*num_space_dim*num_space_dim);
1152
1153 if(cell_dim == 2){
1154 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1155 computeEdgeMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
1156 side,cell_topology,
1157 scratch.get_view());
1158 PHX::Device::execution_space().fence();
1159 } else if(cell_dim == 3){
1160 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1161 computeFaceMeasure(side_weighted_measure, side_jacobian, side_cub_weights,
1162 side,cell_topology,
1163 scratch.get_view());
1164 PHX::Device::execution_space().fence();
1165 }
1166 }
1167
1168
1169 // Copy to the main array
1170 Kokkos::parallel_for("copy surface weighted measure values",policy,KOKKOS_LAMBDA (const int cell,const int point) {
1171 aux(cell,point_offset + point) = side_weighted_measure(cell,point);
1172 });
1173 PHX::Device::execution_space().fence();
1174 }
1175
1176 } else {
1177
1178 auto cell_range = std::make_pair(0,num_evaluate_cells_);
1179 auto s_weighted_measure = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL());
1180 auto cubature_weights = getUniformCubatureWeightsRef(false,force,false);
1181
1182 if (!int_rule->isSide()) {
1183
1184 auto s_jac_det = Kokkos::subview(getJacobianDeterminant(false,force).get_view(),cell_range,Kokkos::ALL());
1185 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1186 computeCellMeasure(s_weighted_measure, s_jac_det, cubature_weights.get_view());
1187
1188 } else if(int_rule->spatial_dimension==3) {
1189
1190 auto s_jac = Kokkos::subview(getJacobian(false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1191 auto scratch = af.template buildStaticArray<Scalar,Point>("scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim);
1192 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1193 computeFaceMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(),
1194 int_rule->side, *int_rule->topology,
1195 scratch.get_view());
1196
1197 } else if(int_rule->spatial_dimension==2) {
1198
1199 auto s_jac = Kokkos::subview(getJacobian(false,force).get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL());
1200 auto scratch = af.template buildStaticArray<Scalar,Point>("scratch_for_compute_measure", num_evaluate_cells_*num_ip*num_space_dim*num_space_dim);
1201 Intrepid2::FunctionSpaceTools<PHX::Device::execution_space>::
1202 computeEdgeMeasure(s_weighted_measure, s_jac, cubature_weights.get_view(),
1203 int_rule->side,*int_rule->topology,
1204 scratch.get_view());
1205
1206 } else {
1207 TEUCHOS_ASSERT(false);
1208 }
1209
1210 }
1211
1212 PHX::Device::execution_space().fence();
1213
1214 // Apply permutations if necessary
1215 if(requires_permutation_)
1216 applyPermutation(aux, permutations_);
1217
1218 if(cache){
1219 weighted_measure = aux;
1220 weighted_measure_evaluated_ = true;
1221 }
1222
1223 return aux;
1224
1225}
1226
1227template <typename Scalar>
1230getWeightedNormals(const bool cache,
1231 const bool force) const
1232{
1233
1234 if(weighted_normals_evaluated_ and not force)
1235 return weighted_normals;
1236
1237 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1238 MDFieldArrayFactory af(prefix_,true);
1239
1240 const int num_space_dim = int_rule->topology->getDimension();
1241 const int num_ip = int_rule->num_points;
1242
1243 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>("weighted_normals",num_cells_,num_ip,num_space_dim);
1244
1245 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() == IntegrationDescriptor::SURFACE,
1246 "IntegrationValues2::getWeightedNormals : Cannot build reference weighted normals for surface integration scheme.");
1247
1248 auto points = af.template buildStaticArray<Scalar,Cell,IP,Dim>("cub_points",num_cells_,num_ip,num_space_dim);
1249
1250 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1251
1252 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1253 auto s_cub_points = Kokkos::subview(points.get_view(),cell_range, Kokkos::ALL(), Kokkos::ALL());
1254 auto s_weighted_normals = Kokkos::subview(aux.get_view(), cell_range, Kokkos::ALL(), Kokkos::ALL());
1255 auto s_node_coord = Kokkos::subview(node_coord, cell_range, Kokkos::ALL(), Kokkos::ALL());
1256
1257 intrepid_cubature->getCubature(s_cub_points,s_weighted_normals,s_node_coord);
1258
1259 PHX::Device::execution_space().fence();
1260
1261 // Apply permutations if necessary
1262 if(requires_permutation_)
1263 applyPermutation(aux, permutations_);
1264
1265 if(cache){
1266 weighted_normals = aux;
1267 weighted_normals_evaluated_ = true;
1268 }
1269
1270 return aux;
1271
1272}
1273
1274template <typename Scalar>
1277getSurfaceNormals(const bool cache,
1278 const bool force) const
1279{
1280
1281 if(surface_normals_evaluated_ and not force)
1282 return surface_normals;
1283
1284 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->isSide(),
1285 "IntegrationValues2::getSurfaceNormals : This call doesn't work with sides (only surfaces).");
1286
1287 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->cv_type != "none",
1288 "IntegrationValues2::getSurfaceNormals : This call does not support control volume integration schemes.");
1289
1290 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->getType() != IntegrationDescriptor::SURFACE,
1291 "IntegrationValues2::getSurfaceNormals : Can only build for surface integrators.");
1292
1293 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1294 MDFieldArrayFactory af(prefix_,true);
1295
1296 const shards::CellTopology & cell_topology = *(int_rule->topology);
1297 const int cell_dim = cell_topology.getDimension();
1298 const int subcell_dim = cell_topology.getDimension()-1;
1299 const int num_subcells = cell_topology.getSubcellCount(subcell_dim);
1300 const int num_space_dim = int_rule->topology->getDimension();
1301 const int num_ip = int_rule->num_points;
1302 const int num_points_on_face = num_ip / num_subcells;
1303
1304 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>("surface_normals",num_cells_,num_ip,num_space_dim);
1305
1306 // We only need the jacobian if we're not in 1D
1307 ConstArray_CellIPDimDim jacobian;
1308 if(cell_dim != 1)
1309 jacobian = getJacobian(false,force);
1310
1311 // We get to build up our surface normals one 'side' at a time
1312 for(int subcell_index=0; subcell_index<num_subcells; ++subcell_index) {
1313
1314 const int point_offset = subcell_index * num_points_on_face;;
1315
1316 // Calculate normals
1317 auto side_normals = Kokkos::DynRankView<Scalar,PHX::Device>("side_normals",num_evaluate_cells_,num_points_on_face,cell_dim);
1318 if(cell_dim == 1){
1319
1320 const int other_subcell_index = (subcell_index==0) ? 1 : 0;
1321 auto in_node_coordinates_k = getNodeCoordinates().get_view();
1322 const auto min = std::numeric_limits<typename Sacado::ScalarType<Scalar>::type>::min();
1323
1324 Kokkos::parallel_for("compute normals 1D",num_evaluate_cells_,KOKKOS_LAMBDA (const int cell) {
1325 Scalar norm = (in_node_coordinates_k(cell,subcell_index,0) - in_node_coordinates_k(cell,other_subcell_index,0));
1326 side_normals(cell,0,0) = norm / fabs(norm + min);
1327 });
1328
1329 } else {
1330
1331 // Iterating over side points
1332 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_face});
1333
1334 auto side_jacobian = Kokkos::DynRankView<Scalar,PHX::Device>("side_jac",num_evaluate_cells_,num_points_on_face,cell_dim,cell_dim);
1335 Kokkos::parallel_for("Copy jacobian to side jacobian",policy,KOKKOS_LAMBDA (const int cell,const int point) {
1336 for(int dim=0;dim<cell_dim;++dim)
1337 for(int dim1=0;dim1<cell_dim;++dim1)
1338 side_jacobian(cell,point,dim,dim1) = jacobian(cell,point_offset+point,dim,dim1);
1339 });
1340
1341 // Get the "physical side normals"
1342 cell_tools.getPhysicalSideNormals(side_normals,side_jacobian,subcell_index,cell_topology);
1343
1344 // Normalize each normal
1345 Kokkos::parallel_for("Normalize the normals",policy,KOKKOS_LAMBDA (const int cell,const int point) {
1346 Scalar n = 0.;
1347 for(int dim=0;dim<cell_dim;++dim){
1348 n += side_normals(cell,point,dim)*side_normals(cell,point,dim);
1349 }
1350 // If n is zero then this is - hopefully - a virtual cell
1351 if(n > 0.){
1352 n = Kokkos::sqrt(n);
1353 for(int dim=0;dim<cell_dim;++dim)
1354 side_normals(cell,point,dim) /= n;
1355 }
1356 });
1357 }
1358
1359 PHX::Device::execution_space().fence();
1360
1361 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_points_on_face});
1362 Kokkos::parallel_for("copy surface normals", policy,KOKKOS_LAMBDA (const int cell,const int point) {
1363 for(int dim=0;dim<cell_dim;++dim)
1364 aux(cell,point_offset + point,dim) = side_normals(cell,point,dim);
1365 });
1366 PHX::Device::execution_space().fence();
1367 }
1368
1369 // Need to correct the virtual cells
1370 {
1371 TEUCHOS_TEST_FOR_EXCEPT_MSG(side_connectivity_ == Teuchos::null,
1372 "IntegrationValues2::getSurfaceNormals : Surface normals require a SubcellConnectivity object pass in through the setupPermutations call");
1373 TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ < 0,
1374 "IntegrationValues2::getSurfaceNormals : Number of virtual cells (see setupPermutations or evaluateValues) must be greater or equal to zero");
1375
1376 // Virtual cell normals need to be reversed
1377 if(num_virtual_cells_ > 0)
1378 correctVirtualNormals(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_);
1379 }
1380
1381 if(cache){
1382 surface_normals = aux;
1383 surface_normals_evaluated_ = true;
1384 }
1385
1386 return aux;
1387
1388}
1389
1390template <typename Scalar>
1393getSurfaceRotationMatrices(const bool cache,
1394 const bool force) const
1395{
1396
1397 if(surface_rotation_matrices_evaluated_ and not force)
1398 return surface_rotation_matrices;
1399
1400 MDFieldArrayFactory af(prefix_,true);
1401
1402 const int num_ip = int_rule->num_points;
1403 const int cell_dim = int_rule->topology->getDimension();
1404
1405 // This call will handle all the error checking
1406 auto normals = getSurfaceNormals(false,force).get_static_view();
1407 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("surface_rotation_matrices",num_cells_, num_ip, 3, 3);
1408
1409 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1410 Kokkos::parallel_for("create surface rotation matrices",policy,KOKKOS_LAMBDA (const int cell,const int point) {
1411 Scalar normal[3];
1412 for(int i=0;i<3;i++)
1413 normal[i]=0.;
1414 for(int dim=0; dim<cell_dim; ++dim)
1415 normal[dim] = normals(cell,point,dim);
1416
1417 Scalar transverse[3];
1418 Scalar binormal[3];
1419 panzer::convertNormalToRotationMatrix(normal,transverse,binormal);
1420
1421 for(int dim=0; dim<3; ++dim){
1422 aux(cell,point,0,dim) = normal[dim];
1423 aux(cell,point,1,dim) = transverse[dim];
1424 aux(cell,point,2,dim) = binormal[dim];
1425 }
1426 });
1427 PHX::Device::execution_space().fence();
1428
1429 // Need to correct the virtual cells
1430 {
1431 TEUCHOS_TEST_FOR_EXCEPT_MSG(side_connectivity_ == Teuchos::null,
1432 "IntegrationValues2::getSurfaceRotationMatrices : Surface normals require a SubcellConnectivity object pass in through the setupPermutations call");
1433 TEUCHOS_TEST_FOR_EXCEPT_MSG(num_virtual_cells_ < 0,
1434 "IntegrationValues2::getSurfaceRotationMatrices : Number of virtual cells (see setupPermutations or evaluateValues) must be greater or equal to zero");
1435
1436 // Virtual cell normals need to be reversed
1437 if(num_virtual_cells_ > 0)
1438 correctVirtualRotationMatrices(aux, num_evaluate_cells_ - num_virtual_cells_, num_virtual_cells_, *int_rule->topology, *side_connectivity_);
1439 }
1440
1441 if(cache){
1442 surface_rotation_matrices = aux;
1443 surface_rotation_matrices_evaluated_ = true;
1444 }
1445
1446 return aux;
1447}
1448
1449template <typename Scalar>
1452getCovarientMatrix(const bool cache,
1453 const bool force) const
1454{
1455
1456 if(covarient_evaluated_ and not force)
1457 return covarient;
1458
1459 MDFieldArrayFactory af(prefix_,true);
1460
1461 const int num_space_dim = int_rule->topology->getDimension();
1462 const int num_ip = int_rule->num_points;
1463
1464 auto jacobian = getJacobian(false,force).get_static_view();
1465 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("covarient",num_cells_, num_ip, num_space_dim,num_space_dim);
1466
1467 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1468 Kokkos::parallel_for("evalaute covarient metric tensor",policy,KOKKOS_LAMBDA (const int cell,const int ip) {
1469 // g^{ij} = \frac{\parital x_i}{\partial \chi_\alpha}\frac{\parital x_j}{\partial \chi_\alpha}
1470 for (int i = 0; i < num_space_dim; ++i) {
1471 for (int j = 0; j < num_space_dim; ++j) {
1472 Scalar value(0.0);
1473 for (int alpha = 0; alpha < num_space_dim; ++alpha)
1474 value += jacobian(cell,ip,i,alpha) * jacobian(cell,ip,j,alpha);
1475 aux(cell,ip,i,j) = value;
1476 }
1477 }
1478 });
1479 PHX::Device::execution_space().fence();
1480
1481 if(cache){
1482 covarient = aux;
1483 covarient_evaluated_ = true;
1484 }
1485
1486 return aux;
1487
1488}
1489
1490template <typename Scalar>
1493getContravarientMatrix(const bool cache,
1494 const bool force) const
1495{
1496
1497 if(contravarient_evaluated_ and not force)
1498 return contravarient;
1499
1500 MDFieldArrayFactory af(prefix_,true);
1501
1502 const int num_space_dim = int_rule->topology->getDimension();
1503 const int num_ip = int_rule->num_points;
1504
1505 auto cov = getCovarientMatrix(false,force);
1506 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim,Dim>("contravarient",num_cells_, num_ip, num_space_dim,num_space_dim);
1507
1508 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1509 auto s_contravarient = Kokkos::subview(aux.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
1510 auto s_covarient = Kokkos::subview(cov.get_view(), cell_range,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL);
1511
1512 Intrepid2::RealSpaceTools<PHX::Device::execution_space>::inverse(s_contravarient, s_covarient);
1513 PHX::Device::execution_space().fence();
1514
1515 if(cache){
1516 contravarient = aux;
1517 contravarient_evaluated_ = true;
1518 }
1519
1520 return aux;
1521
1522}
1523
1524template <typename Scalar>
1527getNormContravarientMatrix(const bool cache,
1528 const bool force) const
1529{
1530
1531 if(norm_contravarient_evaluated_ and not force)
1532 return norm_contravarient;
1533
1534 MDFieldArrayFactory af(prefix_,true);
1535
1536 const int num_space_dim = int_rule->topology->getDimension();
1537 const int num_ip = int_rule->num_points;
1538
1539 auto con = getContravarientMatrix(false,force).get_static_view();
1540 auto aux = af.template buildStaticArray<Scalar,Cell,IP>("norm_contravarient",num_cells_, num_ip);
1541
1542 // norm of g_ij
1543 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<2>> policy({0,0},{num_evaluate_cells_,num_ip});
1544 Kokkos::parallel_for("evaluate norm_contravarient",policy,KOKKOS_LAMBDA (const int cell,const int ip) {
1545 aux(cell,ip) = 0.0;
1546 for (int i = 0; i < num_space_dim; ++i) {
1547 for (int j = 0; j < num_space_dim; ++j) {
1548 aux(cell,ip) += con(cell,ip,i,j) * con(cell,ip,i,j);
1549 }
1550 }
1551 aux(cell,ip) = Kokkos::sqrt(aux(cell,ip));
1552 });
1553 PHX::Device::execution_space().fence();
1554
1555 if(cache){
1556 norm_contravarient = aux;
1557 norm_contravarient_evaluated_ = true;
1558 }
1559
1560 return aux;
1561
1562}
1563
1564template <typename Scalar>
1567getCubaturePoints(const bool cache,
1568 const bool force) const
1569{
1570
1571 if(ip_coordinates_evaluated_ and not force)
1572 return ip_coordinates;
1573
1574 MDFieldArrayFactory af(prefix_,true);
1575
1576 const int num_space_dim = int_rule->topology->getDimension();
1577 const int num_ip = int_rule->num_points;
1578
1579 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ip_coordinates",num_cells_, num_ip, num_space_dim);
1580
1582 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1583
1584 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1585
1586 if(is_cv){
1587
1588 // CV integration uses a single call to map from physical space to the weighted measure - I assume this is slower than what we do with non-cv integration methods
1589 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1590 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1591 auto s_cub_points = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1592
1593 // TODO: We need to pull this apart for control volumes. Right now we calculate both weighted measures/norms and cubature points at the same time
1594 if(int_rule->cv_type == "side"){
1595 auto scratch = af.template buildStaticArray<Scalar,Cell,IP,Dim>("scratch",num_evaluate_cells_,num_ip,num_space_dim);
1596 intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord);
1597 } else {
1598 // I think boundary is included as a weighted measure because it has a side embedded in intrepid_cubature
1599 TEUCHOS_ASSERT((int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_BOUNDARY));
1600 auto scratch = af.template buildStaticArray<Scalar,Cell,IP>("scratch",num_evaluate_cells_,num_ip);
1601 intrepid_cubature->getCubature(s_cub_points,scratch.get_view(),s_node_coord);
1602 }
1603
1604 } else {
1605
1606 // 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
1607 auto const_ref_coord = getCubaturePointsRef(false,force);
1608 auto ref_coord = PHX::getNonConstDynRankViewFromConstMDField(const_ref_coord);
1609
1610 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1611 auto s_ref_coord = Kokkos::subview(ref_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1612 auto s_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1613 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1614
1615 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1616 cell_tools.mapToPhysicalFrame(s_coord, s_ref_coord, s_node_coord, *(int_rule->topology));
1617
1618 }
1619
1620 PHX::Device::execution_space().fence();
1621
1622 if(cache){
1623 ip_coordinates = aux;
1624 ip_coordinates_evaluated_ = true;
1625 }
1626
1627 return aux;
1628
1629}
1630
1631
1632template <typename Scalar>
1635getCubaturePointsRef(const bool cache,
1636 const bool force) const
1637{
1638
1639 if(ref_ip_coordinates_evaluated_)
1640 return ref_ip_coordinates;
1641
1643 const bool is_surface = int_rule->getType() == ID::SURFACE;
1644 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1645
1646 const int num_space_dim = int_rule->topology->getDimension();
1647 const int num_ip = int_rule->num_points;
1648
1649 MDFieldArrayFactory af(prefix_,true);
1650
1651 Intrepid2::CellTools<PHX::Device::execution_space> cell_tools;
1652
1653 auto aux = af.template buildStaticArray<Scalar,Cell,IP,Dim>("ref_ip_coordinates",num_cells_, num_ip, num_space_dim);
1654
1655 if(is_cv){
1656
1657 // Control volume reference points are actually generated from the physical points (i.e. reverse to everything else)
1658
1659 auto node_coord = PHX::getNonConstDynRankViewFromConstMDField(getNodeCoordinates());
1660
1661 // 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
1662 auto const_coord = getCubaturePoints(false,force);
1663 auto coord = PHX::getNonConstDynRankViewFromConstMDField(const_coord);
1664
1665 const auto cell_range = std::make_pair(0,num_evaluate_cells_);
1666 auto s_ref_coord = Kokkos::subview(aux.get_view(),cell_range,Kokkos::ALL(),Kokkos::ALL());
1667 auto s_coord = Kokkos::subview(coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1668 auto s_node_coord = Kokkos::subview(node_coord, cell_range,Kokkos::ALL(),Kokkos::ALL());
1669
1670 cell_tools.mapToReferenceFrame(s_ref_coord, s_coord, s_node_coord, *(int_rule->topology));
1671
1672 } else if(is_surface){
1673
1674 const auto & cell_topology = *int_rule->topology;
1675 const int cell_dim = cell_topology.getDimension();
1676 const int num_sides = (cell_dim==1) ? 2 : cell_topology.getSideCount();
1677 const int subcell_dim = cell_dim-1;
1678 const int num_points_on_face = num_ip / num_sides;
1679 const int order = int_rule->getOrder();
1680
1681 // Scratch space for storing the points for each side of the cell
1682 auto side_cub_points = af.template buildStaticArray<Scalar,IP,Dim>("side_cub_points",num_points_on_face,cell_dim);
1683
1684 Intrepid2::DefaultCubatureFactory cubature_factory;
1685
1686 // We get to build up our cubature one side at a time
1687 for(int side=0; side<num_sides; ++side) {
1688
1689 const int point_offset = side*num_points_on_face;
1690
1691 // Get the cubature for the side
1692 if(cell_dim==1){
1693 // In 1D the surface point is either on the left side of the cell, or the right side
1694 auto side_cub_points_host = Kokkos::create_mirror_view(side_cub_points.get_view());
1695 side_cub_points_host(0,0) = (side==0)? -1. : 1.;
1696 Kokkos::deep_copy(side_cub_points.get_view(),side_cub_points_host);
1697 } else {
1698
1699 // Get the face topology from the cell topology
1700 const shards::CellTopology face_topology(cell_topology.getCellTopologyData(subcell_dim,side));
1701
1702 // Create a cubature for the face of the cell
1703 auto ic = cubature_factory.create<PHX::Device::execution_space,double,double>(face_topology,order);
1704 auto tmp_side_cub_weights = Kokkos::DynRankView<double,PHX::Device>("tmp_side_cub_weights",num_points_on_face);
1705 auto tmp_side_cub_points = Kokkos::DynRankView<double,PHX::Device>("tmp_side_cub_points",num_points_on_face,subcell_dim);
1706
1707 // Get the reference face points
1708 ic->getCubature(tmp_side_cub_points, tmp_side_cub_weights);
1709
1710 // Convert from reference face points to reference cell points
1711 cell_tools.mapToReferenceSubcell(side_cub_points.get_view(), tmp_side_cub_points, subcell_dim, side, cell_topology);
1712 }
1713
1714 PHX::Device::execution_space().fence();
1715
1716 // Copy from the side allocation to the surface allocation
1717 Kokkos::MDRangePolicy<PHX::Device::execution_space,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_points_on_face, num_space_dim});
1718 Kokkos::parallel_for("copy values",policy,KOKKOS_LAMBDA (const int cell,const int point, const int dim) {
1719 aux(cell,point_offset + point,dim) = side_cub_points(point,dim);
1720 });
1721 PHX::Device::execution_space().fence();
1722 }
1723
1724 } else {
1725
1726 // Haven't set this up yet
1727 TEUCHOS_TEST_FOR_EXCEPT_MSG(int_rule->isSide() && num_space_dim==1,
1728 "ERROR: 0-D quadrature rule infrastructure does not exist!!! Will not be able to do "
1729 << "non-natural integration rules.");
1730
1731 auto cub_points = getUniformCubaturePointsRef(false,force,false);
1732
1733 Kokkos::MDRangePolicy<PHX::Device,Kokkos::Rank<3>> policy({0,0,0},{num_evaluate_cells_,num_ip,num_space_dim});
1734 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const int & cell, const int & ip, const int & dim){
1735 aux(cell,ip,dim) = cub_points(ip,dim);
1736 });
1737 }
1738
1739 PHX::Device::execution_space().fence();
1740
1741 if(requires_permutation_)
1742 applyPermutation(aux, permutations_);
1743
1744 if(cache){
1745 ref_ip_coordinates = aux;
1746 ref_ip_coordinates_evaluated_ = true;
1747 }
1748
1749 return aux;
1750
1751}
1752
1753template <typename Scalar>
1754void
1757{
1758
1760 const bool is_surface = int_rule->getType() == ID::SURFACE;
1761 const bool is_cv = (int_rule->getType() == ID::CV_VOLUME) or (int_rule->getType() == ID::CV_SIDE) or (int_rule->getType() == ID::CV_BOUNDARY);
1762 const bool is_side = int_rule->isSide();
1763
1764 resetArrays();
1765
1766 // Base cubature stuff
1767 if(is_cv){
1768 getCubaturePoints(true,true);
1769 getCubaturePointsRef(true,true);
1770 } else {
1771 if(not is_surface){
1772 getUniformCubaturePointsRef(true,true);
1773 getUniformCubatureWeightsRef(true,true);
1774 if(is_side)
1775 getUniformSideCubaturePointsRef(true,true);
1776 }
1777 getCubaturePointsRef(true,true);
1778 getCubaturePoints(true,true);
1779 }
1780
1781 // Measure stuff
1782 getJacobian(true,true);
1783 getJacobianDeterminant(true,true);
1784 getJacobianInverse(true,true);
1785 if(int_rule->cv_type == "side")
1786 getWeightedNormals(true,true);
1787 else
1788 getWeightedMeasure(true,true);
1789
1790 // Surface stuff
1791 if(is_surface){
1792 getSurfaceNormals(true,true);
1793 getSurfaceRotationMatrices(true,true);
1794 }
1795
1796 // Stabilization stuff
1797 if(not (is_surface or is_cv)){
1798 getContravarientMatrix(true,true);
1799 getCovarientMatrix(true,true);
1800 getNormContravarientMatrix(true,true);
1801 }
1802
1803}
1804
1805#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR) \
1806 template class IntegrationValues2<SCALAR>;
1807
1809
1810// Disabled FAD support due to long build times on cuda (in debug mode
1811// it takes multiple hours on some platforms). If we need
1812// sensitivities wrt coordinates, we can reenable.
1813
1814// INTEGRATION_VALUES2_INSTANTIATION(panzer::Traits::FadType)
1815
1816}
#define PANZER_CROSS(a, b, c)
#define PANZER_DOT(a, b)
#define INTEGRATION_VALUES2_INSTANTIATION(SCALAR)
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
const int & getOrder() const
Get order of integrator.
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
const int & getType() const
Get type of integrator.
ConstArray_IPDim getUniformCubaturePointsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature points.
void setupPermutations(const Teuchos::RCP< const SubcellConnectivity > &face_connectivity, const int num_virtual_cells)
Initialize the permutation arrays given a face connectivity.
ConstArray_CellIP getWeightedMeasure(const bool cache=true, const bool force=false) const
Get the weighted measure (integration weights)
ConstArray_CellIPDimDim getJacobian(const bool cache=true, const bool force=false) const
Get the Jacobian matrix evaluated at the cubature points.
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int num_cells=-1, const Teuchos::RCP< const SubcellConnectivity > &face_connectivity=Teuchos::null, const int num_virtual_cells=-1)
Evaluate basis values.
PHX::MDField< const Scalar, Cell, IP > ConstArray_CellIP
ConstArray_CellIPDim getCubaturePointsRef(const bool cache=true, const bool force=false) const
Get the cubature points in the reference space.
ConstArray_CellIPDimDim getSurfaceRotationMatrices(const bool cache=true, const bool force=false) const
Get the surface rotation matrices.
PHX::MDField< const Scalar, IP, Dim > ConstArray_IPDim
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
PHX::MDField< const Scalar, Cell, IP, Dim, Dim > ConstArray_CellIPDimDim
PHX::MDField< const Scalar, IP > ConstArray_IP
IntegrationValues2(const std::string &pre="", const bool allocArrays=false)
Base constructor.
ConstArray_CellIPDim getWeightedNormals(const bool cache=true, const bool force=false) const
Get the weighted normals.
ConstArray_CellIP getNormContravarientMatrix(const bool cache=true, const bool force=false) const
Get the contravarient matrix.
ConstArray_CellIPDimDim getContravarientMatrix(const bool cache=true, const bool force=false) const
Get the contravarient matrix.
ConstArray_CellIPDim getCubaturePoints(const bool cache=true, const bool force=false) const
Get the cubature points in physical space.
PHX::MDField< const Scalar, Cell, IP, Dim > ConstArray_CellIPDim
void setup(const Teuchos::RCP< const panzer::IntegrationRule > &ir, const PHX::MDField< Scalar, Cell, NODE, Dim > &node_coordinates, const int num_cells=-1)
Main setup call for the lazy evaluation interface.
ConstArray_CellIPDimDim getCovarientMatrix(const bool cache=true, const bool force=false) const
Get the covarient matrix.
ConstArray_CellBASISDim getNodeCoordinates() const
Get the node coordinates describing the geometry of the mesh.
ConstArray_IPDim getUniformSideCubaturePointsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature points for a side.
ConstArray_CellIPDim getSurfaceNormals(const bool cache=true, const bool force=false) const
Get the surface normals.
PHX::MDField< const Scalar, Cell, BASIS, Dim > ConstArray_CellBASISDim
ConstArray_IP getUniformCubatureWeightsRef(const bool cache=true, const bool force=false, const bool apply_permutation=true) const
Get the uniform cubature weights.
ConstArray_CellIPDimDim getJacobianInverse(const bool cache=true, const bool force=false) const
Get the inverse of the Jacobian matrix evaluated at the cubature points.
ConstArray_CellIP getJacobianDeterminant(const bool cache=true, const bool force=false) const
Get the determinant of the Jacobian matrix evaluated at the cubature points.
Teuchos::RCP< shards::CellTopology > side_topology
Teuchos::RCP< const shards::CellTopology > topology
KOKKOS_INLINE_FUNCTION void convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])