// ------------------------------------------------------------
// main.cc
//
// ------------------------------------------------------------
#define _USE_MATH_DEFINES
#include <cmath>
#include <complex>
#include <iostream>
#include <string>
#include <vector>
#include "../program2D/source/ugblock2D.h"
#include "source/vectorVar.h"
#include "solver/iterativeSolversSystem.h"


using std::complex;
using namespace ::_COLSAMM_;

#define L2NORM(vector) sqrt(product(vector,vector).real())
#define RESTART 0


double isZero(double x)
{
	return fabs(x) < 1e-10?1.0:0.0;
}



int main(int argc, char** argv) {
	std::ofstream DATEI;

	int n = 10;

	//SectionedRectangle geo(v_x, v_y);  
        //Rectangle geo(1.0, 1.0);
	Disc geo(2.0,1.0);

	Boundary_Marker2D boundary(&geo);
	Unstructured2DGrid_Marker2D inner_points(&geo);
	inner_points.complement(&boundary);

        // boudary integral test 
        Blockgrid2D  blockBoTest(&geo,n);
        Variable2D<double> one(blockBoTest);
        Variable2D<double>  ff(blockBoTest);

        Variable2D<double>  u(blockBoTest);
        Variable2D<double>  interpol(blockBoTest);
        Variable2D<double>  g(blockBoTest);
        Variable2D<double>  p(blockBoTest);	
	
        Local_stiffness_matrix2D<double> boIntegral(blockBoTest);  

	DATEI.open("ug.vtk",std::ios :: out);
	geo.PrintVtk(&DATEI);
	DATEI.close();

	X_coordinate2d X(blockBoTest);
	Y_coordinate2d Y(blockBoTest);

/*
        Local_stiffness_matrix2D<double>  locStiff(blockBoTest);
	//	Test.Calculate(v_()*w_());
        locStiff.Calculate_HelmBoundary();        
	one = 1.0;
	ff = locStiff(one);	
	cout << "inegral Rand: " << product(ff,one) << endl;
*/

        u = X;
        g = Y -0.2;    

        VariableVector<double, Variable2D<double> > vectorU(u,u);
        VariableVector<double, Variable2D<double> > vectorB(g,g);
        VariableVector<double, Variable2D<double> > vectorI(interpol,p);  
        vectorI = vectorU + vectorB;
  
        cout << " Maximum I: " << Maximum(vectorI) << endl;
        cout << " Minimum I: " << Minimum(vectorI) << endl;

        cout << " Maximum U: " << Maximum(vectorU) << endl;
        cout << " Minimum U: " << Minimum(vectorU) << endl;

        cout << " Maximum B: " << Maximum(vectorB) << endl;
        cout << " Minimum B: " << Minimum(vectorB) << endl;
	
	
	DATEI.open("interpol.vtk");
	interpol.Print_VTK(DATEI);
	DATEI.close();	

	DATEI.open("one.vtk");
	one.Print_VTK(DATEI);
	DATEI.close();

	DATEI.open("ff.vtk");
	ff.Print_VTK(DATEI);
	DATEI.close();
}