#ifndef INTEGRATION_H
#define INTEGRATION_H

#include <stdint.h>
#include <stdbool.h>

#include "util/interval.h"
#include "util/region.h"
#include "util/vec3f.h"
#include "util/mat3f.h"
#include "util/plane.h"

typedef struct IntegrationPoint_{
    Vec3f p; //integration point
    float w; //integration weight
    float d; //distance value
    Vec3f grad; //gradient value
} IntegrationPoint;

Vec3f edge_vec(const int e);
float asdf_volume(const ASDF* const asdf);
float asdf_face_area(const ASDF* const asdf, const int f);
Vec3f asdf_vertex_point(const ASDF* const asdf,const int v);
float asdf_cell_dimension(const ASDF* const asdf, const uint8_t n);

float bilinear_interp(float* v, float s1, float s2);
float trilinear_interp(float* v, float s1, float s2, float s3);

Vec3f asdf_interior_point(const ASDF* const asdf, 
	const float wx, const float wy, const float wz);
Vec3f asdf_interior_point_wrt_face(const ASDF* const asdf, 
	const int f, const float w1, const float w2, const float w3);
Vec3f asdf_interior_point_wrt_edge(const ASDF* asdf,const int e, const float w1, const float w2, const float w3);
Vec3f asdf_interior_point_wrt_vertex(const ASDF* asdf,const int v, const float w1, const float w2, const float w3);

Vec3f asdf_internal_ray_wrt_vertex(const ASDF* const asdf, const int v, const float w1, const float w2);
Vec3f adsf_internal_ray_wrt_edge(const ASDF* const asdf, const int e, const float w);

int asdf_find_interior_face(const ASDF* const asdf);
int asdf_find_interior_edge(const ASDF* const asdf);
int asdf_find_interior_vertex(const ASDF* const asdf);

float asdf_vertex_distance(const ASDF* const asdf, const uint8_t v0, const uint8_t v1);


int asdf_integrate_cartesian(const ASDF* const asdf, const ASDF* const boundary, const int N, IntegrationPoint* const result);
int asdf_integrate_wrt_face(const ASDF* const leaf, const ASDF* const boundary, const int N, IntegrationPoint* const result, const int f);
int asdf_integrate_wrt_edge(const ASDF* const leaf, const ASDF* const boundary, const int N, IntegrationPoint* const result, const int e);
int asdf_integrate_wrt_vertex(const ASDF* const leaf, const ASDF* const boundary, const int N, IntegrationPoint* const result, const int vert);

int asdf_integration_points_weights(
	const ASDF* const asdf, const ASDF* const boundary, const int degree,
	IntegrationPoint* const result);

float asdf_integrate_unity(const ASDF* const asdf, const int N);
int asdf_construct_integration_filled(const ASDF* const asdf, const int N, float* const result, int num_pts);
int asdf_construct_integration_leaf(const ASDF* const asdf, const int N, float* const result, int num_pts, int depth, int desired_depth);
int write_integration_points(const ASDF* const asdf, const int N, const int desired_depth);

Vec3f integrate_bodyforce_bspline(const ASDF* const asdf,const ASDF* const boundary, 
 int N, const Vec3f F, const Vec3f X0, const float h);
Vec3f integrate_bodyforce_bspline_focused(const ASDF* const asdf,const ASDF* const boundary, 
 int N, const Vec3f F, const Vec3f X0, const float h, const int n);

Mat3f integrate_stiffness_bspline(const ASDF* const asdf, const ASDF* const boundary,
 int N, const Vec3f X0, const Vec3f X1, const float h, 
 const float mu, const float l);
Mat3f integrate_stiffness_bspline_focused(const ASDF* const asdf, const ASDF* const boundary,
 int N, const Vec3f X0, const Vec3f X1, const float h, 
 const float mu, const float l, const int n);


#endif