epsilon_divdiv_kerngen_v01_initial.hpp

Go to the documentation of this file.

#pragma once
 
#include "../../../quadrature/quadrature.hpp"
#include "communication/shell/communication.hpp"
#include "dense/vec.hpp"
#include "fe/wedge/integrands.hpp"
#include "fe/wedge/kernel_helpers.hpp"
#include "grid/shell/spherical_shell.hpp"
#include "linalg/operator.hpp"
#include "linalg/solvers/gca/local_matrix_storage.hpp"
#include "linalg/trafo/local_basis_trafo_normal_tangential.hpp"
#include "linalg/vector.hpp"
#include "linalg/vector_q1.hpp"
#include "util/timer.hpp"
 
namespace terra::fe::wedge::operators::shell::epsdivdiv_history {
 
using terra::grid::shell::BoundaryConditionFlag;
using terra::grid::shell::BoundaryConditions;
using terra::grid::shell::ShellBoundaryFlag;
 
template < typename ScalarT, int VecDim = 3 >
class EpsilonDivDivKerngenV01Initial
{
  public:
    using SrcVectorType           = linalg::VectorQ1Vec< ScalarT, VecDim >;
    using DstVectorType           = linalg::VectorQ1Vec< ScalarT, VecDim >;
    using ScalarType              = ScalarT;
    using Grid4DDataLocalMatrices = terra::grid::Grid4DDataMatrices< ScalarType, 18, 18, 2 >;
 
  private:
    bool storeLMatrices_ =
        false; // set to let apply_impl() know, that it should store the local matrices after assembling them
    bool applyStoredLMatrices_ =
        false; // set to make apply_impl() load and use the stored LMatrices for the operator application
    Grid4DDataLocalMatrices lmatrices_;
    bool                    single_quadpoint_ = true;
 
    grid::shell::DistributedDomain domain_;
 
    grid::Grid3DDataVec< ScalarT, 3 >                        grid_;
    grid::Grid2DDataScalar< ScalarT >                        radii_;
    grid::Grid4DDataScalar< ScalarType >                     k_;
    grid::Grid4DDataScalar< grid::shell::ShellBoundaryFlag > mask_;
    BoundaryConditions                                       bcs_;
 
    bool treat_boundary_;
    bool diagonal_;
 
    linalg::OperatorApplyMode         operator_apply_mode_;
    linalg::OperatorCommunicationMode operator_communication_mode_;
    linalg::OperatorStoredMatrixMode  operator_stored_matrix_mode_;
 
    communication::shell::SubdomainNeighborhoodSendRecvBuffer< ScalarT, VecDim > send_buffers_;
    communication::shell::SubdomainNeighborhoodSendRecvBuffer< ScalarT, VecDim > recv_buffers_;
 
    grid::Grid4DDataVec< ScalarType, VecDim > src_;
    grid::Grid4DDataVec< ScalarType, VecDim > dst_;
 
  public:
    EpsilonDivDivKerngenV01Initial(
        const grid::shell::DistributedDomain&                           domain,
        const grid::Grid3DDataVec< ScalarT, 3 >&                        grid,
        const grid::Grid2DDataScalar< ScalarT >&                        radii,
        const grid::Grid4DDataScalar< grid::shell::ShellBoundaryFlag >& mask,
        const grid::Grid4DDataScalar< ScalarT >&                        k,
        BoundaryConditions                                              bcs,
        bool                                                            diagonal,
        linalg::OperatorApplyMode         operator_apply_mode = linalg::OperatorApplyMode::Replace,
        linalg::OperatorCommunicationMode operator_communication_mode =
            linalg::OperatorCommunicationMode::CommunicateAdditively,
        linalg::OperatorStoredMatrixMode operator_stored_matrix_mode = linalg::OperatorStoredMatrixMode::Off )
    : domain_( domain )
    , grid_( grid )
    , radii_( radii )
    , k_( k )
    , mask_( mask )
    , treat_boundary_( true )
    , diagonal_( diagonal )
    , operator_apply_mode_( operator_apply_mode )
    , operator_communication_mode_( operator_communication_mode )
    , operator_stored_matrix_mode_( operator_stored_matrix_mode )
    // TODO: we can reuse the send and recv buffers and pass in from the outside somehow
    , send_buffers_( domain )
    , recv_buffers_( domain )
    {
        bcs_[0] = bcs[0];
        bcs_[1] = bcs[1];
    }
 
    const grid::Grid4DDataScalar< ScalarType >& k_grid_data() { return k_; }
 
    /// @brief Getter for domain member
    const grid::shell::DistributedDomain& get_domain() const { return domain_; }
 
    /// @brief Getter for radii member
    grid::Grid2DDataScalar< ScalarT > get_radii() const { return radii_; }
 
    /// @brief Getter for grid member
    grid::Grid3DDataVec< ScalarT, 3 > get_grid() const { return grid_; }
 
    /// @brief Retrives the local matrix stored in the operator
    KOKKOS_INLINE_FUNCTION
    dense::Mat< ScalarT, 6, 6 > get_lmatrix(
        const int local_subdomain_id,
        const int x_cell,
        const int y_cell,
        const int r_cell,
        const int wedge,
        const int dimi,
        const int dimj ) const
    {
        assert( lmatrices_.data() != nullptr );
        dense::Mat< ScalarT, 6, 6 > ijslice;
        for ( int i = 0; i < 6; ++i )
        {
            for ( int j = 0; j < 6; ++j )
            {
                ijslice( i, j ) =
                    lmatrices_( local_subdomain_id, x_cell, y_cell, r_cell, wedge )( i + dimi * 6, j + dimj * 6 );
            }
        }
        return ijslice;
    }
 
    /// @brief Set the local matrix stored in the operator
    KOKKOS_INLINE_FUNCTION
    void set_lmatrix(
        const int                   local_subdomain_id,
        const int                   x_cell,
        const int                   y_cell,
        const int                   r_cell,
        const int                   wedge,
        const int                   dimi,
        const int                   dimj,
        dense::Mat< ScalarT, 6, 6 > mat ) const
    {
        assert( lmatrices_.data() != nullptr );
        for ( int i = 0; i < 6; ++i )
        {
            for ( int j = 0; j < 6; ++j )
            {
                lmatrices_( local_subdomain_id, x_cell, y_cell, r_cell, wedge )( i + dimi * 6, j + dimj * 6 ) =
                    mat( i, j );
            }
        }
    }
 
    /// @brief Setter/Getter for app applyStoredLMatrices_: usage of stored local matrices during apply
    void setApplyStoredLMatrices( bool v ) { applyStoredLMatrices_ = v; }
 
    /// @brief S/Getter for diagonal member
    void set_diagonal( bool v ) { diagonal_ = v; }
 
    /// @brief
    /// allocates memory for the local matrices
    /// calls kernel with storeLMatrices_ = true to assemble and store the local matrices
    /// sets applyStoredLMatrices_, such that future applies use the stored local matrices
    void store_lmatrices()
    {
        storeLMatrices_ = true;
        if ( lmatrices_.data() == nullptr )
        {
            lmatrices_ = Grid4DDataLocalMatrices(
                "LaplaceSimple::lmatrices_",
                domain_.subdomains().size(),
                domain_.domain_info().subdomain_num_nodes_per_side_laterally() - 1,
                domain_.domain_info().subdomain_num_nodes_per_side_laterally() - 1,
                domain_.domain_info().subdomain_num_nodes_radially() - 1 );
            Kokkos::parallel_for(
                "assemble_store_lmatrices", grid::shell::local_domain_md_range_policy_cells( domain_ ), *this );
            Kokkos::fence();
        }
        storeLMatrices_       = false;
        applyStoredLMatrices_ = true;
    }
 
    void set_operator_apply_and_communication_modes(
        const linalg::OperatorApplyMode         operator_apply_mode,
        const linalg::OperatorCommunicationMode operator_communication_mode )
    {
        operator_apply_mode_         = operator_apply_mode;
        operator_communication_mode_ = operator_communication_mode;
    }
 
    void apply_impl( const SrcVectorType& src, DstVectorType& dst )
    {
        util::Timer timer_apply( "vector_laplace_apply" );
        if ( storeLMatrices_ or applyStoredLMatrices_ )
            assert( lmatrices_.data() != nullptr );
 
        if ( operator_apply_mode_ == linalg::OperatorApplyMode::Replace )
        {
            assign( dst, 0 );
        }
 
        src_ = src.grid_data();
        dst_ = dst.grid_data();
 
        if ( src_.extent( 0 ) != dst_.extent( 0 ) || src_.extent( 1 ) != dst_.extent( 1 ) ||
             src_.extent( 2 ) != dst_.extent( 2 ) || src_.extent( 3 ) != dst_.extent( 3 ) )
        {
            throw std::runtime_error( "VectorLaplace: src/dst mismatch" );
        }
 
        if ( src_.extent( 1 ) != grid_.extent( 1 ) || src_.extent( 2 ) != grid_.extent( 2 ) )
        {
            throw std::runtime_error( "VectorLaplace: src/dst mismatch" );
        }
 
        util::Timer timer_kernel( "vector_laplace_kernel" );
        Kokkos::parallel_for( "matvec", grid::shell::local_domain_md_range_policy_cells( domain_ ), *this );
        Kokkos::fence();
        timer_kernel.stop();
 
        if ( operator_communication_mode_ == linalg::OperatorCommunicationMode::CommunicateAdditively )
        {
            util::Timer timer_comm( "vector_laplace_comm" );
 
            communication::shell::pack_send_and_recv_local_subdomain_boundaries(
                domain_, dst_, send_buffers_, recv_buffers_ );
            communication::shell::unpack_and_reduce_local_subdomain_boundaries( domain_, dst_, recv_buffers_ );
        }
    }
 
    KOKKOS_INLINE_FUNCTION void
        operator()( const int local_subdomain_id, const int x_cell, const int y_cell, const int r_cell ) const
    {
        // Compute the local element matrix.
        dense::Mat< ScalarT, 18, 18 > A[num_wedges_per_hex_cell] = {};
 
        if ( !applyStoredLMatrices_ )
        {
            double wedge_surf_phy_coords[2][3][3];
            double quad_surface_coords[2][2][3];
            quad_surface_coords[0][0][0]   = grid_( local_subdomain_id, x_cell, y_cell, 0 );
            quad_surface_coords[0][0][1]   = grid_( local_subdomain_id, x_cell, y_cell, 1 );
            quad_surface_coords[0][0][2]   = grid_( local_subdomain_id, x_cell, y_cell, 2 );
            quad_surface_coords[0][1][0]   = grid_( local_subdomain_id, x_cell, y_cell + 1, 0 );
            quad_surface_coords[0][1][1]   = grid_( local_subdomain_id, x_cell, y_cell + 1, 1 );
            quad_surface_coords[0][1][2]   = grid_( local_subdomain_id, x_cell, y_cell + 1, 2 );
            quad_surface_coords[1][0][0]   = grid_( local_subdomain_id, x_cell + 1, y_cell, 0 );
            quad_surface_coords[1][0][1]   = grid_( local_subdomain_id, x_cell + 1, y_cell, 1 );
            quad_surface_coords[1][0][2]   = grid_( local_subdomain_id, x_cell + 1, y_cell, 2 );
            quad_surface_coords[1][1][0]   = grid_( local_subdomain_id, x_cell + 1, y_cell + 1, 0 );
            quad_surface_coords[1][1][1]   = grid_( local_subdomain_id, x_cell + 1, y_cell + 1, 1 );
            quad_surface_coords[1][1][2]   = grid_( local_subdomain_id, x_cell + 1, y_cell + 1, 2 );
            wedge_surf_phy_coords[0][0][0] = quad_surface_coords[0][0][0];
            wedge_surf_phy_coords[0][0][1] = quad_surface_coords[0][0][1];
            wedge_surf_phy_coords[0][0][2] = quad_surface_coords[0][0][2];
            wedge_surf_phy_coords[0][1][0] = quad_surface_coords[1][0][0];
            wedge_surf_phy_coords[0][1][1] = quad_surface_coords[1][0][1];
            wedge_surf_phy_coords[0][1][2] = quad_surface_coords[1][0][2];
            wedge_surf_phy_coords[0][2][0] = quad_surface_coords[0][1][0];
            wedge_surf_phy_coords[0][2][1] = quad_surface_coords[0][1][1];
            wedge_surf_phy_coords[0][2][2] = quad_surface_coords[0][1][2];
            wedge_surf_phy_coords[1][0][0] = quad_surface_coords[1][1][0];
            wedge_surf_phy_coords[1][0][1] = quad_surface_coords[1][1][1];
            wedge_surf_phy_coords[1][0][2] = quad_surface_coords[1][1][2];
            wedge_surf_phy_coords[1][1][0] = quad_surface_coords[0][1][0];
            wedge_surf_phy_coords[1][1][1] = quad_surface_coords[0][1][1];
            wedge_surf_phy_coords[1][1][2] = quad_surface_coords[0][1][2];
            wedge_surf_phy_coords[1][2][0] = quad_surface_coords[1][0][0];
            wedge_surf_phy_coords[1][2][1] = quad_surface_coords[1][0][1];
            wedge_surf_phy_coords[1][2][2] = quad_surface_coords[1][0][2];
            double r_0                     = radii_( local_subdomain_id, r_cell );
            double r_1                     = radii_( local_subdomain_id, r_cell + 1 );
            double src_local_hex[3][2][6];
            int    dim;
            for ( dim = 0; dim < 3; dim += 1 )
            {
                src_local_hex[dim][0][0] = src_( local_subdomain_id, x_cell, y_cell, r_cell, dim );
                src_local_hex[dim][0][1] = src_( local_subdomain_id, x_cell + 1, y_cell, r_cell, dim );
                src_local_hex[dim][0][2] = src_( local_subdomain_id, x_cell, y_cell + 1, r_cell, dim );
                src_local_hex[dim][0][3] = src_( local_subdomain_id, x_cell, y_cell, r_cell + 1, dim );
                src_local_hex[dim][0][4] = src_( local_subdomain_id, x_cell + 1, y_cell, r_cell + 1, dim );
                src_local_hex[dim][0][5] = src_( local_subdomain_id, x_cell, y_cell + 1, r_cell + 1, dim );
                src_local_hex[dim][1][0] = src_( local_subdomain_id, x_cell + 1, y_cell + 1, r_cell, dim );
                src_local_hex[dim][1][1] = src_( local_subdomain_id, x_cell, y_cell + 1, r_cell, dim );
                src_local_hex[dim][1][2] = src_( local_subdomain_id, x_cell + 1, y_cell, r_cell, dim );
                src_local_hex[dim][1][3] = src_( local_subdomain_id, x_cell + 1, y_cell + 1, r_cell + 1, dim );
                src_local_hex[dim][1][4] = src_( local_subdomain_id, x_cell, y_cell + 1, r_cell + 1, dim );
                src_local_hex[dim][1][5] = src_( local_subdomain_id, x_cell + 1, y_cell, r_cell + 1, dim );
            };
            double k_local_hex[2][6];
            k_local_hex[0][0] = k_( local_subdomain_id, x_cell, y_cell, r_cell );
            k_local_hex[0][1] = k_( local_subdomain_id, x_cell + 1, y_cell, r_cell );
            k_local_hex[0][2] = k_( local_subdomain_id, x_cell, y_cell + 1, r_cell );
            k_local_hex[0][3] = k_( local_subdomain_id, x_cell, y_cell, r_cell + 1 );
            k_local_hex[0][4] = k_( local_subdomain_id, x_cell + 1, y_cell, r_cell + 1 );
            k_local_hex[0][5] = k_( local_subdomain_id, x_cell, y_cell + 1, r_cell + 1 );
            k_local_hex[1][0] = k_( local_subdomain_id, x_cell + 1, y_cell + 1, r_cell );
            k_local_hex[1][1] = k_( local_subdomain_id, x_cell, y_cell + 1, r_cell );
            k_local_hex[1][2] = k_( local_subdomain_id, x_cell + 1, y_cell, r_cell );
            k_local_hex[1][3] = k_( local_subdomain_id, x_cell + 1, y_cell + 1, r_cell + 1 );
            k_local_hex[1][4] = k_( local_subdomain_id, x_cell, y_cell + 1, r_cell + 1 );
            k_local_hex[1][5] = k_( local_subdomain_id, x_cell + 1, y_cell, r_cell + 1 );
            double qp_array[6][3];
            double qw_array[6];
            qp_array[0][0]    = 0.66666666666666663;
            qp_array[1][0]    = 0.16666666666666671;
            qp_array[2][0]    = 0.16666666666666671;
            qp_array[3][0]    = 0.66666666666666663;
            qp_array[4][0]    = 0.16666666666666671;
            qp_array[5][0]    = 0.16666666666666671;
            qp_array[0][1]    = 0.16666666666666671;
            qp_array[1][1]    = 0.66666666666666663;
            qp_array[2][1]    = 0.16666666666666671;
            qp_array[3][1]    = 0.16666666666666671;
            qp_array[4][1]    = 0.66666666666666663;
            qp_array[5][1]    = 0.16666666666666671;
            qp_array[0][2]    = -0.57735026918962573;
            qp_array[1][2]    = -0.57735026918962573;
            qp_array[2][2]    = -0.57735026918962573;
            qp_array[3][2]    = 0.57735026918962573;
            qp_array[4][2]    = 0.57735026918962573;
            qp_array[5][2]    = 0.57735026918962573;
            qw_array[0]       = 0.16666666666666671;
            qw_array[1]       = 0.16666666666666671;
            qw_array[2]       = 0.16666666666666671;
            qw_array[3]       = 0.16666666666666671;
            qw_array[4]       = 0.16666666666666671;
            qw_array[5]       = 0.16666666666666671;
            int cmb_shift     = ( ( treat_boundary_ && diagonal_ == false && r_cell == 0 ) ? ( 3 ) : ( 0 ) );
            int max_rad       = -1 + radii_.extent( 1 );
            int surface_shift = ( ( treat_boundary_ && diagonal_ == false && max_rad == r_cell + 1 ) ? ( 3 ) : ( 0 ) );
            double dst_array[3][2][6] = { 0 };
            int    w                  = 0;
            for ( w = 0; w < 2; w += 1 )
            {
                int q = 0;
                for ( q = 0; q < 6; q += 1 )
                {
                    /* Coefficient evaluation on current wedge w */;
                    double tmpcse_k_eval_0 = ( 1.0 / 2.0 ) * qp_array[q][2];
                    double tmpcse_k_eval_1 = 1.0 / 2.0 - tmpcse_k_eval_0;
                    double tmpcse_k_eval_2 = tmpcse_k_eval_0 + 1.0 / 2.0;
                    double tmpcse_k_eval_3 = -qp_array[q][0] - qp_array[q][1] + 1;
                    double k_eval          = tmpcse_k_eval_1 * tmpcse_k_eval_3 * k_local_hex[w][0] +
                                    tmpcse_k_eval_1 * k_local_hex[w][1] * qp_array[q][0] +
                                    tmpcse_k_eval_1 * k_local_hex[w][2] * qp_array[q][1] +
                                    tmpcse_k_eval_2 * tmpcse_k_eval_3 * k_local_hex[w][3] +
                                    tmpcse_k_eval_2 * k_local_hex[w][4] * qp_array[q][0] +
                                    tmpcse_k_eval_2 * k_local_hex[w][5] * qp_array[q][1];
                    /* Computation + Inversion of the Jacobian */;
                    double tmpcse_J_0 = -1.0 / 2.0 * r_0 + ( 1.0 / 2.0 ) * r_1;
                    double tmpcse_J_1 = r_0 + tmpcse_J_0 * ( qp_array[q][2] + 1 );
                    double tmpcse_J_2 = -qp_array[q][0] - qp_array[q][1] + 1;
                    double J_0_0 = tmpcse_J_1 * ( -wedge_surf_phy_coords[w][0][0] + wedge_surf_phy_coords[w][1][0] );
                    double J_0_1 = tmpcse_J_1 * ( -wedge_surf_phy_coords[w][0][0] + wedge_surf_phy_coords[w][2][0] );
                    double J_0_2 = tmpcse_J_0 * ( tmpcse_J_2 * wedge_surf_phy_coords[w][0][0] +
                                                  qp_array[q][0] * wedge_surf_phy_coords[w][1][0] +
                                                  qp_array[q][1] * wedge_surf_phy_coords[w][2][0] );
                    double J_1_0 = tmpcse_J_1 * ( -wedge_surf_phy_coords[w][0][1] + wedge_surf_phy_coords[w][1][1] );
                    double J_1_1 = tmpcse_J_1 * ( -wedge_surf_phy_coords[w][0][1] + wedge_surf_phy_coords[w][2][1] );
                    double J_1_2 = tmpcse_J_0 * ( tmpcse_J_2 * wedge_surf_phy_coords[w][0][1] +
                                                  qp_array[q][0] * wedge_surf_phy_coords[w][1][1] +
                                                  qp_array[q][1] * wedge_surf_phy_coords[w][2][1] );
                    double J_2_0 = tmpcse_J_1 * ( -wedge_surf_phy_coords[w][0][2] + wedge_surf_phy_coords[w][1][2] );
                    double J_2_1 = tmpcse_J_1 * ( -wedge_surf_phy_coords[w][0][2] + wedge_surf_phy_coords[w][2][2] );
                    double J_2_2 = tmpcse_J_0 * ( tmpcse_J_2 * wedge_surf_phy_coords[w][0][2] +
                                                  qp_array[q][0] * wedge_surf_phy_coords[w][1][2] +
                                                  qp_array[q][1] * wedge_surf_phy_coords[w][2][2] );
                    double J_det = J_0_0 * J_1_1 * J_2_2 - J_0_0 * J_1_2 * J_2_1 - J_0_1 * J_1_0 * J_2_2 +
                                   J_0_1 * J_1_2 * J_2_0 + J_0_2 * J_1_0 * J_2_1 - J_0_2 * J_1_1 * J_2_0;
                    double tmpcse_J_invT_0 = 1.0 / J_det;
                    double J_invT_cse_0_0  = tmpcse_J_invT_0 * ( J_1_1 * J_2_2 - J_1_2 * J_2_1 );
                    double J_invT_cse_0_1  = tmpcse_J_invT_0 * ( -J_1_0 * J_2_2 + J_1_2 * J_2_0 );
                    double J_invT_cse_0_2  = tmpcse_J_invT_0 * ( J_1_0 * J_2_1 - J_1_1 * J_2_0 );
                    double J_invT_cse_1_0  = tmpcse_J_invT_0 * ( -J_0_1 * J_2_2 + J_0_2 * J_2_1 );
                    double J_invT_cse_1_1  = tmpcse_J_invT_0 * ( J_0_0 * J_2_2 - J_0_2 * J_2_0 );
                    double J_invT_cse_1_2  = tmpcse_J_invT_0 * ( -J_0_0 * J_2_1 + J_0_1 * J_2_0 );
                    double J_invT_cse_2_0  = tmpcse_J_invT_0 * ( J_0_1 * J_1_2 - J_0_2 * J_1_1 );
                    double J_invT_cse_2_1  = tmpcse_J_invT_0 * ( -J_0_0 * J_1_2 + J_0_2 * J_1_0 );
                    double J_invT_cse_2_2  = tmpcse_J_invT_0 * ( J_0_0 * J_1_1 - J_0_1 * J_1_0 );
 
                    double scalar_grad[6][3] = { 0 };
                    double tmpcse_grad_i_0   = ( 1.0 / 2.0 ) * qp_array[q][2];
                    double tmpcse_grad_i_1   = tmpcse_grad_i_0 - 1.0 / 2.0;
                    double tmpcse_grad_i_2   = ( 1.0 / 2.0 ) * qp_array[q][0];
                    double tmpcse_grad_i_3   = ( 1.0 / 2.0 ) * qp_array[q][1];
                    double tmpcse_grad_i_4   = tmpcse_grad_i_2 + tmpcse_grad_i_3 - 1.0 / 2.0;
                    double tmpcse_grad_i_5   = J_invT_cse_0_2 * tmpcse_grad_i_2;
                    double tmpcse_grad_i_6   = -tmpcse_grad_i_1;
                    double tmpcse_grad_i_7   = J_invT_cse_1_2 * tmpcse_grad_i_2;
                    double tmpcse_grad_i_8   = J_invT_cse_2_2 * tmpcse_grad_i_2;
                    double tmpcse_grad_i_9   = J_invT_cse_0_2 * tmpcse_grad_i_3;
                    double tmpcse_grad_i_10  = J_invT_cse_1_2 * tmpcse_grad_i_3;
                    double tmpcse_grad_i_11  = J_invT_cse_2_2 * tmpcse_grad_i_3;
                    double tmpcse_grad_i_12  = tmpcse_grad_i_0 + 1.0 / 2.0;
                    double tmpcse_grad_i_13  = -tmpcse_grad_i_12;
                    double tmpcse_grad_i_14  = -tmpcse_grad_i_4;
                    scalar_grad[0][0]        = J_invT_cse_0_0 * tmpcse_grad_i_1 + J_invT_cse_0_1 * tmpcse_grad_i_1 +
                                        J_invT_cse_0_2 * tmpcse_grad_i_4;
                    scalar_grad[0][1] = J_invT_cse_1_0 * tmpcse_grad_i_1 + J_invT_cse_1_1 * tmpcse_grad_i_1 +
                                        J_invT_cse_1_2 * tmpcse_grad_i_4;
                    scalar_grad[0][2] = J_invT_cse_2_0 * tmpcse_grad_i_1 + J_invT_cse_2_1 * tmpcse_grad_i_1 +
                                        J_invT_cse_2_2 * tmpcse_grad_i_4;
                    scalar_grad[1][0] = J_invT_cse_0_0 * tmpcse_grad_i_6 - tmpcse_grad_i_5;
                    scalar_grad[1][1] = J_invT_cse_1_0 * tmpcse_grad_i_6 - tmpcse_grad_i_7;
                    scalar_grad[1][2] = J_invT_cse_2_0 * tmpcse_grad_i_6 - tmpcse_grad_i_8;
                    scalar_grad[2][0] = J_invT_cse_0_1 * tmpcse_grad_i_6 - tmpcse_grad_i_9;
                    scalar_grad[2][1] = J_invT_cse_1_1 * tmpcse_grad_i_6 - tmpcse_grad_i_10;
                    scalar_grad[2][2] = J_invT_cse_2_1 * tmpcse_grad_i_6 - tmpcse_grad_i_11;
                    scalar_grad[3][0] = J_invT_cse_0_0 * tmpcse_grad_i_13 + J_invT_cse_0_1 * tmpcse_grad_i_13 +
                                        J_invT_cse_0_2 * tmpcse_grad_i_14;
                    scalar_grad[3][1] = J_invT_cse_1_0 * tmpcse_grad_i_13 + J_invT_cse_1_1 * tmpcse_grad_i_13 +
                                        J_invT_cse_1_2 * tmpcse_grad_i_14;
                    scalar_grad[3][2] = J_invT_cse_2_0 * tmpcse_grad_i_13 + J_invT_cse_2_1 * tmpcse_grad_i_13 +
                                        J_invT_cse_2_2 * tmpcse_grad_i_14;
                    scalar_grad[4][0] = J_invT_cse_0_0 * tmpcse_grad_i_12 + tmpcse_grad_i_5;
                    scalar_grad[4][1] = J_invT_cse_1_0 * tmpcse_grad_i_12 + tmpcse_grad_i_7;
                    scalar_grad[4][2] = J_invT_cse_2_0 * tmpcse_grad_i_12 + tmpcse_grad_i_8;
                    scalar_grad[5][0] = J_invT_cse_0_1 * tmpcse_grad_i_12 + tmpcse_grad_i_9;
                    scalar_grad[5][1] = J_invT_cse_1_1 * tmpcse_grad_i_12 + tmpcse_grad_i_10;
                    scalar_grad[5][2] = J_invT_cse_2_1 * tmpcse_grad_i_12 + tmpcse_grad_i_11;
                    int dimi;
                    int dimj;
                    for ( dimi = 0; dimi < 3; dimi += 1 )
                    {
                        for ( dimj = 0; dimj < 3; dimj += 1 )
                        {
                            if ( diagonal_ == false )
                            {
                                double grad_u[3][3] = { 0 };
                                double div_u        = 0.0;
                                int    node_idx;
                                for ( node_idx = cmb_shift; node_idx < 6 - surface_shift; node_idx += 1 )
                                {
                                    double E_grad_trial[3][3]     = { 0 };
                                    E_grad_trial[0][dimj]         = scalar_grad[node_idx][0];
                                    E_grad_trial[1][dimj]         = scalar_grad[node_idx][1];
                                    E_grad_trial[2][dimj]         = scalar_grad[node_idx][2];
                                    double tmpcse_symgrad_trial_0 = 0.5 * E_grad_trial[0][1] + 0.5 * E_grad_trial[1][0];
                                    double tmpcse_symgrad_trial_1 = 0.5 * E_grad_trial[0][2] + 0.5 * E_grad_trial[2][0];
                                    double tmpcse_symgrad_trial_2 = 0.5 * E_grad_trial[1][2] + 0.5 * E_grad_trial[2][1];
                                    grad_u[0][0] =
                                        1.0 * E_grad_trial[0][0] * src_local_hex[dimj][w][node_idx] + grad_u[0][0];
                                    grad_u[0][1] =
                                        tmpcse_symgrad_trial_0 * src_local_hex[dimj][w][node_idx] + grad_u[0][1];
                                    grad_u[0][2] =
                                        tmpcse_symgrad_trial_1 * src_local_hex[dimj][w][node_idx] + grad_u[0][2];
                                    grad_u[1][0] =
                                        tmpcse_symgrad_trial_0 * src_local_hex[dimj][w][node_idx] + grad_u[1][0];
                                    grad_u[1][1] =
                                        1.0 * E_grad_trial[1][1] * src_local_hex[dimj][w][node_idx] + grad_u[1][1];
                                    grad_u[1][2] =
                                        tmpcse_symgrad_trial_2 * src_local_hex[dimj][w][node_idx] + grad_u[1][2];
                                    grad_u[2][0] =
                                        tmpcse_symgrad_trial_1 * src_local_hex[dimj][w][node_idx] + grad_u[2][0];
                                    grad_u[2][1] =
                                        tmpcse_symgrad_trial_2 * src_local_hex[dimj][w][node_idx] + grad_u[2][1];
                                    grad_u[2][2] =
                                        1.0 * E_grad_trial[2][2] * src_local_hex[dimj][w][node_idx] + grad_u[2][2];
                                    div_u = div_u + E_grad_trial[dimj][dimj] * src_local_hex[dimj][w][node_idx];
                                };
                                for ( node_idx = cmb_shift; node_idx < 6 - surface_shift; node_idx += 1 )
                                {
                                    double E_grad_test[3][3]     = { 0 };
                                    E_grad_test[0][dimi]         = scalar_grad[node_idx][0];
                                    E_grad_test[1][dimi]         = scalar_grad[node_idx][1];
                                    E_grad_test[2][dimi]         = scalar_grad[node_idx][2];
                                    double tmpcse_symgrad_test_0 = 0.5 * E_grad_test[0][1] + 0.5 * E_grad_test[1][0];
                                    double tmpcse_symgrad_test_1 = 0.5 * E_grad_test[0][2] + 0.5 * E_grad_test[2][0];
                                    double tmpcse_symgrad_test_2 = 0.5 * E_grad_test[1][2] + 0.5 * E_grad_test[2][1];
                                    double tmpcse_pairing_0      = 2 * tmpcse_symgrad_test_0;
                                    double tmpcse_pairing_1      = 2 * tmpcse_symgrad_test_1;
                                    double tmpcse_pairing_2      = 2 * tmpcse_symgrad_test_2;
                                    dst_array[dimi][w][node_idx] =
                                        k_eval *
                                            ( -0.66666666666666663 * div_u * E_grad_test[dimi][dimi] +
                                              tmpcse_pairing_0 * grad_u[0][1] + tmpcse_pairing_0 * grad_u[1][0] +
                                              tmpcse_pairing_1 * grad_u[0][2] + tmpcse_pairing_1 * grad_u[2][0] +
                                              tmpcse_pairing_2 * grad_u[1][2] + tmpcse_pairing_2 * grad_u[2][1] +
                                              2.0 * E_grad_test[0][0] * grad_u[0][0] +
                                              2.0 * E_grad_test[1][1] * grad_u[1][1] +
                                              2.0 * E_grad_test[2][2] * grad_u[2][2] ) *
                                            fabs( J_det ) * qw_array[q] +
                                        dst_array[dimi][w][node_idx];
                                };
                            };
                            if ( dimi == dimj &&
                                 ( diagonal_ || treat_boundary_ && ( r_cell == 0 || r_cell + 1 == max_rad ) ) )
                            {
                                int node_idx;
                                for ( node_idx = surface_shift; node_idx < 6 - cmb_shift; node_idx += 1 )
                                {
                                    double E_grad_test[3][3] = { 0 };
                                    E_grad_test[0][dimj]     = scalar_grad[node_idx][0];
                                    E_grad_test[1][dimj]     = scalar_grad[node_idx][1];
                                    E_grad_test[2][dimj]     = scalar_grad[node_idx][2];
 
                                    double grad_u_diag[3][3]     = { 0 };
                                    double tmpcse_symgrad_test_0 = 0.5 * E_grad_test[0][1] + 0.5 * E_grad_test[1][0];
                                    double tmpcse_symgrad_test_1 = 0.5 * E_grad_test[0][2] + 0.5 * E_grad_test[2][0];
                                    double tmpcse_symgrad_test_2 = 0.5 * E_grad_test[1][2] + 0.5 * E_grad_test[2][1];
                                    grad_u_diag[0][0] = 1.0 * E_grad_test[0][0] * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[0][1] = tmpcse_symgrad_test_0 * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[0][2] = tmpcse_symgrad_test_1 * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[1][0] = tmpcse_symgrad_test_0 * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[1][1] = 1.0 * E_grad_test[1][1] * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[1][2] = tmpcse_symgrad_test_2 * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[2][0] = tmpcse_symgrad_test_1 * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[2][1] = tmpcse_symgrad_test_2 * src_local_hex[dimj][w][node_idx];
                                    grad_u_diag[2][2] = 1.0 * E_grad_test[2][2] * src_local_hex[dimj][w][node_idx];
                                    double tmpcse_pairing_0 = 4 * src_local_hex[dimj][w][node_idx];
                                    double tmpcse_pairing_1 = 2.0 * src_local_hex[dimj][w][node_idx];
                                    dst_array[dimi][w][node_idx] =
                                        k_eval *
                                            ( tmpcse_pairing_0 * pow( tmpcse_symgrad_test_0, 2 ) +
                                              tmpcse_pairing_0 * pow( tmpcse_symgrad_test_1, 2 ) +
                                              tmpcse_pairing_0 * pow( tmpcse_symgrad_test_2, 2 ) +
                                              tmpcse_pairing_1 * pow( E_grad_test[0][0], 2 ) +
                                              tmpcse_pairing_1 * pow( E_grad_test[1][1], 2 ) +
                                              tmpcse_pairing_1 * pow( E_grad_test[2][2], 2 ) -
                                              0.66666666666666663 * E_grad_test[dimi][dimi] * E_grad_test[dimj][dimj] *
                                                  src_local_hex[dimj][w][node_idx] ) *
                                            fabs( J_det ) * qw_array[q] +
                                        dst_array[dimi][w][node_idx];
                                };
                            };
                        };
                    };
                };
            };
            int dim_add;
            for ( dim_add = 0; dim_add < 3; dim_add += 1 )
            {
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell, y_cell, r_cell, dim_add ), dst_array[dim_add][0][0] );
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell + 1, y_cell, r_cell, dim_add ),
                    dst_array[dim_add][0][1] + dst_array[dim_add][1][2] );
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell, y_cell + 1, r_cell, dim_add ),
                    dst_array[dim_add][0][2] + dst_array[dim_add][1][1] );
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell, y_cell, r_cell + 1, dim_add ), dst_array[dim_add][0][3] );
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell + 1, y_cell, r_cell + 1, dim_add ),
                    dst_array[dim_add][0][4] + dst_array[dim_add][1][5] );
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell, y_cell + 1, r_cell + 1, dim_add ),
                    dst_array[dim_add][0][5] + dst_array[dim_add][1][4] );
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell + 1, y_cell + 1, r_cell, dim_add ), dst_array[dim_add][1][0] );
                Kokkos::atomic_add(
                    &dst_( local_subdomain_id, x_cell + 1, y_cell + 1, r_cell + 1, dim_add ),
                    dst_array[dim_add][1][3] );
            };
        }
        else
        {
            // load LMatrix for both local wedges
            A[0] = lmatrices_( local_subdomain_id, x_cell, y_cell, r_cell, 0 );
            A[1] = lmatrices_( local_subdomain_id, x_cell, y_cell, r_cell, 1 );
        }
 
        if ( diagonal_ )
        {
            A[0] = A[0].diagonal();
            A[1] = A[1].diagonal();
        }
 
        if ( storeLMatrices_ )
        {
            // write local matrices to mem
            lmatrices_( local_subdomain_id, x_cell, y_cell, r_cell, 0 ) = A[0];
            lmatrices_( local_subdomain_id, x_cell, y_cell, r_cell, 1 ) = A[1];
        }
        else
        {
            dense::Vec< ScalarT, 18 > src[num_wedges_per_hex_cell];
            for ( int dimj = 0; dimj < 3; dimj++ )
            {
                dense::Vec< ScalarT, 6 > src_d[num_wedges_per_hex_cell];
                extract_local_wedge_vector_coefficients(
                    src_d, local_subdomain_id, x_cell, y_cell, r_cell, dimj, src_ );
 
                for ( int wedge = 0; wedge < num_wedges_per_hex_cell; wedge++ )
                {
                    for ( int i = 0; i < num_nodes_per_wedge; i++ )
                    {
                        src[wedge]( dimj * num_nodes_per_wedge + i ) = src_d[wedge]( i );
                    }
                }
            }
            //extract_local_wedge_vector_coefficients( src, local_subdomain_id, x_cell, y_cell, r_cell, dimj, src_ );
 
            dense::Vec< ScalarT, 18 > dst[num_wedges_per_hex_cell];
 
            dst[0] = A[0] * src[0];
            dst[1] = A[1] * src[1];
 
            //atomically_add_local_wedge_vector_coefficients( dst_, local_subdomain_id, x_cell, y_cell, r_cell, dimi, dst );
            for ( int dimi = 0; dimi < 3; dimi++ )
            {
                dense::Vec< ScalarT, 6 > dst_d[num_wedges_per_hex_cell];
                dst_d[0] = dst[0].template slice< 6 >( dimi * num_nodes_per_wedge );
                dst_d[1] = dst[1].template slice< 6 >( dimi * num_nodes_per_wedge );
 
                atomically_add_local_wedge_vector_coefficients(
                    dst_, local_subdomain_id, x_cell, y_cell, r_cell, dimi, dst_d );
            }
        }
    }
};
 
static_assert( linalg::OperatorLike< EpsilonDivDivKerngenV01Initial< float > > );
static_assert( linalg::OperatorLike< EpsilonDivDivKerngenV01Initial< double > > );
 
} // namespace terra::fe::wedge::operators::shell::epsdivdiv_history