//============================================================================
//                                                                          // 
//           2D  Curvature Driven Grain Growth                              //
//                                                                          //
//                                                                          // 
//               Shlomo Ta'asan,  Eva Eggeling                              // 
//            Department of Mathematical Sciences                           //
//                Carnegie Mellon University                                //
//                   Pittsburgh PA 15213                                    //
//                                                                          //
//                    Ph: 412 268 9560                                      //
//                    email: eggeling@andrew.cmu.edu                        //
//                                                                          //
//                                                                          // 
//    Bug Report:                                                           // 
//                    eggeling@andrew.cmu.edu                               //
//                                                                          //
//                    version 0.08 January 2008, E. Eggeling                // 
//    References:                                                           //
//                                                                          //
//                                                                          //
//============================================================================

#include <list>
#include <vector>
#include <algorithm>
#include <string>
#include <fstream>
#include <iostream>
#include <sstream>
#include <ctime>
#include <iomanip>
#include <cmath>
#include <time.h>
#include <stdio.h>

using namespace std;

int DEBUG;
int WRITE;
int PRINT;
int READ_ORIENTATIONS;
double PI;

double mobilityN;
double mobilityTJ;
double lineMinLength; 
double hMin;
double deltaT;
double mindeltaT;

double current_time;

int triple_ite;                 // Number of triple junctions iterations

int print_times;               // output will be printed print_times module


double epsLine;                // Epsilon to remove shortlines // Later updated to Largest line /20.0
double ratio_epsline;            // update epsLine to Largest line / ratio
double ratio_lines;            // hmin = min(maxLen/ratio_lines,minLen/2.0)
double ratio_split;            // New line is going to be hMin * ratio_split in splitUnstableVerices
//double split_constant = 0.5;           // move nodes in line to exp(-k/split_constant) 

int mesh_factor;
double deltaT_const;      // DeltaT = hmin^2 * DeltaT_cons

vector<int> EmptyLine;    // this is used to keep track of empty loc in the 
                        // the vLine vector so insertion of new lines will 
                        // take place in unsused indices. 
vector<int> EmptyVertex;  // same for vList. 
                        // each time a Line os Vertex is removed, its index 
                        // is inserted in this list. Each time a new one 
                        // is created we take out the most recent one.
vector<int> EmptyGrain;   // similar. 


//===========================================================================
string toString(int x)
{
    ostringstream str;
    str << x << flush;
    return str.str();
}

// ====================== DATA STRUCTURES   ==============================
struct Vec2D
{
    double x;
    double y;

    Vec2D(){x = 0.0; y = 0.0;}
    Vec2D(double z1, double z2) {x=z1; y = z2;}
    Vec2D operator*(double a){Vec2D t; t.x = x*a; t.y = y*a; return t;}
    Vec2D operator/(double a){Vec2D t; t.x = x/a; t.y = y/a; return t;}
    Vec2D operator+(Vec2D a){ 
      Vec2D t; 
      t.x = x + a.x; t.y = y + a.y; 

      if( t.x <= -1.0) t.x += 2.0;
      else if( t.x > 1.0 ) t.x -= 2.0; 
      if( t.y <= -1.0) t.y += 2.0;
      else if( t.y > 1.0 ) t.y -= 2.0; 

      return t;
    }
    Vec2D operator-(Vec2D a){ 
      Vec2D t; 
      t.x = x - a.x; t.y = y - a.y; 

      if( t.x <= -1.0) t.x += 2.0;
      else if( t.x > 1.0 ) t.x -= 2.0; 
      if( t.y <= -1.0) t.y += 2.0;
      else if( t.y > 1.0 ) t.y -= 2.0; 

//       if( t.x <= -1.) t.x += 2.0;
//       else if( t.x > 1. ) t.x -= 2.0; 
//       if( t.y <= -1.0) t.y += 2.0;
//       else if( t.y > 1.0 ) t.y -= 2.0; 
      return t;
    }
    double length(){return sqrt(x*x + y*y);}
    double crossProduct(Vec2D a){ return (x * a.y - y * a.x);}
    double dotProduct(Vec2D a){ return (x * a.x + y * a.y);}
    Vec2D perp(){return Vec2D(-y, x);}
    void print(){cout << x << " " << y << endl; }
};

struct Line
{
    int id;                // id = -1 for dead Lines
    int vertexI;
    int vertexJ;
    int grainI;
    int grainJ;
    double alphaI;
    double alphaJ;
    list<Vec2D> node;      // node position. first note coincide with vertexI
    list<Vec2D> nodeVel;   // the node velocity
    list<Vec2D> normal;    // at mid nodes. 
    double length; 
    double dLdt;  // to determine if short lines shrink or expand
    //   functions
    bool Line::operator<(const Line &a);
    void reverse();
    Line Line::add(Line L);
    void Line::kill();
    void defineLine(int, int, int, int, int, Vec2D, Vec2D, Vec2D);
    
};

struct Vertex
{
    Vec2D pos;
    Vec2D vel;
    int isInside;        // = 1 for inside vertex, 0 outside (in IOqv.o list)
    int id;              // id = -1, for deleted vertices 
    vector<int> line;     // indices of lines
    vector<int> grain;    // indices of grains. grain 0: between line 0 & 1
    void kill();
    Vertex(){pos.x = pos.y = 0.0; vel = pos; isInside = 1; id = -1; }
    void defineVertex(int, int, int, int, int, int, int, int, Vec2D);
    int Vertex::otherVertex( Line l)
	{
	    //int vv=this.id;
	    int oV; // other ending Vertex of this line l 
	    if (l.vertexI == id)
		oV = l.vertexJ;
	    else if (l.vertexJ == id)
		oV= l.vertexI;
	    else
	    {
		cout << " bug in otherVertex " ;
		exit(0);
	    }
	    return oV;
	}
};

struct Grain
{
    int id;                // id = -1 for dead grains
    double alpha;
    vector<int> vertex; 
    vector<int> line;  
    void Grain::kill();
  //void addLineAndVertex(int, int, vector<Line> &);
};

struct Pair
{
  int i;
  int o;
  Pair(int a, int b){i=a; o=b;}
  bool Pair::operator==(const Pair &p);
  Pair flip(){return Pair(o,i);}
};
//=======================================================================
bool Pair::operator==(const Pair &p)
{
// this operator is used to find easier in IOqv list of Pairs
    if (i == p.i)
	return true;
    else
	return false;
}

//========================================================================
vector<Pair> vecFlip( vector<Pair> &p)
{
  vector<Pair> q;
  vector<Pair>::iterator k;
  for(k=p.begin(); k < p.end(); k++) 
    q.push_back( (*k).flip());

  return q;
}
//=======================================================================
bool Line::operator<(const Line &a)
{
    // this is used to sort lines in order to remove redundancy at 
    // initial stage of defining the lines from grains and vertices.

    if( vertexI < a.vertexI ) 
	return true;
    else if(vertexI == a.vertexI && vertexJ <= a.vertexJ) 
	return true;

    return false;
}

//======================================================================
void Line::reverse()
{
    int tmpV = vertexI; 
    vertexI = vertexJ;
    vertexJ = tmpV;
    node.reverse();
    nodeVel.reverse();
    normal.reverse();
}
Line Line::add(Line L)
{
  
}
//=======================================================================
void Line::kill()
{
    EmptyLine.push_back(id);
    id = -1; 
    node.clear();
    nodeVel.clear();
    normal.clear();
    vertexI = -1;
    vertexJ = -1;
    grainI = -1;
    grainJ = -1;
}
//======================================================================
void Grain::kill()
{
    EmptyGrain.push_back(id);
    id = -1;
    vertex.clear();
    line.clear();
    alpha = -1;
}
//=======================================================================
void Vertex::kill()
{
    EmptyVertex.push_back(id);
    id = -1; 
    line.clear();
    grain.clear();
}
//===========================================================================
void Vertex::defineVertex(int vID, int isIn, int g0, int g1, int g2, 
			  int l0, int l1, int l2, Vec2D position)
{
    id = vID;
    if( !grain.empty() )grain.clear();
    grain.push_back(g0);
    grain.push_back(g1);
    grain.push_back(g2);
    if( !line.empty() ) line.clear();
    line.push_back(l0);
    line.push_back(l1);
    line.push_back(l2);
    pos = position;  
    isInside = isIn;
}
//============================================================================
void Line::defineLine(int ID, int vI, int vJ, int gI, int gJ, Vec2D x0, Vec2D x1, Vec2D x2)
{
    id = ID;
    vertexI = vI;
    vertexJ = vJ;
    grainI = gI;
    grainJ = gJ;
    node.clear();
    node.push_back(x0);
    node.push_back(x1);
    node.push_back(x2);

    // Initialize Normals for new line           // Added by Guillermo 
 
    Vec2D v1 = x1 - x0;
    Vec2D n1 = v1.perp();
    n1 = n1 / n1.length();
    normal.push_back(n1);

    Vec2D v2 = x2 - x1;
    Vec2D n2 = v2.perp();
    n2 = n2 / n2.length();
    normal.push_back(n2);

    // Initialize Velocities for new Line      // Added By Guillermo

    nodeVel.push_back(Vec2D(0.0,0.0));
    nodeVel.push_back(Vec2D(0.0,0.0));
    nodeVel.push_back(Vec2D(0.0,0.0));


}


//=====================  functions prototype ==============================

void readData(string, vector<Vertex> &, vector<Grain> &, vector<Pair> &, vector<Pair> &);
void defineDataStructures(vector<Vertex> &, vector<Grain> &, vector<Line> &);
void finalizeGrains(vector<Vertex> &, vector<Grain> &, vector<Line> &);
void finalizeVertices(vector<Vertex> &, vector<Grain> &, vector<Line> &, vector<Pair> &  );
    
void checkConsist(vector<Vertex>, vector<Grain>, vector<Line> );

void periodicStructur(vector<Vertex> &, vector<Grain> &, vector<Line> &, 
		      vector<Pair> &, vector<Pair> &);  
void updateTimeStep(vector<Vertex> &, vector<Line> &);
void adjustNodes(vector<Line> &, Vec2D, int, int);
void recalcVelo(vector<Vertex> &, vector<Line> &);
void printLines(vector<Line>&);
void printLine(Line);
void printLines(list<Line>&);
void printGrains(vector<Grain>&);
void printGrain(Grain);
void printVertices(vector<Vertex>&);
void printVertex(Vertex);

void printGrainsGNU(vector<Grain> &, vector<Vertex> &  ) ;
void printLinesGNU(vector<Line> &, vector<Vertex> & ) ;
void findGrain(vector<Vertex> &, vector<Grain> &,int,int,int &,int &,int,vector<int> :: iterator &);
void splitVertex(Vertex & , vector<Vertex> &, vector<Grain> &, vector<Line> &, 
		 vector<Pair> &, Vertex &, vector<int> &); 
vector<int>::iterator findNextLine(vector<Vertex> &, vector<Grain> &, vector<Line> &,int,int,int, 
				   vector<int> :: iterator &);
//void getMasterAndTwins(int, vector<Vertex>, vector<Pair>, int, vector<int> & );
void remove2SidedGrain(int, vector<Grain> &, vector<Vertex> &, vector<Line> &, vector<int> & );
void Treat2LineVertex(int, vector<Vertex> &, vector<Line> &, vector<Grain> &, vector<int> &);
void spaceNodesEqually(vector<Line> &, int &);
void removeLineAndVertexFromGrain(vector<Grain> &, vector<Line> &, int , int , int );
vector<Pair> vecFlip( vector<Pair> &);

int findPosition(int , vector<int> );
int findLineWithVIVJ(vector<Line> ,int, int) ;
//===================  PHYSICAL PARAMETERS =============================

double sigma(double a1, double a2, Vec2D n)            // a: misorientation  
{
    return 1.0;
}
//======================================================================
double Dsigma(double a1, double a2, Vec2D n)
{
    return 0.0;
}
//======================================================================
double mobility(double a1, double a2) // killed the normal here!!!
{
    return 1.0;
}
//=====================================================================
double TJmobility(double a1, double a2, double a3, 
		  Vec2D n1, Vec2D n2, Vec2D n3)
{
    return 0.0;
}
//=================Boundary Condition================================
// Check that all vectors have the correct direction // 
// created by Guillermo 2/13/06
Vec2D correctDirection ( Vec2D v, int option)
{
    double bound = 1.0;
    double max_len = 0.9;

    vector<Vec2D> vec;
    int done = 0;
    vec.push_back(Vec2D(0,0));               vec.push_back(Vec2D(2*bound,0));  
    vec.push_back(Vec2D(-2*bound,0));        vec.push_back(Vec2D(0,2*bound));  
    vec.push_back(Vec2D(0,-2*bound));        vec.push_back(Vec2D(2*bound,2*bound)); 
    vec.push_back(Vec2D(-2*bound,-2*bound)); vec.push_back(Vec2D(-2*bound,2*bound)); 
    vec.push_back(Vec2D(2*bound,-2*bound));   

    for (int j = 0; j < vec.size(); j++)
    {
	if ( abs((v + vec[j]).x) < bound && abs((v + vec[j]).y) < bound && option == 1 )  
	{
	    v = v + vec[j];
	    done = 1;
	    break;
	}
	if ( (v + vec[j]).length() < max_len && option == 2 ) 
	{
	    v = v + vec[j];
	    done = 1;
	    break;
	}

    }
    if (done == 0)
    {
	cout << "Bug in correctDirection" << endl;
	exit(0);
    }
    return v;
}
//======================== INPUT FUNCTIONS ===========================
void readParameters(string fName)
{
     ifstream inFile(fName.c_str(), ios::in);
    if( inFile == NULL ) 
    {
	cout << " could not open " << fName << endl;
	exit(0);
    }
    string s;    

    inFile >> DEBUG;                getline(inFile,s);
    inFile >> WRITE;                getline(inFile,s);
    inFile >> PRINT;                getline(inFile,s);
    inFile >> READ_ORIENTATIONS;    getline(inFile,s);
    inFile >> PI;                   getline(inFile,s);

    inFile >> mobilityN;           getline(inFile,s);
    inFile >> mobilityTJ;           getline(inFile,s);
    inFile >> lineMinLength;        getline(inFile,s);
    inFile >> hMin;                 getline(inFile,s);

    inFile >> deltaT;               getline(inFile,s);
    inFile >> mindeltaT;            getline(inFile,s);

    inFile >> triple_ite;           getline(inFile,s);     // Number of triple junctions iterations
    
    inFile >> print_times;          getline(inFile,s);     // output will be printed print_times module

    inFile >> epsLine;             getline(inFile,s);     // Epsilon to remove shortlines 
                                                           // Later updated to Largest line /ratio_epsline
    inFile >> ratio_epsline;        getline(inFile,s);         // update epsLine to Largest line / ratio
 
    inFile >> ratio_lines;          getline(inFile,s);     // hmin = min(maxLen/ratio_lines,minLen/2.0)

    inFile >> ratio_split;          getline(inFile,s);     // New line is going to be hMin * ratio_split 
                                                           //in splitUnstableVerices
    //inFile >> splitnnnns_constant;     getline(inFile,s);     // move nodes in line to exp(-k/split_constant) 

    inFile >> mesh_factor;          getline(inFile,s);     // Nodes are distributed by a factor of meshfactor

    inFile >> deltaT_const;         getline(inFile,s);     // DeltaT = hmin^2 * DeltaT_cons

    if (DEBUG)
    {
	cout << DEBUG                 << endl;
	cout << WRITE                 << endl;
	cout << PRINT                 << endl;
	cout << READ_ORIENTATIONS     << endl;
	cout << PI                    << endl;

	cout << mobilityN           << endl;
	cout << lineMinLength         << endl;
	cout << hMin <<  endl;

	cout << deltaT               << endl;
	cout << mindeltaT            << endl; 

	cout << triple_ite           << endl;     // Number of triple junctions iterations
    
	cout << print_times          << endl;     // output will be printed print_times module

	cout << epsLine             << endl;     // Epsilon to remove shortlines 
	// Later updated to Largest line /ratio_epsline
	cout << ratio_epsline        << endl;         // update epsLine to Largest line / ratio
 
	cout << ratio_lines          << endl;     // hmin = min(maxLen/ratio_lines,minLen/2.0)

	cout << ratio_split          << endl;     // New line is going to be hMin * ratio_split 
	//in splitUnstableVerices
	//cout << split_constant;     <<endl;     // move nodes in line to exp(-k/split_constant) 

	cout << mesh_factor          << endl;     // Nodes are distributed by a factor of meshfactor

	cout << deltaT_const         << endl;     // DeltaT = hmin^2 * DeltaT_cons

    }
    inFile.close();
}
//========================================================================
void readData(string fName, vector<Vertex> &vList, vector<Grain> &gList, 
	      vector<Pair> &IOgrCenter, vector<Pair> &IOqv)
{
    int nVertex, nGrain; 
    double x,y;

    ifstream inFile(fName.c_str(), ios::in);
    if( inFile == NULL ) 
    {
	cout << " could not open " << fName << endl;
	exit(0);
    }
    int n, ndim;
    inFile >> ndim ;
    if( ndim != 2 ) 
    {
	cout << "It seems that the data is not for 2D " << endl;
	exit(1);
    }

    inFile >> nVertex >> nGrain >> n;

    // read the vertices

    for(int i=0; i<nVertex; i++)
    {
	Vertex tmpVertex;
	tmpVertex.id = i;
	inFile >> x >> y;
	tmpVertex.pos = Vec2D(x,y);
	vList.push_back(tmpVertex);
    }

    // read the grains (a list of vertices)
    int gInd;
    for(int i=0; i<nGrain; i++)
    {
	Grain tmpGrain;
	tmpGrain.id = i;
	if( !READ_ORIENTATIONS ) 
	{
	    double r = (   (double)rand() / ((double)(RAND_MAX)+(double)(1)) ); // r in [0;1)
	    //tmpGrain.alpha = 2. * PI * ( rand() % 100 ); // randomize orientations
	    tmpGrain.alpha = 0.5 * PI * r ; // random orientations
	    // cout << tmpGrain.alpha << " \t" << (tmpGrain.alpha)*180/PI << endl;
	}
	else
	    inFile >> tmpGrain.alpha; 

	inFile >> n; 
	for(int j=0; j<n; j++)
	{
	    inFile >> gInd;
	    (tmpGrain.vertex).push_back(gInd);
	}
	gList.push_back(tmpGrain);
    }
    
    //---  read IO list of grain centers
    int nGrCenter;
    int IOin, IOout;
    
    inFile >> nGrCenter;
    for (int i=0; i< nGrCenter; i++)
    {
 	inFile >> IOin >> IOout;
 	Pair p = Pair(IOin, IOout) ;
 	IOgrCenter.push_back(p) ;
    }
    //---  read IO list of vertices
    int nqv;
    inFile >> nqv ;
    for (int i=0; i< nqv; i++)
    {
 	inFile >> IOin >> IOout;
 	Pair q = Pair(IOin, IOout) ;
 	IOqv.push_back(q) ;
    }
    
    inFile.close();


    //cout << " -----------check isInside of vertices------ " << endl;
     for (int i=0; i<IOqv.size(); i++)
     	vList[IOqv[i].o].isInside = 0;
        
}
//======================================================================


void checkIniConf(vector<Grain> gList,vector<Vertex> vList)
{
    int marker;
    for (int k=0; k<vList.size(); k++ )
    {
	marker = 0;
	//cout << " check vertex " << vList[k].id << endl;
	for (int i=0; i<gList.size(); i++)
	    for (int j=0; j<gList[i].vertex.size() ; j++)
		if (vList[k].id == vList[gList[i].vertex[j]].id)
		    marker = 1;
	if (marker == 0)
	{
	    cout << " vertex " << vList[k].id << " not found in gList.vertex" << endl;
	    cout << "not used in current configuration, problem." << endl;
	    cout << " Check IniConfig failed!!! " << endl;
	    exit(0);
	}
    }
    cout << "check IniConfig passed " << endl;
}

//======================================================================

void defineDataStructures(vector<Vertex> &vList, vector<Grain> &gList, 
			  vector<Line> &lList, vector<Pair> &IOqv)
{
    vector<Grain>::iterator g;
    vector<Vertex>:: iterator ver;
    vector<int>::iterator v;
    list<Line>::iterator l;
    Line tmpLine;
    cout.precision(16);

    //--- start by asisgning pairs of successive vertices in 
    // gList.vertex to lines, then remove duplications. 

    list<Line> myList;
    for (g=gList.begin(); g != gList.end(); g++)
    {
	int gId = (*g).id;
	double gAlpha = (*g).alpha;
	int gSize = ((*g).vertex).size();
	int *V = new int[gSize];
	int j = 0;
	for(v=((*g).vertex).begin(); v != ((*g).vertex).end(); v++)
	    V[j++]= (*v);

	Line line;
	line.grainI = gId;
	line.alphaI = gAlpha;
	line.grainJ = -1;
	line.alphaJ = -1.0;
	for(int j = 0; j<gSize; j++) 
	{	 
	    line.id = -1;
	    line.vertexI = min(V[j], V[(j+1)%gSize]);
	    line.vertexJ = max(V[j], V[(j+1)%gSize]);
	    //cout << line.vertexI << " " << line.vertexJ << endl;
	    myList.push_back(line);
	}
	delete [] V;
    }
    myList.sort();

    //--- remove duplicates; most lines appears twice
    int I = -1;   int J = -1;
    int lineNo = 0;
    for (l=myList.begin(); l != myList.end(); l++) 
    {
	int vI = (*l).vertexI; 
	int vJ = (*l).vertexJ;
	if( !((vI == I && vJ == J) || (vI == J && vJ == I)) )
	{
	    (*l).id = lineNo++;
	    lList.push_back(*l);
	    I = vI; 
	    J = vJ;
	} 
	else // need to add gId for this line
	{
	    (lList.back()).grainJ = (*l).grainI;  
	    (lList.back()).alphaJ = (*l).alphaI;
	}
    }
    if( lList.size() + 1 != gList.size() + vList.size() ) 
    {
	cout << endl << endl;
	cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
	cout << endl;
	cout << " There seems to be an error in the input file "<< endl;
	cout << " Number of grains, vertices and lines does not satisfy" << endl;
	cout << " The Euler formula:     V + C = L + 1 " << endl;
	cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
	cout << endl << endl;
	exit(0);
    }   
    myList.clear();

    // at this stage we have defined all lines in terms of their 2 end vertices
    //
    // now finalizing connections between data structures
    
    //printVertices(vList);
    //printGrains(gList);
    finalizeGrains(vList, gList, lList); 
    //printGrainsGNU(gList,vList);
    //printLinesGNU(lList,vList);
    finalizeVertices(vList, gList, lList, IOqv);
    
    

}
//========================================================================
int findLineWithVIVJ(vector<Line> lList, int vI, int vJ)
{
    int line=-5;
    int v1, v2;
    for (int i=0; i<lList.size(); i++)
    {
	if (lList[i].id > -1)
	{
	    v1 = lList[i].vertexI;
	    v2 = lList[i].vertexJ;
	    if ((vI == v1 && vJ == v2) || (vI == v2 && vJ == v1) )
	    {// found the right line
		line = lList[i].id;
	    }
	}
    }
    return line;
} 
//========================================================================

void periodicStructur(vector<Vertex> &vList, vector<Grain> &gList, 
		      vector<Line> &lList, vector<Pair> &IOqv, 
		      vector<Pair> &IOgrCenter)  
{       
    // get structure lines, grains, vertices before periodic Structure
 


    //--- fix lList 
    
    // check all vertices in IOqv and replace outside neighbour vertices
    // by corresponding inside partner
    int inVertex, v;
    vector<Pair> OIqv=vecFlip(IOqv);

    for (int j=0; j<IOqv.size(); j++)
    {
	inVertex = IOqv[j].i; //inside vertex 
	//cout << " inVertex is " << vList[inVertex].id << endl;
	//printVertex(vList[inVertex]);
	for (int k=0; k<vList[inVertex].line.size(); k++)
	{
	    int thisLine = vList[inVertex].line[k];
	    //cout << " thisLine " << thisLine ; 
	    //printLine(lList[thisLine]);
            // get ending vertex v of lines of inVertex 
	    if ( vList[inVertex].id == lList[thisLine].vertexI )
		v = lList[thisLine].vertexJ;
	    else if (vList[inVertex].id == lList[thisLine].vertexJ )
		v = lList[thisLine].vertexI;
	    else
	    {
		cout << " BUG 1 in periodic Structur \n";
		exit(0);
	    }
	    
	    vector<Pair>::iterator qq, qBegin=OIqv.begin();
	    
	    if (vList[v].isInside == 0) //replace v by .i partner of IOqv list
	    {
		Pair pp= Pair(vList[v].id,0);
		qq = find(OIqv.begin(), OIqv.end(), pp);
 		if ( (qq) != OIqv.end() )
		{
		    int vNew = (*qq).o; // inside partner of qq
		    //update line
		    if (lList[thisLine].vertexI == v)
			lList[thisLine].vertexI = vNew;
		    else if (lList[thisLine].vertexJ == v )
			lList[thisLine].vertexJ = vNew;
		    else
		    {
			cout << "BUG 2 in periodic Structure \n";
			exit(0);
		    }   
		    
		    
		    if (lList[thisLine].vertexI > lList[thisLine].vertexJ)
		    {// change order of vertices for sorting algorithm
			int temp = lList[thisLine].vertexJ;
			lList[thisLine].vertexJ = lList[thisLine].vertexI;
			lList[thisLine].vertexI = temp ;
		    }
		}
	    }
	}
    }
 
    vector<Pair> trans ;
    
    // check lines with two outside vertices, 
    // find corresponding inside-line 
    // create trans
    
    //--- start checking for outside lines 
    for (int i=0; i<IOqv.size(); i++)
    {
	int vI,vJ;
	int insideLine;
	for (int j=0; j< vList[IOqv[i].o].line.size(); j++)
	{
	    int line = vList[IOqv[i].o].line[j];
	    vI = lList[line].vertexI;
	    if (vI ==  vList[IOqv[i].o].id)
		vJ = lList[line].vertexJ;
	    else
		vJ = lList[line].vertexI; // endpoint
	    // find vJ in IOqv.o
	    vector<Pair>::iterator qq;
	    vector<Pair> OIqv = vecFlip(IOqv);
	    Pair p = Pair(vJ,0);
	    qq = find(OIqv.begin(), OIqv.end(), p);
	    if (qq != OIqv.end() )
	    { 		
		int v1,v2;
		v1 = vList[IOqv[i].i].id; // inside partner 1 
		v2 = (*qq).o;             // inside partner 2
	
		// line has to be replaced later by 'partner line'
		// find line between v1 and v2
		insideLine = -5;
		insideLine = findLineWithVIVJ(lList, v1, v2);
		// add to trans
		if (insideLine != -5 )
		{ 
		    //cout << " found inside line " << insideLine << endl;
		    Pair lineIO = Pair(line, insideLine ); // CHECK ORDER 
		    trans.push_back(lineIO);
		}
		else
		{
		    cout << "could not find inside line, problem " << endl;
		    exit(0);
		}
	    }
	}
	
    }
     
    // create vector tmp, ordered in first component     
    vector<Pair> tmp1; 
    for (int k=0; k<lList.size(); k++)
    {
	Pair p = Pair(-1,-1);
	tmp1.push_back(p);
    }
    
    for (int ll=0; ll<trans.size(); ll++)
    {
	tmp1[trans[ll].i].i = trans[ll].i;
	tmp1[trans[ll].i].o = trans[ll].o;
    }
//     cout << " tmp1: size =  " << tmp1.size() << endl;
//     for (int i=0; i<tmp1.size(); i++)
//     {
// 	cout <<  tmp1[i].i << " " << tmp1[i].o << endl;
//     }

    
    //-- fix grains step 1,  based on tmp1 (only outside lines)
    //cout << "------- update grains, step 1 ------ " << endl;
    for (int k=0; k<gList.size(); k++)
    	if (gList[k].id > -1)
	    for (int j=0; j<(gList[k].line).size() ; j++)
		if ( tmp1[gList[k].line[j]].o != -1 )
		    gList[k].line[j] =  tmp1[gList[k].line[j]].o;
    
     // sort myList
    list<Line> myList;
    vector<Line> newlList;
    list<Line>::iterator l,lprev;
    for (int ll=0; ll<lList.size(); ll++)
      myList.push_back(lList[ll]);
    myList.sort();

    // myList includes only inside vertices
    // it is sorted, but has duplicate lines
        
    //----------------------------------------------------------------------------

    // trans2: transformation to replace old by new lines
    // to be applied to vector.lines and grain.lines 
    vector<Pair> trans2 ;

    int vI, vJ;
    int I = -1;   int J = -1;
    int lineNo = 0;
    
    
    for (l=myList.begin(); l!= myList.end(); l++) 
    {
	vI = (*l).vertexI; 
	vJ = (*l).vertexJ;
	
	if ( (vI == I && vJ == J) || (vI == J && vJ == I)  ) // a duplicate
	{
	    //trans2[lineNo].i = (*l).id;
	    lprev = l;
	    lprev--;
	    Pair p=Pair((*l).id, (*lprev).id);
	    trans2.push_back( p ); 
	    //cout << trans[lineNo].i<<" " << trans[lineNo].o << endl;
	    I = vI; 
	    J = vJ;

            // ---adjust grains of the duplicate lines
	    if ( (*l).grainI > -1 && (*l).grainJ > -1) // both grains are real
	    {
		(*lprev).grainI = (*l).grainI;
		(*lprev).grainJ = (*l).grainJ;
	    }
	    else if ( (*lprev).grainI > -1 && (*lprev).grainJ > -1 )
	    {// one of the grains of line l is -1
		
		(*l).grainI = (*lprev).grainI;
		(*l).grainJ = (*lprev).grainJ;
	    }
	    else
	    {//gather the !=-1 grains and replace 

		
		int g1=-5, g2 = -5;
		if ((*l).grainI > -1 )
		    g1 = (*l).grainI;
		if ( (*l).grainJ > -1 )
		    g1 = (*l).grainJ;
		if ((*lprev).grainI > -1)
		    g2 = (*lprev).grainI;
		else
		    g2 = (*lprev).grainJ;
		(*lprev).grainI = g1;
		(*lprev).grainJ = g2;
		(*l).grainI = g1;
		(*l).grainJ = g2;
	    }
	} 
	else
	{ //add same number for .i and .o 
	    //cout << " not a double line, add same nr " << endl;
	    Pair p = Pair((*l).id,(*l).id);
	    trans2.push_back(p);
	    I = vI; 
	    J = vJ;
	}
	
	lineNo++;			
    }
    //----------------------------------------------------------------------------
     
    // create vector tmp2, ordered in first component
   
    vector<Pair> tmp2; 
    for (int k=0; k<myList.size(); k++)
    {
	Pair p = Pair(-1,-1);
	tmp2.push_back(p);
    }
    
    for (int ll=0; ll<trans2.size(); ll++)
    {
	tmp2[trans2[ll].i].i = trans2[ll].i;
	tmp2[trans2[ll].i].o = trans2[ll].o;
    }

    //----------------------------------------------------------------------------
    //--- remove duplicates; most lines appears twice
    I = -1;   
    J = -1;
    lineNo = 0;
    int largeLineNumber =  myList.size();
    for (l=myList.begin(); l != myList.end(); l++) 
    {
	vI = (*l).vertexI; 
	vJ = (*l).vertexJ;
	
	if (!(vList[vI].isInside == 0 && vList[vJ].isInside == 0) && 
	    ( !((vI == I && vJ == J) || (vI == J && vJ == I))) ) // no duplicate
	{
	    newlList.push_back(*l);
	    I = vI; 
	    J = vJ;
	} 
    }
    
    myList.clear();
    lList.clear();
    for (int ll=0; ll<newlList.size(); ll++)
	lList.push_back(newlList[ll]);
    newlList.clear();


    //------------------------------fix grains ------------------------------------
    //--- update vertices in gList.vertex, based on IOqv list 
    for (int k=0; k<gList.size(); k++)
    	if (gList[k].id > -1) 
	    for (int j=0; j<gList[k].vertex.size(); j++)
	    {
		int v=gList[k].vertex[j];
		if (vList[v].isInside == 0) // outside vertex
		{
		    vector<Pair>::iterator qq;
		    Pair pp= Pair(vList[v].id,0);
		    qq = find(OIqv.begin(), OIqv.end(), pp);
		    if (qq != OIqv.end() )
			gList[k].vertex[j] = (*qq).o; 
		}		
	    }
    
    
    
    //--update lines in gList.lines, based on tmp2   
    for (int k=0; k<gList.size(); k++)
    	if (gList[k].id > -1)
	    for (int j=0; j<(gList[k].line).size(); j++)
	    {
		gList[k].line[j] =  tmp2[gList[k].line[j]].o;
	    }

   
    //------------------------------fix lines ------------------------------------

    //--- update grains in lList  lList[].grainI/J using IOgrCenter
    vector<Pair> OIgrCenter=vecFlip(IOgrCenter);
    for (int k=0; k<lList.size(); k++)
    {
      if (lList[k].id > -1)
	{
      	  int gI = lList[k].grainI;
	  int gJ = lList[k].grainJ;
	  vector<Pair>::iterator q1,q2;
	  
	  // find grain as .o in IOgrCenter list 
	  Pair p1 = Pair(gList[gI].id,0);
	  Pair p2 = Pair(gList[gJ].id,0);
	  
	  q1 = find(OIgrCenter.begin(), OIgrCenter.end(), p1);
	  q2 = find(OIgrCenter.begin(), OIgrCenter.end(), p2);

	  if ( (q1) != OIgrCenter.end() ) //found outside grain q1
	      lList[k].grainI = (*q1).o;
	  
	  if ( (q2) != OIgrCenter.end() ) //found outside grain q2
	      lList[k].grainJ = (*q2).o;
	}
    }    
    //------------------------------fix vertices ------------------------------------

    //--- update grains in vList[].grain using IOgrCenter
    for (int k=0; k<vList.size(); k++)
    	for (int j=0; j<(vList[k].grain).size(); j++)
	{
	    int g=vList[k].grain[j]; //current grain name
	    vector<Pair>::iterator qq;
	    Pair pp = Pair(gList[g].id,0);
	    qq = find(OIgrCenter.begin(), OIgrCenter.end(), pp);
	    if ( qq != OIgrCenter.end() )
		vList[k].grain[j]=(*qq).o ;
	}

    //--- update lines in vList[].line using tmp2 
    vector<Pair>::iterator itL;
    for (int k=0; k<vList.size(); k++)
    	if (vList[k].isInside == 1)
	    for (int j=0; j<vList[k].line.size(); j++)
		vList[k].line[j] =  tmp2[vList[k].line[j]].o;
    
    //------------------- kill unneeded structure elements----------------------------
    //--- kill all grains outside domain
    for (int k=0; k<IOgrCenter.size(); k++)
    	gList[IOgrCenter[k].o].kill();

    //--- kill all lines with vertices outside domain
    for (int k=0; k<lList.size(); k++)
	if (vList[lList[k].vertexI].isInside ==0 || 
	    vList[lList[k].vertexJ].isInside == 0)
	    lList[k].kill();
    
    //--- kill all vertices outside domain
    for (int k=0; k<vList.size(); k++)
	if (vList[k].isInside == 0)
	    vList[k].kill();

    //--- put all the structure elements in ordered vector
    
    //--- lines ---
    Line tmpL = lList[0];
    tmpL.id = -1; 
    tmpL.vertexI = -1;
    tmpL.vertexJ = -1;
    tmpL.grainI = -1;
    tmpL.grainJ = -1;
     
    vector<Line> newLine(largeLineNumber, tmpL);
    
    for (int j=0; j<lList.size(); j++)
    	if( lList[j].id > -1) 
	    newLine[lList[j].id] = lList[j];
    
    // write back on array lList: 
     lList.clear();
     
     for (int j=0; j<newLine.size(); j++)
	 lList.push_back(newLine[j]);
     newLine.clear();

     //cout << "-------- print structure after ordering --- " << endl;
    
          //---vertices
     Vertex tmpV = vList[0];
     tmpV.id = -1;
     vector<Vertex> newVertex(vList.size(), tmpV);
     for (int j=0; j<vList.size(); j++)
    	 if (vList[j].id > -1)
	     newVertex[vList[j].id ] = vList[j];
     
     vList.clear();
     for (int j=0; j<newVertex.size(); j++)
	 vList.push_back(newVertex[j]);

     newVertex.clear();
//      cout << " print at the end of periodicStructure" << endl;
//      for (int k=0; k<vList.size(); k++)
//   	 if (vList[k].id > -1)
// 	     printVertex(vList[k]);
     
//      printGrains(gList); 
//      printLines(lList);
     
//      printLinesGNU(lList,vList); 
//      printGrainsGNU(gList,vList); 
     //exit(0);
     
}

//=========================================================================
void finalizeGrains(vector<Vertex> &vList, vector<Grain> &gList, 
		    vector<Line> &lList)
{

    // update Grains to include additional info
   
    vector<Grain>::iterator g;
    vector<Vertex>:: iterator ver;
    vector<int>::iterator v;
    vector<Line>::iterator l;
    Line tmpLine;
 
    vector<int> *gr = new vector<int> [gList.size()]; 
    vector<int> *vr = new vector<int> [vList.size()];
    vector<Line>::iterator ln; 
    vector<int>::iterator vNext;

    // fix up the line list for each grain. the vertices are ordered 
    // we use the same order for lines. 
    
    for(int i=0; i<gList.size(); i++)
    {
	for(v=(gList[i].vertex).begin()++; v != (gList[i].vertex).end(); v++)
	{
	    int v1 = (*v); int v2; vNext = v; vNext++;
	    if( vNext != (gList[i].vertex).end() ) v2 = (*vNext);
	    else v2 = (gList[i].vertex).front();
	    int k; 
	    for(k=0; k<lList.size(); k++)
	    {
		if( (lList[k].vertexI == v1)&&(lList[k].vertexJ == v2) ||
		    (lList[k].vertexI == v2)&&(lList[k].vertexJ == v1) )
		    break;
	    }
	    if( k == lList.size() )
	    {
		cout << " A BUG: Could not find Line for vertex v1 v2 ";
		cout << v1 << " " << v2 << endl;
		exit(0);
	    }
	    (gList[i].line).push_back(k);
	}

    }
}
//===========================================================================
void checkClockwise(Grain &gr, Vec2D vec1, Vec2D vec2)
{
    
    vector<int>::iterator v;
    double crossPd = vec1.crossProduct(vec2);
    //cout << " cross prd " << gr.id << " " << crossPd << endl;
 
    if (crossPd < 0) // change the orer of points to counterclockwise
    {
	vector<int> tmp; 
	int n = (gr.vertex).size();
	for( v=(gr.vertex).begin(); v!= (gr.vertex).end(); v++)
	    tmp.push_back(*v);
    } 
}
//=======================================================================
void finalizeVertices(vector<Vertex> &vList, vector<Grain> &gList,
		      vector<Line> &lList, vector<Pair> &IOqv )
{
    vector<int>::iterator g,gg;
    vector<int>::iterator lLast, lNext;
    vector<int>::iterator it, v;
    
    for(int i=0; i<lList.size(); i++) // looping over lines
    {
	int vI =  lList[i].vertexI; 
	int vJ =  lList[i].vertexJ; 
 	if( vI < 0 || vJ < 0) 
	{
	    cout << "bug in Lines (negative index for vertex) vI vJ ";
	    cout << vI << '\t' << vJ << endl;
	}
	(vList[vI].line).push_back( lList[i].id );
    	(vList[vJ].line).push_back( lList[i].id );  
    }
    
    for(int i=0; i<gList.size(); i++) // loop over grains
	for(v=(gList[i].vertex).begin(); v!= (gList[i].vertex).end(); v++) 
	    (vList[(*v)].grain).push_back(gList[i].id);
    
    // the lines and grain lists in each Vertex are complete but not ordered 
    // we order them now.
    
    for(int i=0; i<vList.size(); i++) // loop over vertices
    {
	//cout << " ------------------next vertex ----------i =  "  << vList[i].id << endl;
	     //printVertex(vList[i]);

	vector<int> newLList;
	vector<int> newGList;
	
	lNext = (vList[i].line).begin(); // pick first line 
	//cout << " first lNext:  " << *lNext << endl;
	int firstL = (*lNext);
	newLList.push_back(firstL); // add entry in vlist.line of current vertex
	
	int tmpL = -1; // tmpl : index of current line 
	int gPrev = -1;
	
	while( 1 )
	{
	    if( tmpL == -1) tmpL = firstL; // only first round of this loop 
	    int gI = -1; 
	    int gTest = -1;

	    //--- find grain that has this line
	    findGrain( vList, gList, i, tmpL, gI, gTest, gPrev, it);

	    if ( gTest == -1) // found grain in findGrain
	    {
		// reverse order of newLList and newGList 
		//cout << " reverse order of newLList and newGList " << endl;
		std::reverse(newLList.begin(), newLList.end());
		std::reverse(newGList.begin(), newGList.end());

		// switch values for previous grain and current lineL
		gPrev = newGList.back();
		tmpL = newLList.back();
		
                //--- start again grain loop
		findGrain( vList, gList, i, tmpL, gI, gTest, gPrev, it);
	    } // end of gTest loop
	    
	    //cout << " gI pushback in newGList" << gI << endl;
	    newGList.push_back(gI);    
	    
            //--- stopping criterion: 
	    if (newGList.size()==(vList[i].grain).size() && newLList.size()==(vList[i].line).size())
	    {
		//cout << " all grains and lines have been found and ordered, break 1 " << endl;
		break;
	    }       
	    
	    //---  find next line	  	  
	    lNext = it; //  
	    vector<int> :: iterator lNextT = findNextLine(vList,gList,lList,i,tmpL,gI,lNext) ;
	    //cout << "lNextT " << (*lNextT) << endl;
	    lNext = lNextT;
	    newLList.push_back((*lNext));
	    
	    //--- stopping criterion: lines and grains --- 
	    if (newGList.size() == (vList[i].grain).size() && newLList.size() == (vList[i].line).size() )
	    {
		//cout << " all grains and lines have been found and ordered, break 2 " << endl;
		break;
	    }
	    gPrev = gI;
	    gI = -1;
	    tmpL = (*lNext); 
	    	    	    	    
	} // end while

	
      	//---  checking clockwise orientation
	//cout << endl;
	//cout << " check clockwise orientation \n" ;
	
	int v1 = vList[i].otherVertex(lList[ newLList[0] ]);
	int v2 = vList[i].otherVertex(lList[ newLList[1] ]);

	//cout << " v1 = " << v1 << endl;
	//cout << " v2 =" << v2 << endl;
	Vec2D vec1;
	Vec2D vec2;
	vec1 = vList[v1].pos - vList[i].pos;
	vec2 = vList[v2].pos - vList[i].pos;
	
	int nl = newLList.size();
	int ng = newGList.size();
	double crossPd = vec1.crossProduct(vec2); 
	//cout << " crossPd " << crossPd << endl;
	(vList[i].line).clear(); 
	(vList[i].grain).clear();

	vector< int > :: reverse_iterator k2; 

	if (crossPd < 0) // change the order of points to counterclockwise
	{
	    // cout << " change order of points to counterclockwise \n";
	    for ( k2=newLList.rbegin(); k2 != newLList.rend(); k2++)
		(vList[i].line).push_back(*k2);
	    for ( k2=newGList.rbegin(); k2 != newGList.rend(); k2++)
		(vList[i].grain).push_back(*k2);
	    
            // if outside grain, shift grainlist to make sure, 
	    // first grain is between first and second line

	    if (nl == ng ) 
	    {
		int g1 = (vList[i].grain).front(); 
		(vList[i].grain).erase(vList[i].grain.begin());  // get rid of first element
		(vList[i].grain).push_back(g1);
	    }
	} 
	else 
	{	    
	    for(int k=0; k<newLList.size(); k++) 
		(vList[i].line).push_back(newLList[k]);
	    for(int k=0; k<newGList.size(); k++) 
		(vList[i].grain).push_back(newGList[k]);
	}
    }
 }
//=======================================================================

void findGrain(vector<Vertex> &vList, vector<Grain> &gList, int i, int tmpL, 
	       int &gI, int &gTest, int gPrev, vector<int> :: iterator &it)
{
    vector<int>::iterator g;
    for( g = (vList[i].grain).begin(); g !=  (vList[i].grain).end(); g++)
    {
	
	//printGrain(gList[(*g)]);
	it = find((gList[(*g)].line).begin(), (gList[(*g)].line).end(), tmpL); 
	if( it != (gList[(*g)].line).end()  && gI == -1 && (*g) != gPrev ) 
	{
	    gI = (*g);
	    gTest = 1;
	    
	    break;
	}
	
    }   
}
//=======================================================================
vector<int>::iterator findNextLine(vector<Vertex> &vList, vector<Grain> &gList, 
  vector<Line> &lList, int i, int tmpL, int gI, vector<int> :: iterator &lNext)
{
    // check which vertex is current vList[i] (vI or vJ)
    //     check if other vertex is v++ or v--
    //     if v++ ==> lNext--
    //     iv v-- ==> lNext++
    // take care of .end() end .begin(), chose wrt direction 
    //-----------------------------------------------------------------	 
    
    vector<int>::iterator l, lLast;
    vector<int>::iterator v, vLast;
    
    int vv = vList[i].id ; // current vertex vv 
    int vI = lList[tmpL].vertexI; // position in list (*g).vertex;
    int vJ = lList[tmpL].vertexJ; // position in list (*g).vertex;


    v = find( (gList[gI].vertex).begin(), (gList[gI].vertex).end(),vv);
    
    vLast = (gList[gI].vertex).end() ; // 
    vLast--; // this is last vertex in vertex list! 
    lLast = (gList[gI].line).end() ; // 
    lLast--; // this is last line in grain list! 
        
    //--- determine v++ 
    if ( v != vLast)
	v++;
    else   
	v=(gList[gI].vertex).begin();
        
    if (vI == vv) // current vertex is vI
    { 	     
	//--- decide about lNext depending on v++ 
	if ((*v) == vJ) 
	{
	    if( lNext != (gList[gI].line).begin() ) lNext--;
	    else 
	    {
		lNext = (gList[gI].line).end(); 
		lNext--;
	    }
	}
	else
	    if( lNext != lLast )  lNext++;
	    else lNext = (gList[gI].line).begin(); 
    }
    else if (vJ == vv)
    {   
	if ((*v) == vI) 
	{
	    if( lNext != (gList[gI].line).begin() ) lNext--;
	    else 
	    {
		lNext = (gList[gI].line).end()--; 
		lNext--;
	    }
	}
	else
	{ 
	    if( lNext != lLast)  lNext++;
	    else lNext = (gList[gI].line).begin(); 
	}
    }
    else
	cout << " we have a problem, vlist[i] is neither vI nor vJ " << endl;

    //cout << " new lNext:  " << *lNext << endl;
    //cout << " ----------------------------- " << endl;
    return lNext;
}
//=================== Consistency Check ================================
void checkConsist(vector<Vertex> vList, vector<Grain>  gList, vector<Line> lList )
{


    //--- check vertices and lines in gList
    for (int i=0; i<gList.size(); i++)
    {
	if (gList[i].id > -1)
	{
	    int ctrV = 0;
	    for (int j=0; j<gList[i].vertex.size(); j++)
	    {
		int v =gList[i].vertex[j]; // pick a vertex 
		for (int k=0; k<vList[v].grain.size(); k++) 
		    if (vList[v].grain[k] == gList[i].id ) 
			ctrV = 1; // grain shows up in vertex.grain
		if ( ctrV == 0 )
		{
		    cout << " consistency check failed at grain " << gList[i].id 
			 << " and vertex " << vList[v].id << endl;
		    exit(0);
		}
	    }
	    //cout << " Vertices in gList are consistent with grains in vList,  o.k. " << endl;
	    int ctrL = 0;
	    for (int j=0; j<gList[i].line.size(); j++)
	    {
		int l = gList[i].line[j]; // pick a line of grain
		if (lList[l].grainI == gList[i].id || lList[l].grainJ == gList[i].id)
		    ctrL = 1;
		if ( ctrL == 0 )
		{
		    cout << " consistency check failed at grain " << gList[i].id 
			 << " and line "<< lList[l].id << endl;
		    exit(0);
		}
		
	    }
	    //cout << " Lines in gList are consistent with grains in lList, o.k. " << endl;
	}
    }
    //--- check lines and grains in vList
    for (int i=0 ; i< vList.size(); i++)
    {
	if (vList[i].id > -1 )
	{
	    int ctrL = 0; // check lines 
	    for (int j=0; j<vList[i].line.size(); j++)
	    {
		int l = vList[i].line[j] ;
		if (lList[l].vertexI == vList[i].id || lList[l].vertexJ == vList[i].id)
		    ctrL = 1;
		
		if ( ctrL == 0 )
		{
		    cout << " consistency check failed at vertex " << vList[i].id 
			 << " and line "<< lList[l].id << endl;
		    exit(0);
		}
	    }
	    //cout << " lines in vList are consistent with vertices in lList, o.k. " << endl;
	    int ctrG = 0; // check grains
	    for (int j=0; j<vList[i].grain.size(); j++)
	    {
		int g = vList[i].grain[j]; // pick a grain of this vertex 
		for (int k=0; k<gList[g].vertex.size(); k++) 
		    if (gList[g].vertex[k] == vList[i].id )
			ctrG = 1;
		
		if ( ctrG == 0 )
		{
		    cout << " consistency check failed at vertex " << vList[i].id 
			 << " and grain "<< gList[g].id << endl;
		    exit(0);
		}
	    }
	    //cout << " grains in vList are consistent with vertices in gList, o.k. " << endl;
	}
    }
    //---check grains and vertices in lList
    for (int i=0; i<lList.size(); i++)
    {
	//printLine(lList[i]);
	if (lList[i].id > -1) 
	{
	    int ctrG = 0; // check grains
	    int gI= lList[i].grainI; // pick a grain
	    for (int j=0; j<gList[gI].line.size(); j++)
		if ( gList[gI].line[j] == lList[i].id)
		    ctrG = 1;
	    
	    if ( ctrG == 0 )
	    {
		cout << " consistency check failed at line " << lList[i].id 
		     << " and grain "<< gList[gI].id << endl;
		printGrain(gList[gI]);
		printLine(lList[i]);
		exit(0);
	    }	
	    //cout << " grain gI in lList is consistent with lines in gList[gI], o.k. " << endl;
	    ctrG = 0;
	    int gJ= lList[i].grainJ;
	    for (int j=0; j<gList[gJ].line.size(); j++)
		if ( gList[gJ].line[j] == lList[i].id)
		    ctrG = 1;
	    
	    if ( ctrG == 0 ) // here is the current problem !!!!!!!!!!!!!!!!!!!!!
	    {
		cout << " consistency check failed at line " << lList[i].id 
		     << " and grain "<< gList[gJ].id << endl;
		printGrain(gList[gJ]);
		printLine(lList[i]);
		exit(0);
	    }	
	    //cout << " grain gJ in lList is consistent with lines in gList[gJ] " << endl;
	    //------------------------
	    int ctrV = 0; // check vertices
	    int vI = lList[i].vertexI; // pick a vertex
	    for (int j=0; j<vList[vI].line.size(); j++)
	    {
		if ( vList[vI].line[j] == lList[i].id )
		    ctrV = 1;
	    }
	    if ( ctrV == 0 )
	    {
		cout << " consistency check failed at line " << lList[i].id 
		     << " and vertex "<< vList[vI].id << endl;
		exit(0);
	    }
	    //cout << " vertex vI in lList is consistent with lines in vList[vI]" << endl; 
	    ctrV = 0;
	    int vJ = lList[i].vertexJ;
	    for (int j=0; j<vList[vJ].line.size(); j++)
	    {
		if ( vList[vJ].line[j] == lList[i].id )
		    ctrV = 1;
	    }
	    if ( ctrV == 0 )
	    {
		cout << " consistency check failed at line " << lList[i].id 
		     << " and vertex "<< vList[vJ].id << endl;
		exit(0);
	    }
	    //cout << " vertex vJ in lList is consistent with lines in vList[vJ]" << endl; 
	}
    }
    
    cout << " Consistency check passed! " << endl;
}



//====================   PRINTING FUNCTIONS   ===========================
void printLine(Line l)
{
    cout << endl << "Line: " << endl;
    cout << l.id << "   V: " << l.vertexI << " ";
    cout << l.vertexJ;
    cout << "    G:" << l.grainI << " " << l.grainJ << endl;
}
//========================================================================
void printLines(vector<Line> &lList)
{
    cout << " printLines " << endl;
    vector<Line>::iterator l;

    cout << endl << "Lines " << endl << endl;
    for(l=lList.begin(); l != lList.end(); l++) 
    {
	cout << "line: " << (*l).id << "   V: " << (*l).vertexI << " ";
	cout << (*l).vertexJ;
	cout << "    G:" << (*l).grainI << " " << (*l).grainJ << endl;
    }
}
//======================================================================
void printGrains(vector<Grain> &gList)
{
    vector<Grain>::iterator g;
    vector<int>::iterator v;
    vector<int>::iterator l;

    cout << endl << "Grains " << endl << endl ;
    for(g=gList.begin(); g != gList.end(); g++) 
    {
	//if ((*g).id > -1) 
	//{
	    cout << "Grain no :" << (*g).id << " V: ";
	    for(v=((*g).vertex).begin(); v != ((*g).vertex).end(); v++)
	    {
		cout << (*v) << " ";    
	    }
	    cout << "   L: ";
	    for(l=((*g).line).begin(); l != ((*g).line).end(); l++)
		cout << (*l) << " ";
	    cout << endl;
	    //} 
    }
}
//======================================================================
void printGrain(Grain g)
{
    vector<int>::iterator v;
    vector<int>::iterator l;
    cout << "Grain " << endl;
    cout << "Grain no :" << g.id << " V: ";
    for(v=g.vertex.begin(); v != g.vertex.end(); v++)
	cout << (*v) << " ";
    cout << "   L: ";
    for(l=g.line.begin(); l != g.line.end(); l++)
	cout << (*l) << " ";
    cout << endl;
}
//======================================================================
void printVertices(vector<Vertex> &vList)
{
    vector<Vertex>::iterator v;
    vector<int>::iterator l;
    vector<int>::iterator g;

    cout << endl << "Vertices " << endl << endl;
    int m = 0;
    for(v=vList.begin(); v != vList.end(); v++)
    { 
	cout << (*v).id << " " << ((*v).pos).x <<  " " << ((*v).pos).y << endl;	
    }
  
}
//=========================================================================
void printVertex(Vertex v)
{
    vector<int>::iterator l;
    vector<int>::iterator g;

    cout << endl << "Vertex:  " << endl;

    cout << "Id " << v.id << " " << (v.pos).x;
    cout << " " << (v.pos).y << endl;
    cout << " L : ";
    for( l=(v.line).begin(); l != (v.line).end(); l++) 
	cout << (*l) << " ";
    cout << endl;
    cout << " G : ";
    for( g=(v.grain).begin(); g != (v.grain).end(); g++) 
	cout << (*g) << " ";
    cout << endl;
}

//=============================== Get Statistics =========================================

void writeSides_Old(vector<Grain> &gList, string fName)
{
    vector<int> sides;
    int N = 20;

    ofstream oFile;
    oFile.open( fName.c_str(), ios::out);
    if( oFile == NULL ) 
    {
	cout << " could not open output file " << fName << endl;
	exit(0);
    }

    for (int n = 0; n < N; n++)
    {
	sides.push_back(0);
    }


    for(int g=0; g < gList.size(); g++)
    {
	if( gList[g].id != -1 )
	{		  
	    int nsides = (gList[g].line).size();
	    sides[nsides] += 1;
	}
    }
    for (int n = 0; n < sides.size(); n++)
    {
	oFile << n << " " << sides[n] << endl;
    }
  

}
//=====================================================================================
void writeSides(vector<Grain> &gList, string fName,string fName2)
{
    vector<int> sides;
    int N = 20;

    ofstream oFile,oFile2;
    oFile2.open( fName2.c_str(), ios::out);
    if( oFile2 == NULL ) 
    {
	cout << " could not open output file " << fName2 << endl;
	exit(0);
    }

    for (int n = 0; n < N; n++)
    	sides.push_back(0);
    
    for(int g=0; g < gList.size(); g++)
    	if( gList[g].id > -1 )
	{	
	    int nsides = (gList[g].line).size();
	    oFile2 << nsides << endl; 
	    sides[nsides] += 1;
	}
    oFile2.close();

//     oFile.open( fName.c_str(), ios::out);
//     if( oFile == NULL ) 
//     {
// 	cout << " could not open output file " << fName << endl;
// 	exit(0);
//     }
//     for (int n = 0; n < sides.size(); n++)
//     	oFile << n << " " << sides[n] << endl; // gives histogram of sides 
    
}
//=====================================================================================
void writeLineLength(vector<Line> lList, vector<Vertex> vList,  string fName)
{
    ofstream oFile;
    vector<double> Length;
    vector<int> lId;
    Vec2D dir;
    oFile.open( fName.c_str(), ios::out);
    if( oFile == NULL ) 
    {
	cout << " could not open output file " << fName << endl;
	exit(0);
    }

    list<Vec2D>::iterator n, nNext, nLast;
    for (int l=0; l<lList.size(); l++)
    {
	if (lList[l].id > -1)
	{
	  nLast = lList[l].node.end();
	  nLast--;
	  
	  dir = Vec2D(0.0,0.0);
	  double lineLength = 0.;
	  for (n=lList[l].node.begin(); n!= nLast; n++)
	    {
	      nNext = n;
	      nNext++;
	      //dir = vList[lList[l].vertexJ].pos - vList[lList[l].vertexI].pos;
	      dir = ((*n) - (*nNext));
	      lineLength = lineLength + dir.length();
	    }
	  Length.push_back(lineLength);
	  lId.push_back(lList[l].id);
	}
    }
    for (int i=0; i<Length.size(); i++)
    {
	oFile << lId[i] << " " << Length[i] << endl;
    }
    oFile.close();
}

//=====================================================================================
void writeRelativeArea_Eva(vector<Grain> &gList, vector<Line> &lList, vector<Vertex> &vList, 
		       string fName0, string fName1, string fName2, double totalArea, 
			   vector<double> SingArea)
{
    ofstream oFile0, oFile1, oFile2,oFile3;
    
    int x, xNext ;
    int nGrains=0;
    vector<double> Area;
    vector<Vec2D> TdsArea;
    totalArea=0.0;
    vector<int> gId;
    
    for(int g=0; g < gList.size(); g++)
    {
	if( gList[g].id > -1 )
    	{
	    //printGrain(gList[g]);
	    Vec2D tds, xMid;
	    double singArea=0.0;
	    Vec2D sumTds = Vec2D(0.,0.); 
	    
	    x = gList[g].vertex[0];
	    xNext = gList[g].vertex[1];
	    double v1,v2;
	    v1 = x;
	    v2 = xNext;
	    double tmpx = vList[xNext].pos.x - vList[x].pos.x;
	    double tmpy = vList[xNext].pos.y - vList[x].pos.y;

	    sumTds = vList[xNext].pos - vList[x].pos;
	    
	    for (int v=1; v<gList[g].vertex.size(); v++)
	    {
		x = gList[g].vertex[v];
		if (v != gList[g].vertex.size()-1)
		    xNext = gList[g].vertex[v+1];
		else
		    xNext = gList[g].vertex[0];
		tds = vList[xNext].pos - vList[x].pos;
		singArea += sumTds.crossProduct(tds);
		
		sumTds = sumTds + tds;
	    } 
	    
	    singArea = 0.5*abs(singArea);
	    
	    totalArea += singArea;
	    Area.push_back(singArea);
	    
	    SingArea.push_back(singArea);
	    gId.push_back(gList[g].id);

	    nGrains++;
	} 
	else
	{
	    SingArea.push_back(-2);
	    gId.push_back(-1);
	}
    }// end of gList-loop

    // Area contains area of each single grain 

    double avArea=totalArea / nGrains;
    cout <<"total Area " << totalArea << "\t Average area " << avArea << endl;

    if (abs(totalArea - 4.)> 0.1 )
    {
	cout << "totalArea != 4," << endl;
	//exit(0);
    }

    //--- write single Area 
    oFile2.open( fName2.c_str(), ios::out);
    if( oFile2 == NULL ) 
    {
	cout << " could not open output file " << fName2 << endl;
	exit(0);
    }
    int gLineSize;
    for (int ii=0; ii < SingArea.size(); ii++)
    {   
	if (gId[ii] == -1)
	    gLineSize = 0;
	else
	    gLineSize =gList[gId[ii]].line.size();
	
	oFile2 << gId[ii] << " " << SingArea[ii] << " " <<gLineSize << endl; 
    }
    oFile2.close();

    





//     //----------------------------------
//     for ( int i=0; i<Area.size() ; i++ )
//       Area[i] = Area[i]/avArea;     

     // now Area is relative Area of each single grain 
    
    oFile1.open( fName1.c_str(), ios::out);
    if( oFile1 == NULL ) 
    {
	cout << " could not open output file " << fName1 << endl;
	exit(0);
    }
    for (int i=0; i<Area.size(); i++)
	oFile1 << Area[i] << "\n " ;
    oFile1.close(); 


     //---- for histogram 
    double minArea = 1000.;
    double maxArea = 1.;
    for (int j=0; j<Area.size(); j++)
    {
	if ( minArea > Area[j] )
	    minArea = Area[j];
	if ( maxArea < Area[j] ) 
	    maxArea = Area[j];
    }
    
    double histArea[30];
    for(int i=0; i< 30; i++) histArea[i] = 0.0; 
    
    for (int j = 0; j< Area.size() ; j++)
    {
	int k = (int)( (Area[j] - minArea)/(maxArea - minArea) * 30);
	histArea[k]++;
    }
    
    //--- write relative Area for Histogram 
    oFile0.open( fName0.c_str(), ios::out);
    if( oFile0 == NULL ) 
    {
	cout << " could not open output file " << fName0 << endl;
	exit(0);
    }
    
    for (int i=0; i<30; i++)
	oFile0 << histArea[i] << "\t " ;
      
}



//=====================================================================================

void writeRelativeArea_Nodes(vector<Grain> &gList, vector<Line> &lList, vector<Vertex> &vList, 
		       string fName0, string fName1, string fName2, double totalArea, 
			   vector<double> SingArea)
{
    ofstream  oFile2;
    
    int x, xNext ;
    int nGrains= 0;
    int nLines = 0;
    
    string GrainArea;
    
    vector<int> gId;
    vector<double> Area;
    vector<Vec2D> TdsArea;
 
    totalArea=0.0;
    //printGrains(gList);

    for(int g=0; g < gList.size(); g++)
    {
	if( gList[g].id > -1 )
	{ 
	    vector<Vec2D> z; // temporary nodes
	    list<Vec2D>:: iterator node, n, nFirst, nLast, nNext;
	    nLines = (gList[g].line).size();
	    
	    int firstL = lList[gList[g].line[0]].id;
	    int vJ_old = vList[gList[g].vertex[0]].id;
	    int vI = lList[firstL].vertexI;
	    int vJ = lList[firstL].vertexJ;

	    Vec2D zLast;
	    zLast = vList[vJ_old].pos; 
	    z.push_back(zLast);
	    

	    for (int ll = 0; ll < (gList[g].line).size();  ll++)
	    {
		int reverse = 0;
		int l = gList[g].line[ll];
		
		vI = lList[l].vertexI;
		vJ = lList[l].vertexJ;
		
		if (vList[vI].id == vList[vJ_old].id)
		    vJ_old = vJ;
		
		else // reverse current nodelist 
		{
		    reverse = 1;
		    (lList[l].node).reverse();
		    vJ_old = vI;
		}
		
		nFirst = lList[l].node.begin();
		nLast = lList[l].node.end();
		nLast--;
		
		for (n = lList[l].node.begin(); n != nLast; n++)
		{
		    nNext = n;
		    nNext++;
		    
		    double tmpx = (*nNext).x - (zLast).x; // distance switched nodes
		    double tmpy = (*nNext).y - (zLast).y;
		    
		    if (abs(tmpx) < 1. && abs(tmpy) < 1.)
		    {
			zLast = (*nNext);
			z.push_back(*nNext);
			
		    }
		    else //--- switch nodes
		    {
			if (abs(tmpx) > 1. && abs(tmpy) <= 1. )
			{  
			    double new_nNext_x;
			    
			    if ((*nNext).x < 0)
				new_nNext_x = (*nNext).x + 2.;
			    else
				new_nNext_x = (*nNext).x - 2.;
			    zLast.x = new_nNext_x ;
			    zLast.y = (*nNext).y ;
			    
			}
			if (abs(tmpx) <= 1 && abs(tmpy) > 1 )
			{   
			    
			    double new_nNext_y;
			    if ((*nNext).y < 0)
				new_nNext_y = (*nNext).y + 2.;
			    else
				new_nNext_y = (*nNext).y - 2.;
			    zLast.x = (*nNext).x ;
			    zLast.y = new_nNext_y ;
			}
			
			if (abs(tmpx) > 1 && abs(tmpy) > 1 )
			{   
			    double new_nNext_x, new_nNext_y;
			    if ((*nNext).x < 0)
				new_nNext_x = (*nNext).x + 2.;
			    else
				new_nNext_x = (*nNext).x - 2.;
			    
			    if ((*nNext).y < 0)
				new_nNext_y = (*nNext).y + 2.;
			    else
				new_nNext_y = (*nNext).y - 2.;
			    zLast.x = new_nNext_x;
			    zLast.y = new_nNext_y;
			}
			
			z.push_back(zLast);
			
		    } // end else
		   
		} // end of node-loop (of line l)
		
		
		if (reverse == 1)
		    lList[l].node.reverse();
	    
	    } // end of lines (of grain g)  
	    //----------------------------------------------
	    // calculate the area void calcArea(z);
	    
	    Vec2D tds;
	    double singArea=0.0;
	    double sumTds = 0.;
	    
	    Vec2D z1, z2, normal;
	    vector<Vec2D>::iterator it,zit, zEnd;
	    
	    zEnd = z.end();
	    zEnd--;

	    for (it = z.begin(); it !=zEnd; it++)
	    {
		zit = it;
		zit++;
		
		z1 = (*it);
		z2 = (*zit);

		tds.x = z1.x + (z2.x - z1.x)*0.5;
		tds.y = z1.y + (z2.y - z1.y)*0.5;
		normal = (z2 - z1).perp(); 
		sumTds = sumTds + tds.dotProduct(normal) ;
	    }
	    
	    singArea = sumTds;
	    singArea = 0.5*abs(singArea); ;
	    
	    totalArea += singArea;
	    
	    Area.push_back(singArea);
	    
	    SingArea.push_back(singArea);
	    gId.push_back(gList[g].id);
	    
	    nGrains++;
	   
	    z.empty();
	    
	} // end of g.Id > -1
	
    }// end of gList-loop
    
    // Area contains area of each single grain 
    
    double avArea=totalArea / nGrains;
    cout <<"total Area " << totalArea << "\t Average area " << avArea << endl;
    
    if (abs(totalArea - 4.)> 0.1 )
    {
	cout << "totalArea != 4," << endl;
	//exit(0);
    }
    
    //--- write single Area 
    oFile2.open( fName2.c_str(), ios::out);
    if( oFile2 == NULL ) 
    {
	cout << " could not open output file " << fName2 << endl;
	exit(0);
    }
    int gLineSize;
    for (int ii=0; ii < SingArea.size(); ii++)
    {   
	if (gId[ii] == -1)
	    gLineSize = 0;
	else
	    gLineSize =gList[gId[ii]].line.size();
	
	oFile2 << gId[ii] << " " << SingArea[ii] << " " <<gLineSize << endl; 
    }
    oFile2.close();
 
}


// //=========================================================================
void neighbors(vector<Grain> &gList, vector<Line> &lList, string &fName )
{
    ofstream oFile;
    oFile.open( fName.c_str(), ios::out);
    if( oFile == NULL ) 
    {
	cout << " could not open output file " << fName << endl;
	exit(0);
    }
    
    int gg;
    for (int g=0; g<gList.size(); g++)
    {
	vector<int> nrSides;
	if (gList[g].id > -1)
	{
	    oFile << gList[g].id << " " << gList[g].line.size() << " " ;
	    
	    for (int l=0; l<gList[g].line.size(); l++)
	    {
		
		int ll = gList[g].line[l];
		if (lList[ll].grainI == gList[g].id )
		    gg = lList[ll].grainJ;
		else 
		    gg = lList[ll].grainI;
		int sides = gList[gg].line.size();
		oFile << sides << " " ;
		nrSides.push_back(sides);
	    }
	    oFile << endl;

	}
	//else
	//  oFile << -1 << " " << -2 <<  endl; 
    }

    oFile.close();   

}

//=====================  DISCRETIZATION FUNCTIONS ==========================
void initializeNodes(vector<Vertex> &vList, vector<Line> &lList)
{
    double maxLen = -1.0;
    double minLen = 10000000.0;

    double aveLen=0.;
    int counter = 0;
    Vec2D zero = Vec2D(0.0,0.0);


    // get average length line 
    for(int l=0; l<lList.size(); l++)
    {
	if ( lList[l].id > -1)
	{
	    counter++;
	    int vI = lList[l].vertexI;
	    int vJ = lList[l].vertexJ;

	    Vec2D dir = vList[vJ].pos - vList[vI].pos;
	    
	    double vecLen = dir.length();
	    aveLen += vecLen;
	    lList[l].length = vecLen;
	    if( maxLen < vecLen ) maxLen = vecLen;
	    if( minLen > vecLen ) minLen = vecLen;
	}
	
 
	(lList[l].node).clear();
	(lList[l].nodeVel).clear();
	(lList[l].normal).clear();
    }
    
    aveLen = aveLen/counter;
    cout <<"minLen " << minLen << " maxLen " << maxLen <<  " aveLen " << aveLen << endl;
    hMin = aveLen/3;		//

    epsLine = hMin * 2.;// *1.5; // THINK ABOUT THIS !!! 
    //if (DEBUG) 
    cout <<   "hmin = \t" << hMin << endl;  
    cout << " epsLine = \t" << epsLine << endl;

    //-------

    for(int l=0; l<lList.size(); l++)
    {
	if (lList[l].id > -1)
	{
	    Vec2D X = vList[lList[l].vertexI].pos; 
	    Vec2D Y = vList[lList[l].vertexJ].pos; 
	    
	    Vec2D segment = vList[lList[l].vertexJ].pos - vList[lList[l].vertexI].pos;
 
	    // Calculate mesh factor

	    //--- find proper nPts (2^n +1 points)
	    double  dn = ( log(segment.length()/(2.*hMin))  / log(2.)) + 0.5;
	    
	    int n = (int) (dn);
	    n = max(1,n);
	 
	    
	    dn = double(n);
	    dn  = pow(2., dn );
	    
	    int nPts = (int) (dn) + 1;

	    //nPts = 1;  // this makes it a 3 node code!
 	    Vec2D DX = segment / ((double) (nPts-1));
	    
	    for ( int i=0; i <  nPts-1; i++) 
	    {
		(lList[l].node).push_back(X);
		(lList[l].nodeVel).push_back(zero);
		X = X +  DX; 
	    }
	    
	    (lList[l].node).push_back(Y);
	    (lList[l].nodeVel).push_back(zero);
	    
	}
	if (lList[l].node.size() != lList[l].nodeVel.size())
	{
	    cout << " nodes != nVel " << endl;
	    exit(0);
	}
	
    }
}

//==============================================================================
void removeNodes(int l, vector<Line> &lList)
{
    vector<Vec2D> temp;
    Vec2D newNode ;
    list<Vec2D>::iterator n, nNext, nLast;
    list<Vec2D>::iterator normal, nVel, nVelNext;
    
    if ((lList[l].node.size()/2) * 2 == lList[l].node.size() )
    {	
	cout << " Attention, even number of nodes " << endl;
	cout << "line " << lList[l].id << " nodes " << lList[l].node.size() << endl;
	printLine(lList[l]);
	exit(0);
    }
 
    
    n    = lList[l].node.begin();
    nVel= lList[l].nodeVel.begin();
    
    //----
    //normal = lList[l].normal.begin();
    //----
    
    nLast = lList[l].node.end();
    nLast--;
    int nNodes=lList[l].node.size();

    for (n = lList[l].node.begin(); n!= nLast; n++)
    {
	nNext = n;
	nNext++;
	nVelNext = nVel ;
	nVelNext++; 


	if (nNext != nLast )
	{
	    //cout << " erase a node " << endl;
	    nNext  = (lList[l].node).erase(nNext); 
	    nVelNext = (lList[l].nodeVel).erase(nVelNext);
	    nVel++;
	    if (nNodes == 4)
		break;
	     
	}
	else if (nNext == nLast && nNodes != 4)
	{
	    cout << " erase the first node of this segment, even # nodes, then break" << endl;
	    (lList[l].node).erase(n);	
	    (lList[l].nodeVel).erase(nVel);
	    break;
	}
	
	//----
	//normal = (lList[l].normal).erase(normal);
	//nVel   = (lList[l].nVel).erase(nVel);
	//normal++;
	//nVel++;
	//----
    }
    //---- (lList[l].normal).erase(normal);
    
}
//==============================================================================
void addNodes(int l, vector<Line> &lList)
{
    Vec2D newNode, newVel ;
    list<Vec2D>::iterator n, next, nPrev, nLast, m, nn;
    list<Vec2D>::iterator nVel;
 
    n=lList[l].node.begin();
    nVel = lList[l].nodeVel.begin();
    newVel = (*nVel);
    
    nPrev = n;
        
    nLast = lList[l].node.end();
    nLast--;
    
    while( n!= nLast)
    {
	next = n;
	next++;
	newNode = (*n) + ((*next) - (*n))*0.5;
	n++;
	m = (lList[l].node).insert(next,newNode);

	nVel++;
	newVel = (*nVel);
	(lList[l].nodeVel).insert(nVel,newVel);

	
	
	// what about normal
    }
    
//}
}
//==============================================================================
void redistributeNodes(vector<Vertex> &vList, vector<Line> &lList)
{
    cout << "=======redistribute  hMin " << hMin <<  endl;
    for (int i=0; i<lList.size(); i++)
    	if (lList[i].id > -1)
	{
	    Vec2D segment = Vec2D(0.0,0.0);
	    double segLen = 0.;
	    double hLoc = 0.;
	    list<Vec2D>::iterator nLast, next, n;
	    nLast = lList[i].node.end();
	    nLast--;

	    for (n=lList[i].node.begin(); n!= nLast; n++)
	    {
		next = n;
		next++;
		segment =  ((*next) - (*n));
		segLen = segLen + segment.length();

	    }
	    hLoc = segLen/(lList[i].node.size()-1);
	    
	    if (hLoc > 1.5 * hMin )
	    {
		addNodes(i, lList);
	    }
	    else if (hLoc < 0.75 * hMin && lList[i].node.size() > 3) 
	    {
		removeNodes(i, lList);
	    }
	    
	}
    
}


//==============================================================================
void printNodes(vector<Line> &lList,string &fName)
{
    ofstream oFile;
    oFile.open( fName.c_str(), ios::out);
    if( oFile == NULL ) 
    {
	cout << " could not open output file " << fName << endl;
	exit(0);
    }

    list<Vec2D>::iterator  no;
    for (int i=0; i<lList.size(); i++)
    {
	if (lList[i].id > -1 )
	{	
	    //cout << " line " << lList[i].id << endl;
	    //printLine(lList[i]);
	    for (no =lList[i].node.begin(); no!= lList[i].node.end(); no++)
		oFile << (*no).x << " " << (*no).y << endl;
		    
	    oFile << endl;
	}
    }
    oFile.close();
    
}

//=========================================================================
void initializeNormals(vector<Line> &lList, vector<Vertex> &vList )
{
    list<Vec2D>::iterator n, nNext, normal;
    cout << " in initializeNormals " << endl;

    for(int l = 0; l < lList.size(); l++) 
    {
	if (lList[l].id > -1)
	{
	    int vI, vJ;
	    
	    vI = lList[l].vertexI;
	    vJ = lList[l].vertexJ;
	    Vec2D dx = vList[vI].pos - vList[vJ].pos;
	    dx = dx / dx.length();
	    (lList[l].normal).push_back(dx.perp());
	    //(lList[l].nodeVel).push_back(Vec2D(0.,0.));
	    
	    
	    normal = (lList[l].normal).begin();
	    Vec2D tempNorm = *normal;
	    
	    // check normal
	    double dotPd = -2;
	    if (l == 8) 
	    {
		Vec2D vec1= vList[lList[l].vertexI].pos- vList[lList[l].vertexJ].pos;
		double dotPd= vec1.x * tempNorm.x + vec1.y * tempNorm.y;
	    }
	}
    }
}


//========================================================================
    void defNewVertIndex(int &newV, vector<Vertex> &vList, vector<int> &EmptyVertex  )
{
    if( EmptyVertex.empty() ) 
	newV = vList.size(); 
    else
    {
	newV = EmptyVertex.back();
	EmptyVertex.pop_back();
    }
    
}
//========================================================================
void defNewLineIndex(int &newL, vector<Line> &lList, vector<int> &EmtpyLine)
{
    if( EmptyLine.empty() ) 
    	newL = lList.size();
    else 
    {
	newL = EmptyLine.back(); 
	EmptyLine.pop_back(); 
    }
}

//================ CRITICAL EVENTS FUNCTIONS ============================


void findSplitDir( int v, vector<Line> lList, vector<Vertex> vList, 
		   Vec2D &splitVec, int &splitInd)
{
   
    vector<int>::iterator l;
    vector<Vec2D> T; 
    
    for (int iLine=0; iLine < (vList[v].line).size(); iLine++) // number of lines
    {   // iLine is current index of the line e.g. (0,1,.. 3) 

        //--- find iterator l of VV.line at position iLine  
	l = vList[v].line.begin(); 
	for(int jj=0; jj < iLine; jj++)
	    l++;
		
	if ( (*l) > -1 )
	{
	    //--- calculate T for this line (*l) of Vertex VV
	    double aI =  (lList[(*l)]).alphaI;
	    double aJ =  (lList[(*l)]).alphaJ;

	    // define Vec1 and Vec2 depending on position
	    
	    int v2 = vList[v].otherVertex(lList[(*l)]);
	    
	    Vec2D Vec1, Vec2;
	    
            //--- for nodes only 
	    // if (NODES) 
	    //{
	    list<Vec2D>::iterator nLast, n1, n2;
	    nLast = lList[(*l)].node.end();
	    nLast--;
	    
	    if ( v2 == lList[(*l)].vertexJ )
	    {
		//cout << " take first node, print next one  " << endl;
		n1 = lList[(*l)].node.begin(); // first node of this line
		n2 = n1;
		n2++;
		//cout << (*n2).x << " " << (*n2).y << endl;
		
	    }
	    else if (v2 == lList[(*l)].vertexI )
	    {
		//cout << " take last node, print previous one " << endl;
		n1 = nLast; // last node of this line 
		n2 = n1 ;
		n2--;
		//cout << (*n2).x << " " << (*n2).y << endl;
	    }	
	    else
	    {
		cout << " problem here, dont find other vertex " << endl;
		exit(0);
	    }
	    
	    Vec1 = (*n1);
	    Vec2 = (*n2);
	    
            //}
	    //----

	    Vec2D dx = Vec2 - Vec1; 
	    //cout << " dx splitDir" << dx.x << " " << dx.y << endl;
	    
	    double h = dx.length() + 0.0000001;
	    
	    Vec2D bVec = dx * 0.2 / h; 
	    
	    Vec2D nVec = bVec.perp(); // normal to this line l 
	    double sig = sigma(aI, aJ, nVec);
	    double Dsig = Dsigma(aI, aJ, nVec);
	    Vec2D Ttemp = nVec * Dsig + bVec * sig;          
	    //cout << "Ttemp " << Ttemp.x << " " << Ttemp.y << endl;
	    
	    if (Ttemp.x != 0 && Ttemp.y != 0)
	    { // found good T, go to next line
		T.push_back(Ttemp);
		//cout << " length of Ttemp " << Ttemp.length()<< endl;
	    }
	    else // else T is 0, fake line! take next vertex of listOfV
	    {
		cout << " Ttemp is zero, we have a problem in findSplitDir " << endl;
		exit(0);
	    }    
	}
    }
    
    //--- now we have a proper T, get splitting direction
    splitInd = -5;
    double maxLength = -1.;
    double minLength = 1000.;
    
    // ofstream oFile;
//     oFile.open("../GNUPLOT/Tvec.dat", ios::out);
//     oFile.precision(8 );
    
    for(int k=0; k<T.size(); k++) 
    { 
	Vec2D tVec ; 
	double tAngle;
	// if (k == T.size()-1)
	// 	{
	// 	  tVec = T[k] + T[0] ;
	// 	}
	//       else
	// 	tVec = T[k+1] + T[k];
	
	if (k== T.size() -1)
	    tAngle = T[k].dotProduct(T[0]);
	else
	    tAngle = T[k+1].dotProduct(T[k]);
	
	double tLen = tVec.length();
	if( minLength > tAngle ) 
	{
	    if ( k == T.size()-1 )
		splitInd = k+1-T.size();
	    else
		splitInd = k+1;
	    
	    minLength = tAngle;
	    
	    // *0.1 * 10 because of fancy addition (to handle periodicity)
	    if (splitInd == T.size()-1)
	    {
		splitVec = T[splitInd]*0.1 + T[splitInd+1-T.size() ]*0.1;
		splitVec = splitVec*10.;
	    }
	    else
	    {
		splitVec = T[splitInd]*0.1 + T[splitInd+1]*0.1;
		splitVec = splitVec * 10.;
	    }
	}
    }
   
    // line splitInd and splitInd + 1  have smallest angle = largest T
    splitVec = splitVec / splitVec.length(); 
    splitVec = splitVec* (hMin/2.); 
      
} // Output: splitVec,  splitInd

//======================================================================

void defineNewVertexAndLine(Vertex VV, vector<Vertex> &vList, vector<Line> &lList, 
			    int gprev, int gmax, int gnext, int lmax, int lnext, 
			     Vec2D splitVec, int newV, int newL) 
{
    
    Vec2D newpos = VV.pos + splitVec; // *0.5;
     
    int vID=VV.id;
    //vList[vID].pos = vList[vID].pos - splitVec*0.5;

    int isIn = VV.isInside;
    if (newV != vList.size() )
	vList[newV].defineVertex(newV,isIn, gmax, gnext, gprev, lmax, lnext, newL,  newpos); 
    else
    {
	Vertex tempVertex;
	tempVertex.id = newV;
	tempVertex.grain.push_back(gmax);
	tempVertex.grain.push_back(gnext);
	tempVertex.grain.push_back(gprev);
	
	tempVertex.line.push_back(lmax);
	tempVertex.line.push_back(lnext);
	tempVertex.line.push_back(newL);
	tempVertex.pos = newpos;  
	tempVertex.isInside = isIn; 
	vList.push_back(tempVertex);
    }      

    //printVertex(vList[newV]);
    // define the new line
    // a new line is added from the new vertex to vertex number vID. 	

    Vec2D midpos =  VV.pos + splitVec*0.5;		
    	
    if (newL != lList.size())
	lList[newL].defineLine(newL, vID, newV, gnext, gprev, VV.pos, midpos, vList[newV].pos);
	//lList[newL].defineLine(newL, vID, newV, gnext, gprev, VV.pos, vList[newV].pos);
    else
    {
	Line tempLine;
	tempLine.id = newL;
	tempLine.vertexI = vID;
	tempLine.vertexJ = newV;
	tempLine.grainI = gnext;
	tempLine.grainJ = gprev;
	tempLine.node.clear();
	tempLine.node.push_back(VV.pos);
	tempLine.node.push_back(midpos);
	tempLine.node.push_back( vList[newV].pos);
		    
	//Initialize Normals for new line  	    
	Vec2D v1 = midpos - VV.pos;
	Vec2D n1 = v1.perp();
	n1 = n1 / n1.length();
	tempLine.normal.push_back(n1);
	
	Vec2D v2 = vList[newV].pos - midpos;
	//Vec2D v2 = vList[newV].pos - VV.pos;
	Vec2D n2 = v2.perp();
	n2 = n2 / n2.length();
	tempLine.normal.push_back(n2);	
	    
	// Initialize Velocities for new Line  		    
	tempLine.nodeVel.push_back(Vec2D(0.0,0.0));
	tempLine.nodeVel.push_back(Vec2D(0.0,0.0));
	tempLine.nodeVel.push_back(Vec2D(0.0,0.0));
	
	lList.push_back(tempLine);

	
    }
}
//======================================================================
void splitUnstableVertices(vector<Vertex> &vList, vector<Grain> &gList, 
			   vector<Line> &lList, vector<Pair> &IOqv)
{
    string PATH = "./results/";
    string confFileName;
    string underline = "_";
    string fName; 
    confFileName = "Vertex"; 

    int dummy2, dummy1;
    

    for(int i=0; i<vList.size(); i++) 
    {
		
	if( (vList[i].line).size() == 2 ) 
	{
	    cout << " we have an inner 2 line vertex in splitUnstables \n";
	    exit(0);
	}
	
	while( (vList[i].line).size() > 3 ) // split quadruple junctions or more 
	{
	    //cout << " split vertex " << endl;
	    //--- find splitting direction
	    Vec2D splitVec;  
	    int splitInd;
	    findSplitDir(vList[i].id, lList, vList, splitVec, splitInd);
	    
	    int lprevInd, lprevprevInd, lnextInd, lnextnextInd;                  
	    //lnextInd = (splitInd + 1)%(vList[i].line).size();
	    
	    if (splitInd != vList[i].line.size()-1)
		lnextInd = splitInd+1;
	    else
		lnextInd = splitInd-vList[i].line.size()+1;

	    //lnextnextInd = (splitInd + 2)%(vList[i].line).size();
	    if (splitInd != vList[i].line.size()-2)
		lnextnextInd = splitInd + 2;
	    else
		lnextnextInd = splitInd-vList[i].line.size()+2;

	    //lprevInd = (splitInd - 1)%(vList[i].line).size();
	    if (splitInd != 0)
		lprevInd = splitInd-1;
	    else
		lprevInd = splitInd + vList[i].line.size()-1;
	    
	    //lprevprevInd = (splitInd - 2)%(vList[i].line).size(); 	
	    if (splitInd > 1)
		lprevprevInd = splitInd - 2;
	    else
		lprevprevInd = splitInd + vList[i].line.size()-2;

	    vector<int>::iterator gcenterIt, splitIt;
	    
	    int thisV = vList[i].id;
	    
	    //---get indices (names) for important lines,grains
	    
	    int gprev, gcenter, gnext;
	    int lnext,lnextnext, lprev, lprevprev, splitL;
	    
	    gcenter = vList[thisV].grain[splitInd];
	    gnext = vList[thisV].grain[lnextInd];
	    gprev = vList[thisV].grain[lprevInd];
	    
	    //cout<< " gcenter, gnext, gprev " 
	    //	<< gcenter << " " << gnext << " " << gprev << endl;

	    splitL = vList[thisV].line[splitInd]; // the name of the line
	    lnext = vList[thisV].line[lnextInd];
	    lnextnext = vList[thisV].line[lnextnextInd];
	    lprev = vList[thisV].line[lprevInd];
	    lprevprev = vList[thisV].line[lprevprevInd];
	    
	    // get number of fake lines 
	    //int counter=0;
	    //int totalNLines = vList[thisV].line.size();
	    // for (int k=0; k<totalNLines; k++) /// Eva 10-1-2008 not needed anywhere 
// 	      if (vList[thisV].line[k]<=-1)
// 		counter++; 
	    
	    int newV, newL ;
	    
	    defNewVertIndex(newV, vList, EmptyVertex );
	    int newVInd = newV;
	    defNewLineIndex(newL, lList, EmptyLine );
	    int newLInd=newL;
	    
	    
	    defineNewVertexAndLine(vList[thisV], vList, lList, gprev, gcenter, gnext, 
				   splitL, lnext, splitVec, newV,newL);
	    	    
	    // thisV was moved by move = splitVec to create newV
	    // adjust splitL and lnext by this move

	    //-- update vertex thisV: thisV.line and thisV.grain 
	    gcenterIt = vList[thisV].grain.begin();
	    splitIt = vList[thisV].line.begin();
	    for(int ii=0; ii < splitInd; ii++) 
	      {
		gcenterIt++;
		splitIt++;
	      }
	    
	    vList[thisV].line[lnextInd] = newL;
	    (vList[thisV].line).erase(splitIt);
	    
	    (vList[thisV].grain).erase(gcenterIt);
	    
	    //---update lines splitL, lnext 
	    //if (NODE)
	    list<Vec2D>::iterator ll, lf;
	    //---
	    if ( splitL > -1 )
	    {
		if( thisV == lList[splitL].vertexI ) 
		{
		    lList[splitL].vertexI = newV;
		    // if (NODE) 
		    ll = lList[splitL].node.begin();
		    (*ll) = vList[newV].pos;
		    //---
		}
		else if( thisV == lList[splitL].vertexJ) 
		{
		    lList[splitL].vertexJ = newV;
		    //if (NODE)
		    lf = lList[splitL].node.end();
		    lf--;
		    (*lf) = vList[newV].pos;
		    //---
		}
		//cout << " splitL " << splitL << endl; 
		//printLine(lList[splitL]);
		
		//cout << " in SplitUnstable " ;
		spaceNodesEqually(lList, splitL);

		// or adjust? 
		int vv = vList[newV].otherVertex(lList[splitL]);
		//cout << " newV " << newV << " vv " << vv << endl;

		//adjustNodes(lList, splitVec, splitL, vv);
		
	    }
	    if ( lnext > -1 )
	    {
		//cout << " lnext " << lnext << endl;
		
		if( thisV == lList[lnext].vertexI ) 
		{
		    lList[lnext].vertexI = newV;
		    // if (NODE) 
		    ll = lList[lnext].node.begin();
		    (*ll) = vList[newV].pos;
		    //---
		}
		else if( thisV == lList[lnext].vertexJ) 
		{
		    lList[lnext].vertexJ = newV;
		     //if (NODE)
		    lf = lList[lnext].node.end();
		    lf--;
		    (*lf) = vList[newV].pos;
		    //---
		}
		// printLine(lList[lnext]);
		spaceNodesEqually(lList, lnext);
		
                // or adjust? 
		int vv = vList[newV].otherVertex(lList[lnext]);
		//cout << " newV " << newV << " vv " << vv << endl;

		//adjustNodes(lList, splitVec, lnext, vv);
	    }
	    //Update grain gcenter: replace VV.id with new vertex 
	    
	    vector<int>::iterator gv,gl;
	    if (gList[gcenter].id > -1) 
	    {
		gv = find((gList[gcenter].vertex).begin(),(gList[gcenter].vertex).end(),thisV);
		if (gv != (gList[gcenter].vertex).end() )
		    (*gv) = newV; //!!!!
		else
		{
		    cout << " could not find gv in gLIst[gcenter].vertex " << endl;
		    exit(0);
		}
	    }
	    // Update gprev and gnext 
	    if ( gList[gprev].id >-1 ) 
	    {
		//printGrain(gList[gprev]);
		gv = find((gList[gprev].vertex).begin(),(gList[gprev].vertex).end(),thisV);
		if (gv != gList[gprev].vertex.end())
		    (gList[gprev].vertex).insert(gv,newV); // insert newV before position gv
		else
		{
		    cout << " could not find thisV " <<thisV << " in vertex list of gprev " 
			 << gList[gprev].id << endl;
		    exit(0);
		}  
		gl = find((gList[gprev].line).begin(),(gList[gprev].line).end(), lprev);
		if (gl != (gList[gprev].line).end() )
		    (gList[gprev].line).insert(gl ,newL);
		else
		{
		    cout << " could not find lprev " << lprev << " in line list of gprev " 
			 << gList[gprev].id << endl;
		    exit(0);
		}  
		//	printGrain(gList[gprev]);
	      }
	    
	    if ( gList[gnext].id > -1 ) 
	    {
		//cout << " thisV " << thisV << endl;
		//printGrain(gList[gnext]);
		gv = find((gList[gnext].vertex).begin(),(gList[gnext].vertex).end(),thisV);
		if (gv != gList[gnext].vertex.end())
		{
		      
		    gv++;
		    //cout << " gv++ " << (*gv) << endl;
		    (gList[gnext].vertex).insert(gv,newV); // insert newV before position gv
		}
		else
		{
		    cout << " could not find thisV in vertex list of grain " 
			 << gList[gnext].id << endl;
		    exit(0);
		}
		gl = find((gList[gnext].line).begin(),(gList[gnext].line).end(), lnext);
		if (gl != gList[gnext].line.end())
		{
		    (gList[gnext].line).insert(gl ,newL);
		    //printGrain(gList[gnext]);
		}
		else
		{
		    cout << " could not find lnext in line list of grain " 
			 << gList[gnext].id << endl;
		    exit(0);
		}
	    }
	    //printVertex(vList[i]);
	} // end of while  
	
    }  
   
} 

//============== FUNCTIONS IN remove a.o.
 
int findPosition(int l, vector<int> vec)
{
  int res = -1; 
  for(int i=0; i<vec.size(); i++) 
    if( vec[i] == l) 
    {
      res = i; 
      break; 
    }
  return res;
}
//================================================================================
void removeLineAndVertexFromGrain(vector<Grain> &gList, vector<Line> &lList,int g,
				  int l, int v)
{

  int shift =0;
    if( gList[g].id > -1 )
    { 
	//printGrain(gList[g]);

	vector<int>::iterator ll = gList[g].line.end(); // last line of g
	ll--;
	vector<int>::iterator vv = gList[g].vertex.begin(); // 1. vertex of g
    
	vector<int>::iterator itL=find((gList[g].line).begin(),(gList[g].line).end(),l);
	vector<int>::iterator itV=find((gList[g].vertex).begin(),(gList[g].vertex).end(),v);

	if (itL == ll && itV == vv)
	  shift = 1; // we need to shift lines by one to keep the correct order

	if (itL != gList[g].line.end() )
	    (gList[g].line).erase(itL);
	else
	{
	    cout << " could not find this line " << l << " in gList.line \n" << endl;
	    exit(0);
	}    
	
	if (itV != (gList[g].vertex).end())
	    (gList[g].vertex).erase(itV);
	else
	{
	    cout << " could not find this vertex " << v << " in gList.vertex  \n" << endl;
	    printGrain(gList[g]);
	    exit(0);
	}    
    }

//    printGrain(gList[g]);

//     // check order of lines and vertices in grain

    
    int lFirst = gList[g].line[0];

//     printLine(lList[lFirst]);

//     int vFirst = gList[g].vertex[0];
//     int vSecond = gList[g].vertex[1];

//     int vI = lList[lFirst].vertexI;
//     int vJ = lList[lFirst].vertexJ;
    
//     //if ( l == (*ll) && v == (*vv)   )
//     if (!((vFirst == vI && vSecond == vJ) || (vFirst == vJ && vSecond == vI)))

    if (shift == 1)
    {
	(gList[g].line).push_back(lFirst);
	gList[g].line.erase(gList[g].line.begin());
    }
    
}
//=========================================================================================== 
void replaceVertexInGrains(vector<Grain> &gList, vector<int> tempG, int vI, int vJ)
{
    vector<int>:: iterator vv;
    for(int k=0; k<tempG.size(); k++) 
	
	if( tempG[k] > -1) 
	{
	    vv=find(gList[tempG[k]].vertex.begin(), gList[tempG[k]].vertex.end(), vJ);
	    (*vv)=vI;
	    
	}
}
//=========================================================================================== 
void replaceVertexInLines(vector<Line> &lList, vector<Vertex> &vList, vector<int> tempL, int vI, int vJ)
{ 
    
    for(int k=0; k<tempL.size(); k++) 
    	if( tempL[k] > -1) 
	{
	    //if (NODE)
	    list<Vec2D>::iterator lf, ll;
	    //--- 
	    //cout << "replaceVertexInLines " << tempL[k] << endl;
	    //cout << " lList[tempL[k]])" << endl;
	    //printLine(lList[tempL[k]]);

	    if (lList[tempL[k]].vertexI == vJ)
	    {
		    lList[tempL[k]].vertexI = vI;
		    // if (NODE)
                    //--- update first node
		    lf = lList[tempL[k]].node.begin();
		    (*lf) = vList[vI].pos;
		    //---
	    }
	    else if ( lList[tempL[k]].vertexJ == vJ)
	    {
		lList[tempL[k]].vertexJ = vI;
		// if (NODE)
		//--- update last node
		ll = lList[tempL[k]].node.end();
		ll--;
		(*ll) = vList[vI].pos;
		//---
	    }
	    else
	    {
		cout << "error in removeShortLines, could not update line "<< tempL[k] <<endl ;
		exit(0);
	    } 

	    
	}
}
//===========================================================================================
void remove(int l, vector<Vertex> &vList, vector<Grain> &gList,  vector<Line> &lList, 
	    vector<Pair> &IOqv, int execute, vector<int> &gCrit, vector<int> &vCrit)
{
    
//     cout << " function remove for line " << l << endl;
    //printLine(lList[l]);

    int idL = lList[l].id;
    int gI = lList[idL].grainI; 
    int gJ = lList[idL].grainJ;
    
    
    int vI = lList[idL].vertexI; 
    int vJ = lList[idL].vertexJ;

    //printVertex(vList[vI]);
  //  printVertex(vList[vJ]);

    // information of critical grains, vertices for further processing
    //  2-sided grain, 2-line vertex :  gI, gJ, vI
    
    int posLineI = findPosition(idL, vList[vI].line); 
    int posLineJ = findPosition(idL, vList[vJ].line); 

    

    //---------------------
    if (gList[gI].id > -1)
	removeLineAndVertexFromGrain(gList, lList, gI, lList[l].id, vJ);
    if (gList[gJ].id > -1) 
	removeLineAndVertexFromGrain(gList, lList, gJ, lList[l].id, vJ);
    
    // insert lines of vJ into vI.line
    vector<int> tempL;
    int vJLineSize = vList[vJ].line.size();
    int ll = posLineJ + 1; 
    
    for(int k = 0; k < vJLineSize - 1; k++) 
	tempL.push_back(vList[vJ].line[ (ll+k) % vJLineSize]); 
    
    vector<int>::iterator itL;
    itL = find((vList[vI].line).begin(), (vList[vI].line).end(), lList[l].id); 
    (vList[vI].line).insert(itL,tempL.begin(), tempL.end());
    
    itL = find((vList[vI].line).begin(), (vList[vI].line).end(), lList[l].id); 
    (vList[vI].line).erase(itL);
    
    // insert grains of vJ into vI.grain 
    
    vector<int> tempG;
    
    int vJGrainSize = vList[vJ].grain.size();
    ll = posLineJ + 1;   
    for(int k = 0; k < vJGrainSize - 2; k++)  
    {
	//if ( (ll+k) < vJGrainSize )
	tempG.push_back(vList[vJ].grain[ (ll+k) % vJGrainSize]);  
    }
    int idG = vList[vI].grain[posLineI];
    itL = find((vList[vI].grain).begin(), (vList[vI].grain).end(), idG ); 
    (vList[vI].grain).insert(itL,tempG.begin(),tempG.end());
    
    replaceVertexInGrains(gList, tempG, vI, vJ);
    replaceVertexInLines(lList, vList, tempL, vI, vJ); 


 //    // adjust nodes of lines in tempL equally
//     Vec2D move = vList[vI].pos - vList[vJ].pos;

    //cout << " tempL.size() " << tempL.size() << endl;
    for (int k=0; k<tempL.size(); k++)
    { 

// 	int vv = vList[vI].otherVertex(lList[tempL[k]]);
// 	cout << " adjustNodes in remove " << endl;
// 	adjustNodes(lList, move, tempL[k], vv);


	spaceNodesEqually(lList,tempL[k]);
    }
    
    
    if (lList[l].id > -1)
	lList[l].kill();
    
    if (vList[vJ].id > -1) 
	vList[vJ].kill();
    
    
    if (gList[gI].line.size() == 2)
	gCrit.push_back(gI);
    
    
    if (gList[gJ].line.size() == 2)
	gCrit.push_back(gJ);
    
    
    if (vList[vI].line.size() == 2)
	vCrit.push_back(vI);
    
}


//==================================================================================== 
void removeShortLines(vector<Vertex> &vList,vector<Grain> &gList,vector<Line> &lList, 
		      vector<Pair> &IOqv, int execute, double deltaTl)
{
    int idL;
    Line thisLine;
    vector<int>::iterator itL;
    //cout << " hMin in removeShort " << hMin << endl; 
   
    for(int l=0; l < lList.size(); l++)
    {
	vector<int> IlistOfV, JlistOfV, tempJlist;
	if ( lList[l].id != -1 ) // line exists  
	{
	    //--- calculate length of line 
	    list<Vec2D>::iterator  nnNext, nn;
	    list<Vec2D>::iterator nnLast=lList[l].node.end();
	    nnLast--;
	    	    
	    
	    double dx=0.;
	    Vec2D dxVec=Vec2D(0.0,0.0);
	    for ( nn = lList[l].node.begin(); nn!= nnLast; nn++)
	    {
		nnNext = nn;
		nnNext++; 
		dx += ((*nnNext)-(*nn)).length();
	    }
	    
	    double hLoc = dx/(lList[l].node.size()-1);
	  	
	    dxVec   = vList[lList[l].vertexI].pos - vList[lList[l].vertexJ].pos;
	    Vec2D lVel = vList[lList[l].vertexI].vel - vList[lList[l].vertexJ].vel;

	    // check if its a too small 3 node line
	    int kill = 0;
	    
	    if (lList[l].node.size() == 3 && hLoc < 0.75*hMin)
	    {
		kill = 1;
		//	cout << " need to kill this small 3 node line " << l << endl;
		
	    }

	                 	                  
	    if( dxVec.length() < epsLine && lVel.dotProduct(dxVec) < 0) // || kill == 1) 
	    //	if( (dxVec.length() < epsLine && lVel.dotProduct(dxVec) < 0) || kill == 1) 
	    {
		//cout << " remove l " << l <<  endl;
		vector<int> gCrit;
		vector<int> vCrit;
		int gMarker = -1;
	
		remove(l, vList,gList, lList, IOqv, execute, gCrit, vCrit);
		
		if ( gCrit.size()!= 0 || vCrit.size()!= 0)
		  gMarker = 1;
				
				
		while( gMarker == 1 ) 
		{ //--- handle 2-sided grains and 2-line vertices
		   
		      for (int j=0; j<gCrit.size(); j++)
		      {
			if ( gList[gCrit[j] ].id > -1  ) 
			    remove2SidedGrain(gCrit[j], gList, vList, lList, vCrit );
		      }
		    
		    gCrit.clear();
		    
		    //--- handle 2-line vertices
		    for (int jj=0; jj<vCrit.size(); jj++)
		    	Treat2LineVertex(vCrit[jj], vList, lList, gList, gCrit);
		    
		    vCrit.clear();
		    //cout << " size of vCrit after clear " << vCrit.size() << endl;
		    //cout << " size of gCrit after Treat2LineVertex " << gCrit.size() << endl;
		    
		    gMarker = -1;
		    if (gCrit.size() != 0)
		    {
			gMarker = 1; // continue while-loop, check again for 2-sided grain
		    }
		    
		} // end of while
	    }  // end of dx
	   
	} // end of lList.id > -1
    } // end of loop over lList 
    
}

//======================================================================================
void Treat2LineVertex(int i, vector<Vertex> &vList, vector<Line> &lList, vector<Grain> &gList, 
			vector<int> &gCritList)
{
    Vertex VV = vList[i];
    int vID = vList[i].id;
    
   
    int l1 = (VV.line).front();
    int l2 = (VV.line).back(); // line to be erased 
    int v11 = lList[l1].vertexI; int v12 = lList[l1].vertexJ;
    int v21 = lList[l2].vertexI; int v22 = lList[l2].vertexJ;
    int g11 = lList[l1].grainI;  int g12 = lList[l1].grainJ;
    int g21 = lList[l2].grainI;  int g22 = lList[l2].grainJ;
    
	
    vector<int>::iterator iLJ;
    list<Vec2D>::iterator n;
    
   // cout << " work on line " << lList[l1].id << endl;
   // printLine(lList[l1]);

   // Extend L1, delete L2, 
    
    if(!(v11 == v21 || v11== v22 || v12 == v21 || v12 == v22 ) )			
    {
	cout << " we have a problem, vertex/line dont coincide, think! " << endl;
	exit(0);
    }
    
    //replace VV in l1 with correct ending point of l2   
    // and update/ redistribute its nodes

    list<Vec2D>::iterator no;
    // case a 
    if( v11 == vID && v21 == vID)
    {
	lList[l1].vertexI = v22; // new other node, replace first node
	// update nodes of l1: insert nodes of l2 at right pos in l1.node
	// same for nodeVel (even if recalculated, the dimension has to be right)
	// reverse list lList[l2].node and insert before first node of l1
		

	//--- update nodes
	lList[l2].node.pop_front();
	(lList[l2].node).reverse();
	
	(lList[l1].node).splice(lList[l1].node.begin(),lList[l2].node ); //update nodes

       
	
	//--- update nodeVel
	lList[l2].nodeVel.pop_front();
	(lList[l2].nodeVel).reverse();
	(lList[l1].nodeVel).splice(lList[l1].nodeVel.begin(),lList[l2].nodeVel );

	
	//--- update lines of vertex 
	iLJ = find((vList[v22].line).begin(),(vList[v22].line).end(),l2);
	if (iLJ != vList[v22].line.end() ) 
	    (*iLJ) = l1;
	else
	{
	    cout << " could not find l2 in vList[v22].line " << endl;
	    exit(0);
	}
    }
    // case b 
    else if( v11 == vID && v22 == vID)
    {	 
	lList[l1].vertexI = v21;
		       
	
	//--- update nodes 
	(lList[l2].node).pop_back();
	(lList[l1].node).splice(lList[l1].node.begin(),lList[l2].node );

       
	//--- update nodeVel
	(lList[l2].nodeVel).pop_back();
	(lList[l1].nodeVel).splice(lList[l1].nodeVel.begin(),lList[l2].nodeVel );

	iLJ = find((vList[v21].line).begin(),(vList[v21].line).end(),l2);
	if (iLJ != vList[v21].line.end() ) 
	    (*iLJ) = l1;
	else
	{
	    cout << " could not find l2 in vList[v21].line " << endl;
	    exit(0);
	}

	
    }
    // case c 
    else if( v12 == vID && v21 == vID)
    {
	lList[l1].vertexJ = v22;


	lList[l2].node.pop_front();
	(lList[l1].node).splice(lList[l1].node.end(),lList[l2].node ); //update nodes

	//--- update nodeVel
        lList[l2].nodeVel.pop_front();
	(lList[l1].nodeVel).splice(lList[l1].nodeVel.end(),lList[l2].nodeVel ); 
		
	//--- update lines of vertex
	iLJ = find((vList[v22].line).begin(),(vList[v22].line).end(),l2);
	if (iLJ != vList[v22].line.end() )
	    (*iLJ) = l1;
	else
	{
	    cout << " could not find l2 in vList[v22].line " << endl;
	    exit(0);
	}
    }
    // case d 
    else if( v12 == vID && v22 == vID)
    {	
	
	lList[l1].vertexJ = v21;
	
       	//--- update node 
	lList[l2].node.pop_back();
	(lList[l2].node).reverse();
	(lList[l1].node).splice(lList[l1].node.end(),lList[l2].node );

	//--- update nodeVel
	lList[l2].nodeVel.pop_back();
	(lList[l2].nodeVel).reverse();
	(lList[l1].nodeVel).splice(lList[l1].nodeVel.end(),lList[l2].nodeVel );

	//--- update lines of vertex 
	iLJ = find((vList[v21].line).begin(),(vList[v21].line).end(),l2);
	if (iLJ != vList[v21].line.end() )
	    (*iLJ) = l1;
	else
	{
	    cout << " could not find l2 in vList[v21].line " << endl;
	    exit(0);
	}
    }
    
    
    //--- erase l2 and vertex vID from grain l1.gI and l1.gJ
    removeLineAndVertexFromGrain(gList, lList, g11, l2, vID);
    removeLineAndVertexFromGrain(gList, lList, g12, l2, vID);
    
    

    // check here if we have a 2-sided grain now.!!!!!
    if (gList[g11].line.size() == 2)
	gCritList.push_back(g11);

    if (gList[g12].line.size() == 2)
	gCritList.push_back(g12);
   
    vList[i].kill();
    lList[l2].kill();
 
    
    spaceNodesEqually(lList, l1);

}
//===========================================================================
void spaceNodesEqually(vector<Line> &lList, int &l)
{
    
    list<Vec2D> newNodes, newVelos;
    Vec2D zero = Vec2D(0.0,0.0);
    
    // Calculate mesh factor 
    // if # nodes != 2^n, get new nPts and new h for this line
    
    double totalLen = 0.;
    double hMin = 1.;
    double hMax = 0.;
    double h=0.;
    list<Vec2D>::iterator no, nPrev, nLast, nLastLast, nNext;
    nLast = lList[l].node.end();
    nLast--;

    nLastLast = nLast;
    nLastLast--;           // node before! last one
    for (no=lList[l].node.begin(); no != nLast; no++)
    {
	nNext = no;
	nNext++;
	h = ((*nNext)-(*no)).length();
	
	totalLen = totalLen + h; // total Length of line l
	
	if (h < hMin)
	    hMin = h;
	if (h > hMax)
	    hMax = h;
    }
    
    h = (hMax+ hMin)*0.5;
    
    //--- find proper nPts (2^n +1 points)
    double  dn = ( log(totalLen/(h))  / log(2.)) + 0.5;
    int n = (int) (dn);
    n = max(1,n);
    dn = double(n);
    dn  = pow(2., dn );
    
    int nPts = (int) (dn) + 1;
    
    // if (nPts != lList[l].node.size())
    //	cout << " change nPts for this line " << endl;

    double hNew = totalLen/ ((double) (nPts-1));
   
    newNodes.push_back(*(lList[l].node.begin())); // first node is the same
    newVelos.push_back(zero);

    list<Vec2D>::iterator n1=lList[l].node.begin();
    list<Vec2D>::iterator n2;
    double sum=0.;
    double prevSum, step;

    no=lList[l].node.begin();
    for (int i=0; i<nPts-1; i++)
    {
	while (sum < hNew*(i+1) ) // keep on adding h's
	{
	    prevSum = sum;
	    
	    if (no != nLast)
	    {
		nPrev = no;
		no++;
		sum = sum + ((*no)-(*nPrev)).length();
	    }
	    else
		break;
	}
	if ((sum < totalLen || i < nPts -2)  )
	{
	    step = hNew*(i+1) - prevSum ;
	    Vec2D nNew = (*no)-(*nPrev);
	    nNew = nNew/nNew.length();
	    nNew = (*nPrev) + nNew * step;
	    
	    newNodes.push_back(nNew); 
	    newVelos.push_back(zero);
	}
	else
	{
	    newNodes.push_back((*nLast));
	    newVelos.push_back(zero);
	    break;
	}
    }
        
    (lList[l].node).assign(newNodes.begin(),newNodes.end());
    (lList[l].nodeVel).assign(newVelos.begin(),newVelos.end());

}

//======================================================================================
void remove2SidedGrain(int g0, vector<Grain> &gList, vector<Vertex> &vList, 
		       vector<Line> &lList, vector<int> &vCritList )
{
    int l1,l2;
    vector<int> vCritTemp;
    
    for (int kk=0; kk<gList[g0].vertex.size(); kk++)
    	vCritTemp.push_back(gList[g0].vertex[kk]);
    
    l1 = (gList[g0].line).front(); 
    l2 = (gList[g0].line).back(); //line to remove 
    
    int g1, g2;
    if (g0 ==lList[l1].grainI)
	g1 = lList[l1].grainJ;
    else if (g0 == lList[l1].grainJ)
	g1 = lList[l1].grainI;
    else
    {
	cout << " BUG, could not find g1 in L1.grain\n";
	exit(0);
    }
    if (g0 == lList[l2].grainI )
	g2 = lList[l2].grainJ;
    else if (g0 == lList[l2].grainJ)
	g2 = lList[l2].grainI;
    else
    {
	cout << " BUG, could not find g2 in l2.grain\n";
	exit(0);
    } 
    
    // determine vertex v2 (other ending vertex of l1 than vI)
    int vI = lList[l1].vertexI;
    int v2 = -5;
    v2 = vList[vI].otherVertex(lList[l1]);
    
    if ( v2 == -5) 
    {
 	cout << " BUG, could not find v2 in l1.vertexI/J \n";
 	exit(0);
    } 
    
    if (vI != (gList[g0].vertex).front() && vI != (gList[g0].vertex).back() )
    {
	cout << " BUG, could not find vI in g0.vertex \n";
 	exit(0);
    }
    if (v2 != gList[g0].vertex.front() && v2 != gList[g0].vertex.back() )
    {
	cout << " BUG, could not find v2 in g0.vertex \n";
 	exit(0);
    }
    
    //--- replace line l2 of 2sided grain by line l1 in grain g2
    if (gList[g2].id > -1)
    {
	//printGrain(gList[g2]);
	vector<int>::iterator itL;
	itL = find( (gList[g2].line).begin(),(gList[g2].line).end(),l2);
	if (itL != gList[g2].line.end() )
	{
	    int posl2_g2 = findPosition((*itL), gList[g2].line);
	    gList[g2].line[posl2_g2] = l1;
	}
	else
	{
	    cout << " could not find this line " <<l2<< " in gList[g2].line \n" ;
	    exit(0);
	} 
	
    }
    
    //--- remove grain g0 and line l2 from vertex vI
    vector<int>::iterator eG,eL,it;  
    //printVertex(vList[vI]);
    eG= find((vList[vI].grain).begin(),(vList[vI].grain).end(),g0);
    eL= find((vList[vI].line).begin(),(vList[vI].line).end(),l2);
    
    int posG, posL;
    posG = findPosition((*eG), vList[vI].grain);
    posL = findPosition((*eL), vList[vI].line);
    
    if (posG != posL)
    {
	int size = vList[vI].grain.size();
	
	if (posL == 0 && posG == size-1 )
	{   // grain position is one before line position
	    // l2 > l1, grain g0 between them
	    
	    it = (vList[vI].grain).begin(); 
	    int grain=(*it);
	    (vList[vI].grain).erase(it);
	    (vList[vI].grain).push_back(grain);
	}
	
    }
    
    eG= find((vList[vI].grain).begin(),(vList[vI].grain).end(),g0);
    eL= find((vList[vI].line).begin(),(vList[vI].line).end(),l2);
    if ( eG != vList[vI].grain.end() )
	(vList[vI].grain).erase(eG); 
    else
    {
	cout << " grain g0 " << g0 
	     << " does not exit in vList[v1].grain, cannot be erased "<< endl;
	printVertex(vList[vI]);
	exit(0);
    }
    if (eL != vList[vI].line.end() )
	(vList[vI].line).erase(eL); 
    else
    {
	cout << " line eL " << l2 
	     << " does not exit in vList[vI].line, cannot be erased "<< endl;
	printVertex(vList[vI]);
	exit(0);
    }
    
    
    //--- remove grain g0 and line l2 from vertex v2
    
    eG= find((vList[v2].grain).begin(),(vList[v2].grain).end(),g0);
    eL= find((vList[v2].line).begin(),(vList[v2].line).end(),l2);
    posG = findPosition((*eG), vList[v2].grain);
    posL = findPosition((*eL), vList[v2].line);
    
    if (posG != posL)
    {
	int size = vList[v2].grain.size();
	
	if (posL == 0 && posG == size-1 )
	{   // grain position is one before line position
	    // l2 > l1, grain g0 between them
	    
	    it = (vList[v2].grain).begin(); 
	    int grain=(*it);
	    (vList[v2].grain).erase(it);
	    (vList[v2].grain).push_back(grain);
	    
	}
    }
    eG= find((vList[v2].grain).begin(),(vList[v2].grain).end(),g0);
    eL= find((vList[v2].line).begin(),(vList[v2].line).end(),l2);
    if (eG != vList[v2].grain.end())
	(vList[v2].grain).erase(eG); 
    else
    {
	cout << " grain g0 " << g0 << " does not exist in v2" << endl;
	exit(0);
    }
    
    if (eL != vList[v2].line.end())
	(vList[v2].line).erase(eL); 
    else
    {
	cout << " line eL " << l2 
	     << " does not exist in vList[v2].line, cannot be erased "<< endl;
	printVertex(vList[v2]);
	exit(0);
    }
    
    //--- update line l1, replace g0 by g2
    if (g0 ==lList[l1].grainI)
	lList[l1].grainI = g2;
    else if (g0 == lList[l1].grainJ)
	lList[l1].grainJ = g2;
    else
    {
	cout << " could not find g0 in grainlist of l1 " << endl;
	exit(0);
    }
    
    //---- check vertices in  critList, if 2-line vertex
    for (int i=0; i<vCritTemp.size(); i++)
    	if (vList[vCritTemp[i]].line.size() == 2) 
	    vCritList.push_back(vCritTemp[i]); //2-line vertex
    
    if (lList[l2].id > -1)
	lList[l2].kill();
    
    if (gList[g0].id > -1)
	gList[g0].kill();
}

//============================================================================
void getVeloVert(vector<Vertex> &vList , vector<Line> &lList)
{
    //--- calc velocities on vertices  
    for( int v=0; v < vList.size(); v++) 
    {
	Vec2D vel    = Vec2D(0.0, 0.0);
	vList[v].vel = Vec2D(0.0, 0.0);

	if (vList[v].id > -1)
	{
	    list<Vec2D>::iterator nodeVel, nLast, nTemp;
	   
	    for (int j=0; j<vList[v].line.size(); j++)
	    {
		int l = vList[v].line[j];
		int vv = vList[v].otherVertex(lList[l]);
		
		// get velocity of first/last node - depending on the direction
		if ( vv == lList[l].vertexJ)
		    nodeVel = lList[l].nodeVel.begin(); // first node of this line
		else if (vv == lList[l].vertexI)
		{
		    nLast = lList[l].nodeVel.end();
		    nLast--;
		    nodeVel = nLast;
		}
		else
		{
		    cout << "Bug finding other vertex" << endl;
		    exit(0);
		}
		vel = vel + (*nodeVel); // add all vel's of lines of this triple junction
	    }
	    vList[v].vel = vel;
	    
	}
    }

}



//====================== EVOLUTION ===========================================

void calcVelocities(vector<Line> &lList)
{
    //--- calculate velocities of the nodes 
    Vec2D zero = Vec2D(0.0,0.0);        
    list<Vec2D>::iterator l;
    for (int i=0; i<lList.size(); i++)
    	for (l =lList[i].nodeVel.begin(); l != lList[i].nodeVel.end(); l++)
	    (*l) = zero; 
  
    for(int l=0; l<lList.size(); l++)
     	if (lList[l].id > -1) 
	{
	    double aI = lList[l].alphaI; 
	    double aJ = lList[l].alphaJ;

	    list<Vec2D>::iterator n, n2, nLast; 
	    list<Vec2D>::iterator nVel, nVelNext;	    
	    nLast = lList[l].node.end();
	    nLast--;
	    nVel = lList[l].nodeVel.begin();	   
 
	    for(n= (lList[l].node).begin(); n != nLast; n++)
	    {
		n2 = n;
		n2++;
		if (n2 == NULL)
		    cout << " not defined " << endl;
		Vec2D dx = (*n2) - (*n); 		
		nVelNext = nVel;
		nVelNext++;		
		double h = dx.length() + 0.0000001;
		if( h > 0 ) 
		{
		    Vec2D bVec = dx * 0.2 / h;  
		    //normalized vector between node1 and node2
		    Vec2D nVec = bVec.perp();
		    nVec = nVec * 0.2 / nVec.length(); // normal  

		    double sig = sigma(aI, aJ, nVec);
		    double Dsig = Dsigma(aI, aJ, nVec);

		    Vec2D Tvec = nVec * Dsig + bVec * sig;  // main formula
		    (*nVel)      = (*nVel)    + (Tvec*0.1);
		    (*nVelNext) = (*nVelNext) - (Tvec*0.1);
		}
		else
		{
		    cout << "h < =0 " << endl;
		    exit(0);
		}
		nVel++;
	    }  

	    //--- multiply by 10 to neutralize division by 0.1 above 
	    nLast = lList[l].nodeVel.end();
	    for ( nVel= (lList[l].nodeVel).begin(); nVel != nLast; nVel++ )
		(*nVel) = (*nVel) * 10.;
	}
        
}
// =======================================================================

Vec2D getMiddlePoint(Line &ll)
{
    Vec2D nMid;
    list<Vec2D>::iterator nLast, nFirst;
    nFirst = ll.node.begin();
    nLast = ll.node.end();
    nLast--;
    
    double tmpx = (*nLast).x - (*nFirst).x; // distance switched nodes
    double tmpy = (*nLast).y - (*nFirst).y;
    
    if (abs(tmpx) < 1. && abs(tmpy) < 1.)
    {
	nMid = ((*nFirst) + (*nLast))*0.5;
    }
    else //--- find correct middle points switch nodes
    {
	if (abs(tmpx) > 1. && abs(tmpy) <= 1. )
	{  
	    double new_nLast_x;
	    
	    if ((*nLast).x < 0)
		new_nLast_x = (*nLast).x + 2.;
	    else
		new_nLast_x = (*nLast).x - 2.;
	    nMid.x = ( new_nLast_x + (*nFirst).x ) *0.5 ;
	    nMid.y = ((*nLast).y   + (*nFirst).y ) *0.5 ;
	}
	if (abs(tmpx) <= 1 && abs(tmpy) > 1 )
	{   
	    
	    double new_nLast_y;
	    if ((*nLast).y < 0)
		new_nLast_y = (*nLast).y + 2.;
	    else
		new_nLast_y = (*nLast).y - 2.;
	    nMid.x = ((*nLast).x   + (*nFirst).x ) *0.5 ;
	    nMid.y = ( new_nLast_y + (*nFirst).y ) *0.5 ;
	    
	}
	
	if (abs(tmpx) > 1 && abs(tmpy) > 1 )
	{   
	    double new_nLast_x, new_nLast_y;
	    if ((*nLast).x < 0)
		new_nLast_x = (*nLast).x + 2.;
	    else
		new_nLast_x = (*nLast).x - 2.;
	    
	    if ((*nLast).y < 0)
		new_nLast_y = (*nLast).y + 2.;
	    else
		new_nLast_y = (*nLast).y - 2.;
	    nMid.x  =( new_nLast_x + (*nFirst).x) *0.5;
	    nMid.y  =( new_nLast_y + (*nFirst).y) *0.5;
	}
	cout << " first " << (*nFirst).x << " " << (*nFirst).y << endl;
	cout << " nMid " << nMid.x << " " << nMid.y << endl;
	cout << "last  "  << (*nLast).x << " " << (*nLast).y << endl;

	cout << endl;
    } // end of else 

    return nMid;
}

//=================================================================================
    
void evolveInterior(vector<Line> &lList, vector<Vertex> &vList, double &deltaTl, double extime)
{
    calcVelocities(lList);             //--- calculate 

    // deltaTl = deltaTl/100.;   // reduce timestep
    // deltaTl = .3*(hMin*hMin);
    current_time += deltaTl;

    //--- evolve inner points
    list<Vec2D>::iterator nn, nnLast, nnNext;
    
    double linetotal=0;
    double counter=0; 
    
    for(int l=0; l<lList.size(); l++)
     	if (lList[l].id > -1) 
	{
	    linetotal++;
	    list<Vec2D>::iterator n, nFirst, nLast;
	    nFirst = lList[l].node.begin();
	    nLast = lList[l].node.end();
	    nLast--;
	    double len = ((*nLast) - (*nFirst)).length();
	    nn = lList[l].node.begin();

	    //------------- Eva Feb 7
	    if (lList[l].node.size() == 3  && len < epsLine )
	    {
		//move inner point by averaging
		//cout << " move inner point by averaging" <<endl;
		// Vec2D nMid = getMiddlePoint(lList[l]);
		
 		nn++; // is middle node of this line  

		Vec2D nMid = (*nFirst) + ((*nLast) - (*nFirst))*0.5; 

 		(*nn) = nMid;
	    }
	    else
	    {
		//---------------------------------------
		list<Vec2D>::iterator nBeg, nNext, nPrev, nVel, nVel2;
		double aI = lList[l].alphaI; 
		double aJ = lList[l].alphaJ;
		nBeg = lList[l].node.begin();
		nBeg++;
		//nLast = lList[l].node.end();
		//nLast--;
		
		nVel = lList[l].nodeVel.begin();
		nVel++;
		
		for(n= nBeg; n != nLast; n++)
		{
		    double h, hPrev; 
		    
		    nPrev = n;
		    nPrev--;
		    nNext = n;
		    nNext++;
		    hPrev = ((*n)-(*nPrev)).length();
		    h     = ((*nNext)-(*n)).length();
		    
		    (*n) = (*n) + (*nVel) * (mobilityN * deltaTl) * (2./((h+hPrev)/2.) ); // main
		    
		    nVel++;
		}
		 }
	    
	}  
    cout << "total # lines " << linetotal << " fraction > 0.5 " << counter 
	 << "\t " << counter/linetotal << endl;
    
    
}   



//============================================================================
void evolveTripleJunctions(vector<Line> &lList, vector<Vertex> &vList, 
			   double &deltaTl, int g_left)
{
    // for each triple junction find the 3 lines that connect to it, 
    // evolution using velocity according to generalized Herring. Note the signs
    
    ofstream oFile1, oFile2;
    string fName_Nodes = "../Results/Nodes_Before ";  
    
    //oFile1.open(fName_Nodes.c_str(), ios::out);
    //printNodes(lList, fName_Nodes);
    //oFile1.close();
    
    

    double deltaTl_TJ = deltaTl;
    
    //cout << " deltaTl in evolve TJs " << deltaTl _TJ<< endl;

    vector<Vec2D> move;
    Vec2D zero  = Vec2D(0.0,0.0);
    for (int i=0; i<vList.size(); i++)
	move.push_back(zero);

    double res; 
    int nVert;
    //--- move TJs
    for (int iter=0; iter < triple_ite; iter++)
    {
	res = 0.;
	//cout << iter << endl;
	recalcVelo(vList, lList);  // recalculate velocity on first and last node
	getVeloVert(vList, lList); //--- get velocities on vertices
	nVert=0;
	for( int v=0; v < vList.size(); v++) 
	    if (vList[v].id > -1)
	    {   nVert++;
		vList[v].pos = vList[v].pos + vList[v].vel * mobilityTJ * deltaTl_TJ; // main formula
		
		move[v] = move[v] +  vList[v].vel * mobilityTJ * deltaTl_TJ; // for function adjust
		
		res += vList[v].vel.length();

		//--- update three 'ending' nodes of lines, that coincide with vertex vList[v] 
		list<Vec2D>::iterator first, last, n ;
		
		for (int j=0; j<vList[v].line.size(); j++)
		{
		    int l = vList[v].line[j];
		    int vv = vList[v].otherVertex(lList[l]);
		    
		    first = lList[l].node.begin();
		    last  = lList[l].node.end();
		    last--;
		    
		    if ( vv == lList[l].vertexJ)
			(*first) = vList[v].pos;
		    else if ( vv == lList[l].vertexI)
			(*last) = vList[v].pos;
		    else
		    {
			cout << " Bug here, cannot update node " << endl;
			exit(0);
		    }
		    
		    
		}
	    } // end of loop over vList  
	
        
    } // end of iter

    double mill = 1000000. ;
    cout << "res = " << res << " nVert= " << nVert <<  " " << (mill * res)/nVert   << endl; // if this is < 1 Hering fullfilled 
    //exit(0);

    int adjust = 1; // 1 = adjust nodes, 0 = don't adjust 
    if (adjust == 1)
    {
	for( int v=0; v < vList.size(); v++) 
	    if (vList[v].id > -1)
	    {
		for (int j=0; j<vList[v].line.size(); j++)
		{
		    int l = vList[v].line[j];
		    int vv = vList[v].otherVertex(lList[l]);

		    adjustNodes(lList, move[v] , l, vv);
		       
		}
	    }
    }
    //string fName_Nodes2 = "../Results/Nodes_After ";  
    
   // oFile2.open(fName_Nodes2.c_str(), ios::out);
    //printNodes(lList, fName_Nodes2);
    //oFile2.close();
    
}  

//===============================================================================

void recalcVelo(vector<Vertex> &vList, vector<Line> &lList)
{// calculates only the velocity of the first and the last node

    Vec2D zero = Vec2D(0.0,0.0);        
            
    for(int l=0; l<lList.size(); l++)
    { 
	if (lList[l].id > -1) 
	{
	    double aI = lList[l].alphaI; 
	    double aJ = lList[l].alphaJ;

	    list<Vec2D>::iterator n, n2, nVel; 
	    	    
	    //--- update first node velocity nVel
	    n= (lList[l].node).begin(); 
	    nVel = lList[l].nodeVel.begin();
	    //--------

	    (*nVel) = zero;
	    n2 = n;
	    n2++;
	    if (n2 == NULL)
		cout << " not defined " << endl;
	    Vec2D dx = (*n2) - (*n); 		
	    double h = dx.length() + 0.0000001;
	    
	    if( h > 0 ) 
	    {
		Vec2D bVec = dx * 0.2 / h;  
		
		//normalized vector between node1 and node2
		Vec2D nVec = bVec.perp();
		nVec = nVec * 0.2 / nVec.length(); // normal  
		
		double sig = sigma(aI, aJ, nVec);
		double Dsig = Dsigma(aI, aJ, nVec);
		
		Vec2D Tvec = nVec * Dsig + bVec * sig;  // main formula
		(*nVel)      = (*nVel)    + (Tvec*0.1);
		
	    }
	    else
	    {
		cout << "h < =0 " << endl;
		exit(0);
	    }
	      	    
	    (*nVel) = (*nVel) * 10.; // multiply by 10 to neutralize division by 0.1 above 
	    

	    //--- update last velocity nVel

	    n = lList[l].node.end(); 
	    n--;
	    n--;

	    nVel = lList[l].nodeVel.end();
	    nVel--; 
	    
	    //--------
	    (*nVel) = zero;
	    n2 = n;
	    n2++;

	    if (n2 == NULL)
		cout << " not defined " << endl;
	    
	    dx = (*n2) - (*n); 		
	    h = dx.length() + 0.0000001;
	    
	    if( h > 0 ) 
	    {
		Vec2D bVec = dx * 0.2 / h;  
		
		//normalized vector between node1 and node2
		Vec2D nVec = bVec.perp();
		nVec = nVec * 0.2 / nVec.length(); // normal  
		
		double sig = sigma(aI, aJ, nVec);
		double Dsig = Dsigma(aI, aJ, nVec);
		
		Vec2D Tvec = nVec * Dsig + bVec * sig;  // main formula
		
		(*nVel) = (*nVel) - (Tvec*0.1);
	    }
	    else
	    {
		cout << "h < =0 " << endl;
		exit(0);
	    }
	    	  
	    	    
	    (*nVel) = (*nVel)* 10.; //--- multiply by 10 to neutralize division by 0.1 above 
	}
    }  
}
//================================================================================
    void adjustNodes(vector<Line> &lList, Vec2D move, int l, int vv)
{

     //--- adjust the nodes, according to TJ-movement, to keep equally spacenodes
    //    if move  << h , dont do anything
    list<Vec2D>::iterator nBeg, nLast, n, nNext;
    nBeg = lList[l].node.begin();
    nBeg++;
    nLast = lList[l].node.end();
    nLast--;
    
    double hSum=0.;
    double hMax = 0.;
    double hMin = 1.;
    double h=0.;
    for(n= nBeg; n != nLast; n++)
    {
	nNext = n;
	nNext++;
	h = ((*nNext)-(*n)).length();

	if (h > hMax)
	    hMax = h;
	if (h < hMin)
	    hMin = h;

	hSum     = hSum + ((*nNext)-(*n)).length();
    }
    hSum = hSum/(lList[l].node.size()-1);
    
    
// get the modified move for each line
    Vec2D update;
    list<Vec2D>::iterator first, last, nn;
    first = lList[l].node.begin();
    last  = lList[l].node.end();
    last--;
    
    Vec2D bVec;
    if ( vv == lList[l].vertexI) // line pointing towards TJ
    {
	Vec2D x1 = (*last);
	nn = last;
	nn--;
	Vec2D x2 = (*nn);
	bVec = x1-x2; //tangent
    } 
    else if ( vv == lList[l].vertexJ) //line pointing outwards
    {
	Vec2D x1 = (*first);
	nn = first;
	nn++;
	Vec2D x2 = (*nn);
	bVec = x1-x2; //tangent
    }
    else
    {
	cout << " Bug here, cannot calc tangent in adjust " << endl;
	exit(0);
    }
    
    double fact = move.dotProduct(bVec);
    fact = fact/(bVec.length()*bVec.length());
    update = bVec * fact;
    
//==========
    
    
    //if (hMax/hMin > 2.)
    //{
	int N = lList[l].node.size(); 
	list<Vec2D>::iterator node;
	if ( vv == lList[l].vertexJ)
	{
	    node  = lList[l].node.begin();
	    node++;
	    
	    for (int i=1; i<N; i++)
	    {
		(*node) = (*node) + update*(1-(double)(i)/(double)(N-1));
		node++; 
	    }
	}
	else if ( vv == lList[l].vertexI)
	{
	    node  = lList[l].node.end();
	    node--;
	    node--;
	    for (int i=1; i<N; i++)
	    {
		(*node) = (*node) + update*(1-(double)(i)/(double)(N-1));
		node--;
	    }
	}
	else
	{
	    cout << " Bug here, cannot space nodes equally " << endl;
	    exit(0);
	}
	//}
    
    
}



//=========================================================================
void evolveInterior_Old(vector<Line> &lList)
{
    // evolve all interior points
    vector<Line>::iterator l;
    vector<Line>::iterator lp;
    list<Vec2D>::iterator n, np;

    // calculate nodes velocity
    for(l=lList.begin(); l != lList.end(); l++)
    {
	double aI = (*l).alphaI; 
	double aJ = (*l).alphaJ;
     
	for (n=((*l).nodeVel).begin(); n != ((*l).nodeVel).end(); n++)
	    (*n) = Vec2D(0.,0.);
    
	list<Vec2D>::iterator nVel = ((*l).nodeVel).begin();
	list<Vec2D>::iterator nVelNext;
      
	double hPrev = 0.0;
	double h;
	list<Vec2D>::iterator N = ((*l).normal).begin();
	list<Vec2D>::iterator prevN; 

	for(n= (((*l).node).begin()); n != --((*l).node).end(); n++)
	{
	    np = n;		 
	    np++;
     
	    nVelNext = nVel; nVelNext++;  
	    Vec2D dx = (*np) - (*n); 
	    dx = correctDirection(dx,2);

	    h = dx.length();
	    Vec2D bVec = dx / h;
	    Vec2D nVec = bVec.perp();
	    nVec = nVec / nVec.length();
	    (*N) = nVec; 
	  
	    if ( h == 0)
	    {
	      
		cout << "DIVISION BY ZERO 2 " << endl;
	      
		cout << "line " << (*l).id << " n " << (*n).x << " " << (*n).y << endl;
		cout << "line " << (*l).id << " np " << (*np).x << " " << (*np).y << endl;

	    }

	    double sig = sigma(aI, aJ, nVec);
	    double Dsig = Dsigma(aI, aJ, nVec);
	  
	    Vec2D Tvec = nVec * Dsig + bVec * sig;  // main formula
	  
	    (*nVel) = (*nVel) + Tvec;   
	    (*nVelNext) = (*nVelNext) - Tvec; 

	    //if ( n == ((*l).node).begin() )
	    //cout << "Node I Velocity in line " << (*l).id << " " <<
	    // (*nVel).x << " " << (*nVel).y << endl;

	    //if ( np == --((*l).node).end() )
	    //cout << "Node J Velocity in line " << (*l).id << " " << 
	    // (*nVelNext).x << " " << (*nVelNext).y << endl;

	    if( hPrev > 0.0 )
	    {
		prevN = N; prevN--;
		Vec2D nMid = ((*N)+(*prevN))/2.0;
		nMid = nMid / nMid.length();
		(*n) = (*n) + (*nVel) * (deltaT * mobility(aI, aJ) /  ((h + hPrev)/2.0) ) ;
		(*n) = correctDirection((*n),1);

	    }
	    nVel++;  N++;
	    hPrev = h;
	}
    }
}



//=========================================================================================== 
void printVertexEvo(vector<Vertex> &vList,vector<Line> &lList, vector<Grain> &gList, 
		    string fName )
{
    vector<Line> :: iterator l;
    
    //int firstV, lastV;
    int v1, v2;
    vector<int>::iterator v ;
    vector<Grain>::iterator g;
    ofstream oFile;
    oFile.open( fName.c_str(), ios::out);
    
    if( oFile == NULL ) 
    {
	cout << " could not open output file " << fName << endl;
	exit(0);
    }
    
//    int gg=557;
//    for (int l=0; l<  (gList[gg].line).size(); l++)
    for( int l=0; l < lList.size(); l++) 
    {
        //if (lList[gList[gg].line[l]].id > -1)
	if (lList[l].id > -1)
	{
	    v1 = lList[l].vertexI;
	    v2 = lList[l].vertexJ;
	    double domainSize = 1.;
	    double tmpx = vList[v1].pos.x -vList[v2].pos.x;
	    double tmpy = vList[v1].pos.y -vList[v2].pos.y;
	    double lineLength = tmpx*tmpx + tmpy*tmpy;
	    lineLength = sqrt(lineLength);
	    // if( lineLength < 1.2 * domainSize  ) 
// 		//    if( lList[l].id > -1 ) 
// 	    {
// 		firstV= lList[l].vertexI;
// 		lastV = lList[l].vertexJ;
// 		if (vList[firstV].id > -1)
// 		    oFile << vList[firstV].pos.x << '\t' <<  vList[firstV].pos.y << endl;
// 		if (vList[lastV].id > -1)
// 		    oFile << vList[lastV].pos.x << '\t' << vList[lastV].pos.y  << endl << endl;
// 	    }
	    
	    if (vList[v1].id > -1  && vList[v2].id > -1)
	    {
		if (abs(tmpx) > 1 && abs(tmpy) <= 1 )
		{// take x mirror
		    // cout << " print Evo using x mirrors" << endl;
// 		    cout << "#v1 " << " " << vList[v1].id << " v2 " << vList[v2].id << endl;
		    double new_v1_x;
		    if (vList[v1].pos.x < 0)
			new_v1_x = vList[v1].pos.x + 2.;
		    else
			new_v1_x = vList[v1].pos.x - 2.;
		    oFile << new_v1_x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl << endl; 
		    
		    double new_v2_x;
		    if (vList[v2].pos.x < 0)
			new_v2_x = vList[v2].pos.x + 2.;
		    else
			new_v2_x = vList[v2].pos.x - 2.;
		    
		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << new_v2_x << '\t' << vList[v2].pos.y << endl << endl;
		}
		if (abs(tmpx) <= 1 && abs(tmpy) > 1 )
		{//take y mirror
		    // cout << "#v1 " << " " << vList[v1].id << " v2 " << vList[v2].id << endl;
// 		    cout << "# print Evo using y mirrors" << endl;
		    double new_v1_y;
		    if (vList[v1].pos.y < 0)
			new_v1_y = vList[v1].pos.y + 2.;
		    else
			new_v1_y = vList[v1].pos.y - 2.;

		    oFile << vList[v1].pos.x << '\t' << new_v1_y << endl;
		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl << endl;
		    		    
		    double new_v2_y;
		    if (vList[v2].pos.y < 0)
			new_v2_y = vList[v2].pos.y +2.;
		    else
			new_v2_y = vList[v2].pos.y -2.;
		    
		    
		    oFile << vList[v2].pos.x << '\t' << new_v2_y << endl;
		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl<< endl;
		    		}
		if (abs(tmpx) > 1 && abs(tmpy) > 1 )
		{// take diagonal mirror
		    // cout << "#v1 " << " " << vList[v1].id << " v2 " << vList[v2].id << endl;
// 		    cout << "# print Evo using x and y mirrors" << endl;
		    double new_v1_x, new_v2_x;
		    double new_v1_y, new_v2_y;
		    if (vList[v1].pos.x < 0)
			new_v1_x = vList[v1].pos.x + 2.;
		    else
			new_v1_x = vList[v1].pos.x - 2.;

		    if (vList[v1].pos.y < 0)
			new_v1_y = vList[v1].pos.y + 2.;
		    else
			new_v1_y = vList[v1].pos.y - 2.;

		    if (vList[v2].pos.x < 0)
			new_v2_x = vList[v2].pos.x + 2.;
		    else
			new_v2_x = vList[v2].pos.x - 2.;
		    if (vList[v2].pos.y < 0)
			new_v2_y = vList[v2].pos.y +2.;
		    else
			new_v2_y = vList[v2].pos.y -2.;

		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << new_v2_x << '\t' << new_v2_y << endl << endl;

		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl;
		    oFile << new_v1_x << '\t' << new_v1_y << endl << endl;
		}
		if (abs(tmpx) <= 1 && abs(tmpy) <= 1)
		{
		    //cout << "# normal line " << endl;
		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl << endl;
		}
	    }
	    
	}
    }
}	    
//=========================================================================================== 
void printVertexLabelEvo(vector<Vertex> &vList,vector<Line> &lList, string fName )
{
    vector<Line> :: iterator l;
    
    int firstV, lastV;
    vector<int>::iterator v ;
    vector<Grain>::iterator g;
    ofstream oFile;
    oFile.open( fName.c_str(), ios::out);
    
    if( oFile == NULL ) 
    {
	cout << " could not open output file " << fName << endl;
	exit(0);
    }
	
    oFile.precision(16);
    cout << " print to file ../GNUPLOT/evoVertexLabel.dat " << endl;
    
    for ( int vv=0; vv < vList.size(); vv++ )
    {
	if (vList[vv].id != -1) 
	    oFile << " set label  ' " << vv  << " '\t  at \t" << vList[vv].pos.x << "\t," << vList[vv].pos.y << endl;
    }
    oFile.close();
    
}

//=====================print for GNUPLOT grains and lines ================ 
void printGrainsGNU(vector<Grain> &gList,  vector<Vertex> &vList ) 
{

    ofstream oFile;
        
    vector<Grain> :: iterator g;
    vector<int> :: iterator v, firstV;
    //vector<Vertex> :: iterator vv;
    
    oFile.open("../GNUPLOT/evolvedGrains.dat", ios::out);
    oFile.precision(8 );
    cout << " print to file ../GNUPLOT/evolvedGrains.dat " << endl;
    
    for( g=gList.begin(); g!=gList.end(); g++) 
    {
	int gId = (*g).id;
	firstV= ((*g).vertex).begin();
	
        for( v= ((*g).vertex).begin(); v!= ((*g).vertex).end(); v++)
	{
	    if (vList[*v].id != -1 ) 
		oFile << vList[(*v)].pos.x << '\t' <<  vList[(*v)].pos.y << endl;
	}
	if (vList[*firstV].id != -1) 
	    oFile << vList[(*firstV)].pos.x << '\t' << vList[(*firstV)].pos.y  << endl;
	oFile << endl;
    }
    
    oFile.close();
    oFile.open("../GNUPLOT/evoVertexLabel.dat", ios::out);
    oFile.precision(16);
    cout << " print to file ../GNUPLOT/evoVertexLabel.dat " << endl;
    
    for ( int vv=0; vv < vList.size(); vv++ )
    {
	if (vList[vv].id != -1) 
	//cout << " vv = " << vv << endl ;
	    oFile << " set label  ' " << vv  << " '\t  at \t" << vList[vv].pos.x << "\t," << vList[vv].pos.y 	      << endl;
    }
    oFile.close();
}
//=========================================================================================== 
void printLinesGNU(vector<Line> &lList,  vector<Vertex> &vList ) 
{

    ofstream oFile;
        
    vector<Line> :: iterator l;
    vector<int> :: iterator v; 
    int firstV, lastV;
    int v1, v2;

    //vector<Vertex> :: iterator vv;
    
    oFile.open("../GNUPLOT/evolvedLines.dat", ios::out);
    oFile.precision(8 );
    cout << " print to file ../GNUPLOT/evolvedLines.dat " << endl;
    double domainSize = 1.;
    double lineLength;
    
    for( int l=0; l < lList.size(); l++) 
    {
	if (lList[l].id > -1)
	{
	    double tmpx = vList[lList[l].vertexI].pos.x -vList[lList[l].vertexJ].pos.x;
	    double tmpy = vList[lList[l].vertexI].pos.y -vList[lList[l].vertexJ].pos.y;
	    lineLength = tmpx*tmpx + tmpy*tmpy;
	    lineLength = sqrt(lineLength);
	
// 	    //if( lineLength < 1.5 * domainSize  ) 
// 	    //if( lList[l].id > -1  ) 
// 	    {
// 		firstV= lList[l].vertexI;
// 		lastV = lList[l].vertexJ;
// 		if (vList[firstV].id > -1)
// 		    oFile << vList[firstV].pos.x << '\t' <<  vList[firstV].pos.y << endl;
// 		if (vList[lastV].id > -1)
// 		    oFile << vList[lastV].pos.x << '\t' << vList[lastV].pos.y  << endl << endl;
// 	    }
//	}
	    //  }
    
//============
	    v1 = lList[l].vertexI;
	    v2 = lList[l].vertexJ;
	    if (vList[v1].id > -1  && vList[v2].id > -1)
	    {
		if (abs(tmpx) > 1 && abs(tmpy) <= 1 )
		{// take x mirror
		    
		    double new_v1_x;
		    if (vList[v1].pos.x < 0)
			new_v1_x = vList[v1].pos.x + 2.;
		    else
			new_v1_x = vList[v1].pos.x - 2.;
		    oFile << new_v1_x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl << endl; 
		    
		    double new_v2_x;
		    if (vList[v2].pos.x < 0)
			new_v2_x = vList[v2].pos.x + 2.;
		    else
			new_v2_x = vList[v2].pos.x - 2.;
		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << new_v2_x << '\t' << vList[v2].pos.y << endl << endl;
		    
		}
		if (abs(tmpx) <= 1 && abs(tmpy) > 1 )
		{//take y mirror
		   
		    double new_v1_y;
		    if (vList[v1].pos.y < 0)
			new_v1_y = vList[v1].pos.y + 2.;
		    else
			new_v1_y = vList[v1].pos.y - 2.;

		    oFile << vList[v1].pos.x << '\t' << new_v1_y << endl;
		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl << endl;
		    
		    double new_v2_y;
		    if (vList[v2].pos.y < 0)
			new_v2_y = vList[v2].pos.y +2.;
		    else
			new_v2_y = vList[v2].pos.y -2.;
		    
		    
		    oFile << vList[v2].pos.x << '\t' << new_v2_y << endl;
		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl<< endl;
		    
		}
		if (abs(tmpx) > 1 && abs(tmpy) > 1 )
		{// take diagonal mirror
		    
		    double new_v1_x, new_v2_x;
		    double new_v1_y, new_v2_y;
		    if (vList[v1].pos.x < 0)
			new_v1_x = vList[v1].pos.x + 2.;
		    else
			new_v1_x = vList[v1].pos.x - 2.;

		    if (vList[v1].pos.y < 0)
			new_v1_y = vList[v1].pos.y + 2.;
		    else
			new_v1_y = vList[v1].pos.y - 2.;

		    if (vList[v2].pos.x < 0)
			new_v2_x = vList[v2].pos.x + 2.;
		    else
			new_v2_x = vList[v2].pos.x - 2.;
		    if (vList[v2].pos.y < 0)
			new_v2_y = vList[v2].pos.y +2.;
		    else
			new_v2_y = vList[v2].pos.y -2.;

		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << new_v2_x << '\t' << new_v2_y << endl << endl;
		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl;
		    oFile << new_v1_x << '\t' << new_v1_y << endl << endl;

		}
		if (abs(tmpx) <= 1 && abs(tmpy) <= 1)
		{
		    
		    oFile << vList[v1].pos.x << '\t' <<  vList[v1].pos.y << endl;
		    oFile << vList[v2].pos.x << '\t' << vList[v2].pos.y  << endl << endl;
		}
	    }
	}
    }
//==============

    oFile.close();
    oFile.open("../GNUPLOT/evoLinesLabel.dat", ios::out);
    oFile.precision(16);
    cout << " print to file ../GNUPLOT/evoLinesLabel.dat " << endl;
    
    for ( int vv=0; vv < vList.size(); vv++ )
    	if (vList[vv].id > -1 ) 
	    oFile << " set label  ' " <<vv<< " '\t  at \t" <<vList[vv].pos.x<<"\t,"
		  << vList[vv].pos.y 	      << endl;
    oFile.close();
}
//====================================================================================
void updateTimeStep(vector<Vertex> &vList, vector<Line> &lList)
{
    // Update time step to order hmin^2
    // Done by Guillermo 11/15/05

    double minLen = 10000000.0;
    double h;

    for(int l=0; l<lList.size(); l++)
    {
	if ( lList[l].id != -1)
	{
	    int vI = lList[l].vertexI;
	    int vJ = lList[l].vertexJ;
	    Vec2D dir = vList[vJ].pos - vList[vI].pos;
	    

	    double vecLen = dir.length();
	    if( minLen > vecLen ) minLen = vecLen;
	}
    }
   
    h = minLen;
    deltaT = min( (h*h) * deltaT_const , mindeltaT);
}
//====================================================================================
double updatehMin(vector<Vertex> &vList, vector<Line> &lList, double &factor)
{
    double aveLen = 0.0;
    double minLen = 1.;
    double h;
    double counter=0.0;
    list<Vec2D>::iterator n, nLast,nNext;

    // calculate average line length
    for(int l=0; l<lList.size(); l++)
    	if ( lList[l].id != -1)
	{
	    double len=0.0;
	    nLast = lList[l].node.end();
	    nLast--;
	    
	    for (n=lList[l].node.begin(); n!= nLast; n++ ) 
	    {
		nNext=n;
		nNext++;
		Vec2D dir = (*nNext)-(*n);
		len += dir.length();
	    }
	    aveLen += len;
	    if (len < minLen)
		minLen = len;
	    counter++;
	}
    aveLen=aveLen/counter; // average Linelength at time t
    hMin = aveLen * factor ; // factor = hMinIni/aveLenIni (at time 0)
    //hMin = minLen * factor ;
    //cout << " new hMin " << hMin << endl; 
    return hMin; 
}




//===============================================================================================
  int main(int argc, char **argv)
{
    if (!DEBUG) srand(static_cast<unsigned>(time(0)));
    string inputFile;
    int extime = 5;

    current_time = 0;

    if( argc > 1 ) 
	inputFile.append(argv[1]);
    else
    {
	cout << "Usage: grainGrowth fileName [time] " << endl;
	exit(1);
    }

    if( argc > 2 ) 
	extime = atoi(argv[2]);

    // main data structure
    vector<Vertex> vList;
    vector<Grain> gList;
    vector<Line> lList;
    vector<Pair> IOgrCenter;
    vector<Pair> IOqv;
    
    // Output Files
    string PATH = "../Results/";
    string confFileName;
    string dot = ".";
    string fName; 
    string zero1 = "0";
    string zero2 = "00";
    string zero3 = "000";
    string zero4 = "0000";
    string zero5 = "00000";
    string jpeg = ".jpeg";
    string fName_jpeg;
    string fName_sides;
    string fName_Neigh;
    string fName_HistSides;
    string fName_RelArea;
    string fName_HistRelArea;
    string fName_SingArea;
    string fName_Length;
    string fName2;
    string fName_Vert;
    string fName_Label;
    string fName_Nodes;
    
    
    ofstream oFile, oFile2;

    readParameters("parameters.graingrowth");
    readData(inputFile, vList, gList,IOgrCenter,IOqv);

    checkIniConf(gList, vList);
    

    defineDataStructures(vList, gList, lList, IOqv); 
    periodicStructur(vList, gList, lList, IOqv, IOgrCenter);

    checkConsist(vList, gList, lList );

    cout << " triple_ite " << triple_ite << endl;

    if( PRINT ) 
    {
	printLines(lList);
	printGrains(gList);
	printVertices(vList);
	confFileName = "initintitnf"; 
	fName = PATH + confFileName;
	printVertexEvo(vList,lList, gList, fName);
    }

    double small_l = 100000.0;
    double big_l = -1.0;
    int lines = 0;

    double aveLen0=0.0;
    for (int l =0; l < lList.size(); l++)
    {
	if (lList[l].id > -1)
	{
	    Vec2D actual = vList[lList[l].vertexI].pos - vList[lList[l].vertexJ].pos;
	    double len = actual.length();
	    aveLen0 += len;
	    if ( len < small_l )
		small_l = len;
	    if ( len > big_l )
		big_l = len;
	    
	    lines++;
	}
    }
    aveLen0 = aveLen0/lines;
   
    cout << " big_l " << big_l << " small_l " << small_l << endl;
    cout << "ratio lines biggest/smallest " << big_l/small_l << " number of lines " 
	 << lines << endl;
    
    
    cout << " epsLine = " << epsLine << endl;

    //updateTimestep(vList, lList);

    
    cout << " in main hMin" << hMin <<  endl;

    cout << "assigning nodes" << endl;
    initializeNodes(vList, lList);
   
    double factor=hMin/aveLen0;  // factor  = hMinIni/aveLenIni (at time 0)
 
    cout << "factor " << factor << " hMin=" << hMin << " aveLen0 " << aveLen0 << endl;
    
    //---
    int nodes = 0;
    int maxnodeslines = -1;
    
    for (int l = 0; l < lList.size(); l++ )
    	if (lList[l].id != -1 )
	{
	    int nodel = (lList[l].node).size();
	    nodes += nodel;
	    if ( nodel > maxnodeslines)
		maxnodeslines = nodel;
	}
        
    cout << "after initializing nodes, number of nodes: " << nodes 
	 << " Max number of nodes in lines: " << maxnodeslines << endl;
    //---

    fName_Nodes = PATH + "Nodes";
    string fName_Nodes1 = PATH + "Nodes1";

    time_t start, end,step;
    time (&start);
    int g_left, g_prev;

    int nGrains = 0;
    for (int k=0; k<gList.size(); k++)
	if (gList[k].id > -1)
	    nGrains++;
    
    cout << "Initialnumber of grains " << nGrains << endl << endl;
    g_left = nGrains;
    
    double Area, aveArea;
    vector<double> writeAveArea;
    vector<double> timeVec;
    double totalArea;
    vector<double> SingArea;

    vector<double> NGrainsLeft;
    double deltaTl = mindeltaT; // first deltaTl
   
    

//----------------------------- begin of main loop --------------------------	
    for(int t=0; t<extime+1; t++) 
    {
	g_prev=g_left;

	//--- to print average area:------ 
	double aveArea_0;
	aveArea = Area/g_left;
	if (t==0 )
	   aveArea_0 = aveArea; 
	
//	printNodes(lList, fName_Nodes1);
//	printLinesGNU(lList, vList);
	
	aveArea=aveArea/aveArea_0;
	
	writeAveArea.push_back(aveArea);
	timeVec.push_back(current_time);
	//---------------------------------

	//if ( t % extime == 0)
	//{
	    cout << "======================================================" << endl;
	    cout << "Executing Time  " << t << " of " << extime << endl;
	    g_prev=g_left;
	    g_left=0;
	    for (int m=0; m<gList.size();m++)
		if (gList[m].id > -1)
		    g_left++;

	    cout << endl << "Grains left " << g_left << " Deleted in this interval " 
		         << g_prev-g_left << endl;
	    time(&step);

	    cout << "Accumulated Time  " <<  difftime(step,start)/60.0 << " min" << endl;
	    cout << " gleft , nGrains " << g_left << " " << nGrains <<" " << 0.25*nGrains << endl;
	    
	    
	    if (g_left < 0.25*nGrains) //  || g_left < 1000)
	    {
		cout << "========================================================" << endl;
		cout << " DONE, more than 75% grains remov " << endl;

		cout << "Writing File at time " << t << endl;
		// To make the files in order
		if ( t <= 9)      fName = confFileName + zero5 + toString(t);
		else if ( t <= 99)  fName = confFileName + zero4 + toString(t);
		else if ( t <= 999)  fName = confFileName + zero3 + toString(t);
		else if ( t <= 9999)  fName = confFileName + zero2 + toString(t);
		else if ( t <= 99999)  fName = confFileName + zero1 + toString(t);
		else if ( t <= 999999)  fName = confFileName + toString(t);

		//fName_sides = PATH + "sides_" + fName;	  
		fName_sides = "../Results/Sides/" + toString(t);		
		fName_HistSides = PATH + "HistSides_" + fName;
		fName_HistRelArea = PATH + "HistRelArea_" + fName;
		fName_RelArea = PATH + "RelArea_" + fName;
		//fName_SingArea = PATH + "SingArea_" + fName;
		fName_SingArea = "../Results/Area/" + toString(t);
		//fName_Length = PATH + "Length_" + fName;
		fName_Length = "../Results/Length/"  + toString(t);
		fName_Neigh = "../Results/Neighbor/" + toString(t);

		fName = PATH + fName;
		
		//if (DEBUG == 3)
		//{
		//---  Get statistics
		
		writeRelativeArea_Nodes(gList, lList, vList, fName_HistRelArea, fName_RelArea, 
				    fName_SingArea, totalArea, SingArea);	
		//writeLineLength(lList, vList, fName_Length);
		neighbors(gList, lList, fName_Neigh);

		// --- print Average Area --- 
		fName2 = PATH + "AveArea_" + toString(nGrains);
		oFile.open(fName2.c_str(),ios::out);
		for (int k=0; k<writeAveArea.size(); k++)
		    oFile << timeVec[k] << " " << writeAveArea[k] << endl;
 		oFile.close();
		//}
		
                // --- print number of grains left
		string fName3 = PATH + "NGrains_" + toString(nGrains);
		oFile2.open(fName3.c_str(),ios::out);
		for (int k=0; k<NGrainsLeft.size(); k++)
		    oFile2 << timeVec[k] << "  " << NGrainsLeft[k]<< " \t  "<<nGrains - NGrainsLeft[k] <<  endl;
		oFile2.close();

		cout << "Done writing statistics" << endl;
		cout << "exiting ..." << endl;
		exit(0);
	    }
	//}

	
	hMin = updatehMin(vList, lList, factor); // --- update hMin 
	epsLine = hMin  * 1.5 ;                    //--- update epsLine 

	

//-------------------  Critical Events --------------------------------------
	cout << "====critical Events ================= \n";
	
	//cout << " removeShorts " << endl;
	removeShortLines(vList, gList, lList, IOqv, 1, deltaTl ); 
	
	//cout << " splitUnst " << endl;
	splitUnstableVertices(vList, gList, lList, IOqv);
	
//-------------------- Evolution ---------------------------------------------
	cout << "====Evolution ================= \n";

	//--- get number of grains 
	g_left=0;
	for (int m=0; m<gList.size();m++)
	    if (gList[m].id > -1)
		g_left++;
	//---

	deltaTl = hMin*hMin * 0.1 ; // * 0.25;   // update deltaTl
	
	evolveInterior(lList,vList, deltaTl, extime);
	evolveTripleJunctions(lList, vList, deltaTl, g_left);

	redistributeNodes(vList, lList);

//----------------------------------------------------------------------------
	if (PRINT) 
	{
	    printLines(lList);
	    printGrains(gList);
	    printVertices(vList);
	    fName_Vert = PATH + "Vert0_" + zero1;
 	    printVertexEvo(vList,lList, gList, fName_Vert);

	}
	
        // here create function void Printing(WRITE, t , print_times) 

	if ( WRITE && t % print_times == 0 )
	{
	    cout << "Writing File at time " << t << endl;
	    // To make the files in order
	    if ( t <= 9)      fName = confFileName + zero5 + toString(t);
	    else if ( t <= 99)  fName = confFileName + zero4 + toString(t);
	    else if ( t <= 999)  fName = confFileName + zero3 + toString(t);
	    else if ( t <= 9999)  fName = confFileName + zero2 + toString(t);
	    else if ( t <= 99999)  fName = confFileName + zero1 + toString(t);
	    else if ( t <= 999999)  fName = confFileName + toString(t);

	    
	    fName_sides = "../Results/Sides/" + toString(t);
	    fName_HistSides = PATH + "HistSides_" + fName;
	    fName_HistRelArea = PATH + "HistRelArea_" + fName;
	    fName_RelArea = PATH + "RelArea_" + fName;	
	    fName_SingArea = "../Results/Area/" + toString(t);
	    fName_Vert = "../Results/Evolution/" + toString(t);
	    fName_Label = PATH + "Label_" + fName;
	    fName_Length = "../Results/Length/" + toString(t);
	    fName_Neigh = "../Results/Neighbor/" + toString(t);
   	    fName = PATH + fName;

	    //printVertexEvo(vList,lList, gList, fName_Vert);
	    printNodes(lList,fName_Vert);

            //if (DEBUG == 3)
            //---  Get statistics
	    //{		
	    
	    writeRelativeArea_Nodes(gList, lList, vList, fName_HistRelArea, fName_RelArea, 
	    			    fName_SingArea, totalArea, SingArea);
	    //writeLineLength(lList, vList, fName_Length);
	    neighbors(gList, lList, fName_Neigh);

	    cout << "Done writing statistics" << endl;
	    //}	
	}

	// g_left=0;
// 	for (int m=0; m<gList.size();m++)
// 	    if (gList[m].id > -1)
// 		g_left++;
	
	cout << endl << "Grains left " << g_left << " Deleted in this interval " 
	     << g_prev-g_left << endl;
	NGrainsLeft.push_back(g_left);
	
    }
    
    time (&end);

    cout << "Total Time for  " << extime << " : "                            
	 <<  difftime(end,start)/60.0 << " min" << endl;
        
// --- print Average Area --- 
    fName2 = PATH + "AveArea_" + toString(nGrains);
    oFile.open(fName2.c_str(),ios::out);
    for (int k=0; k<writeAveArea.size(); k++)
	oFile << timeVec[k] << " " << writeAveArea[k] << endl;
    oFile.close();
    
// --- print number of grains left
    string fName3 = PATH + "NGrains_" + toString(nGrains);
    oFile2.open(fName3.c_str(),ios::out);
    for (int k=0; k<NGrainsLeft.size(); k++)
	oFile2 << timeVec[k] << "\t" << NGrainsLeft[k] << " \t" << nGrains - NGrainsLeft[k] << endl;
    oFile2.close();



    return 0;
}