integratorRHS.h

/**********************************************************************************
* Author: Christoph Pflaum, Riccarda Scherner-Griesshammer, Tim Rheinfels
 *                Department Informatik Lehrstuhl 10 - Systemsimulation
 *                Friedrich-Alexander Universität Erlangen-Nürnberg
 *
*********************************************/
#ifndef INTERATORRHS_H
#define INTERATORRHS_H
 
#include <functional>
#include <limits>
 
#include "integral.hpp"
#include "../interfaceMatrices.h"
 
template <typename F, size_t N>
    class IntegratorRHS : public InterfaceLocalStiffnessMatrices<N>
{
 
    private:
 
        class TermRHS
        {
                
            private:
                
                std::array<double,N> h_u;              
                std::array<double,N> h_v;              
                std::array<double,N> p_u;              
                std::array<double,N> p_v;              
                std::array<int,N> d_u;              
                std::array<int,N> d_v;              
                F m_variable_coefficient;              
                double m_gamma;                        
 
                
            public:
 
                inline TermRHS( const std::array<double, N>& h_u_, const std::array<double, N>&  h_v_, const std::array<double, N>& p_u_, const std::array<double, N>&  p_v_,const std::array<int, N>& d_u_, const std::array<int, N>&  d_v_, const F& variable_coefficient, double gamma )
                    : h_u(h_u_),
                      h_v(h_v_),
                      p_u(p_u_),
                      p_v(p_v_),
                      d_u(d_u_),
                      d_v(d_v_),
                      m_variable_coefficient(variable_coefficient),
                      m_gamma(gamma)
                {
                    
                }
 
                inline double operator() (std::array<double, N>& x ) const
                {
 
 
                    
                    double value = 1.0;
                    
                    for( size_t i = 0; i < N; i++ )
                    {
                        // Multiply by the first base function
 
                        if(d_u[i]!=0){
                            double suppleft=p_u[i]-h_u[i];
                            double suppright=p_u[i]+h_u[i];
 
                                if(suppleft <= x[i] && x[i] <= suppright){
                                    double slope = 1.0/h_u[i];
                                    if(x[i]<p_u[i]){
                                        double c = -(p_u[i]-h_u[i])/h_u[i];
                                        value *= (slope*x[i]+c);
                                    } else{
                                        double c = + (p_u[i]+h_u[i])/h_u[i];
                                        value *= (-1.0*slope*x[i]+c);
                                    }
                                }else{
                                    value = 0.0;
                                }
                        }else{
                            if(p_u[i]==0.0)
                                value *= -x[i]+1.0;
                            if(p_u[i]==1.0)
                                value *= x[i];
                        }
                        if(d_v[i]!=0) {
                            double suppleft = p_v[i] - h_v[i];
                            double suppright = p_v[i] + h_v[i];
                            if (suppleft <= x[i] && x[i] <= suppright) {
                                double slope = 1.0 / h_v[i];
                                if (x[i] < p_v[i]) {
 
                                    double c = -(p_v[i] - h_v[i]) / h_v[i];
                                    value *= (slope * x[i] + c);
                                } else {
 
                                    double c = +(p_v[i] + h_v[i]) / h_v[i];
                                    value *= (-1.0 * slope * x[i] + c);
                                }
                            } else {
                                value = 0.0;
                            }
                        }else{
 
                            if(p_v[i]==0.0){
                                value *= -x[i]+1.0;
                            }
 
                            if(p_v[i]==1.0){
                                value *= x[i];
                            }
 
                        }
 
                    }
 
 
                    
                    // 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:
        
        IntegratorRHS( 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 integration(std::array<double, N> l_b, std::array<double, N> r_b, std::array<double, N> p_u, std::array<double, N> p_v, std::array<double, N> h_u, std::array<double, N> h_v,
                           array<int, N> d_u, std::array<int, N> d_v) const
        {
 
            TermRHS ht(h_u, h_v, p_u,p_v,d_u,d_v,m_variable_coefficient, m_gamma );
 
 
 
 
 
            // Integrate
            return integral<double, TermRHS, N>( ht, l_b, r_b, m_epsilon, m_depth_min, m_depth_max, m_gamma );
 
        }
 
        double stencil_integration(double p_left[N], double p_right[N],
                                   BasisFunctionType u[N], BasisFunctionType v[N]) const {return 0;};
 
 
 
};
 
#endif