helmIntegrator.h

/**********************************************************************************
* Author: Christoph Pflaum, Riccarda Scherner-Griesshammer, Tim Rheinfels
 *                Department Informatik Lehrstuhl 10 - Systemsimulation
 *                Friedrich-Alexander Universität Erlangen-Nürnberg
 *
*********************************************/
 
#ifndef INTERATORHELM_H
#define INTERATORHELM_H
 
#include <cmath>
#include <functional>
#include <limits>
 
#include "integral.hpp"
 
template <typename F, size_t N>
    class IntegratorHelm : public InterfaceLocalStiffnessMatrices<N>
{
 
    private:
 
        class HelmholtzTerm
        {
                
            private:
                
                std::array<double, N> m_p_left;        
                std::array<double, N> m_p_right;       
                std::array<BasisFunctionType, N> m_u;  
                std::array<BasisFunctionType, N> m_v;  
                F m_variable_coefficient;              
                double m_gamma;                        
 
                
            public:
 
                inline HelmholtzTerm( const std::array<double, N>& p_left, const std::array<double, N>&  p_right,
                                      const std::array<BasisFunctionType, N>& u, const std::array<BasisFunctionType, N>& v,
                                      const F& variable_coefficient, double gamma )
                    : m_p_left(p_left),
                      m_p_right(p_right),
                      m_u(u),
                      m_v(v),
                      m_variable_coefficient(variable_coefficient),
                      m_gamma(gamma)
                {
                    
                }
 
                inline double operator() ( const std::array<double, N>& x ) const
                {
                    
                    double value = 1.0;
                    
                    for( size_t i = 0; i < N; ++i )
                    {
                        
                        // === Muliply by the first base function ===
                        switch( m_u[i] )
                        {
                            case BasisFunctionType::leftBasis:
                                // p(x) = (x-b)/(a-b)
                                value *= (x[i] - m_p_right[i]) / (m_p_left[i] - m_p_right[i]);
                                break;
                                
                            case BasisFunctionType::rightBasis:
                                // p(x) = (x-a)/(b-a)
                                value *= (x[i] - m_p_left[i]) / (m_p_right[i] - m_p_left[i]);
                                break;
                                
                            case BasisFunctionType::gradLeftBasis:
                                // p'(x) =  1/(b-a)
                                value /= (m_p_left[i] - m_p_right[i]);
                                break;
                                
                            case BasisFunctionType::gradRightBasis:
                                // p'(x) = 1/(a-b)
                                value /= (m_p_right[i] - m_p_left[i]);
                                break;
                        }
                        
                        // === Muliply by the second base function ===
                        switch( m_v[i] )
                        {
                            case BasisFunctionType::leftBasis:
                                // p(x) = (x-b)/(a-b)
                                value *= (x[i] - m_p_right[i]) / (m_p_left[i] - m_p_right[i]);
                                break;
                                
                            case BasisFunctionType::rightBasis:
                                // p(x) = (x-a)/(b-a)
                                value *= (x[i] - m_p_left[i]) / (m_p_right[i] - m_p_left[i]);
                                break;
                                
                            case BasisFunctionType::gradLeftBasis:
                                // p'(x) =  1/(b-a)
                                value /= (m_p_left[i] - m_p_right[i]);
                                break;
                                
                            case BasisFunctionType::gradRightBasis:
                                // p'(x) = 1/(a-b)
                                value /= (m_p_right[i] - m_p_left[i]);
                                break;
                        }
 
                    }
                    
                    // Evaluate the variable coefficient replacing floating point inf by gamma
                    double c_eval = m_variable_coefficient(x);
                    
                    if( (std::abs(c_eval) > m_gamma) || (std::fpclassify(c_eval) == FP_INFINITE ))
                    {
                        c_eval = (std::signbit(c_eval) ? (-m_gamma) : m_gamma);
                    }
                    
                    // Mulitply by the variable coefficient
                    value *= c_eval;
                    
                    // Return Helmholtz function with variable coefficient evaluated at x
                    return value;
 
                }
                    
        };
        
        
        F m_variable_coefficient;  
        double m_epsilon;          
        size_t m_depth_min;        
        size_t m_depth_max;        
        double m_gamma;            
 
        
    public:
        
        IntegratorHelm( const F& c, double epsilon, size_t depth_min = 0, size_t depth_max = std::numeric_limits<size_t>::max(), double gamma = std::numeric_limits<double>::max() )
            : m_variable_coefficient(c),
              m_epsilon(epsilon),
              m_depth_min(depth_min),
              m_depth_max(depth_max),
              m_gamma(gamma)
        {
 
            /* Nothing to do here... */
 
        }
 
        double stencil_integration(double p_left[N], double p_right[N], BasisFunctionType u[N], BasisFunctionType v[N]) const
        {
 
            // Duplicate input data as std::array
            std::array<double, N> p_left_array;
            std::array<double, N> p_right_array;
 
            std::array<BasisFunctionType, N> u_array;
            std::array<BasisFunctionType, N> v_array;
            
            for( size_t i = 0; i < N; ++i )
            {
 
                p_left_array[i] = p_left[i];
                p_right_array[i] = p_right[i];
 
                u_array[i] = u[i];
                v_array[i] = v[i];
                
            }
            
            // Create local helmholtz term
            HelmholtzTerm ht(p_left_array, p_right_array, u_array, v_array, m_variable_coefficient, m_gamma);
 
            
            // Integrate
            return integral<double, HelmholtzTerm, N>( ht, p_left_array, p_right_array, m_epsilon, m_depth_min, m_depth_max, m_gamma );
 
        }
 
};
 
#endif