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