#include "inmost.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "../../Source/Misc/utils.h"


using namespace INMOST;

#ifndef M_PI
#define M_PI 3.141592653589
#endif

#if defined(USE_MPI)
#define BARRIER MPI_Barrier(MPI_COMM_WORLD);
#else
#define BARRIER
#endif

//shortcuts
typedef Storage::bulk bulk;
typedef Storage::real real;
typedef Storage::integer integer;
typedef Storage::enumerator enumerator;
typedef Storage::real_array real_array;
typedef Storage::var_array var_array;

bool print_niter = false; //save file on nonlinear iterations

//If true, then tag ELASTIC_TENSOR contains stiffness tensor C,
//determining strain to stress relation C*\varepsilon = \sigma, otherwise
//the tag contains compliance tensor S, determining stress to strain relation
//S*\sigma = \varepsilon. Here \varepsilon = \frac{G+G^T}{2},
//with G - gradient of displacement matrix
bool inverse_tensor = true;
bool reference_solution = false;
bool print_matrix = false;


template<typename T>
static void SVD2Eigen(const Matrix<T> & U, Matrix<T> & S, Matrix<T> & V)
{
	for (unsigned i = 0; i < V.Cols(); ++i)
	{
		double dot = 0.0;
		for (unsigned j = 0; j < V.Rows(); ++j)
			dot += get_value(U(j, i))*get_value(V(j, i));
		if (dot < 0.0)
		{
			S(i, i) *= -1;
			for (unsigned j = 0; j < V.Rows(); ++j)
				V(j, i) *= -1;
		}
	}
	//check
	//if ((U - V).FrobeniusNorm() > 1.0e-8)
	//	(U - V).Print();
}


class BlockRow
{
	std::map< unsigned, rMatrix > entries;
	unsigned height, width;
	rMatrix & GetAdd(unsigned col)
	{
		std::map< unsigned, rMatrix >::iterator s = entries.find(col/width);
		if( s == entries.end() )
			return (entries[col/width] = rMatrix(height,width,0.0));
		else
			return s->second;
	}
public:
	BlockRow(unsigned h, unsigned w) : height(h), width(w) {}
	rMatrix & operator [] (unsigned col) { return GetAdd(col); }
	double & operator () (unsigned i, unsigned j) { return GetAdd(j)(i%height,j%width);}
	double operator () (unsigned i, unsigned j) const
	{
		std::map< unsigned, rMatrix >::const_iterator s = entries.find(j/width);
		if( s != entries.end() )
			return s->second(i%height,j%width);
		throw -1;
	}
	std::map< unsigned, rMatrix >::iterator Begin() {return entries.begin();}
	std::map< unsigned, rMatrix >::iterator End() {return entries.end();}
};

//indices for tensor
const int IC11 = 0;
const int IC12 = 1;
const int IC13 = 2;
const int IC14 = 3;
const int IC15 = 4;
const int IC16 = 5;
const int IC22 = 6;
const int IC23 = 7;
const int IC24 = 8;
const int IC25 = 9;
const int IC26 = 10;
const int IC33 = 11;
const int IC34 = 12;
const int IC35 = 13;
const int IC36 = 14;
const int IC44 = 15;
const int IC45 = 16;
const int IC46 = 17;
const int IC55 = 18;
const int IC56 = 19;
const int IC66 = 20;


/*
 
 stored tensor
 C_1  C_2  C_3  C_4  C_5  C_6
      C_7  C_8  C_9  C_10 C_11
           C_12 C_13 C_14 C_15
                C_16 C_17 C_18
                     C_19 C_20
                          C_21
 
 
 voigt notation tensor
 
 c_11 c_12 c_13 c_14 c_15 c_16
 c_12 c_22 c_23 c_24 c_25 c_26
 c_13 c_23 c_33 c_34 c_35 c_36
 c_14 c_24 c_34 c_44 c_45 c_46
 c_15 c_25 c_35 c_45 c_55 c_56
 c_16 c_26 c_36 c_46 c_56 c_66
 
 co-normal tensors
 
     | K_1   \vec{n} \cdot \nabla u + K_4   \vec{n} \cdot \nabla v + K_5 \vec{n} \cdot \nabla w |
 F = | K_4^T \vec{n} \cdot \nabla u + K_2   \vec{n} \cdot \nabla v + K_6 \vec{n} \cdot \nabla w |
     | K_5^T \vec{n} \cdot \nabla u + K_6^T \vec{n} \cdot \nabla v + K_3 \vec{n} \cdot \nabla w |
 
 tensors representations
 
       | c_11 c_16 c_15 |
 K_1 = | c_16 c_66 c_56 |
       | c_15 c_56 c_55 |
 
       | c_66 c_26 c_46 |
 K_2 = | c_26 c_22 c_24 |
       | c_46 c_24 c_44 |
 
       | c_55 c_45 c_35 |
 K_3 = | c_45 c_44 c_34 |
       | c_35 c_34 c_33 |
 
       | c_16 c_66 c_56 |
 K_4 = | c_12 c_26 c_25 |
       | c_14 c_46 c_45 |
 
       | c_15 c_56 c_55 |
 K_5 = | c_14 c_46 c_45 |
       | c_13 c_36 c_35 |
 
       | c_56 c_25 c_45 |
 K_6 = | c_46 c_24 c_44 |
       | c_36 c_23 c_34 |
 
 two-point flux approximation
 
 F = T (\vec{u}_2 - \vec{u}_1)
 
 here the transimissibility matrix is
 
 T = C'_1 C'_2 (d_1 C'_2 + d_2 C'_1)^{-1}
 
      | k_1 k_4 k_5 |
 C' = | k_4 k_2 k_6 |
      | k_5 k_6 k_3 |
 
 with
 
 k_i = \vec{n} \cdot K_i \vec{n} = \vec{n} \cdot K_i^T \vec{n}
 
 */


void CTensor(const real_array & Cv, rMatrix & C)
{
	C = rMatrix::FromTensor(Cv.data(),Cv.size(),6);
}

void KTensor(const real_array & Cv, rMatrix & K)
{
	int si = 0, sj = 0;
	rMatrix C = rMatrix::FromTensor(Cv.data(),Cv.size(),6);
	
	K(0+si,0+sj) = C(0,0); K(0+si,1+sj) = C(0,5); K(0+si,2+sj) = C(0,4);
	K(1+si,0+sj) = C(0,5); K(1+si,1+sj) = C(5,5); K(1+si,2+sj) = C(4,5);
	K(2+si,0+sj) = C(0,4); K(2+si,1+sj) = C(4,5); K(2+si,2+sj) = C(4,4);
	
	si = sj = 3;
	
	K(0+si,0+sj) = C(5,5); K(0+si,1+sj) = C(1,5); K(0+si,2+sj) = C(3,5);
	K(1+si,0+sj) = C(1,5); K(1+si,1+sj) = C(1,1); K(1+si,2+sj) = C(1,3);
	K(2+si,0+sj) = C(3,5); K(2+si,1+sj) = C(1,3); K(2+si,2+sj) = C(3,3);
	
	si = sj = 6;
	
	K(0+si,0+sj) = C(4,4); K(0+si,1+sj) = C(3,4); K(0+si,2+sj) = C(2,4);
	K(1+si,0+sj) = C(3,4); K(1+si,1+sj) = C(3,3); K(1+si,2+sj) = C(2,3);
	K(2+si,0+sj) = C(2,4); K(2+si,1+sj) = C(2,3); K(2+si,2+sj) = C(2,2);
	
	si = 0, sj = 3;
	
	K(0+si,0+sj) = C(0,5); K(0+si,1+sj) = C(0,1); K(0+si,2+sj) = C(0,3);
	K(1+si,0+sj) = C(5,5); K(1+si,1+sj) = C(1,5); K(1+si,2+sj) = C(3,5);
	K(2+si,0+sj) = C(4,5); K(2+si,1+sj) = C(1,4); K(2+si,2+sj) = C(3,4);
	
	si = 3, sj = 0;
	
	K(0+si,0+sj) = C(0,5); K(0+si,1+sj) = C(5,5); K(0+si,2+sj) = C(4,5);
	K(1+si,0+sj) = C(0,1); K(1+si,1+sj) = C(1,5); K(1+si,2+sj) = C(1,4);
	K(2+si,0+sj) = C(0,3); K(2+si,1+sj) = C(3,5); K(2+si,2+sj) = C(3,4);
	
	si = 0, sj = 6;
	
	K(0+si,0+sj) = C(0,4); K(0+si,1+sj) = C(0,3); K(0+si,2+sj) = C(0,2);
	K(1+si,0+sj) = C(4,5); K(1+si,1+sj) = C(3,5); K(1+si,2+sj) = C(2,5);
	K(2+si,0+sj) = C(4,4); K(2+si,1+sj) = C(3,4); K(2+si,2+sj) = C(2,4);
	
	si = 6, sj = 0;
	
	K(0+si,0+sj) = C(0,4); K(0+si,1+sj) = C(4,5); K(0+si,2+sj) = C(4,4);
	K(1+si,0+sj) = C(0,3); K(1+si,1+sj) = C(3,5); K(1+si,2+sj) = C(3,4);
	K(2+si,0+sj) = C(0,2); K(2+si,1+sj) = C(2,5); K(2+si,2+sj) = C(2,4);
	
	si = 3, sj = 6;
	
	K(0+si,0+sj) = C(4,5); K(0+si,1+sj) = C(3,5); K(0+si,2+sj) = C(2,5);
	K(1+si,0+sj) = C(1,4); K(1+si,1+sj) = C(1,3); K(1+si,2+sj) = C(1,2);
	K(2+si,0+sj) = C(3,4); K(2+si,1+sj) = C(3,3); K(2+si,2+sj) = C(2,3);
	
	si = 6, sj = 3;
	
	K(0+si,0+sj) = C(4,5); K(0+si,1+sj) = C(1,4); K(0+si,2+sj) = C(3,4);
	K(1+si,0+sj) = C(3,5); K(1+si,1+sj) = C(1,3); K(1+si,2+sj) = C(3,3);
	K(2+si,0+sj) = C(2,5); K(2+si,1+sj) = C(1,2); K(2+si,2+sj) = C(2,3);
}

void GetBC(const real_array & bc, const rMatrix & n, rMatrix & Ra, rMatrix & Rb, rMatrix & r)
{
	const rMatrix I = rMatrix::Unit(3);
	double alpha_perp, alpha_parallel, beta_perp, beta_parallel;

	if( bc.size() == 6 )
	{
		double alpha, beta, proj;
		alpha = bc[0];
		beta = bc[1];
		proj = bc[2];
		r(0,0) = bc[3];
		r(1,0) = bc[4];
		r(2,0) = bc[5];
		if( proj )
		{
			alpha_perp = alpha;
			beta_parallel = beta;
			alpha_parallel = beta_perp = 0;
		}
		else
		{
			alpha_perp = alpha_parallel = alpha;
			beta_perp = beta_parallel = beta;
		}
	}
	else if( bc.size() == 7 )
	{
		alpha_perp = bc[0];
		beta_perp = bc[1];
		alpha_parallel = bc[2];
		beta_parallel = bc[3];
		r(0,0) = bc[4];
		r(1,0) = bc[5];
		r(2,0) = bc[6];
	}
	Ra = alpha_parallel*I + (alpha_perp-alpha_parallel)*n*n.Transpose();
	Rb = beta_parallel*I + (beta_perp-beta_parallel)*n*n.Transpose();
}

void PrintSV(const rMatrix & A)
{
	rMatrix U,S,V;
	A.SVD(U,S,V);
	std::cout << "singular values:";
	int cnt = 0;
	for(unsigned k = 0; k < A.Cols(); ++k) if( fabs(S(k,k)) > 1.0e-10 )
	{
		std::cout << " " << S(k,k);
		cnt++;
	}
	else
		std::cout << " 0";
	std::cout << " count: " << cnt << "/" << A.Cols() << std::endl;
}

void PrintRS(const rMatrix & A)
{
	double sum;
	std::cout << "row sum:";
	for(unsigned i = 0; i < A.Rows(); ++i)
	{
		sum = 0;
		for(unsigned j = 0; j < A.Cols(); ++j)
			sum+= A(i,j);
		std::cout << " " << sum;
	}
	std::cout << std::endl;
}

int main(int argc,char ** argv)
{
    Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
#if defined(USE_PARTITIONER)
    Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
#endif
    if( argc > 1 )
    {

        double ttt; // Variable used to measure timing
        bool repartition = false; // Is it required to redistribute the mesh?
        Mesh * m = new Mesh(); // Create an empty mesh
        { // Load the mesh
            ttt = Timer();
            m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
            if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
            if( m->isParallelFileFormat(argv[1]) ) //The format is
            {
                m->Load(argv[1]); // Load mesh from the parallel file format
                repartition = true; // Ask to repartition the mesh
            }
            else if( m->GetProcessorRank() == 0 ) m->Load(argv[1]); // Load mesh from the serial file format
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Load the mesh: " << Timer()-ttt << std::endl;
        }

		
#if defined(USE_PARTITIONER)
        if (m->GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner
        {
            { // Compute mesh partitioning
                ttt = Timer();
                Partitioner p(m); //Create Partitioning object
                p.SetMethod(Partitioner::INNER_KMEANS,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner
                p.Evaluate(); // Compute the partitioner and store new processor ID in the mesh
                BARRIER
                if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;
            }

            { //Distribute the mesh
                ttt = Timer();
                m->Redistribute(); // Redistribute the mesh data
                m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
                BARRIER
                if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
            }
        }
#endif

        { // prepare geometrical data on the mesh
            ttt = Timer();
            Mesh::GeomParam table;
            table[BARYCENTER]    = CELL | FACE; //Compute averaged center of mass
            table[NORMAL]      = FACE;        //Compute normals
            table[ORIENTATION] = FACE;        //Check and fix normal orientation
            table[MEASURE]     = CELL | FACE; //Compute volumes and areas
            //table[BARYCENTER]  = CELL | FACE; //Compute volumetric center of mass
			m->RemoveGeometricData(table);
            m->PrepareGeometricData(table); //Ask to precompute the data
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
        }

        // data tags for
        TagRealArray tag_UVW;// Displacement
        TagRealArray tag_S;  // Stress (for error)
        TagRealArray tag_C;  // Elasticity tensor in Voigt format
        TagRealArray tag_F;  // Forcing term
        TagRealArray tag_BC; // Boundary conditions
        TagRealArray tag_W;  // Gradient matrix
        TagRealArray tag_Ws; // Matrix to reconstruct stress from fluxes
		TagRealArray tag_FLUX; // Flux (for error)
		TagRealArray tag_f_coefs, tag_i_coefs, tag_f_rhs, tag_i_rhs;
		TagReferenceArray tag_f_elems, tag_i_elems;
		TagIntegerArray tag_Row;

        if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
        { // Exchange ghost cells
            ttt = Timer();
            m->ExchangeGhost(1,FACE); // Produce layer of ghost cells
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Exchange ghost: " << Timer()-ttt << std::endl;
        }
		
		std::cout << "Cells " << m->TotalNumberOf(CELL) << std::endl;
		std::cout << "Faces " << m->TotalNumberOf(FACE) << std::endl;
		std::cout << "Edges " << m->TotalNumberOf(EDGE) << std::endl;
		std::cout << "Nodes " << m->TotalNumberOf(NODE) << std::endl;
		
		const double vB[] =
		{
			1,0,0,0,0,0,
			0,0,0,0,0,1,
			0,0,0,0,1,0,
			0,0,0,0,0,1,
			0,1,0,0,0,0,
			0,0,0,1,0,0,
			0,0,0,0,1,0,
			0,0,0,1,0,0,
			0,0,1,0,0,0
		};
		const double viBtB[] =
		{
			1,0,0,0.0,0.0,0.0,
			0,1,0,0.0,0.0,0.0,
			0,0,1,0.0,0.0,0.0,
			0,0,0,0.5,0.0,0.0,
			0,0,0,0.0,0.5,0.0,
			0,0,0,0.0,0.0,0.5
		};
		// |  0  -z   y | u   | w_y - v_z |
		// |  z   0  -x | v = | u_z - w_x |
		// | -y   x   0 | w   | v_x - u_y |
		//
		// 0  0  0  0  0 -z  0 +y  0
		// 0  0 +z  0  0  0 -x  0  0
		// 0 -y  0 +x  0  0  0  0  0
		//
		// 0  0  0  0  0  0  0  0  0
		// 0  0  0  0  0  0  0  1  0
		// 0  0  0  0  0 -1  0  0  0
		// 0  0  0  0  0  0 -1  0  0
		// 0  0  0  0  0  0  0  0  0
		// 0  0  1  0  0  0  0  0  0
		// 0  0  0  1  0  0  0  0  0
		// 0 -1  0  0  0  0  0  0  0
		// 0  0  0  0  0  0  0  0  0
		//
		const double vCurl[] =
		{
			0,  0,  0,  0,  0,  0,  0,  0,  0,
			0,  0,  0,  0,  0,  0,  0,  1,  0,
			0,  0,  0,  0,  0, -1,  0,  0,  0,
			0,  0,  0,  0,  0,  0, -1,  0,  0,
			0,  0,  0,  0,  0,  0,  0,  0,  0,
			0,  0,  1,  0,  0,  0,  0,  0,  0,
			0,  0,  0,  1,  0,  0,  0,  0,  0,
			0, -1,  0,  0,  0,  0,  0,  0,  0,
			0,  0,  0,  0,  0,  0,  0,  0,  0
		};
		const rMatrix B(vB,9,6);
		const rMatrix iBtB(viBtB,6,6);
		const rMatrix Curl(vCurl,9,9);
		const rMatrix I = rMatrix::Unit(3);
		const rMatrix I9 = rMatrix::Unit(9);
		
		//PrintSV(B);
		std::cout << "B^T*B" << std::endl;
		(B.Transpose()*B).Print();
		
		std::cout << "B*iBtB*B^T" << std::endl;
		(B*iBtB*B.Transpose()).Print();
		
		std::cout << "B^T*B*iBtB" << std::endl;
		(B.Transpose()*B*iBtB).Print();
		
		std::cout << "B*B^T" << std::endl;
		(B*B.Transpose()).Print();

        { //initialize data
            if( m->HaveTag("ELASTIC_TENSOR") ) // is elasticity tensor already defined on the mesh?
                tag_C = m->GetTag("ELASTIC_TENSOR"); // get the elasticity tensor
			else
			{
				std::cout << "No ELASTIC_TENSOR on mesh" << std::endl;
				return -1;
			}

			if( m->HaveTag("FORCE") ) //Is there force on the mesh?
			{
				tag_F = m->GetTag("FORCE"); //initial force
				assert(tag_F.isDefined(CELL)); //assuming it was defined on cells
			} // end of force

            if( m->HaveTag("BOUNDARY_CONDITION") ) //Is there boundary condition on the mesh?
            {
                tag_BC = m->GetTag("BOUNDARY_CONDITION");
                //initialize unknowns at boundary
            }
			tag_i_elems = m->CreateTag("i_elems",DATA_REFERENCE,FACE,NONE);
			tag_f_elems = m->CreateTag("f_elems",DATA_REFERENCE,FACE,NONE);
			tag_i_coefs = m->CreateTag("i_coefs",DATA_REAL,FACE,NONE);
			tag_f_coefs = m->CreateTag("f_coefs",DATA_REAL,FACE,NONE);
			tag_i_rhs = m->CreateTag("i_rhs",DATA_REAL,FACE,NONE,3);
			tag_f_rhs = m->CreateTag("f_rhs",DATA_REAL,FACE,NONE,3);
            tag_W = m->CreateTag("W",DATA_REAL,CELL,NONE);
            tag_Ws = m->CreateTag("Ws",DATA_REAL,CELL,NONE);
			
            ttt = Timer();
            //Assemble gradient matrix W on cells
			
			tag_Row = m->CreateTag("row_index",DATA_INTEGER,FACE,NONE,2);
#if defined(USE_OMP)
//#pragma omp parallel
#endif
			{
				rMatrix N, R, L, M(9,9), T, K(9,9), iK(9,9), C(6,6), Q, W1, W2, W3, W3s, U,S,V, w, u, v;
				rMatrix x(3,1), xf(3,1), n(3,1);
				double area; //area of the face
				//double volume; //volume of the cell
				double dist; //distance from cell center to face
#if defined(USE_OMP)
#pragma omp for
#endif
				for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
				{
					Cell cell = m->CellByLocalID(q);
					ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
					int NF = (int)faces.size(); //number of faces;
					//~ volume = cell->Volume(); //volume of the cell
					cell->Barycenter(x.data());;
					//get permeability for the cell
					KTensor(tag_C[cell],K);
					CTensor(tag_C[cell],C);
					iK = K.PseudoInvert(1.0e-11);
					if( !inverse_tensor ) K = K.Invert();
					//K += rMatrix::Unit(9)*1.0e-6*K.FrobeniusNorm();
					//PrintSV(K);
					N.Resize(3*NF,9); //co-normals
					//NQ.Resize(3*NF,9);
					//T.Resize(3*NF,9); //transversals
					R.Resize(3*NF,9); //directions
					L.Resize(3*NF,3*NF);
					//M.Resize(3*NF,3*NF);
					L.Zero();
					
					for(int k = 0; k < NF; ++k) //loop over faces
					{
						area = faces[k].Area();
						faces[k].Barycenter(xf.data());
						faces[k].OrientedUnitNormal(cell->self(),n.data());
						dist = n.DotProduct(xf-x);
						// assemble matrix of directions
						R(3*k,3*(k+1),0,9) = I.Kronecker((xf-x).Transpose());
						// assemble matrix of co-normals
						
						L(3*k,3*(k+1),3*k,3*(k+1)) = (area/dist)*(I.Kronecker(n.Transpose())*K*I.Kronecker(n));
						
						N(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose());
#if defined(USE_OMP)
#pragma omp critical
#endif
						{
							if( faces[k].BackCell() == cell )
								tag_Row[faces[k]][0] = k; //detect my row
							else
								tag_Row[faces[k]][1] = k; //detect my row
						}
					} //end of loop over faces
					//R += N*Curl;
					tag_Ws[cell].resize(6*3*NF);
					tag_Ws(cell,6,3*NF) = ((R*B*iBtB).Transpose()*N*B).Invert()*(R*B*iBtB).Transpose();
					tag_W[cell].resize(9*NF*NF);
					
					if( true )
					{
						double alpha = 0;
						double beta = alpha;
						
						//M = B*iBtB*B.Transpose();
						//R = R*(I9 + M)*0.5;
						W1 = (N*K+alpha*L*R)*((N*K+alpha*L*R).Transpose()*R).PseudoInvert(1.0e-11)*(N*K+alpha*L*R).Transpose();
						W2 = L - (1+beta)*(L*R)*((L*R).Transpose()*R).PseudoInvert(1.0e-11)*(L*R).Transpose();
					}
					else
					{
						M = B*iBtB*B.Transpose();
						W1 = (N*B*C)*((N*B*C).Transpose()*R*B*iBtB).Invert()*(N*B*C).Transpose();
						W2 = L - (L*R*B*iBtB)*((L*R*B*iBtB).Transpose()*R*B*iBtB).Invert()*(L*R*B*iBtB).Transpose();
						
					}
					
					tag_W(cell,3*NF,3*NF) = W1+W2;// + rMatrix::Unit(3*NF)*volume;
					
					if(false)if(!K.isSymmetric() )
					{
						std::cout << "K nonsymmetric: " << std::endl;
						(K-K.Transpose()).Print();
					}
					
					//if( false )
					if( !tag_W(cell,3*NF,3*NF).isSymmetric(1.0e-5) )
					{
						std::cout << "W nonsymmetric: " << std::endl;
						(tag_W(cell,3*NF,3*NF)-tag_W(cell,3*NF,3*NF).Transpose()).Print();
						//std::cout << "(N*C)^T*R" << std::endl;
						//((N*C).Transpose()*R).Print();
					}
					
					if( false )
						if( (N*C - tag_W(cell,3*NF,3*NF)*R).FrobeniusNorm() > 1.0e-3 )
						{
							std::cout << "error: " << std::endl;
							(N*C - tag_W(cell,3*NF,3*NF)*R).Print();
						}
					
					
					
					if( false )
					{
						std::cout << "total       ";
						PrintSV(tag_W(cell,3*NF,3*NF));
					}
					
					
				} //end of loop over cells
				rMatrix C1(3,3),C2(3,3), iC12sum, C1sum(3,3), C2sum(3,3), CC(3,3), CCI(3,3), CCF(3,3), B_D(3,3), B_N(3,3), B_R(3,1);
				rMatrix W[2];
#if defined(USE_OMP)
#pragma omp for
#endif
				for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
				{
					Face f = m->FaceByLocalID(q);
					
					Storage::integer_array rows = f->IntegerArray(tag_Row);
					ElementArray<Face> faces[2];
					ElementArray<Cell> cells = f->getCells();
					//std::cout << "face " << f->LocalID() << " NC " << cells.size() << std::endl;
					for(int q = 0; q < (int)cells.size(); ++q)
					{
						Cell cell = cells[q];
						faces[q] = cell->getFaces(); //obtain faces of the cell
						W[q] = raMatrixMake(cell->RealArray(tag_W).data(),3*faces[q].size(),3*faces[q].size());
					}
					Storage::reference_array f_elems = f->ReferenceArray(tag_f_elems);
					Storage::reference_array i_elems = f->ReferenceArray(tag_i_elems);
					Storage::real_array f_coefs = f->RealArray(tag_f_coefs);
					Storage::real_array i_coefs = f->RealArray(tag_i_coefs);
					tag_f_rhs(f,3,1).Zero();
					tag_i_rhs(f,3,1).Zero();
					real aF = 1;
					//f_rhs = i_rhs = 0;
					real multi = 1, multf = -aF;//, multfbnd = -aF;
					if( cells.size() == 2 ) //internal face
					{
						C1 = W[0](rows[0]*3,rows[0]*3+3,rows[0]*3,rows[0]*3+3)*aF;
						C2 = W[1](rows[1]*3,rows[1]*3+3,rows[1]*3,rows[1]*3+3)*aF;
						iC12sum = (C1+C2).Invert();
						C1sum.Zero();
						C2sum.Zero();
						//add back cell faces and compute sum
						for(int q = 0; q < (int)faces[0].size(); ++q)
						{
							CC = W[0](rows[0]*3,rows[0]*3+3,q*3,q*3+3);
							if( CC.FrobeniusNorm() > 1.0e-12 )
							{
								C1sum += CC;
								CC = iC12sum*CC;
								if( faces[0][q] != f )
								{
									CCI = -CC*multi;
									CCF = -CC*C2*multf;
									i_coefs.insert(i_coefs.end(),CCI.data(),CCI.data()+9);
									f_coefs.insert(f_coefs.end(),CCF.data(),CCF.data()+9);
									i_elems.push_back(faces[0][q]);
									f_elems.push_back(faces[0][q]);
									
								}
							}
							
						}
						//add front cell faces and compute sum
						for(int q = 0; q < (int)faces[1].size(); ++q)
						{
							CC = W[1](rows[1]*3,rows[1]*3+3,q*3,q*3+3);
							if( CC.FrobeniusNorm() > 1.0e-12 )
							{
								C2sum += CC;
								CC = iC12sum*CC;
								if( faces[1][q] != f )
								{
									CCI = -CC*multi;
									CCF = CC*C1*multf;
									i_coefs.insert(i_coefs.end(),CCI.data(),CCI.data()+9);
									f_coefs.insert(f_coefs.end(),CCF.data(),CCF.data()+9);
									i_elems.push_back(faces[1][q]);
									f_elems.push_back(faces[1][q]);
								}
							}
						}
						//already have back cell and front cell
						raMatrixMake(f_coefs.data()+0,3,3) = iC12sum*C1sum*C2*multf; //BackCell
						assert(f_elems[0] == f->BackCell());
						raMatrixMake(f_coefs.data()+9,3,3) = -iC12sum*C1*C2sum*multf; //FrontCell
						assert(f_elems[1] == f->FrontCell());
						//i_coefs[0] = iC12sum*C1sum*multi;
						//i_coefs[1] = iC12sum*C2sum*multi;
					}
					else //Boundary face
					{
						B_D.Zero();
						B_N = I;
						B_R.Zero();
						if( tag_BC.isValid() && f.HaveData(tag_BC) )
						{
							f.UnitNormal(n.data());
							GetBC(f.RealArray(tag_BC),n,B_D,B_N,B_R);
						}
						C1 = W[0](rows[0]*3,rows[0]*3+3,rows[0]*3,rows[0]*3+3)*aF;
						C1sum.Zero();
						//TODO
						//iC12sum = (alpha*I+beta*C1).Invert();
						
						tag_f_rhs(f,3,1).Zero();
						tag_i_rhs(f,3,1) = iC12sum*B_R;
						for(int q = 0; q < (int)faces[0].size(); ++q) if( fabs(W[0](rows[0],q)) > 1.0e-12 )
						{
							CC = W[0](rows[0]*3,rows[0]*3+3,q*3,q*3+3);
							if( CC.FrobeniusNorm() > 1.0e-12 )
							{
								//TODO
								//f_elems.push_back(faces[0][q]);
								//f_coefs.push_back(-CC*multfbnd);
								if( faces[0][q] != f )
								{
									//i_elems.push_back(faces[0][q]);
									//i_coefs.push_back(-beta*iC12sum*CC);
								}
								C1sum += CC;
							}
						}
						//TODO
						//f_coefs[0] = C1sum*multfbnd;
						//assert(f_elems[0] == f->BackCell());
						//i_coefs[0] = beta*C1sum*iC12sum;
					}
				} //end of loop over faces
				
			}
            std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
			
			/*
			ttt = Timer();
#if defined(USE_OMP)
#pragma omp parallel
#endif
			{
				
#if defined(USE_OMP)
#pragma omp for
#endif
				for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
				{
					Face face = m->FaceByLocalID(q);
				}
			}
			std::cout << "Construct fluxes and interpolations: " << Timer() - ttt << std::endl;
			*/
			if( m->HaveTag("UVW") ) //Is there a displacement on the mesh?
				tag_UVW = m->GetTag("UVW"); //Get the pressure
			
			if( !tag_UVW.isValid() || !tag_UVW.isDefined(CELL) ) // Pressure was not initialized or was not defined on nodes
			{
				//srand(1); // Randomization
				tag_UVW = m->CreateTag("UVW",DATA_REAL,CELL|FACE,NONE,3); // Create a new tag for the displacement
				for(Mesh::iteratorElement e = m->BeginElement(CELL|FACE); e != m->EndElement(); ++e) //Loop over mesh cells
				{
					//tag_UVW(*e,3,1).Zero();//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
					tag_UVW[*e][0] = (rand()*1.0)/(RAND_MAX*1.0);
					tag_UVW[*e][1] = (rand()*1.0)/(RAND_MAX*1.0);
					tag_UVW[*e][2] = (rand()*1.0)/(RAND_MAX*1.0);
				}
			}
			
			if( !tag_UVW.isDefined(FACE) )
			{
				tag_UVW = m->CreateTag("UVW",DATA_REAL,FACE,NONE,3);
				for(Mesh::iteratorElement e = m->BeginElement(FACE); e != m->EndElement(); ++e) //Loop over mesh cells
				{
					//tag_UVW(*e,3,1).Zero();//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
					tag_UVW[*e][0] = (rand()*1.0)/(RAND_MAX*1.0);
					tag_UVW[*e][1] = (rand()*1.0)/(RAND_MAX*1.0);
					tag_UVW[*e][2] = (rand()*1.0)/(RAND_MAX*1.0);
				}
			}
			m->ExchangeData(tag_UVW,CELL|FACE,0); //Synchronize initial solution with boundary unknowns
			
			if( reference_solution && m->HaveTag("REFERENCE_SOLUTION") )
			{
				TagRealArray tag_R = m->GetTag("REFERENCE_SOLUTION");
				for(Mesh::iteratorElement it = m->BeginElement(CELL|FACE); it != m->EndElement(); ++it)
					tag_UVW(*it,3,1) = tag_R(*it,3,1);
			}
			
			tag_FLUX = m->CreateTag("FLUX",DATA_REAL,FACE,NONE,3);
			tag_S = m->CreateTag("STRESS",DATA_REAL,CELL,NONE,6);
        } //end of initialize data
		

        std::cout << "Initialization done" << std::endl;


        integer nit = 0;
        ttt = Timer();

        { //Main loop for problem solution
            Automatizator aut; // declare class to help manage unknowns
            Automatizator::MakeCurrent(&aut);
			BlockEntry UVW(CELL|FACE);
			UVW.AddTag(tag_UVW);
			aut.RegisterEntry(UVW);
            aut.EnumerateEntries(); //enumerate all primary variables
            std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;

            Residual R("",aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::LockService Locks(aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::AnnotationService Text(aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::Vector Update  ("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for update
            {//Annotate matrix
                for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                {
                    Cell cell = m->CellByLocalID(q);
                    if( cell.GetStatus() != Element::Ghost )
					{
                        Text.SetAnnotation(UVW.Index(cell,0),"Divergence, U " + to_string(UVW.Index(cell,0)));
						Text.SetAnnotation(UVW.Index(cell,1),"Divergence, V " + to_string(UVW.Index(cell,1)));
						Text.SetAnnotation(UVW.Index(cell,2),"Divergence, W " + to_string(UVW.Index(cell,2)));
					}
                }
                for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                {
                    Face face = m->FaceByLocalID(q);
                    if( face.GetStatus() != Element::Ghost )
                    {
                        if( tag_BC.isValid() && face.HaveData(tag_BC) )
						{
                            Text.SetAnnotation(UVW.Index(face,0),"Boundary condition, U " + to_string(UVW.Index(face,0)));
							Text.SetAnnotation(UVW.Index(face,1),"Boundary condition, V " + to_string(UVW.Index(face,1)));
							Text.SetAnnotation(UVW.Index(face,2),"Boundary condition, W " + to_string(UVW.Index(face,2)));
						}
                        else
						{
                            Text.SetAnnotation(UVW.Index(face,0),"Flux continuity, U " + to_string(UVW.Index(face,0)));
							Text.SetAnnotation(UVW.Index(face,1),"Flux continuity, V " + to_string(UVW.Index(face,1)));
							Text.SetAnnotation(UVW.Index(face,2),"Flux continuity, W " + to_string(UVW.Index(face,2)));
						}
                    }
                }
            }

            std::cout << "Matrix was annotated" << std::endl;
			double condest = 0;
            do
			{
                R.Clear(); //clean up the residual
                double tttt = Timer();
                //~ int total = 0, dmp = 0;
#if defined(USE_OMP)
#pragma omp parallel
#endif
				{
					
					vMatrix U,T; //vector of pressure differences and fluxes on faces
					rMatrix Ra(3,3), Rb(3,3), r(3,1), n(3,1);
#if defined(USE_OMP)
#pragma omp for
#endif
					for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
					{
						Face face = m->FaceByLocalID(q);
						tag_FLUX(face,3,1).Zero();
					}
#if defined(USE_OMP)
#pragma omp for
#endif
					for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells
					{
						Cell cell = m->CellByLocalID(q);
						ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
						int NF = (int)faces.size();
						
						raMatrix W = raMatrixMake(tag_W[cell].data(),3*NF,3*NF); //Matrix for gradient
						
						U.Resize(3*NF,1);
						T.Resize(3*NF,1);
						
						for(int k = 0; k < NF; ++k)
							U(3*k,3*(k+1),0,1) = (UVW[faces[k]] - UVW[cell]);
						T = W*U; //fluxes on faces
						tag_S(cell,6,1) = tag_Ws(cell,6,3*NF)*T;
						if( cell.GetStatus() != Element::Ghost )
						{
							for(int k = 0; k < NF; ++k) //loop over faces of current cell
								R[UVW.Index(cell)] += T(3*k,3*(k+1),0,1);
						}
						for(int k = 0; k < NF; ++k) //loop over faces of current cell
						{
							if( faces[k].GetStatus() == Element::Ghost ) continue;
							int index = UVW.Index(faces[k],0);
							Locks.Lock(index);
							if( tag_BC.isValid() && faces[k].HaveData(tag_BC) )
							{
								faces[k]->UnitNormal(n.data());
								GetBC(faces[k].RealArray(tag_BC),n,Ra,Rb,r);
								if( Rb.FrobeniusNorm() ) Rb.Print();
								R[UVW.Index(faces[k])] -= Ra*UVW[faces[k]] + Rb*T(3*k,3*(k+1),0,1)/faces[k].Area() - r;
								
								
								tag_FLUX(faces[k],3,1) += rMatrix(T(3*k,3*(k+1),0,1))/faces[k].Area();
							}
							else
							{
								R[UVW.Index(faces[k])] -= T(3*k,3*(k+1),0,1);
								
								tag_FLUX(faces[k],3,1) += 0.5*rMatrix(T(3*k,3*(k+1),0,1))/faces[k].Area()*(faces[k].FaceOrientedOutside(cell)?1.:-1.);
							}
							Locks.UnLock(index);
						}
					} //end of loop over cells
					
					
					if( print_matrix )
					{
						std::fstream ff;
						ff.open("mat.txt",std::ios::out);
						for(Mesh::iteratorElement it = m->BeginElement(CELL|FACE); it != m->EndElement(); ++it)
						{
							
							ff << "Element " << ElementTypeName(it->GetElementType()) << ":" << it->LocalID();
							if( it->GetElementType() == FACE && it->getAsFace()->Boundary() )
								ff << " boundary";
							ff << std::endl;
							const int bs = 3;
							std::vector<double> VR(3);
							BlockRow BR(bs,bs);
							
							for(int q = 0; q < bs; ++q)
							{
								VR[q] = R[UVW.Index(it->self(),q)].GetValue();
								Sparse::Row & r = R[UVW.Index(it->self(),q)].GetRow();
								//r.Print();
								for(INMOST_DATA_ENUM_TYPE l = 0; l < r.Size(); ++l) if( r.GetValue(l) )
									BR(q,r.GetIndex(l)) = r.GetValue(l);
							}
							
							ff << std::setprecision(5);
							ff << std::scientific;
							std::map< unsigned, rMatrix >::iterator jt;
							for(int q = 0; q < bs; ++q)
							{
								
								ff << std::setw(13) << VR[q] << ": ";
								for(jt = BR.Begin(); jt != BR.End(); ++jt)
								{
									for(int l = 0; l < bs; ++l)
										ff << std::setw(13) << jt->second(q,l) << " ";
									ff << "|";
								}
								ff << std::endl;
							}
							
							rMatrix Sum(bs,bs);
							Sum.Zero();
							for(jt = BR.Begin(); jt != BR.End(); ++jt)
								Sum += jt->second;
							ff << "Sum: " << std::endl;
							for(int q = 0; q < bs; ++q)
							{
								for(int l = 0; l < bs; ++l)
									ff << std::setw(13) << Sum(q,l) << " ";
								ff << std::endl;
							}
							ff << "diagonal block: " << UVW.Index(it->self(),0)/bs << std::endl;
							rMatrix Diag = BR[UVW.Index(it->self(),0)/bs*bs];
							rMatrix iDiag = Diag.PseudoInvert(1.0e-13);
							ff << "Diag: " << std::endl;
							for(int q = 0; q < bs; ++q)
							{
								for(int l = 0; l < bs; ++l)
									ff << std::setw(13) << Diag(q,l) << " ";
								ff << std::endl;
							}
							ff << "iDiag: " << std::endl;
							for(int q = 0; q < bs; ++q)
							{
								for(int l = 0; l < bs; ++l)
									ff << std::setw(13) << iDiag(q,l) << " ";
								ff << std::endl;
							}
							
							rMatrix W(bs,bs), U(bs,bs), S(bs,bs), V(bs,bs);
							Sum.Zero();
							bool diag = false;
							for(jt = BR.Begin(); jt != BR.End(); ++jt)
							{
								int nz = 0;
								bool curdiag = false;
								W = jt->second;
								W.SVD(U,S,V);
								SVD2Eigen(U,S,V);
								ff << "Block " << jt->first << std::endl;
								if( jt->first == UVW.Index(it->self(),0)/bs )
								{
									ff << "diagonal" << std::endl;
									diag = true;
									curdiag = true;
								}
								for(int k = 0; k < bs; ++k) if( fabs(S(k,k)) < 1.0e-4 )
									nz++;
								if( curdiag && nz ) ff << "problem: " << bs-nz << "/" << bs << std::endl;
								ff << "A | U | S | V" << std::endl;
								for(int k = 0; k < bs; ++k)
								{
									for(int l = 0; l < bs; ++l)
										ff << std::setw(13) << W(k,l) << " ";
									ff << "|";
									for(int l = 0; l < bs; ++l)
										ff << std::setw(13) << U(k,l) << " ";
									ff << "|";
									for(int l = 0; l < bs; ++l)
										ff << std::setw(13) << S(k,l) << " ";
									ff << "|";
									for(int l = 0; l < bs; ++l)
										ff << std::setw(13) << V(k,l) << " ";
									ff << std::endl;
								}
							}
							if( !diag ) std::cout << "no diagonal block!" << std::endl;
						}
						ff.close();
					}


					 if( tag_F.isValid() )
					 {
#if defined(USE_OMP)
#pragma omp for
#endif
						 for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
						 {
							 Cell cell = m->CellByLocalID(q);
							 if( cell.GetStatus() == Element::Ghost ) continue;
							 if( cell->HaveData(tag_F) ) R[UVW.Index(cell)] += tag_F(cell,3,1)*cell->Volume();
						 }
					 }
				}
				std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;

                std::cout << "Nonlinear residual: " << R.Norm() << "\t\t" << std::endl;

                if( R.Norm() < 1.0e-4 ) break;
				tttt = Timer();
				
				//R.GetJacobian().Save("A.mtx",&Text);

				//Solver S(Solver::INNER_MPTILU2);
                Solver S(Solver::INNER_MPTILUC);
				//Solver S("superlu");
                S.SetParameter("relative_tolerance", "1.0e-14");
                S.SetParameter("absolute_tolerance", "1.0e-12");
                S.SetParameter("drop_tolerance", "1.0e-4");
                S.SetParameter("reuse_tolerance", "1.0e-6");

                S.SetMatrix(R.GetJacobian());
				
                if( S.Solve(R.GetResidual(),Update) )
                {
#if defined(USE_OMP)
#pragma omp parallel for
#endif
                    for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                    {
                        Cell cell = m->CellByLocalID(q);
						if( cell->GetStatus() == Element::Ghost ) continue;
                        tag_UVW(cell,3,1) -= Update[UVW.Index(cell)];
                    }
#if defined(USE_OMP)
#pragma omp parallel for
#endif
                    for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                    {
                        Face face = m->FaceByLocalID(q);
						if( face->GetStatus() == Element::Ghost ) continue;
                        tag_UVW(face,3,1) -= Update[UVW.Index(face)];
                    }
                    m->ExchangeData(tag_UVW, CELL|FACE, 0);
					
					if( print_niter )
                    {
                        std::stringstream str;
                        str << "iter" << nit;
                        if( m->GetProcessorsNumber() == 1 )
                            str << ".vtk";
                        else
                            str << ".pvtk";
                        m->Save(str.str());
                    }
                    condest = 0;//S.Condest(1.0e-12,500);
                    std::cout << "condition number " << condest << "            " << std::endl;
                }
                else
                {
                    std::cout << "Unable to solve: " << S.ReturnReason() << std::endl;
                    break;
                }
				std::cout << "solved in " << Timer() - tttt << "\t\t\t\t" << std::endl;
                ++nit;
            } while( R.Norm() > 1.0e-4 && nit < 10); //check the residual norm
        }
        std::cout << "Solved problem in " << Timer() - ttt << " seconds with " << nit << " iterations " << std::endl;
        
        
        std::fstream stat;
		
		if( m->GetProcessorRank() == 0 )
		{
			stat.open("stat.csv",std::ios::in);
			if( stat.fail() )
			{
				stat.clear();
				stat.open("stat.csv",std::ios::out);
				stat << "cells; faces; ";
				stat << "uvw_C; uvw_L2; ";
				stat << "uvw_f_C; uvw_f_L2; ";
				stat << "F_C; F_L2; ";
				stat << "S_C; S_L2; ";
				stat << "energy; grid; problem;" << std::endl;
				stat.close();
			}
			else stat.close();
			
			stat.open("stat.csv",std::ios::app);
			
			stat << m->TotalNumberOf(CELL) << "; ";
			stat << m->TotalNumberOf(FACE) << "; ";
		}
        if( m->HaveTag("REFERENCE_SOLUTION") )
        {
            TagReal tag_E = m->CreateTag("ERROR",DATA_REAL,CELL,NONE,1);
            TagRealArray tag_R = m->GetTag("REFERENCE_SOLUTION");
            real C, L2, volume;
            C = L2 = volume = 0.0;
            for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
            {
                Cell cell = m->CellByLocalID(q);
                real err = (tag_UVW(cell,3,1) - tag_R(cell,3,1)).FrobeniusNorm();
                real vol = cell->Volume();
                if( C < fabs(err) ) C = fabs(err);
                L2 += err*err*vol;
                volume += vol;
                tag_E[cell] = err;
            }
			C = m->AggregateMax(C);
			L2 = m->Integrate(L2);
			volume = m->Integrate(volume);
            L2 = sqrt(L2/volume);
            if( m->GetProcessorRank() == 0 ) 
            {
				std::cout << "Error on cells, C-norm " << C << " L2-norm " << L2 << std::endl;
				stat << C << "; " << L2 << "; ";
			}
            C = L2 = volume = 0.0;
            if( tag_R.isDefined(FACE) )
            {
                tag_E = m->CreateTag("ERROR",DATA_REAL,FACE,NONE,1);
                for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                {
                    Face face = m->FaceByLocalID(q);
                    real err = (tag_UVW(face,3,1) - tag_R(face,3,1)).FrobeniusNorm();
                    real vol = (face->BackCell()->Volume() + (face->FrontCell().isValid() ? face->FrontCell()->Volume() : 0))*0.5;
                    if( C < fabs(err) ) C = fabs(err);
                    L2 += err*err*vol;
                    volume += vol;
                    tag_E[face] = err;
                }
				C = m->AggregateMax(C);
				L2 = m->Integrate(L2);
				volume = m->Integrate(volume);
                L2 = sqrt(L2/volume);
                if( m->GetProcessorRank() == 0 ) 
                {
					std::cout << "Error on faces, C-norm " << C << " L2-norm " << L2 << std::endl;
					stat << C << "; " << L2 << "; ";
				}
            }
            else 
            {
				std::cout << "Reference solution was not defined on faces" << std::endl;
				stat << "NA; NA; ";
			}
        }
        else stat << "NA; NA; NA; NA; ";
		
		if( m->HaveTag("REFERENCE_FLUX") )
		{
			TagReal tag_E = m->CreateTag("ERROR_FLUX",DATA_REAL,FACE,NONE,1);
			TagRealArray tag_R = m->GetTag("REFERENCE_FLUX");
			real C, L2, volume;
			C = L2 = volume = 0.0;
			for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
			{
				Face face = m->FaceByLocalID(q);
				real err = (tag_FLUX(face,3,1) - tag_R(face,3,1)).FrobeniusNorm();
				real vol = face->BackCell()->Volume()/(real)face->BackCell().nbAdjElements(FACE);
				if( face->FrontCell().isValid() )
					vol += face->FrontCell()->Volume()/(real)face->FrontCell().nbAdjElements(FACE);
				if( C < fabs(err) ) C = fabs(err);
				L2 += err*err*vol;
				volume += vol;
				tag_E[face] = err;
			}
			C = m->AggregateMax(C);
			L2 = m->Integrate(L2);
			volume = m->Integrate(volume);
			L2 = sqrt(L2/volume);
			if( m->GetProcessorRank() == 0 ) 
			{
				std::cout << "Error of fluxes, C-norm " << C << " L2-norm " << L2 << std::endl;
				stat << C << "; " << L2 << "; ";
			}
		}
		else stat << "NA; NA; ";
		if( m->HaveTag("REFERENCE_STRESS") )
		{
			TagReal tag_E = m->CreateTag("ERROR_STRESS",DATA_REAL,CELL,NONE,1);
			TagRealArray tag_R = m->GetTag("REFERENCE_STRESS");
			real C, L2, volume;
			C = L2 = volume = 0.0;
			for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
            {
                Cell cell = m->CellByLocalID(q);
                real err = (tag_S(cell,6,1) - tag_R(cell,6,1)).FrobeniusNorm();
                real vol = cell->Volume();
				if( C < fabs(err) ) C = fabs(err);
				L2 += err*err*vol;
				volume += vol;
				tag_E[cell] = err;
			}
			C = m->AggregateMax(C);
			L2 = m->Integrate(L2);
			volume = m->Integrate(volume);
			L2 = sqrt(L2/volume);
			if( m->GetProcessorRank() == 0 ) 
			{
				std::cout << "Error of stress C-norm " << C << " L2-norm " << L2 << std::endl;
				stat << C << "; " << L2 << "; ";
			}
		}
		else stat << "NA; NA; ";
		
		//compute energy by applying inverse of elastic tensor on stress to get strain
		{
			real energy = 0;
#if defined(USE_OMP)
#pragma omp parallel reduction(+:energy)
#endif
			{
				rMatrix s(6,1), e(6,1);
#if defined(USE_OMP)
#pragma omp for
#endif
				for(int i = 0; i < m->CellLastLocalID(); ++i) if( m->isValidCell(i) )
				{
					Cell c1 = m->CellByLocalID(i);
					//extract stress
					s = rMatrix::FromVector(tag_S[c1].data(),6);
					e = rMatrix::FromTensor(tag_C[c1].data(),21,6).Solve(s);
					//calculate energy
					energy += e.DotProduct(s)*c1.Volume();
					
				}
			}
			energy *= 0.5;
			std::cout << "Energy: " << energy << std::endl;
			stat << energy << "; ";
		}
		
		if( m->HaveTag("GRIDNAME") )
		{
			Tag tagname = m->GetTag("GRIDNAME");
			Storage::bulk_array name = m->self().BulkArray(tagname);
			stat << std::string(name.begin(),name.end()) << ";";
		}
		else stat << "NA;";
		if( m->HaveTag("PROBLEMNAME") )
		{
			Tag tagname = m->GetTag("PROBLEMNAME");
			Storage::bulk_array name = m->self().BulkArray(tagname);
			stat << std::string(name.begin(),name.end()) << ";";
		}
		else stat << "NA;";
		stat << std::endl;
		
		stat.close();
		

        if( m->GetProcessorsNumber() == 1 )
            m->Save("out.vtk");
        else
            m->Save("out.pvtk");
		m->Save("out.pmf");

        delete m; //clean up the mesh
    }
    else
    {
        std::cout << argv[0] << " mesh_file" << std::endl;
    }

#if defined(USE_PARTITIONER)
    Partitioner::Finalize(); // Finalize the partitioner activity
#endif
    Solver::Finalize(); // Finalize solver and close MPI activity
    return 0;
}