/* 
 * File:   MaterialModel.cpp
 * Author: David
 * 
 * Created on 24. March 2014
 */

#include "MaterialModelImpl.h"

#include <cmath>
#include <cstdio>
using namespace std;

#include <openvdb/openvdb.h>
#include <openvdb/tools/Interpolation.h>
namespace vdb = openvdb::v2_2_0;

namespace FractureSim{
	MaterialModel* createMaterialModel(string spec, double youngsMod, double poissonsRatio, double density, double strength, double toughness, double compress){
		if(spec.compare("?")==0){ // print options
			printf(
				"\nthe following material model specifications are implemented:\n"
				"default ... homogeneous material\n"
				"'lay(nx,ny,nz,f,p,K)' ... layered material with normal (nx,ny,nz), frequency f, phase p and secondary toughness K\n"
				"'test' ... dummy material model (used for testing)\n" );
		}
		if(spec.compare(0,3,"lay")==0){ // format: "lay(nx,ny,nz,f,p,Kc_hard)" or "lay(nx,ny,nz,f,p,K_hard/snx,sny,snz,sf,sp,Sc_hard)"
			double nx,ny,nz,f,p,k;
			char c;
			int done = sscanf(spec.substr(3).c_str(), "(%lg,%lg,%lg,%lg,%lg,%lg%c", &nx,&ny,&nz,&f,&p,&k,&c);
			if(done == 7 && c == ')'){
				Eigen::Vector3d n(nx,ny,nz); n.normalize();
				return new LayeredMaterialModel( youngsMod, poissonsRatio,density, strength, compress, toughness, k, 2.0*M_PI*f, p, n);
			}else if(done == 7 && c == '/'){
                // also use varying strength
                double snx,sny,snz,sf,sp,s;
                int done = sscanf(spec.substr( spec.find_first_of('/') ).c_str(), 
                    "/%lg,%lg,%lg,%lg,%lg,%lg%c", &snx,&sny,&snz,&sf,&sp,&s,&c);
                if(done == 7 && c == ')'){
    				Eigen::Vector3d  n(nx,ny,nz);     n.normalize();
                    Eigen::Vector3d sn(snx,sny,snz); sn.normalize();
                    return new LayeredMaterialModel(
                        youngsMod, poissonsRatio,density, strength, compress, toughness,
                        k, 2.0*M_PI*f, p, n,    s, 2.0*M_PI*sf, sp, sn );
                }
            }
		}
		if(spec.compare("tube")==0){
			return new TubedMaterialModel(youngsMod, poissonsRatio,density, strength, compress, toughness, 2.0*toughness);
		}
		if(spec.compare("test")==0){
			return new DummyMaterialModel(youngsMod, poissonsRatio,density, strength,toughness, compress, 3.0);
		}
        if(spec.compare(0,3,"vdb")==0){ // format: "vdb(filename/toughnessGridName/strengthGridName,scaleFactor,voxelSize,useUDF)"
			std::string file, toughname="toughness", strengthname="strength";
			file=spec.substr(4, spec.find_first_of(',')-4);                     //printf("\n1)     file=\"%s\"",file.c_str());
			if(file.find_first_of('/')!=file.npos){
				toughname=file.substr(file.find_first_of('/')+1);               //printf("\n2)    tough=\"%s\"", toughname.c_str());
				file=file.substr(0,file.find_first_of('/'));                    //printf("\n2)     file=\"%s\"",file.c_str());
				strengthname=toughname;
			}
			if(toughname.find_first_of('/')!=toughname.npos){
				strengthname=toughname.substr(toughname.find_first_of('/')+1);  //printf("\n3) strength=\"%s\"", strengthname.c_str());
				toughname=toughname.substr(0,toughname.find_first_of('/'));     //printf("\n3)    tough=\"%s\"",toughname.c_str());
			}                                                                   //printf("\n4) remainder = \"%s\"", spec.substr(spec.find_first_of(',')+1).c_str());
			double scale, voxelSize; int readUDF; char check=0;
			sscanf( spec.substr(spec.find_first_of(',')+1).c_str(), "%lf,%lf,%d%c", &scale, &voxelSize, &readUDF, &check);
			if(check!=')'){
				printf("\nillegal definition of VDB material: %s, use "
					   "\"vdb(filename/toughnessGridName/strengthGridName,scaleFactor,voxelSize,useUDF)\"", spec.c_str());
			}else{
				printf("\n%% ... VDB material model: %s / %s / %s scale %.3lg, voxel size %.3lg, udf %d",
					file.c_str(), toughname.c_str(), strengthname.c_str(), scale, voxelSize, (readUDF!=0));
				return new VDBMaterialModel(
					youngsMod, poissonsRatio,density, strength,toughness, compress,
					file,toughname,strengthname,readUDF!=0, scale ,voxelSize);
			}
        }
        printf("\n%% ... using default material model");
		return new HomogeneousMaterialModel( youngsMod, poissonsRatio,density, strength,toughness, compress);
	}

	// DummyMaterialModel ==================================================================================
	//// "y-split" - doesn't have a proper gradient, testing how badly this affects results
	//double DummyMaterialModel::fractureToughness(const Eigen::Vector3d& x) const{
	//	return Kc* (1.0 + KcDelta* (x[1]>0.0?1.0:0.0));
	//}
	//void DummyMaterialModel::fractureToughnessGradient(Eigen::Vector3d& dKc, const Eigen::Vector3d& x) const{
	//	dKc.setZero(); double h=1e-2;
	//	dKc[1]=2.0*h*(fractureToughness(x+h*Eigen::Vector3d::UnitY())-fractureToughness(x-h*Eigen::Vector3d::UnitY()));
	//}

	// cos-"bubbles"
	double DummyMaterialModel::fractureToughness(const Eigen::Vector3d& x) const{
        const double w=15*M_PI;
		return Kc* (1.0 + KcDelta* (0.5 + 0.5*cos(w*x[0])*cos(w*x[1])*cos(w*x[2])));
	}
	void DummyMaterialModel::fractureToughnessGradient(Eigen::Vector3d& dKc, const Eigen::Vector3d& x) const{
    	const double w=15*M_PI;
		dKc[0]= -Kc*KcDelta* 0.5*w*sin(w*x[0])*cos(w*x[1])*cos(w*x[2]);
		dKc[1]= -Kc*KcDelta* 0.5*w*cos(w*x[0])*sin(w*x[1])*cos(w*x[2]);
		dKc[2]= -Kc*KcDelta* 0.5*w*cos(w*x[0])*cos(w*x[1])*sin(w*x[2]);
	}
	// =====================================================================================================

	// LayeredMaterialModel ==================================================================================
	double LayeredMaterialModel::tensileStrength(const Eigen::Vector3d& x) const {
		double w = 0.5 + 0.5*cos(sf * x.dot(sn) + sp); // w is in [0, 1]
//        printf("\n strength eval at (%+.2lf, %+.2lf, %+.2lf) is %.3lg (w=%.3lf)",
//            x[0],x[1],x[2] , (1.0-w)*Sc + w*ScHard ,w);
		return (1.0-w)*Sc + w*ScHard;
    }
    double LayeredMaterialModel::fractureToughness(const Eigen::Vector3d& x) const{
		double w = 0.5 + 0.5*cos(f * x.dot(n) + p); // w is in [0, 1]
		return (1.0-w)*Kc + w*KcHard;
	}
	void LayeredMaterialModel::fractureToughnessGradient(Eigen::Vector3d& dKc, const Eigen::Vector3d& x) const{
		dKc= (KcHard-Kc)* (-0.5*f)*n*sin( f * x.dot(n) + p );
	}
	// =====================================================================================================

	// TubedMaterialModel ==================================================================================
	double TubedMaterialModel::fractureToughness(const Eigen::Vector3d& x) const{
		double f1=15*M_PI, f2=30*M_PI, p1=0.0, p2=0.0;
		Eigen::Vector3d n1(1.0, 0.0, 0.0), n2(0.0, 0.0, 1.0);
		n1.normalize(); n2.normalize();

		double w = 0.5+0.5*cos(f1*x.dot(n1)+p1)*cos(f2*x.dot(n2)+p2); // w is in [0, 1]
		//printf("\n %.3le ", (1-w)*Kc + w*KcHard);
		return (1-w)*Kc + w*KcHard;
	}
	void TubedMaterialModel::fractureToughnessGradient(Eigen::Vector3d& dKc, const Eigen::Vector3d& x) const{
		double f1=15*M_PI, f2=30*M_PI, p1=0.0, p2=0.0;
		Eigen::Vector3d n1(1.0, 0.0, 0.0), n2(0.0, 0.0, 1.0);
		n1.normalize(); n2.normalize();

		dKc= (KcHard-Kc)*(
			-0.5*f1*n1*sin(f1*x.dot(n1)+p1)*cos(f2*x.dot(n2)+p2)
			-0.5*f2*n2*cos(f1*x.dot(n1)+p1)*sin(f2*x.dot(n2)+p2)
		);
		//printf(" grad( %.3le %.3le %.3le )",dKc[0],dKc[1],dKc[2]);
	}
	// =====================================================================================================

	// VDBMaterialModel ==================================================================================
	class VDBMaterialModel::vdbData{
    public:
        vdbData() : ux(1.0,0.0,0.0), uy(0.0,1.0,0.0), uz(0.0,0.0,1.0) {};
        vdb::FloatGrid::Ptr toughnessGrid, strengthGrid;
        vdb::math::Transform::Ptr xform;
        vdb::Vec3d ux,uy,uz;
		vdb::math::CoordBBox bbox;
    };
    void VDBMaterialModel::init() const{
        haveToughnessMap=false; haveStrengthMap=false;
        if(data==NULL) data = new vdbData();
        vdb::initialize();
        vdb::io::File inFile(vdbFile);
        if(!inFile.open()){
			printf("\n%% Failed to open file %s", vdbFile.c_str());
		}else{
			vdb::GridBase::Ptr baseGrid;
			if( inFile.hasGrid(kName) ){
				baseGrid = inFile.readGrid(kName);
				data->toughnessGrid=vdb::gridPtrCast<vdb::FloatGrid>(baseGrid);
				data->bbox = data->toughnessGrid->evalActiveVoxelBoundingBox();
				haveToughnessMap=true;
			}
			if( kName.compare(sName)) data->strengthGrid = data->toughnessGrid; // share the grid ptr when using the same data for both strength and toughness
			else if( inFile.hasGrid(sName) ){
				baseGrid = inFile.readGrid(sName);
				data->strengthGrid=vdb::gridPtrCast<vdb::FloatGrid>(baseGrid);
				data->bbox.expand( data->strengthGrid->evalActiveVoxelBoundingBox() );
				haveStrengthMap=true;
			}
			inFile.close();
			data->xform = vdb::math::Transform::createLinearTransform(voxelSize);
			//printf("\n%% VDB material model initialized%s%s\n",
			//	haveToughnessMap?", toughness map loaded":"",
			//	haveStrengthMap ? ", strength map loaded":""
			//);
		}
    }
    double VDBMaterialModel::tensileStrength(const Eigen::Vector3d& x) const{
        if(data==NULL) init();
        if(!haveStrengthMap) return Sc;
        vdb::tools::GridSampler<vdb::FloatGrid::ConstAccessor, vdb::tools::BoxSampler> sampler(
            data->strengthGrid->getConstAccessor(), *(data->xform)
        );
		vdb::Vec3d xi = toPeriodicIndexSpace<vdb::Vec3d>(x);
        float gridValue = sampler.isSample(xi) / data->strengthGrid->background();
        if(useUDF) gridValue=std::abs(gridValue);
        return Sc + scale*gridValue;

    }
	double VDBMaterialModel::fractureToughness(const Eigen::Vector3d& x) const{
        if(data==NULL) init();
        if(!haveToughnessMap) return Kc;
        vdb::tools::GridSampler<vdb::FloatGrid::ConstAccessor, vdb::tools::BoxSampler> sampler(
            data->toughnessGrid->getConstAccessor(), *(data->xform)
        );
		vdb::Vec3d xi = toPeriodicIndexSpace<vdb::Vec3d>(x);
        float gridValue = sampler.isSample(xi) / data->toughnessGrid->background();
        if(useUDF) gridValue=std::abs(gridValue);
		//printf("\n eval Kc at (%.3lf,%.3lf,%.3lf) df %.3lf Kc %.3le",x[0],x[1],x[2],gridValue,Kc + scale*gridValue);
        return Kc + scale*gridValue;
    }
	void VDBMaterialModel::fractureToughnessGradient(Eigen::Vector3d& dKc, const Eigen::Vector3d& x) const{
		if(data==NULL) init();
        if(!haveToughnessMap){
            dKc.setZero();
            return;
        }
        vdb::tools::GridSampler<vdb::FloatGrid::ConstAccessor, vdb::tools::BoxSampler> sampler(
            data->toughnessGrid->getConstAccessor(), *(data->xform)
        );
        float x0,x1,y0,y1,z0,z1, bg=data->toughnessGrid->background();
		vdb::Vec3d xi = toPeriodicIndexSpace<vdb::Vec3d>(x);
        x0=sampler.isSample(xi-0.25*data->ux) / bg;
        x1=sampler.isSample(xi+0.25*data->ux) / bg;
        y0=sampler.isSample(xi-0.25*data->uy) / bg;
        y1=sampler.isSample(xi+0.25*data->uy) / bg;
        z0=sampler.isSample(xi-0.25*data->uz) / bg;
        z1=sampler.isSample(xi+0.25*data->uz) / bg;
        if(useUDF){
            x0=std::abs(x0); x1=std::abs(x1);
            y0=std::abs(y0); y1=std::abs(y1);
            z0=std::abs(z0); z1=std::abs(z1);
		}
        dKc[0] = scale*(x1-x0)*2.0/voxelSize;
        dKc[1] = scale*(y1-y0)*2.0/voxelSize;
        dKc[2] = scale*(z1-z0)*2.0/voxelSize;
		//if(dKc.norm() > FLT_EPSILON){
		//	printf("\n x(%.3lf, %.3lf) grad_x %.3le",x0,x1,dKc[0]);
		//	printf("\n y(%.3lf, %.3lf) grad_y %.3le",y0,y1,dKc[1]);
		//	printf("\n z(%.3lf, %.3lf) grad_z %.3le",z0,z1,dKc[2]);
		//}
    }
	template <class Vec3>
	Vec3 VDBMaterialModel::toPeriodicIndexSpace(const Eigen::Vector3d& x) const{
		Vec3 xi = data->xform->worldToIndex(Vec3(x[0],x[1],x[2]));
		//printf("\n pb_map w(%.3lf,%.3lf,%.3lf) -> i(%.1lf,%.1lf,%.1lf)", x[0],x[1],x[2], xi[0],xi[1],xi[2]);
		vdb::math::BBox<Vec3> box ( data->bbox.min().asVec3d(), data->bbox.max().asVec3d() );
		if( box.isInside(xi) ) return xi;
		xi-=box.min(); // shift min to (0,0,0)
		//printf("\n--> s(%.1lf,%.1lf,%.1lf)", xi[0],xi[1],xi[2]);
		xi[0] = std::fmod( xi[0], box.extents()[0] ); if(xi[0]<0.0) xi[0]+=box.extents()[0];
		xi[1] = std::fmod( xi[1], box.extents()[1] ); if(xi[1]<0.0) xi[1]+=box.extents()[1];
		xi[2] = std::fmod( xi[2], box.extents()[2] ); if(xi[2]<0.0) xi[2]+=box.extents()[2];
		//printf("\n--> M(%.1lf,%.1lf,%.1lf), extents (%.1lf,%.1lf,%.1lf)", xi[0],xi[1],xi[2], box.extents()[0],box.extents()[1],box.extents()[2]);
		xi+=box.min(); // undo shift
		//printf("\n--> b(%.1lf,%.1lf,%.1lf) in box (%.0lf,%.0lf,%.0lf)-(%.0lf,%.0lf,%.0lf)", xi[0],xi[1],xi[2], box.min()[0],box.min()[1],box.min()[2], box.max()[0],box.max()[1],box.max()[2]);
		return xi;
	}
	// =====================================================================================================
}