// DATE : MAY 2018
// SUMMARY 	:
// USAGE 	: LGPL
// ORG		: LGAD, Université de Nice Sophia-Antipolis
// AUTHOR	: Olivier Pantz
// EMAIL	: olivier.pantz@polytechnique.org

/*
 This file is part of Freefem++
 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/

/*
 Large chunks of the code has been inspired by iovtk.cpp written by Jacques Morice 
*/

#include  <iostream>
#include <vector>
#include <fstream>
#include <sstream>
#include <boost/algorithm/string.hpp>
using namespace std;

#include "error.hpp"
#include "AFunction.hpp"
#include "rgraph.hpp"
#include "fem.hpp"
#include "FESpace.hpp" 
#include "MeshPoint.hpp"
#include "lgfem.hpp"
#include "lgmesh3.hpp"

class SVG_WriteMesh_Op : public E_F0mps
{
	private:
		Expression eTh;			// Mesh
		Expression filename;		// Name of output SVG File
		Expression u;			// Finite element (if given)
		static const int n_name_param=6; 			// Nb of optional parameters
	public:
		static basicAC_F0::name_and_type name_param[];
		Expression nargs[n_name_param];				// to store the optional parameters

	public:
		typedef long  Result;
		static ArrayOfaType  typeargs() { return  ArrayOfaType( atype<string *>(), atype<pmesh>(), true); }
		static  E_F0 * f(const basicAC_F0 & args) { return new SVG_WriteMesh_Op(args);}
		SVG_WriteMesh_Op(const basicAC_F0 &  args)
  		{
			args.SetNameParam(n_name_param,name_param,nargs);
			if ((args.size()!=2) and (args.size()!=3)) {
				cerr<<"The number of arguments for savesvg should be two or three."<<endl;
				ExecError("Incorect number of arguments for savesvg");
			}
			if(verbosity > 2) cout << "Write Mesh and Solutions in SVG Formats" << endl;
			if (BCastTo<string *>(args[0])) {
				filename = CastTo<string *>(args[0]);
			} else {
				cerr<<"The First argument of savesvg should be a string"<<endl;
				ExecError("Incorect type in savesvg. First argument should be a string");
			}
 			if (BCastTo<pmesh>(args[1])) {
				eTh= CastTo<pmesh>(args[1]);
			} else {
				cerr<<"The Second argument of savesvg should be a mesh."<<endl;
				ExecError("Incorect type in savesvg. Second argument should be a mash");
			}
			u=0;
			if(args.size()==3)
			{// A Finite element function is given in argument...
				if (BCastTo<double>(args[2])){
					u=to<double>(args[2]);
				} else {
					cerr<<"Incorrect third argument. Finite element expected."<<endl;
					ExecError("Incorect type in savesvg. Third argument could only be a Finite element function");				
				}
			} 
		
		}
		AnyType operator()(Stack stack)  const
		{
			int width=500;			// Width of the SV
			string stroke="black";		// stroke color
			string strokeEdge="red";		// stroke color
			int strokeWidth=0;		// stroke width
			int strokeWidthEdge=0;		// stroke width
			KN<double> default_rgba(8);	// default value for fill color... There is probably a smarter way to initialize it
			default_rgba[0]=0; default_rgba[1]=0; default_rgba[2]=0; default_rgba[3]=0;
			default_rgba[4]=0; default_rgba[5]=0; default_rgba[6]=0; default_rgba[7]=1;
			KN<double>* rgba=&default_rgba;	

			if(nargs[0])  width= GetAny<long>((*nargs[0])(stack));
			if(nargs[1])  strokeWidth= GetAny<long>((*nargs[1])(stack));
			if(nargs[2])  stroke= *(GetAny<string*>((*nargs[2])(stack)));
			if(nargs[3])  rgba= (GetAny<KN<double>*>((*nargs[3])(stack)));
			if(nargs[4])  strokeWidthEdge= GetAny<long>((*nargs[4])(stack));
			if(nargs[5])  strokeEdge= *(GetAny<string*>((*nargs[5])(stack)));
			if((*rgba).size()!=8){
				cerr << "rgba should be an array of 8 numbers"<< endl;
				ExecError("rgba should be an array of 8 numbers");
			};
			string * pffname= GetAny<string *>((*filename)(stack));
			Mesh * pTh= GetAny<Mesh *>((*eTh)(stack));
			ffassert(pTh);
			Mesh &Th=*pTh;
			FILE *fp = fopen( (*pffname).c_str(), "wb");
			if(!fp){
				cerr << "Unable to write in file " << (*pffname).c_str() << endl;
				ExecError("error in saving svg file");
			}
			// We compute the min/max of X and Y
			double minX,maxX,minY,maxY;
			{
				minX=Th.vertices[0].x;minY=Th.vertices[0].y;
				maxX=minX;maxY=minY;
				for(int iv=0;iv<Th.nv;iv++){
					double xv=Th.vertices[iv].x;				
					double yv=Th.vertices[iv].y;
					if(xv<minX) minX=xv;if(yv<minY) minY=yv;
					if(xv>maxX) maxX=xv;if(yv>maxY) maxY=yv;
				}
			}
			double scale=width/(maxX-minX);
			int height=(maxY-minY)*scale;
			MeshPoint *mp3(MeshPointStack(stack)); 

			// We compute the min and max of u
			double minu,maxu;
			if(u!=0){
				int it=0,k=0;
				mp3->setP(&Th,it,k);
				maxu=GetAny<double>((*u)(stack));
				minu=maxu;
				for(it=0;it<Th.nt;it++)
				for(k=0;k<3;k++)
				{
					mp3->setP(&Th,it,k);
					double val=GetAny<double>((*u)(stack));
					if(val>maxu) maxu=val;
					if(val<minu) minu=val;	
				}
			}
			// We write the mesh in SVG format
			// At this stage, it is basic...
			if(strokeWidth==0){ 
				fprintf(fp,"<svg version='1.1' xmlns='http://www.w3.org/2000/svg' height='%d' width='%d' shape-rendering='crispEdges'>\n",height,width);
			} else 
				fprintf(fp,"<svg version='1.1' xmlns='http://www.w3.org/2000/svg' height='%d' width='%d'>\n",height,width);
			// Display the function and the mesh
			if(strokeWidth==0){
				fprintf(fp,"<g stroke='none'>\n");
			} else 
				fprintf(fp,"<g stroke='%s' stroke-width='%d' stroke-linejoin='round'>>\n",stroke.c_str(),strokeWidth);
			for(int it=0;it<Th.nt;it++)			// Display of the triangles
			{
				string style;
				const Mesh::Triangle  & K(Th.t(it));
				if(u!=0){
					// a Finit element function is given, we define the local gradient
					// First, we get the value of u on the vertices
					double value[3];
					for(int k=0;k<3;k++){
						mp3->setP(&Th,it,k);
						value[k]=GetAny<double>((*u)(stack));
					}
					// computation of the gradient up to a multiplicative constant 
					double dxu=0,dyu=0;
					for(int k=0;k<3;k++){
						int k1=Th.operator()(K[(k+1)%3]);
						int k2=Th.operator()(K[(k+2)%3]);
						const  Mesh::Vertex & P1 = Th.vertices[k1];
						const  Mesh::Vertex & P2 = Th.vertices[k2];
						double Tx=(P2.x-P1.x)/(2.*K.area);
						double Ty=(P2.y-P1.y)/(2.*K.area);
						dxu+=-Ty*value[k];	
						dyu+=Tx*value[k];	
					}
					// We compute the extremal points on K of [dxu,dyu] . (x,y) and the value of u
					double Xmin,Ymin,Xmax,Ymax,uAtMin,uAtMax;
					{
						int k=0;
						int iv=Th.operator()(K[k]);
	                                        const  Mesh::Vertex & P = Th.vertices[iv];
						Xmin=P.x;Xmax=P.x;Ymin=P.y;Ymax=P.y;
						uAtMin=value[0];uAtMax=value[0];
						double minVal=dxu*P.x+dyu*P.y;
						double maxVal=minVal;	
						for(k=1;k<3;k++){
							int iv=Th.operator()(K[k]);
	                                        	const  Mesh::Vertex & P = Th.vertices[iv];
							double Val=dxu*P.x+dyu*P.y;
							if(Val<minVal){minVal=Val;Xmin=P.x;Ymin=P.y;uAtMin=value[k];}
				                        if(Val>maxVal){maxVal=Val;Xmax=P.x;Ymax=P.y;uAtMax=value[k];}
						}
						
					}
					if(uAtMin==uAtMax){// u is constant on the triangle, no gradient is applied
						ostringstream strs;
						strs<<" opacity='"<<(uAtMin-minu)/(maxu-minu)*((*rgba)[7]-(*rgba)[3])+(*rgba)[3]<<"'";
						strs<<" fill='rgb("<<(int)((uAtMin-minu)/(maxu-minu)*((*rgba)[4]-(*rgba)[0])+(*rgba)[0])<<","<<(int)((uAtMin-minu)/(maxu-minu)*((*rgba)[5]-(*rgba)[1])+(*rgba)[1])<<","<<(int)((uAtMin-minu)/(maxu-minu)*((*rgba)[6]-(*rgba)[2])+(*rgba)[2])<<")'";
						style=strs.str();
					} else {
						// We define a linear Gradient
						double x1=Xmin,y1=Ymin;
						double x2=Xmax,y2=Ymax;
						double s=(-(x1-x2)*dyu+(y1-y2)*dxu)/(dxu*dxu+dyu*dyu);
						x2=x2-s*dyu;y2=y2+s*dxu;
						fprintf(fp,"<linearGradient id='Grad-%d' gradientUnits='userSpaceOnUse' x1='%f' y1='%f' x2='%f' y2='%f'>\n",it,scale*(x1-minX),scale*(maxY-y1),scale*(x2-minX),scale*(maxY-y2));
						fprintf(fp,"\t<stop offset='0' stop-opacity='%f' stop-color='rgb(%d,%d,%d)'/>\n",(uAtMin-minu)/(maxu-minu)*((*rgba)[7]-(*rgba)[3])+(*rgba)[3], (int)((uAtMin-minu)/(maxu-minu)*((*rgba)[4]-(*rgba)[0])+(*rgba)[0]), (int)((uAtMin-minu)/(maxu-minu)*((*rgba)[5]-(*rgba)[1])+(*rgba)[1]), (int)((uAtMin-minu)/(maxu-minu)*((*rgba)[6]-(*rgba)[2])+(*rgba)[2]));
						
						fprintf(fp,"\t<stop offset='1' stop-opacity='%f' stop-color='rgb(%d,%d,%d)'/>\n",(uAtMax-minu)/(maxu-minu)*((*rgba)[7]-(*rgba)[3])+(*rgba)[3],(int)((uAtMax-minu)/(maxu-minu)*((*rgba)[4]-(*rgba)[0])+(*rgba)[0]), (int)((uAtMax-minu)/(maxu-minu)*((*rgba)[5]-(*rgba)[1])+(*rgba)[1]), (int)((uAtMax-minu)/(maxu-minu)*((*rgba)[6]-(*rgba)[2])+(*rgba)[2]));
						fprintf(fp,"</linearGradient>\n");
						ostringstream strs;
						strs<<"fill='url(#Grad-"<<it<<")' ";
						style=strs.str();
					}
				}
				fprintf(fp,"<path %s d='M",style.c_str());
				for(int k=0;k<3;k++){
					int iv=Th.operator()(K[k]);
					const  Mesh::Vertex & P = Th.vertices[iv];
					fprintf(fp," %f %f ",scale*(P.x-minX),scale*(maxY-P.y));
					if(k==2) {fprintf(fp,"z'");} else fprintf(fp,"L");
				}
				fprintf(fp,"/>\n");
			}
			fprintf(fp,"</g>\n");
			// Display the Boundary
			fprintf(fp,"<g stroke='%s' stroke-width='%d' shape-rendering='crispEdges'>\n",strokeEdge.c_str(),strokeWidthEdge);
			for(long ibe=0;ibe<Th.neb;ibe++){
				const Mesh::BorderElement &K( Th.be(ibe) );
				fprintf(fp,"<path d='M ");
				for(int ii=0; ii<2; ii++){
					int ien = Th.operator()(K[ii]);
					const  Mesh::Vertex & P = Th.vertices[ien];
					fprintf(fp," %f %f ",scale*(P.x-minX),scale*(maxY-P.y));
					if(ii==0) fprintf(fp,"L ");
				}
				fprintf(fp,"'/>\n");
			}
			fprintf(fp,"</g>\n");
			fprintf(fp,"</svg>");
		}
};
		basicAC_F0::name_and_type SVG_WriteMesh_Op::name_param[]= {	// Optional parameters
			{"width", &typeid(long)},				// width of the SVG (default value = 500)
			{"strokeWidth", &typeid(long)},				// width of the stroke for the mesh (default value = 0)
			{"stroke",&typeid(string*)},				// color of the stroke for the mesh (default value = black)
			{"rgba",&typeid(KN<double>*)},				// rgba fill  = [rm,gm,bm,am,RM,GM,BM,AM]
										// (rm,gm,bm,am) = rgba color for the min value
										// (RM,GM,BM,AM) = rgba color for the max value
			{"strokeWidthEdge", &typeid(long)},			// width of the stroke for the boundary edges (default value = 0)
			{"strokeEdge",&typeid(string*)}				// color of the stroke for the boundary edges (default value = red)
		};

static void finit()
{
    
    typedef const Mesh *pmesh;
    
    Global.Add("savesvg","(",new OneOperatorCode<SVG_WriteMesh_Op>);
}

LOADFUNC(finit);