Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
view/fenl_functors.hpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Kokkos: Manycore Performance-Portable Multidimensional Arrays
6 // Copyright (2012) 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 H. Carter Edwards (hcedwar@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 */
43 
44 #ifndef KOKKOS_EXAMPLE_FENLFUNCTORS_HPP
45 #define KOKKOS_EXAMPLE_FENLFUNCTORS_HPP
46 
47 #include <stdio.h>
48 
49 #include <iostream>
50 #include <fstream>
51 #include <iomanip>
52 #include <cstdlib>
53 #include <cmath>
54 #include <limits>
55 
56 #include <Kokkos_Core.hpp>
57 #include <Kokkos_Pair.hpp>
58 #include <Kokkos_UnorderedMap.hpp>
59 #include <Kokkos_StaticCrsGraph.hpp>
60 
61 #include <impl/Kokkos_Timer.hpp>
62 
63 #include <BoxElemFixture.hpp>
64 #include <HexElement.hpp>
65 
66 #include "Sacado.hpp"
67 
68 //----------------------------------------------------------------------------
69 //----------------------------------------------------------------------------
70 
71 namespace Kokkos {
72 namespace Example {
73 namespace FENL {
74 
75 template< typename ValueType , class Space >
76 struct CrsMatrix {
77 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE // Don't remove this until Kokkos has removed the deprecated code path probably around September 2018
78  typedef Kokkos::StaticCrsGraph< unsigned , Space , void , unsigned > StaticCrsGraphType ;
79 #else
80  typedef Kokkos::StaticCrsGraph< unsigned , Space , void , void , unsigned > StaticCrsGraphType ;
81 #endif
82  typedef View< ValueType * , Space > coeff_type ;
83 
86 
87  CrsMatrix() : graph(), coeff() {}
88 
89  CrsMatrix( const StaticCrsGraphType & arg_graph )
90  : graph( arg_graph )
91  , coeff( "crs_matrix_coeff" , arg_graph.entries.extent(0) )
92  {}
93 };
94 
95 template< class ElemNodeIdView , class CrsGraphType , unsigned ElemNode >
96 class NodeNodeGraph {
97 public:
98 
99  typedef typename ElemNodeIdView::execution_space execution_space ;
100  typedef pair<unsigned,unsigned> key_type ;
101 
102  typedef Kokkos::UnorderedMap< key_type, void , execution_space > SetType ;
103  typedef typename CrsGraphType::row_map_type::non_const_type RowMapType ;
104  typedef Kokkos::View< unsigned , execution_space > UnsignedValue ;
105 
106  // Static dimensions of 0 generate compiler warnings or errors.
107  typedef Kokkos::View< unsigned*[ElemNode][ElemNode] , execution_space >
109 
110  struct TagFillNodeSet {};
111  struct TagScanNodeCount {};
112  struct TagFillGraphEntries {};
113  struct TagSortGraphEntries {};
114  struct TagFillElementGraph {};
115 
116 private:
117 
123 
124  const unsigned node_count ;
125  const ElemNodeIdView elem_node_id ;
130  PhaseType phase ;
131 
132 public:
133 
134  CrsGraphType graph ;
136 
137  struct Times
138  {
139  double ratio;
140  double fill_node_set;
141  double scan_node_count;
142  double fill_graph_entries;
143  double sort_graph_entries;
144  double fill_element_graph;
145  };
146 
147  NodeNodeGraph( const ElemNodeIdView & arg_elem_node_id ,
148  const unsigned arg_node_count,
149  Times & results
150  )
151  : node_count(arg_node_count)
152  , elem_node_id( arg_elem_node_id )
153  , row_total( "row_total" )
154  , row_count(Kokkos::ViewAllocateWithoutInitializing("row_count") , node_count ) // will deep_copy to 0 inside loop
155  , row_map( "graph_row_map" , node_count + 1 )
156  , node_node_set()
157  , phase( FILL_NODE_SET )
158  , graph()
159  , elem_graph()
160  {
161  //--------------------------------
162  // Guess at capacity required for the map:
163 
164  Kokkos::Impl::Timer wall_clock ;
165 
166  wall_clock.reset();
167  phase = FILL_NODE_SET ;
168 
169  // upper bound on the capacity
170  size_t set_capacity = (28ull * node_count) / 2;
171  unsigned failed_insert_count = 0 ;
172 
173  do {
174  // Zero the row count to restart the fill
175  Kokkos::deep_copy( row_count , 0u );
176 
177  node_node_set = SetType( ( set_capacity += failed_insert_count ) );
178 
179  // May be larger that requested:
180  set_capacity = node_node_set.capacity();
181 
182  Kokkos::parallel_reduce( Kokkos::RangePolicy<execution_space,TagFillNodeSet>(0,elem_node_id.extent(0))
183  , *this
184  , failed_insert_count );
185 
186  } while ( failed_insert_count );
187 
188  execution_space().fence();
189  results.ratio = (double)node_node_set.size() / (double)node_node_set.capacity();
190  results.fill_node_set = wall_clock.seconds();
191  //--------------------------------
192 
193  wall_clock.reset();
195 
196  // Exclusive scan of row_count into row_map
197  // including the final total in the 'node_count + 1' position.
198  // Zero the 'row_count' values.
199  Kokkos::parallel_scan( node_count , *this );
200 
201  // Zero the row count for the fill:
202  Kokkos::deep_copy( row_count , 0u );
203 
204  unsigned graph_entry_count = 0 ;
205 
206  Kokkos::deep_copy( graph_entry_count , row_total );
207 
208  // Assign graph's row_map and allocate graph's entries
209  graph.row_map = row_map ;
210  graph.entries = typename CrsGraphType::entries_type( "graph_entries" , graph_entry_count );
211 
212  //--------------------------------
213  // Fill graph's entries from the (node,node) set.
214 
215  execution_space().fence();
216  results.scan_node_count = wall_clock.seconds();
217 
218  wall_clock.reset();
220  Kokkos::parallel_for( node_node_set.capacity() , *this );
221 
222  execution_space().fence();
223  results.fill_graph_entries = wall_clock.seconds();
224 
225  //--------------------------------
226  // Done with the temporary sets and arrays
227  wall_clock.reset();
229 
231  row_count = RowMapType();
232  row_map = RowMapType();
233  node_node_set.clear();
234 
235  //--------------------------------
236 
237  Kokkos::parallel_for( node_count , *this );
238 
239  execution_space().fence();
240  results.sort_graph_entries = wall_clock.seconds();
241 
242  //--------------------------------
243  // Element-to-graph mapping:
244  wall_clock.reset();
246  elem_graph = ElemGraphType("elem_graph", elem_node_id.extent(0) );
247  Kokkos::parallel_for( elem_node_id.extent(0) , *this );
248 
249  execution_space().fence();
250  results.fill_element_graph = wall_clock.seconds();
251  }
252 
253  //------------------------------------
254  // parallel_for: create map and count row length
255 
257  void operator()( const TagFillNodeSet & , unsigned ielem , unsigned & count ) const
258  {
259  // Loop over element's (row_local_node,col_local_node) pairs:
260  for ( unsigned row_local_node = 0 ; row_local_node < elem_node_id.extent(1) ; ++row_local_node ) {
261 
262  const unsigned row_node = elem_node_id( ielem , row_local_node );
263 
264  for ( unsigned col_local_node = row_local_node ; col_local_node < elem_node_id.extent(1) ; ++col_local_node ) {
265 
266  const unsigned col_node = elem_node_id( ielem , col_local_node );
267 
268  // If either node is locally owned then insert the pair into the unordered map:
269 
270  if ( row_node < row_count.extent(0) || col_node < row_count.extent(0) ) {
271 
272  const key_type key = (row_node < col_node) ? make_pair( row_node, col_node ) : make_pair( col_node, row_node ) ;
273 
274  const typename SetType::insert_result result = node_node_set.insert( key );
275 
276  // A successfull insert: the first time this pair was added
277  if ( result.success() ) {
278 
279  // If row node is owned then increment count
280  if ( row_node < row_count.extent(0) ) { atomic_increment( & row_count( row_node ) ); }
281 
282  // If column node is owned and not equal to row node then increment count
283  if ( col_node < row_count.extent(0) && col_node != row_node ) { atomic_increment( & row_count( col_node ) ); }
284  }
285  else if ( result.failed() ) {
286  ++count ;
287  }
288  }
289  }
290  }
291  }
292 
294  void fill_graph_entries( const unsigned iset ) const
295  {
296  if ( node_node_set.valid_at(iset) ) {
297  // Add each entry to the graph entries.
298 
299  const key_type key = node_node_set.key_at(iset) ;
300  const unsigned row_node = key.first ;
301  const unsigned col_node = key.second ;
302 
303  if ( row_node < row_count.extent(0) ) {
304  typedef typename std::remove_reference< decltype( row_count(0) ) >::type atomic_incr_type;
305  const unsigned offset = graph.row_map( row_node ) + atomic_fetch_add( & row_count( row_node ) , atomic_incr_type(1) );
306  graph.entries( offset ) = col_node ;
307  }
308 
309  if ( col_node < row_count.extent(0) && col_node != row_node ) {
310  typedef typename std::remove_reference< decltype( row_count(0) ) >::type atomic_incr_type;
311  const unsigned offset = graph.row_map( col_node ) + atomic_fetch_add( & row_count( col_node ) , atomic_incr_type(1) );
312  graph.entries( offset ) = row_node ;
313  }
314  }
315  }
316 
318  void sort_graph_entries( const unsigned irow ) const
319  {
320  const unsigned row_beg = graph.row_map( irow );
321  const unsigned row_end = graph.row_map( irow + 1 );
322  for ( unsigned i = row_beg + 1 ; i < row_end ; ++i ) {
323  const unsigned col = graph.entries(i);
324  unsigned j = i ;
325  for ( ; row_beg < j && col < graph.entries(j-1) ; --j ) {
326  graph.entries(j) = graph.entries(j-1);
327  }
328  graph.entries(j) = col ;
329  }
330  }
331 
333  void fill_elem_graph_map( const unsigned ielem ) const
334  {
335  for ( unsigned row_local_node = 0 ; row_local_node < elem_node_id.extent(1) ; ++row_local_node ) {
336 
337  const unsigned row_node = elem_node_id( ielem , row_local_node );
338 
339  for ( unsigned col_local_node = 0 ; col_local_node < elem_node_id.extent(1) ; ++col_local_node ) {
340 
341  const unsigned col_node = elem_node_id( ielem , col_local_node );
342 
343  unsigned entry = ~0u ;
344 
345  if ( row_node + 1 < graph.row_map.extent(0) ) {
346 
347  const unsigned entry_end = graph.row_map( row_node + 1 );
348 
349  entry = graph.row_map( row_node );
350 
351  for ( ; entry < entry_end && graph.entries(entry) != col_node ; ++entry );
352 
353  if ( entry == entry_end ) entry = ~0u ;
354  }
355 
356  elem_graph( ielem , row_local_node , col_local_node ) = entry ;
357  }
358  }
359  }
360 
362  void operator()( const unsigned iwork ) const
363  {
364 /*
365  if ( phase == FILL_NODE_SET ) {
366  operator()( TagFillNodeSet() , iwork );
367  }
368  else */
369  if ( phase == FILL_GRAPH_ENTRIES ) {
370  fill_graph_entries( iwork );
371  }
372  else if ( phase == SORT_GRAPH_ENTRIES ) {
373  sort_graph_entries( iwork );
374  }
375  else if ( phase == FILL_ELEMENT_GRAPH ) {
376  fill_elem_graph_map( iwork );
377  }
378  }
379 
380  //------------------------------------
381  // parallel_scan: row offsets
382 
383  typedef unsigned value_type ;
384 
386  void operator()( const unsigned irow , unsigned & update , const bool final ) const
387  {
388  // exclusive scan
389  if ( final ) { row_map( irow ) = update ; }
390 
391  update += row_count( irow );
392 
393  if ( final ) {
394  if ( irow + 1 == row_count.extent(0) ) {
395  row_map( irow + 1 ) = update ;
396  row_total() = update ;
397  }
398  }
399  }
400 
401  // For the reduce phase:
403  void init( const TagFillNodeSet & , unsigned & update ) const { update = 0 ; }
404 
406  void join( const TagFillNodeSet &
407  , volatile unsigned & update
408  , volatile const unsigned & input ) const { update += input ; }
409 
410  // For the scan phase::
412  void init( unsigned & update ) const { update = 0 ; }
413 
415  void join( volatile unsigned & update
416  , volatile const unsigned & input ) const { update += input ; }
417 
418  //------------------------------------
419 };
420 
421 } /* namespace FENL */
422 } /* namespace Example */
423 } /* namespace Kokkos */
424 
425 //----------------------------------------------------------------------------
426 //----------------------------------------------------------------------------
427 
428 namespace Kokkos {
429 namespace Example {
430 namespace FENL {
431 
432 template< class ExecutionSpace , BoxElemPart::ElemOrder Order ,
433  class CoordinateMap , typename ScalarType >
434 class ElementComputationBase
435 {
436 public:
437 
440 
441  //------------------------------------
442 
443  typedef ExecutionSpace execution_space ;
444  typedef ScalarType scalar_type ;
445 
448  typedef Kokkos::View< scalar_type* , Kokkos::LayoutLeft, execution_space > vector_type ;
449 
450  //------------------------------------
451 
452  static const unsigned SpatialDim = element_data_type::spatial_dimension ;
453  static const unsigned TensorDim = SpatialDim * SpatialDim ;
454  static const unsigned ElemNodeCount = element_data_type::element_node_count ;
455  static const unsigned FunctionCount = element_data_type::function_count ;
457 
458  //------------------------------------
459 
462  typedef Kokkos::View< scalar_type*[FunctionCount][FunctionCount] , execution_space > elem_matrices_type ;
463  typedef Kokkos::View< scalar_type*[FunctionCount] , execution_space > elem_vectors_type ;
464 
466 
467  //------------------------------------
468 
469 
470  //------------------------------------
471  // Computational data:
472 
479  const vector_type solution ;
480  const vector_type residual ;
482 
484  : elem_data()
486  , node_coords( rhs.node_coords )
487  , elem_graph( rhs.elem_graph )
490  , solution( rhs.solution )
491  , residual( rhs.residual )
492  , jacobian( rhs.jacobian )
493  {}
494 
495  ElementComputationBase( const mesh_type & arg_mesh ,
496  const vector_type & arg_solution ,
497  const elem_graph_type & arg_elem_graph ,
498  const sparse_matrix_type & arg_jacobian ,
499  const vector_type & arg_residual )
500  : elem_data()
501  , elem_node_ids( arg_mesh.elem_node() )
502  , node_coords( arg_mesh.node_coord() )
503  , elem_graph( arg_elem_graph )
504  , elem_jacobians()
505  , elem_residuals()
506  , solution( arg_solution )
507  , residual( arg_residual )
508  , jacobian( arg_jacobian )
509  {}
510 
511  //------------------------------------
512 
515  const double grad[][ FunctionCount ] , // Gradient of bases master element
516  const double x[] ,
517  const double y[] ,
518  const double z[] ,
519  double dpsidx[] ,
520  double dpsidy[] ,
521  double dpsidz[] ) const
522  {
523  enum { j11 = 0 , j12 = 1 , j13 = 2 ,
524  j21 = 3 , j22 = 4 , j23 = 5 ,
525  j31 = 6 , j32 = 7 , j33 = 8 };
526 
527  // Jacobian accumulation:
528 
529  double J[ TensorDim ] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
530 
531  for( unsigned i = 0; i < FunctionCount ; ++i ) {
532  const double x1 = x[i] ;
533  const double x2 = y[i] ;
534  const double x3 = z[i] ;
535 
536  const double g1 = grad[0][i] ;
537  const double g2 = grad[1][i] ;
538  const double g3 = grad[2][i] ;
539 
540  J[j11] += g1 * x1 ;
541  J[j12] += g1 * x2 ;
542  J[j13] += g1 * x3 ;
543 
544  J[j21] += g2 * x1 ;
545  J[j22] += g2 * x2 ;
546  J[j23] += g2 * x3 ;
547 
548  J[j31] += g3 * x1 ;
549  J[j32] += g3 * x2 ;
550  J[j33] += g3 * x3 ;
551  }
552 
553  // Inverse jacobian:
554 
555  double invJ[ TensorDim ] = {
556  static_cast<double>( J[j22] * J[j33] - J[j23] * J[j32] ) ,
557  static_cast<double>( J[j13] * J[j32] - J[j12] * J[j33] ) ,
558  static_cast<double>( J[j12] * J[j23] - J[j13] * J[j22] ) ,
559 
560  static_cast<double>( J[j23] * J[j31] - J[j21] * J[j33] ) ,
561  static_cast<double>( J[j11] * J[j33] - J[j13] * J[j31] ) ,
562  static_cast<double>( J[j13] * J[j21] - J[j11] * J[j23] ) ,
563 
564  static_cast<double>( J[j21] * J[j32] - J[j22] * J[j31] ) ,
565  static_cast<double>( J[j12] * J[j31] - J[j11] * J[j32] ) ,
566  static_cast<double>( J[j11] * J[j22] - J[j12] * J[j21] ) };
567 
568  const double detJ = J[j11] * invJ[j11] +
569  J[j21] * invJ[j12] +
570  J[j31] * invJ[j13] ;
571 
572  const double detJinv = 1.0 / detJ ;
573 
574  for ( unsigned i = 0 ; i < TensorDim ; ++i ) { invJ[i] *= detJinv ; }
575 
576  // Transform gradients:
577 
578  for( unsigned i = 0; i < FunctionCount ; ++i ) {
579  const double g0 = grad[0][i];
580  const double g1 = grad[1][i];
581  const double g2 = grad[2][i];
582 
583  dpsidx[i] = g0 * invJ[j11] + g1 * invJ[j12] + g2 * invJ[j13];
584  dpsidy[i] = g0 * invJ[j21] + g1 * invJ[j22] + g2 * invJ[j23];
585  dpsidz[i] = g0 * invJ[j31] + g1 * invJ[j32] + g2 * invJ[j33];
586  }
587 
588  return detJ ;
589  }
590 
591 };
592 
594  Analytic,
595  FadElement,
598 };
599 
600 template< class FiniteElementMeshType ,
601  class SparseMatrixType ,
602  AssemblyMethod Method
603  >
604 class ElementComputation ;
605 
606 template< class ExecutionSpace , BoxElemPart::ElemOrder Order ,
607  class CoordinateMap , typename ScalarType >
608 class ElementComputation
609  < Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
610  CrsMatrix< ScalarType , ExecutionSpace > ,
611  Analytic > :
612  public ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
613  ScalarType> {
614 public:
615 
616  typedef ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
617  ScalarType> base_type;
618 
621 
622  static const unsigned FunctionCount = base_type::FunctionCount;
623  static const unsigned IntegrationCount = base_type::IntegrationCount;
624  static const unsigned ElemNodeCount = base_type::ElemNodeCount;
625 
626  typedef Kokkos::View<scalar_type[FunctionCount],Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged> elem_vec_type;
627  typedef Kokkos::View<scalar_type[FunctionCount][FunctionCount],Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged> elem_mat_type;
628 
630 
632  const typename base_type::mesh_type & arg_mesh ,
633  const typename base_type::vector_type & arg_solution ,
634  const typename base_type::elem_graph_type & arg_elem_graph ,
635  const typename base_type::sparse_matrix_type & arg_jacobian ,
636  const typename base_type::vector_type & arg_residual ) :
637  base_type(arg_mesh, arg_solution, arg_elem_graph,
638  arg_jacobian, arg_residual) {}
639 
640  //------------------------------------
641 
642  void apply() const
643  {
644  const size_t nelem = this->elem_node_ids.extent(0);
645  parallel_for( nelem , *this );
646  }
647 
649  void gatherSolution(const unsigned ielem,
650  const elem_vec_type& val,
651  unsigned node_index[],
652  double x[], double y[], double z[],
653  const elem_vec_type& res,
654  const elem_mat_type& mat) const
655  {
656  for ( unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
657  const unsigned ni = this->elem_node_ids( ielem , i );
658 
659  node_index[i] = ni ;
660 
661  x[i] = this->node_coords( ni , 0 );
662  y[i] = this->node_coords( ni , 1 );
663  z[i] = this->node_coords( ni , 2 );
664 
665  val(i) = this->solution( ni ) ;
666  res(i) = 0 ;
667 
668  for( unsigned j = 0; j < FunctionCount ; j++){
669  mat(i,j) = 0 ;
670  }
671  }
672  }
673 
675  void scatterResidual(const unsigned ielem,
676  const unsigned node_index[],
677  const elem_vec_type& res,
678  const elem_mat_type& mat) const
679  {
680  for( unsigned i = 0 ; i < FunctionCount ; i++ ) {
681  const unsigned row = node_index[i] ;
682  if ( row < this->residual.extent(0) ) {
683  atomic_add( & this->residual( row ) , res(i) );
684 
685  for( unsigned j = 0 ; j < FunctionCount ; j++ ) {
686  const unsigned entry = this->elem_graph( ielem , i , j );
687  if ( entry != ~0u ) {
688  atomic_add( & this->jacobian.coeff( entry ) , mat(i,j) );
689  }
690  }
691  }
692  }
693  }
694 
697  const elem_vec_type& dof_values ,
698  const double x[],
699  const double y[],
700  const double z[],
701  const elem_vec_type& elem_res ,
702  const elem_mat_type& elem_mat ) const
703  {
704  double coeff_k = 3.456;
705  double coeff_src = 1.234;
706  double advection[] = { 1.1, 1.2, 1.3 };
707  double dpsidx[ FunctionCount ] ;
708  double dpsidy[ FunctionCount ] ;
709  double dpsidz[ FunctionCount ] ;
710  for ( unsigned i = 0 ; i < IntegrationCount ; ++i ) {
711 
712  const double integ_weight = this->elem_data.weights[i];
713  const double* bases_vals = this->elem_data.values[i];
714  const double detJ =
715  this->transform_gradients( this->elem_data.gradients[i] ,
716  x , y , z ,
717  dpsidx , dpsidy , dpsidz );
718  const double detJ_weight = detJ * integ_weight;
719  const double detJ_weight_coeff_k = detJ_weight * coeff_k;
720 
721  scalar_type value_at_pt = 0 ;
722  scalar_type gradx_at_pt = 0 ;
723  scalar_type grady_at_pt = 0 ;
724  scalar_type gradz_at_pt = 0 ;
725  for ( unsigned m = 0 ; m < FunctionCount ; m++ ) {
726  value_at_pt += dof_values(m) * bases_vals[m] ;
727  gradx_at_pt += dof_values(m) * dpsidx[m] ;
728  grady_at_pt += dof_values(m) * dpsidy[m] ;
729  gradz_at_pt += dof_values(m) * dpsidz[m] ;
730  }
731 
732  const scalar_type source_term =
733  coeff_src * value_at_pt * value_at_pt ;
734  const scalar_type source_deriv =
735  2.0 * coeff_src * value_at_pt ;
736 
737  const scalar_type advection_x = advection[0];
738  const scalar_type advection_y = advection[1];
739  const scalar_type advection_z = advection[2];
740 
741  const scalar_type advection_term =
742  advection_x*gradx_at_pt +
743  advection_y*grady_at_pt +
744  advection_z*gradz_at_pt ;
745 
746  for ( unsigned m = 0; m < FunctionCount; ++m) {
747  const double bases_val_m = bases_vals[m] * detJ_weight ;
748  const double dpsidx_m = dpsidx[m] ;
749  const double dpsidy_m = dpsidy[m] ;
750  const double dpsidz_m = dpsidz[m] ;
751 
752  elem_res(m) +=
753  detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
754  dpsidy_m * grady_at_pt +
755  dpsidz_m * gradz_at_pt ) +
756  bases_val_m * ( advection_term + source_term ) ;
757 
758  for( unsigned n = 0; n < FunctionCount; n++) {
759  const double dpsidx_n = dpsidx[n] ;
760  const double dpsidy_n = dpsidy[n] ;
761  const double dpsidz_n = dpsidz[n] ;
762  elem_mat(m,n) +=
763  detJ_weight_coeff_k * ( dpsidx_m * dpsidx_n +
764  dpsidy_m * dpsidy_n +
765  dpsidz_m * dpsidz_n ) +
766  bases_val_m * ( advection_x * dpsidx_n +
767  advection_y * dpsidy_n +
768  advection_z * dpsidz_n +
769  source_deriv * bases_vals[n] ) ;
770  }
771  }
772  }
773  }
774 
776  void operator()( const unsigned ielem ) const
777  {
778  double x[ FunctionCount ] ;
779  double y[ FunctionCount ] ;
780  double z[ FunctionCount ] ;
781  unsigned node_index[ ElemNodeCount ];
782 
783  scalar_type local_val[FunctionCount], local_res[FunctionCount], local_mat[FunctionCount*FunctionCount];
784  elem_vec_type val(local_val, FunctionCount) ;
785  elem_vec_type elem_res(local_res, FunctionCount ) ;
786  elem_mat_type elem_mat(local_mat, FunctionCount, FunctionCount ) ;
787 
788  // Gather nodal coordinates and solution vector:
789  gatherSolution(ielem, val, node_index, x, y, z, elem_res, elem_mat);
790 
791  // Compute nodal element residual vector and Jacobian matrix
792  computeElementResidualJacobian( val, x, y, z, elem_res , elem_mat );
793 
794  // Scatter nodal element residual and Jacobian in global vector and matrix:
795  scatterResidual( ielem, node_index, elem_res, elem_mat );
796  }
797 }; /* ElementComputation */
798 
799 template< class ExecutionSpace , BoxElemPart::ElemOrder Order ,
800  class CoordinateMap , typename ScalarType >
801 class ElementComputation
802  < Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
803  CrsMatrix< ScalarType , ExecutionSpace > ,
804  FadElement > : public ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
805  ScalarType> {
806 public:
807 
808  typedef ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
809  ScalarType> base_type;
810 
813 
814  static const unsigned FunctionCount = base_type::FunctionCount;
815  static const unsigned IntegrationCount = base_type::IntegrationCount;
816  static const unsigned ElemNodeCount = base_type::ElemNodeCount;
817 
819 
820  //typedef Kokkos::View<fad_scalar_type[FunctionCount],Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged> elem_vec_type;
821  typedef Kokkos::View<fad_scalar_type*,Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged> elem_vec_type; // possibly fix warning and internal compiler error with gcc 4.7.2????
822 
824 
826  const typename base_type::mesh_type & arg_mesh ,
827  const typename base_type::vector_type & arg_solution ,
828  const typename base_type::elem_graph_type & arg_elem_graph ,
829  const typename base_type::sparse_matrix_type & arg_jacobian ,
830  const typename base_type::vector_type & arg_residual ) :
831  base_type(arg_mesh, arg_solution, arg_elem_graph,
832  arg_jacobian, arg_residual) {}
833 
834  //------------------------------------
835 
836  void apply() const
837  {
838  const size_t nelem = this->elem_node_ids.extent(0);
839  parallel_for( nelem , *this );
840  }
841 
843  void gatherSolution(const unsigned ielem,
844  const elem_vec_type& val,
845  unsigned node_index[],
846  double x[], double y[], double z[],
847  const elem_vec_type& res) const
848  {
849  for ( unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
850  const unsigned ni = this->elem_node_ids( ielem , i );
851 
852  node_index[i] = ni ;
853 
854  x[i] = this->node_coords( ni , 0 );
855  y[i] = this->node_coords( ni , 1 );
856  z[i] = this->node_coords( ni , 2 );
857 
858  val(i).val() = this->solution( ni );
859  val(i).diff( i, FunctionCount );
860  res(i) = 0.0;
861  }
862  }
863 
865  void scatterResidual(const unsigned ielem,
866  const unsigned node_index[],
867  const elem_vec_type& res) const
868  {
869  for( unsigned i = 0 ; i < FunctionCount ; i++ ) {
870  const unsigned row = node_index[i] ;
871  if ( row < this->residual.extent(0) ) {
872  atomic_add( & this->residual( row ) , res(i).val() );
873 
874  for( unsigned j = 0 ; j < FunctionCount ; j++ ) {
875  const unsigned entry = this->elem_graph( ielem , i , j );
876  if ( entry != ~0u ) {
877  atomic_add( & this->jacobian.coeff( entry ) ,
878  res(i).fastAccessDx(j) );
879  }
880  }
881  }
882  }
883  }
884 
886  void computeElementResidual(const elem_vec_type& dof_values ,
887  const double x[],
888  const double y[],
889  const double z[],
890  const elem_vec_type& elem_res ) const
891  {
892  double coeff_k = 3.456;
893  double coeff_src = 1.234;
894  double advection[] = { 1.1, 1.2, 1.3 };
895  double dpsidx[ FunctionCount ] ;
896  double dpsidy[ FunctionCount ] ;
897  double dpsidz[ FunctionCount ] ;
898  for ( unsigned i = 0 ; i < IntegrationCount ; ++i ) {
899 
900  const double integ_weight = this->elem_data.weights[i];
901  const double* bases_vals = this->elem_data.values[i];
902  const double detJ =
903  this->transform_gradients( this->elem_data.gradients[i] ,
904  x , y , z ,
905  dpsidx , dpsidy , dpsidz );
906  const double detJ_weight = detJ * integ_weight;
907  const double detJ_weight_coeff_k = detJ_weight * coeff_k;
908 
909  fad_scalar_type value_at_pt = 0 ;
910  fad_scalar_type gradx_at_pt = 0 ;
911  fad_scalar_type grady_at_pt = 0 ;
912  fad_scalar_type gradz_at_pt = 0 ;
913  for ( unsigned m = 0 ; m < FunctionCount ; m++ ) {
914  value_at_pt += dof_values(m) * bases_vals[m] ;
915  gradx_at_pt += dof_values(m) * dpsidx[m] ;
916  grady_at_pt += dof_values(m) * dpsidy[m] ;
917  gradz_at_pt += dof_values(m) * dpsidz[m] ;
918  }
919 
920  const fad_scalar_type source_term =
921  coeff_src * value_at_pt * value_at_pt ;
922 
923  const fad_scalar_type advection_term =
924  advection[0]*gradx_at_pt +
925  advection[1]*grady_at_pt +
926  advection[2]*gradz_at_pt;
927 
928  for ( unsigned m = 0; m < FunctionCount; ++m) {
929  const double bases_val_m = bases_vals[m] * detJ_weight ;
930  const double dpsidx_m = dpsidx[m] ;
931  const double dpsidy_m = dpsidy[m] ;
932  const double dpsidz_m = dpsidz[m] ;
933 
934  elem_res(m) +=
935  detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
936  dpsidy_m * grady_at_pt +
937  dpsidz_m * gradz_at_pt ) +
938  bases_val_m * ( advection_term + source_term ) ;
939  }
940  }
941  }
942 
944  void operator()( const unsigned ielem ) const
945  {
946  double x[ FunctionCount ] ;
947  double y[ FunctionCount ] ;
948  double z[ FunctionCount ] ;
949  unsigned node_index[ ElemNodeCount ];
950 
951  scalar_type local_val[FunctionCount*(FunctionCount+1)], local_res[FunctionCount*(FunctionCount+1)];
953  elem_vec_type elem_res(local_res, FunctionCount, FunctionCount+1) ;
954 
955  // Gather nodal coordinates and solution vector:
956  gatherSolution( ielem, val, node_index, x, y, z, elem_res );
957 
958  // Compute nodal element residual vector:
959  computeElementResidual( val, x, y, z, elem_res );
960 
961  // Scatter nodal element residual in global vector:
962  scatterResidual( ielem, node_index, elem_res );
963  }
964 }; /* ElementComputation */
965 
966 template< class ExecutionSpace , BoxElemPart::ElemOrder Order ,
967  class CoordinateMap , typename ScalarType >
968 class ElementComputation
969  < Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
970  CrsMatrix< ScalarType , ExecutionSpace > ,
972  public ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
973  CrsMatrix< ScalarType , ExecutionSpace > ,
974  FadElement > {
975 public:
976 
977  typedef ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
978  CrsMatrix< ScalarType , ExecutionSpace > ,
980 
983 
984  static const unsigned FunctionCount = base_type::FunctionCount;
985  static const unsigned IntegrationCount = base_type::IntegrationCount;
986  static const unsigned ElemNodeCount = base_type::ElemNodeCount;
987 
988  typedef typename base_type::fad_scalar_type fad_scalar_type;
989 
990  typedef Kokkos::View<scalar_type[FunctionCount],Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged> scalar_elem_vec_type;
991  typedef typename base_type::elem_vec_type fad_elem_vec_type;
992 
994 
996  const typename base_type::mesh_type & arg_mesh ,
997  const typename base_type::vector_type & arg_solution ,
998  const typename base_type::elem_graph_type & arg_elem_graph ,
999  const typename base_type::sparse_matrix_type & arg_jacobian ,
1000  const typename base_type::vector_type & arg_residual ) :
1001  base_type(arg_mesh, arg_solution, arg_elem_graph,
1002  arg_jacobian, arg_residual) {}
1003 
1004  //------------------------------------
1005 
1006  void apply() const
1007  {
1008  const size_t nelem = this->elem_node_ids.extent(0);
1009  parallel_for( nelem , *this );
1010  }
1011 
1013  void gatherSolution(const unsigned ielem,
1014  const scalar_elem_vec_type& val,
1015  unsigned node_index[],
1016  double x[], double y[], double z[],
1017  const fad_elem_vec_type& res) const
1018  {
1019  for ( unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
1020  const unsigned ni = this->elem_node_ids( ielem , i );
1021 
1022  node_index[i] = ni ;
1023 
1024  x[i] = this->node_coords( ni , 0 );
1025  y[i] = this->node_coords( ni , 1 );
1026  z[i] = this->node_coords( ni , 2 );
1027 
1028  val(i) = this->solution( ni );
1029  res(i) = 0.0;
1030  }
1031  }
1032 
1035  const double x[],
1036  const double y[],
1037  const double z[],
1038  const fad_elem_vec_type& elem_res ) const
1039  {
1040  double coeff_k = 3.456;
1041  double coeff_src = 1.234;
1042  double advection[] = { 1.1, 1.2, 1.3 };
1043  double dpsidx[ FunctionCount ] ;
1044  double dpsidy[ FunctionCount ] ;
1045  double dpsidz[ FunctionCount ] ;
1046  for ( unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1047 
1048  const double integ_weight = this->elem_data.weights[i];
1049  const double* bases_vals = this->elem_data.values[i];
1050  const double detJ =
1051  this->transform_gradients( this->elem_data.gradients[i] ,
1052  x , y , z ,
1053  dpsidx , dpsidy , dpsidz );
1054  const double detJ_weight = detJ * integ_weight;
1055  const double detJ_weight_coeff_k = detJ_weight * coeff_k;
1056 
1061  for ( unsigned m = 0 ; m < FunctionCount ; m++ ) {
1062  value_at_pt.val() += dof_values(m) * bases_vals[m] ;
1063  value_at_pt.fastAccessDx(m) = bases_vals[m] ;
1064 
1065  gradx_at_pt.val() += dof_values(m) * dpsidx[m] ;
1066  gradx_at_pt.fastAccessDx(m) = dpsidx[m] ;
1067 
1068  grady_at_pt.val() += dof_values(m) * dpsidy[m] ;
1069  grady_at_pt.fastAccessDx(m) = dpsidy[m] ;
1070 
1071  gradz_at_pt.val() += dof_values(m) * dpsidz[m] ;
1072  gradz_at_pt.fastAccessDx(m) = dpsidz[m] ;
1073  }
1074 
1075  const fad_scalar_type source_term =
1076  coeff_src * value_at_pt * value_at_pt ;
1077 
1078  const fad_scalar_type advection_term =
1079  advection[0]*gradx_at_pt +
1080  advection[1]*grady_at_pt +
1081  advection[2]*gradz_at_pt;
1082 
1083  for ( unsigned m = 0; m < FunctionCount; ++m) {
1084  const double bases_val_m = bases_vals[m] * detJ_weight ;
1085  const double dpsidx_m = dpsidx[m] ;
1086  const double dpsidy_m = dpsidy[m] ;
1087  const double dpsidz_m = dpsidz[m] ;
1088 
1089  elem_res(m) +=
1090  detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
1091  dpsidy_m * grady_at_pt +
1092  dpsidz_m * gradz_at_pt ) +
1093  bases_val_m * ( advection_term + source_term ) ;
1094  }
1095  }
1096  }
1097 
1099  void operator()( const unsigned ielem ) const
1100  {
1101  double x[ FunctionCount ] ;
1102  double y[ FunctionCount ] ;
1103  double z[ FunctionCount ] ;
1104  unsigned node_index[ ElemNodeCount ];
1105 
1106  scalar_type local_val[FunctionCount], local_res[FunctionCount*(FunctionCount+1)];
1107  scalar_elem_vec_type val(local_val, FunctionCount) ;
1108  fad_elem_vec_type elem_res(local_res, FunctionCount, FunctionCount+1) ;
1109 
1110  // Gather nodal coordinates and solution vector:
1111  gatherSolution( ielem, val, node_index, x, y, z, elem_res );
1112 
1113  // Compute nodal element residual vector:
1114  computeElementResidual( val, x, y, z, elem_res );
1115 
1116  // Scatter nodal element residual in global vector:
1117  this->scatterResidual( ielem, node_index, elem_res );
1118  }
1119 }; /* ElementComputation */
1120 
1121 template< class ExecutionSpace , BoxElemPart::ElemOrder Order ,
1122  class CoordinateMap , typename ScalarType >
1123 class ElementComputation
1124  < Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
1125  CrsMatrix< ScalarType , ExecutionSpace > ,
1126  FadQuadPoint > :
1127  public ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
1128  CrsMatrix< ScalarType , ExecutionSpace > ,
1129  Analytic > {
1130 public:
1131 
1132  typedef ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
1133  CrsMatrix< ScalarType , ExecutionSpace > ,
1135 
1138 
1139  static const unsigned FunctionCount = base_type::FunctionCount;
1140  static const unsigned IntegrationCount = base_type::IntegrationCount;
1141  static const unsigned ElemNodeCount = base_type::ElemNodeCount;
1142 
1143  typedef typename base_type::elem_vec_type elem_vec_type;
1144  typedef typename base_type::elem_mat_type elem_mat_type;
1145 
1147 
1149 
1151  const typename base_type::mesh_type & arg_mesh ,
1152  const typename base_type::vector_type & arg_solution ,
1153  const typename base_type::elem_graph_type & arg_elem_graph ,
1154  const typename base_type::sparse_matrix_type & arg_jacobian ,
1155  const typename base_type::vector_type & arg_residual ) :
1156  base_type(arg_mesh, arg_solution, arg_elem_graph,
1157  arg_jacobian, arg_residual) {}
1158 
1159  //------------------------------------
1160 
1161  void apply() const
1162  {
1163  const size_t nelem = this->elem_node_ids.extent(0);
1164  parallel_for( nelem , *this );
1165  }
1166 
1169  const elem_vec_type& dof_values ,
1170  const double x[],
1171  const double y[],
1172  const double z[],
1173  const elem_vec_type& elem_res ,
1174  const elem_mat_type& elem_mat ) const
1175  {
1176  double coeff_k = 3.456;
1177  double coeff_src = 1.234;
1178  double advection[] = { 1.1, 1.2, 1.3 };
1179  double dpsidx[ FunctionCount ] ;
1180  double dpsidy[ FunctionCount ] ;
1181  double dpsidz[ FunctionCount ] ;
1182 
1183  fad_scalar_type value_at_pt(4, 0, 0.0) ;
1184  fad_scalar_type gradx_at_pt(4, 1, 0.0) ;
1185  fad_scalar_type grady_at_pt(4, 2, 0.0) ;
1186  fad_scalar_type gradz_at_pt(4, 3, 0.0) ;
1187  for ( unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1188 
1189  const double integ_weight = this->elem_data.weights[i];
1190  const double* bases_vals = this->elem_data.values[i];
1191  const double detJ =
1192  this->transform_gradients( this->elem_data.gradients[i] ,
1193  x , y , z ,
1194  dpsidx , dpsidy , dpsidz );
1195  const double detJ_weight = detJ * integ_weight;
1196  const double detJ_weight_coeff_k = detJ_weight * coeff_k;
1197 
1198  value_at_pt.val() = 0.0 ;
1199  gradx_at_pt.val() = 0.0 ;
1200  grady_at_pt.val() = 0.0 ;
1201  gradz_at_pt.val() = 0.0 ;
1202  for ( unsigned m = 0 ; m < FunctionCount ; m++ ) {
1203  value_at_pt.val() += dof_values(m) * bases_vals[m] ;
1204  gradx_at_pt.val() += dof_values(m) * dpsidx[m] ;
1205  grady_at_pt.val() += dof_values(m) * dpsidy[m] ;
1206  gradz_at_pt.val() += dof_values(m) * dpsidz[m] ;
1207  }
1208 
1209  const fad_scalar_type source_term =
1210  coeff_src * value_at_pt * value_at_pt ;
1211 
1212  const fad_scalar_type advection_term =
1213  advection[0]*gradx_at_pt +
1214  advection[1]*grady_at_pt +
1215  advection[2]*gradz_at_pt;
1216 
1217  for ( unsigned m = 0; m < FunctionCount; ++m) {
1218  const double bases_val_m = bases_vals[m] * detJ_weight ;
1219  fad_scalar_type res =
1220  detJ_weight_coeff_k * ( dpsidx[m] * gradx_at_pt +
1221  dpsidy[m] * grady_at_pt +
1222  dpsidz[m] * gradz_at_pt ) +
1223  bases_val_m * ( advection_term + source_term ) ;
1224 
1225  elem_res(m) += res.val();
1226 
1227  for( unsigned n = 0; n < FunctionCount; n++) {
1228  elem_mat(m,n) += res.fastAccessDx(0) * bases_vals[n] +
1229  res.fastAccessDx(1) * dpsidx[n] +
1230  res.fastAccessDx(2) * dpsidy[n] +
1231  res.fastAccessDx(3) * dpsidz[n];
1232  }
1233  }
1234  }
1235  }
1236 
1238  void operator()( const unsigned ielem ) const
1239  {
1240  double x[ FunctionCount ] ;
1241  double y[ FunctionCount ] ;
1242  double z[ FunctionCount ] ;
1243  unsigned node_index[ ElemNodeCount ];
1244 
1245  scalar_type local_val[FunctionCount], local_res[FunctionCount], local_mat[FunctionCount*FunctionCount];
1246  elem_vec_type val(local_val, FunctionCount) ;
1247  elem_vec_type elem_res(local_res, FunctionCount) ;
1248  elem_mat_type elem_mat(local_mat, FunctionCount, FunctionCount) ;
1249 
1250  // Gather nodal coordinates and solution vector:
1251  this->gatherSolution( ielem, val, node_index, x, y, z, elem_res, elem_mat );
1252 
1253  // Compute nodal element residual vector and Jacobian matrix:
1254  computeElementResidualJacobian( val, x, y, z, elem_res, elem_mat );
1255 
1256  // Scatter nodal element residual and Jacobian in global vector and matrix:
1257  this->scatterResidual( ielem, node_index, elem_res, elem_mat );
1258  }
1259 }; /* ElementComputation */
1260 
1261 } /* namespace FENL */
1262 } /* namespace Example */
1263 } /* namespace Kokkos */
1264 
1265 //----------------------------------------------------------------------------
1266 
1267 #endif /* #ifndef KOKKOS_EXAMPLE_FENLFUNCTORS_HPP */
double weights[integration_count]
Definition: HexElement.hpp:219
ElementComputation(const typename base_type::mesh_type &arg_mesh, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual)
KOKKOS_INLINE_FUNCTION void computeElementResidualJacobian(const elem_vec_type &dof_values, const double x[], const double y[], const double z[], const elem_vec_type &elem_res, const elem_mat_type &elem_mat) const
KOKKOS_INLINE_FUNCTION void fill_graph_entries(const unsigned iset) const
Kokkos::View< scalar_type *[FunctionCount], execution_space > elem_vectors_type
static const unsigned element_node_count
Definition: HexElement.hpp:215
CrsMatrix< ScalarType, ExecutionSpace > sparse_matrix_type
KOKKOS_INLINE_FUNCTION void computeElementResidualJacobian(const elem_vec_type &dof_values, const double x[], const double y[], const double z[], const elem_vec_type &elem_res, const elem_mat_type &elem_mat) const
static const unsigned function_count
Definition: HexElement.hpp:217
ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap >, CrsMatrix< ScalarType, ExecutionSpace >, Analytic > base_type
Kokkos::View< const double *[SpaceDim], Device > node_coord_type
double gradients[integration_count][spatial_dimension][function_count]
Definition: HexElement.hpp:221
Kokkos::View< const unsigned *[ElemNode], Device > elem_node_type
KOKKOS_INLINE_FUNCTION void join(volatile unsigned &update, volatile const unsigned &input) const
Kokkos::View< scalar_type *, Kokkos::LayoutLeft, execution_space > vector_type
KOKKOS_INLINE_FUNCTION void computeElementResidual(const scalar_elem_vec_type dof_values, const double x[], const double y[], const double z[], const fad_elem_vec_type &elem_res) const
NodeNodeGraph< elem_node_type, sparse_graph_type, ElemNodeCount >::ElemGraphType elem_graph_type
expr val()
#define KOKKOS_INLINE_FUNCTION
ElemNodeIdView::execution_space execution_space
KOKKOS_INLINE_FUNCTION void operator()(const unsigned irow, unsigned &update, const bool final) const
KOKKOS_INLINE_FUNCTION void init(unsigned &update) const
ElementComputationBase(const ElementComputationBase &rhs)
KOKKOS_INLINE_FUNCTION void fill_elem_graph_map(const unsigned ielem) const
Kokkos::Example::HexElement_Data< mesh_type::ElemNode > element_data_type
CrsGraphType::row_map_type::non_const_type RowMapType
ElementComputationBase(const mesh_type &arg_mesh, const vector_type &arg_solution, const elem_graph_type &arg_elem_graph, const sparse_matrix_type &arg_jacobian, const vector_type &arg_residual)
KOKKOS_INLINE_FUNCTION void operator()(const TagFillNodeSet &, unsigned ielem, unsigned &count) const
CrsMatrix(const StaticCrsGraphType &arg_graph)
Kokkos::StaticCrsGraph< unsigned, Space, void, void, unsigned > StaticCrsGraphType
KOKKOS_INLINE_FUNCTION double transform_gradients(const double grad[][FunctionCount], const double x[], const double y[], const double z[], double dpsidx[], double dpsidy[], double dpsidz[]) const
double values[integration_count][function_count]
Definition: HexElement.hpp:220
Do not initialize the derivative array.
KOKKOS_INLINE_FUNCTION void operator()(const unsigned iwork) const
static const unsigned integration_count
Definition: HexElement.hpp:216
Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap > mesh_type
Kokkos::View< scalar_type[FunctionCount][FunctionCount], Kokkos::LayoutRight, execution_space, Kokkos::MemoryUnmanaged > elem_mat_type
NodeNodeGraph(const ElemNodeIdView &arg_elem_node_id, const unsigned arg_node_count, Times &results)
KOKKOS_INLINE_FUNCTION void init(const TagFillNodeSet &, unsigned &update) const
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, const scalar_elem_vec_type &val, unsigned node_index[], double x[], double y[], double z[], const fad_elem_vec_type &res) const
View< ValueType *, Space > coeff_type
sparse_matrix_type::StaticCrsGraphType sparse_graph_type
Kokkos::View< unsigned, execution_space > UnsignedValue
ElementComputation(const typename base_type::mesh_type &arg_mesh, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual)
ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap >, CrsMatrix< ScalarType, ExecutionSpace >, FadElement > base_type
ElementComputation(const typename base_type::mesh_type &arg_mesh, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual)
ElementComputation(const typename base_type::mesh_type &arg_mesh, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual)
KOKKOS_INLINE_FUNCTION void scatterResidual(const unsigned ielem, const unsigned node_index[], const elem_vec_type &res, const elem_mat_type &mat) const
KOKKOS_INLINE_FUNCTION void sort_graph_entries(const unsigned irow) const
Kokkos::View< scalar_type *[FunctionCount][FunctionCount], execution_space > elem_matrices_type
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, const elem_vec_type &val, unsigned node_index[], double x[], double y[], double z[], const elem_vec_type &res) const
KOKKOS_INLINE_FUNCTION void join(const TagFillNodeSet &, volatile unsigned &update, volatile const unsigned &input) const
Kokkos::View< unsigned *[ElemNode][ElemNode], execution_space > ElemGraphType
int n
Generate a distributed unstructured finite element mesh from a partitioned NX*NY*NZ box of elements...
KOKKOS_INLINE_FUNCTION void computeElementResidual(const elem_vec_type &dof_values, const double x[], const double y[], const double z[], const elem_vec_type &elem_res) const
Kokkos::UnorderedMap< key_type, void, execution_space > SetType
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, const elem_vec_type &val, unsigned node_index[], double x[], double y[], double z[], const elem_vec_type &res, const elem_mat_type &mat) const
static const unsigned spatial_dimension
Definition: HexElement.hpp:214