//#####################################################################
// Copyright 2012, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#ifndef __Quadrature_Integration__
#define __Quadrature_Integration__
#include <Core/Vectors/VECTOR.h>
#include <functional>
namespace PhysBAM{

struct QUAD_RULE_1 {double w,u;};
struct QUAD_RULE_2 {double w,u,v;};

struct QUADRATURE_RULE_SEG
{
    bool use_center;
    double center_weight;
    int n;
    QUAD_RULE_1 r[3];
};

extern const QUADRATURE_RULE_SEG quadrature_rule_seg[12];

struct QUADRATURE_RULE_TRI
{
    bool use_center;
    double center_weight;
    int n3;
    QUAD_RULE_1 r3[3];
    int n6;
    QUAD_RULE_2 r6[1];
};

extern const QUADRATURE_RULE_TRI quadrature_rule_tri[9];

struct QUADRATURE_RULE_TET
{
    bool use_center;
    double center_weight;
    int n4;
    QUAD_RULE_1 r4[3];
    int n6;
    QUAD_RULE_1 r6[3];
    int n12;
    QUAD_RULE_2 r12[1];
};

extern const QUADRATURE_RULE_TET quadrature_rule_tet[7];

struct QUADRATURE_RULE_SQ
{
    bool use_center;
    double center_weight;
    int na;
    QUAD_RULE_1 ra[2];
    int nb;
    QUAD_RULE_1 rb[2];
    int n8;
    QUAD_RULE_2 r8[2];
};

extern const QUADRATURE_RULE_SQ quadrature_rule_sq[12];

struct QUADRATURE_RULE_CUBE
{
    bool use_center;
    double center_weight;
    int n6;
    QUAD_RULE_1 r6[2];
    int n8;
    QUAD_RULE_1 r8[2];
    int n12;
    QUAD_RULE_1 r12[2];
    int n24;
    QUAD_RULE_2 r24[2];
};

extern const QUADRATURE_RULE_CUBE quadrature_rule_cube[8];

// NOTE: order = total polynomial order that is exactly integrated.

// Returns integrated average over a segment, triangle, or tetrahedron.  The
// function f is called with the barycentric coordinates of the sample points.
// To get the integral, multiple the return value by the size of the geometrical
// element.
template<class R,class T> R Quadrature_Integration(int order, std::function<R(VECTOR<T,2> w)> f)
{
    assert(order<sizeof(quadrature_rule_seg)/sizeof(*quadrature_rule_seg));
    R r{};
    const QUADRATURE_RULE_SEG& rule=quadrature_rule_seg[order];
    if(rule.use_center) r+=(T)rule.center_weight*f({(T).5,(T).5});
    for(int i=0;i<rule.n && i<3;i++){
        auto [w_,a_]=rule.r[i];
        T w=w_,a=a_;
        r+=w*(f({1-a,a})+f({a,1-a}));}
    return r;
}

template<class R,class T> R Quadrature_Integration(int order, std::function<R(VECTOR<T,3> w)> f)
{
    assert(order<sizeof(quadrature_rule_tri)/sizeof(*quadrature_rule_tri));
    R r=R();
    const QUADRATURE_RULE_TRI& rule=quadrature_rule_tri[order];
    if(rule.use_center) r+=(T)rule.center_weight*f({(T)1/3,(T)1/3,(T)1/3});
    for(int i=0;i<rule.n3;i++){
        auto [w_,a_]=rule.r3[i];
        T w=w_,a=a_;
        r+=w*(f({1-2*a,a,a})+f({a,1-2*a,a})+f({a,a,1-2*a}));}
    for(int i=0;i<rule.n6;i++){
        auto [w_,a_,b_]=rule.r6[i];
        T w=w_,a=a_,b=b_;
        R s{};
        s+=f({1-a-b,a,b});
        s+=f({1-a-b,b,a});
        s+=f({a,1-a-b,b});
        s+=f({b,1-a-b,a});
        s+=f({a,b,1-a-b});
        s+=f({b,a,1-a-b});
        r+=w*s;}
    return r;
}

template<class R,class T> R Quadrature_Integration(int order, std::function<R(VECTOR<T,4> w)> f)
{
    assert(order<sizeof(quadrature_rule_tet)/sizeof(*quadrature_rule_tet));
    R r=R();
    const QUADRATURE_RULE_TET& rule=quadrature_rule_tet[order];
    if(rule.use_center) r+=(T)rule.center_weight*f({(T)1/4,(T)1/4,(T)1/4,(T)1/4});
    for(int i=0;i<rule.n4;i++){
        auto [w_,a_]=rule.r4[i];
        T w=w_,a=a_,b=1-3*a;
        r+=w*(f({b,a,a,a})+f({a,b,a,a})+f({a,a,b,a})+f({a,a,a,b}));}
    for(int i=0;i<rule.n6;i++){
        auto [w_,a_]=rule.r6[i];
        T w=w_,a=a_,b=(T).5-a;
        R s{};
        s+=f({a,a,b,b});
        s+=f({a,b,a,b});
        s+=f({a,b,b,a});
        s+=f({b,a,a,b});
        s+=f({b,a,b,a});
        s+=f({b,b,a,a});
        r+=w*s;}
    for(int i=0;i<rule.n12;i++){
        auto [w_,a_,b_]=rule.r12[i];
        T w=w_,a=a_,b=b_,c=1-2*a-b;
        R s{};
        s+=f({a,a,b,c});
        s+=f({a,a,c,b});
        s+=f({a,b,a,c});
        s+=f({a,c,a,b});
        s+=f({a,b,c,a});
        s+=f({a,c,b,a});
        s+=f({b,a,a,c});
        s+=f({c,a,a,b});
        s+=f({b,a,c,a});
        s+=f({c,a,b,a});
        s+=f({b,c,a,a});
        s+=f({c,b,a,a});
        r+=w*s;}
    return r;
}

// Returns integrated average over a segment, square, or cube.  The function f
// is called with the coordinates (0..1) of the unit box.  To get the integral,
// multiple the return value by the size of the geometrical element.
template<class R,class T> R Quadrature_Integration_Box(int order, std::function<R(VECTOR<T,1> w)> f)
{
    assert(order<sizeof(quadrature_rule_seg)/sizeof(*quadrature_rule_seg));
    R r{};
    const QUADRATURE_RULE_SEG& rule=quadrature_rule_seg[order];
    if(rule.use_center) r+=(T)rule.center_weight*f({(T).5,(T).5});
    for(int i=0;i<rule.n && i<3;i++){
        auto [w_,a_]=rule.r[i];
        T w=w_,a=a_;
        r+=w*(f({a})+f({1-a}));}
    return r;
}
template<class R,class T> R Quadrature_Integration_Box(int order, std::function<R(VECTOR<T,2> w)> f)
{
    assert(order<sizeof(quadrature_rule_sq)/sizeof(*quadrature_rule_sq));
    R r=R();
    const QUADRATURE_RULE_SQ& rule=quadrature_rule_sq[order];
    T h=(T).5;
    if(rule.use_center) r+=(T)rule.center_weight*f({h,h});
    for(int i=0;i<rule.na;i++){
        auto [w_,a_]=rule.ra[i];
        T w=w_,a=a_,b=1-a;
        r+=w*(f({a,a})+f({b,a})+f({a,b})+f({b,b}));}
    for(int i=0;i<rule.nb;i++){
        auto [w_,a_]=rule.rb[i];
        T w=w_,a=a_,b=1-a;
        r+=w*(f({h,a})+f({h,b})+f({a,h})+f({b,h}));}
    for(int i=0;i<rule.n8;i++){
        auto [w_,a_,b_]=rule.r8[i];
        T w=w_,a=a_,b=b_,c=1-a,d=1-b;
        R s{};
        s+=f({a,b});
        s+=f({c,b});
        s+=f({a,d});
        s+=f({c,d});
        s+=f({b,a});
        s+=f({b,c});
        s+=f({d,a});
        s+=f({d,c});
        r+=w*s;}
    return r;
}
template<class R,class T> R Quadrature_Integration_Box(int order, std::function<R(VECTOR<T,3> w)> f)
{
    assert(order<sizeof(quadrature_rule_cube)/sizeof(*quadrature_rule_cube));
    R r=R();
    const QUADRATURE_RULE_CUBE& rule=quadrature_rule_cube[order];
    T h=(T).5;
    if(rule.use_center) r+=(T)rule.center_weight*f({h,h,h});
    for(int i=0;i<rule.n6;i++){
        auto [w_,a_]=rule.r6[i];
        T w=w_,a=a_,b=1-a;
        R s{};
        s+=f({a,h,h});
        s+=f({h,a,h});
        s+=f({h,h,a});
        s+=f({b,h,h});
        s+=f({h,b,h});
        s+=f({h,h,b});
        r+=w*s;}
    for(int i=0;i<rule.n8;i++){
        auto [w_,a_]=rule.r8[i];
        T w=w_,a=a_,b=1-a;
        R s{};
        s+=f({a,a,a});
        s+=f({b,a,a});
        s+=f({a,b,a});
        s+=f({b,b,a});
        s+=f({a,a,b});
        s+=f({b,a,b});
        s+=f({a,b,b});
        s+=f({b,b,b});
        r+=w*s;}
    for(int i=0;i<rule.n12;i++){
        auto [w_,a_]=rule.r12[i];
        T w=w_,a=a_,b=1-a;
        R s{};
        s+=f({a,a,h});
        s+=f({a,h,a});
        s+=f({h,a,a});
        s+=f({b,a,h});
        s+=f({b,h,a});
        s+=f({h,b,a});
        s+=f({a,b,h});
        s+=f({a,h,b});
        s+=f({h,a,b});
        s+=f({b,b,h});
        s+=f({b,h,b});
        s+=f({h,b,b});
        r+=w*s;}
    for(int i=0;i<rule.n24;i++){
        auto [w_,a_,b_]=rule.r24[i];
        T w=w_,a=a_,b=b_,c=1-a,d=1-b;
        R s{};
        s+=f({a,a,b});
        s+=f({a,b,a});
        s+=f({b,a,a});
        s+=f({c,a,b});
        s+=f({c,b,a});
        s+=f({b,c,a});
        s+=f({a,c,b});
        s+=f({a,b,c});
        s+=f({b,a,c});
        s+=f({c,c,b});
        s+=f({c,b,c});
        s+=f({b,c,c});
        s+=f({a,a,d});
        s+=f({a,d,a});
        s+=f({d,a,a});
        s+=f({c,a,d});
        s+=f({c,d,a});
        s+=f({d,c,a});
        s+=f({a,c,d});
        s+=f({a,d,c});
        s+=f({d,a,c});
        s+=f({c,c,d});
        s+=f({c,d,c});
        s+=f({d,c,c});
        r+=w*s;}
    return r;
}

}
#endif
