#include "sparselizardbase.h"

#include "string.h"

using namespace std;
using namespace mathop;


struct Mesher {
	mesh model;

	vector<string> files;
	string tempName;
	bool __initialized;
	int segmentation_number;


	Mesher() {
		this->__initialized = false;
		this->tempName = "temp.msh";
		this->segmentation_number = 1000000;//997;
		// 16127
	}

	int add(const char* filename) {
		this->files.push_back(filename);
		return this->getfileregion(this->files.size()); 
	}

	void load() {
		this->model.load(false, this->files, 5, false, true, this->segmentation_number);
	}

	void addit(const char* filename) {
		vector<string> temp;
		if(this->__initialized) {
			this->write(this->tempName);
			temp.push_back(this->tempName);
		} else {
			this->__initialized = true;
		}
		temp.push_back(filename);
		this->model.load(true, temp, 5, false);
	}

	void write(string filename) {
		this->model.write(filename);
	}

	void rotate	  (double x,double y,double z) { this->model.rotate(x,y,z); }
	void translate(double x,double y,double z) { this->model.shift(x,y,z); }
	void scale	  (double x,double y,double z) { this->model.scale(x,y,z); }

	void rotate	  (int r, double x,double y,double z) { this->model.rotate(r,x,y,z); }
	void translate(int r, double x,double y,double z) { this->model.shift(r,x,y,z); }
	void scale	  (int r, double x,double y,double z) { this->model.scale(r,x,y,z); }

	void debug(){ this->model.printphysicalregions(); }

	int getphysicalregion(int file, int region, int dimension) {
		int result = (this->segmentation_number + 1000 * file)*(dimension + 1) + region;
		return result;
	}

	int getfileregion(int file) {
		return 999000000 + file;
	}

	void printphysicalregions() {
		vector<int> regions = this->model.getphysicalregions()->getallnumbers();
		for (int i=0;i<regions.size();i++) {
			int temp 	= regions[i];
			int region 	= temp % (this->segmentation_number);
			int dim 	= temp / (this->segmentation_number) -1;
			int pd 		= (1000*(dim+1));
			int file 	= region / pd;
			int model 	= region % pd;
			cout << temp << "\tregion: " << region << "\tmodel: " << model << "\tdim: " << dim << "\tfile: " << file << endl;
		}
	}
};





int resolve(Mesher mesh) {

	int Stator = mesh.getfileregion(1);
	int Magnet = mesh.getfileregion(2);
	int All = regionunion({Stator, Magnet}); 
	
	int StatorFace = mesh.getphysicalregion(1, 102, 2);
	int StatorFaceBack = mesh.getphysicalregion(1, 103, 2);
	int MagnetFace = mesh.getphysicalregion(2, 5, 	2);

    // cout << endl << "[!]Continuitycondition: " << StatorFace << " , " << MagnetFace << endl << endl;
	// cout << "a";
	// int Intersection = regionintersection({StatorFace, MagnetFace});
	// cout << "a";


    // Nodal shape functions 'h1' for the z component of the vector potential.
    // field a("h1xyz");
	spanningtree spantree({Stator, Magnet});
    field a("hcurl", spantree);
    a.setgauge(All);
    a.setorder(All, 2);

    field x("x"), y("y"), z("z");

    // Vacuum magnetic permeability [H/m]: 
    double mu0 = 4.0*getpi()*1e-7;
    // Define the permeability in all regions.
    // Taking into account saturation and measured B-H curves can be easily done
    // by defining an expression based on a 'spline' object (see documentation).
    parameter mu;

    mu|All = 5000.0*mu0;
    // Overwrite on non-magnetic regions:
    mu|Magnet = 1.05*mu0;

    // The remanent induction field in the magnet is 0.5 Tesla perpendicular to the magnet:
    expression magnetdirection = array3x1(0,0,1);
    expression bremanent = 0.5 * magnetdirection;     


	formulation magnetostatics;
    // The strong form of the magnetostatic formulation is curl( 1/mu * curl(a) ) = j, with b = curl(a):
    magnetostatics += integral(All, 1/mu* curl(dof(a)) * curl(tf(a)) );
    // Add the remanent magnetization of the rotor magnet:
    magnetostatics += integral(Magnet, -1/mu* bremanent * curl(tf(a)));
    // Add the current density source js [A/m2] in the z direction in the stator windings.
    // A three-phased actuation is used. The currents are dephased by the mechanical angle
    // times the number of pole pairs. This gives a stator field rotating at synchronous speed.

    // Change the phase (degrees) to tune the electric angle: 
    // double phase = 0;

    // parameter jsz;

    // jsz|winda = Imax/windarea * sin( (phase + 4.0*alpha - 0.0) * getpi()/180.0);
    // jsz|windb = Imax/windarea * sin( (phase + 4.0*alpha - 60.0) * getpi()/180.0);
    // jsz|windc = Imax/windarea * sin( (phase + 4.0*alpha - 120.0) * getpi()/180.0);

    // magnetostatics += integral(windings, -array3x1(0,0,jsz) * tf(a));


    magnetostatics += continuitycondition(StatorFace, MagnetFace, a, a, 1, false);

    // Rotor-stator continuity condition (including antiperiodicity settings with factor '-1'):
    // magnetostatics += continuitycondition(gammastat, gammarot, az, az, {0,0,0}, alpha, 45.0, -1.0);
    // Rotor and stator antiperiodicity condition:
    // magnetostatics += periodicitycondition(gamma1, gamma2, az, {0,0,0}, {0,0,45.0}, -1.0);


    solve(magnetostatics);

    // int Intersection = regionintersection({Stator, Magnet});
    // curl(a).write(all, "results/MagneticField"+std::to_string((int)alpha)+".vtu", 2);
    curl(a).write(All, "results/MagneticField.vtu", 2);



	return 0;
}















int main(void)
{
    SlepcInitialize(0,{},0,0);

    Mesher testmesh;
    int Stator	= testmesh.add("models/Stator.msh");
    int Magnet	= testmesh.add("models/Magnet45.msh");
    testmesh.load();
    testmesh.printphysicalregions();

    // testmesh.rotate		(Magnet, 0, 0, 0);
    // testmesh.translate	(Magnet, 0, 0, 2);
	testmesh.write("models/combined.msh");

    resolve(testmesh);

    SlepcFinalize();

    return 0;
}