/*
 $Id: intA.cc,v 1.3 1996/10/11 15:15:29 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "intA.h"

#include "utils.h"
#include "familyA.h"
#include "nodeco.h"
#include "dirichlet.h"
#include "triang.h"
#include "triangA.h"

//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// --    store node Numbers to vector node for change local <-> global:


void LinElemInt:: getGlobalNodes(const PATCH* t, Vector<int>& node) const
{
    int i, n, nGlob, nLoc = 1;
    const Vector<PT*>& pt = t->getPoints();

    FORALL(pt,i)
    {
	nGlob = pt[i]->getNode();
	if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
	else	   FORNCOMP(n) node[nLoc++] = 0; 
    }  

    //node[1] = e->p1->getNode();
    //node[2] = e->p2->getNode(); ...
}
//-------------------------------------------------------------------------


void LinElemInt:: getGlobalMLNodes(const PATCH* t, Vector<int>& node, int depth)
									   const
{
    int i, n, nGlob, nLoc = 1;

    const Vector<PT*>& pt = t->getPoints();

    FORALL(pt,i)
    {
	nGlob = pt[i]->MLNode(depth);
	if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++; 
	else	   FORNCOMP(n) node[nLoc++] = 0; 
    }	

    //node[1] = e->p1->MLNode(depth);
    //node[2] = e->p2->MLNode(depth); ...
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LQuadElemInt:: getGlobalNodes(const PATCH* t, Vector<int>& node) const
{
    int i, n, nGlob, nLoc = 1;
    const Vector<PT*>& pt = t->getPoints();

    FORALL(pt,i)
    {
	nGlob = pt[i]->getNode();
	if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++; 
	else	   FORNCOMP(n) node[nLoc++] = 0; 
    }	

    const Vector<EDG*>& edg = t->getEdges();

    FORALL(edg,i)
    {
	nGlob = edg[i]->getNode();
	if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++; 
	else	   FORNCOMP(n) node[nLoc++] = 0; 
    }
}
//-------------------------------------------------------------------------


void LQuadElemInt:: getGlobalMLNodes(const PATCH* t, Vector<int>& node,
				     int depth) const
{
    int i, n, nGlob, nLoc = 1;
    const Vector<PT*>& pt = t->getPoints();

    FORALL(pt,i)
    {
	nGlob = pt[i]->MLNode(depth);
	FORNCOMP(n) node[nLoc++] = nGlob++;
    }	

    const Vector<EDG*>& edg = t->getEdges();

    FORALL(edg,i)
    {
	nGlob = edg[i]->getNode();
	FORNCOMP(n) node[nLoc++] = nGlob++;
    }

    // node[1] = e->p1->MLNode(depth);
    // node[2] = e->p2->MLNode(depth);
    // node[3] = e->getNode(); ...
}
//-------------------------------------------------------------------------


void EFTetraInt:: getGlobalNodes(const PATCH* t, Vector<int>& node) const
{
    int i, n, nGlob, nLoc = 1;
    const Vector<PT*>& pt = t->getPoints();

    FORALL(pt,i)
    {
	nGlob = pt[i]->getNode();
	if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++; 
	else	   FORNCOMP(n) node[nLoc++] = 0; 
    }	

    const Vector<EDG*>& edg = t->getEdges();

    FORALL(edg,i)
    {
	nGlob = edg[i]->getNode();
	if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++; 
	else	   FORNCOMP(n) node[nLoc++] = 0; 
    }

    Vector<int> faceNodes(4);
    t->getFaceNodes(faceNodes);

    FORALL(faceNodes,i)
    {
	nGlob = faceNodes[i];
	if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++; 
	else	   FORNCOMP(n) node[nLoc++] = 0; 
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LinElemInt:: setNodeNumbers()
{
    int    node;
    PATCH* p;
    
    mesh->resetPtIter();
    while ((p=mesh->ptIterAll())) p->unMark();

    mesh->resetPtIter();

    if (noOfNodes.high() < 0)
    {
	node = 1;
	while ((p=mesh->ptIterAll())) 
	{ 
	    p->setNode(node);  
	    node += nComp; 
	}
    }
    else
    {
	node = noOfNodes.Top() + 1;

	while ((p=mesh->ptIterAll())) 
	  if (p->getNode() == 0) 
	  {
	      p->setNode(node);	
	      p->setMark();			// new nodes are marked
	      node += nComp;
	  }
    }

    node -= 1;
    noOfNodes.push(node);
}
//-------------------------------------------------------------------------


void LQuadElemInt:: setNodeNumbers()
{
    int    node;
    PATCH* p;
    
    mesh->resetPtIter();
    while ((p=mesh->ptIterAll())) p->unMark();

    mesh->resetPtIter();
    mesh->resetEdgIter();

    if (noOfNodes.high() < 0)
    {
	node = 1;
	while ((p=mesh->ptIterAll()))  
	{
	    p->setNode(node);
	    p->setMark();
	    node += nComp;
	}
	node -= 1;
	pointNodes.push(node);
    }
    else
    {
	node = pointNodes.Top() + 1;

	while ((p=mesh->ptIterAll())) 
	  if (p->getNode() == 0) 		// new node
	  {
	      p->setNode(node);
	      p->setMark();
	      node += nComp;
	  }
	
	node -= 1;
	pointNodes.push(node);
    }
    
    node += 1;
    while ((p=mesh->edgIterAll())) 
    {
	p->setNode(node);
	node += nComp;
    }

    node -= 1;
    noOfNodes.push(node);
}
//-------------------------------------------------------------------------


void HQuadElemInt:: setNodeNumbers()
{
    int    node;
    PATCH* p;
    
    mesh->resetPtIter();
    while ((p=mesh->ptIterAll())) p->unMark();

    mesh->resetPtIter();
    mesh->resetEdgIter();

    if (noOfNodes.high() < 0)
    {
	node = 1;
	while ((p=mesh->ptIterAll()))  
	{ 
	    p->setNode(node);
	    p->setMark();
	    node += nComp;
	}
	node -= 1;
	pointNodes.push(node);
    }
    else
    {
	node = pointNodes.Top() + 1;

	while ((p=mesh->ptIterAll())) 
	  if (p->getNode() == 0) 		// new node
	  {
	      p->setNode(node);
	      p->setMark();
	      node += nComp;
	  }
	node -= 1;
	pointNodes.push(node);
     }

    node += 1;

    if (dirichletBCs->constant())
    {
	while ((p=mesh->edgIterAll())) 
	  if (p->isDirichlet()) p->setNode(0);
	  else 		
	  {
	      p->setNode(node);
	      node += nComp;
	  }
    }
    else 
    {
	while ((p=mesh->edgIterAll())) 
	{
	    p->setNode(node);
	    node += nComp;
	}
    }

    node -= 1;
    noOfNodes.push(node);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LinElemInt:: interpolateSolution(Vector<Num>& u) const
{
    PATCH* ed, *pm;
    int    n, node1, node2, pmNode;

    mesh->resetEdgIter();
    while ((ed=mesh->edgIterAllOfHistory())) 
    {
	if (ed->refined())	
	{
	    pm = ed->getMidPoint();

	    if (pm->marked())		// new node here
	    {
    		const Vector<PT*>& p = ed->getPoints();

		pmNode = pm->getNode();
		node1 = p[1]->getNode();
		node2 = p[2]->getNode();

		FORNCOMP(n) u[pmNode++] = 0.5*(u[node1++] + u[node2++]);
	    }
	}
    }
}
//-------------------------------------------------------------------------


void LQuadElemInt:: interpolateSolution(Vector<Num>& u) const
{
    int i, n;
    PATCH* ed, *pm;
    int edNode, node1, node2, node;

    Vector<Num> uOld(u.high());
    FORALL(u,i) uOld[i] = u[i];

    mesh->resetEdgIter();
    while ((ed=mesh->edgIterAllOfHistory())) 	     // set the new point nodes
    {
	if (ed->refined())	
	{
	    pm = ed->getMidPoint();
	    if (pm->marked()) 
	    {
		node = pm->getNode();
		edNode = ed->getNode();
		FORNCOMP(n)  u[node++] = uOld[edNode++];
	    }
	}
    }

    mesh->resetEdgIter();
    while ((ed=mesh->edgIterAll()))  			// the edge nodes
    {
	const Vector<PT*>& p = ed->getPoints();

	edNode = ed->getNode();
	node1 = p[1]->getNode();
	node2 = p[2]->getNode();

	FORNCOMP(n) u[edNode++] = 0.5*(u[node1++] + u[node2++]);
    }
}
//-------------------------------------------------------------------------


void HQuadElemInt:: interpolateSolution(Vector<Num>& u) const
{
    PATCH* ed, *pm;
    int    i, n, node, node1, node2, edNode;

    Vector<Num> uOld(u.high());
    FORALL(u,i) uOld[i] = u[i];


    mesh->resetEdgIter();
    while ((ed=mesh->edgIterAllOfHistory())) 		// the new point nodes
    {
	if (ed->refined())	
	{
	    pm = ed->getMidPoint();

	    if (pm->marked())			// new point node here
	    {
		edNode = ed->getNode();
		node   = pm->getNode();

		if (edNode) FORNCOMP(n) u[node++] = uOld[edNode++];
		else	    FORNCOMP(n) u[node++] = 0.0;
  
		const Vector<PT*>& p = ed->getPoints();

		node  = pm->getNode();
		node1 = p[1]->getNode();
		node2 = p[2]->getNode();

		FORNCOMP(n) u[node++] += 0.5 * (uOld[node1++] + uOld[node2++]);
	    }
	}
    }

    mesh->resetEdgIter();
    while ((ed=mesh->edgIterAll()))   			// the edge nodes
    {
	edNode=ed->getNode();
	if (edNode)  FORNCOMP(n) u[edNode++] = 0.0;
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LinElemInt:: updateDirichletBCs(Real time) const
{
    int    n, node;
    PATCH* p;
    Vector<Real> x(spaceDim);
    const Bool constDBCs = dirichletBCs->constant();
    
    dirichletBCs->extend(noOfNodes.Top());

    mesh->resetPtIter();
    while ((p=mesh->ptIterAll())) 
    {
	FORNCOMP(n)
	{
	    if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
	    {
		if (!constDBCs) p->getNodeCoord(x);
		node = p->getNode() + n-1;
		dirichletBCs->setBC(node, p->Class(), x, n, time); 
	    }	
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LQuadElemInt:: updateDirichletBCs(Real time) const
{
    int    n, node;
    PATCH* p;
    Vector<Real> x(spaceDim);
    const Bool constDBCs = dirichletBCs->constant();
    
    dirichletBCs->extend(noOfNodes.Top());

    mesh->resetPtIter();
    while ((p=mesh->ptIterAll())) 
    {
	FORNCOMP(n)
	{
	    if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
	    {
		if (!constDBCs) p->getNodeCoord(x);
		node = p->getNode() + n-1;
		dirichletBCs->setBC(node, p->Class(), x, n, time); 
	    }	
	}
    }

    mesh->resetEdgIter();
    while ((p=mesh->edgIterAll())) 
    {
	FORNCOMP(n)
	{
	    if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
	    {
		if (!constDBCs) p->getNodeCoord(x);
		node = p->getNode() + n-1;
		dirichletBCs->setBC(node, p->Class(), x, n, time); 
	    }	
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void HQuadElemInt:: updateDirichletBCs(Real time) const
{
    int    n, node;
    PATCH* p;
    Real val, valP1, valP2;
    Vector<Real>  x(spaceDim);

    const Bool constDBCs = dirichletBCs->constant();
    
    dirichletBCs->extend(noOfNodes.Top());

    mesh->resetPtIter();
    while ((p=mesh->ptIterAll())) 
    {
	FORNCOMP(n)
	{
	    if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
	    {
		if (!constDBCs) p->getNodeCoord(x);
		node = p->getNode() + n-1;
		dirichletBCs->setBC(node, p->Class(), x, n, time); 
	    }	
	}
    }

    mesh->resetEdgIter();
    while ((p=mesh->edgIterAll())) 
    {
	FORNCOMP(n)
	{
	    if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
	    {
		if (!constDBCs) p->getNodeCoord(x);

		node = p->getNode() + n-1;
		dirichletBCs->setBC(node, p->Class(), x, n, time); 

		const Vector<PT*>& pt = p->getPoints();

		valP1 = dirichletBCs->value(pt[1]->getNode() + n-1); 
		valP2 = dirichletBCs->value(pt[2]->getNode() + n-1); 
		val   = dirichletBCs->value(node); 

		val -= 0.5 * (valP1 + valP2);   	// hier. basis
		dirichletBCs->setValue(val, node);
	    }	
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LinElemInt:: getNodeCoordinates(NodeCoordinates& nc) const
{
    PATCH* p;
    Vector<Real> x(spaceDim);

    mesh->resetPtIter();
    while ((p=mesh->ptIterAll())) 
    {
	p->getNodeCoord(x);
	nc.setCoordinate(p->getNode(), x);
    }
}
//-------------------------------------------------------------------------


void LQuadElemInt:: getNodeCoordinates(NodeCoordinates& nc) const
{
    PATCH* p;
    Vector<Real> x(spaceDim);

    mesh->resetPtIter();
    mesh->resetEdgIter();

    while ((p=mesh->ptIterAll())) 
    {
	p->getNodeCoord(x);
	nc.setCoordinate(p->getNode(), x);
    }
    while ((p=mesh->edgIterAll()))
    {
	if (p->getNode())
	{
	    p->getNodeCoord(x);
	    nc.setCoordinate(p->getNode(), x);
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LinElemInt:: setMLNodeNumbers(int maxDepth, int maxDepthM1, int targetDepth)
{
    int  depth;
    PATCH* p;

    if (maxDepth > maxDepthM1)  noOfNodes.push(0);  

    mesh->resetPtIter();

    if (maxDepth == 0)  
    {
	 while ((p=mesh->ptIterAll()))  p->setNode(0);
	 mesh->resetPtIter();
    }
    
    while ((p=mesh->ptIterAll())) 
    {
	PT* pt = p->castToPT();

	if (pt->getNode() == 0)			// a new node
	{
	    depth = pt->Depth();
	    mesh->newNodeStack(pt, depth, depth+targetDepth);

	    pt->pushNode(noOfNodes[depth] + 1);
	    noOfNodes[depth] += nComp;

	    pt->setMark();			// new nodes are marked
	}
	else  pt->unMark();
	
	for (depth=pt->maxDepth()+1; depth<=maxDepth; ++depth) 
	{
	    pt->pushNode(noOfNodes[depth] + 1);
	    noOfNodes[depth] += nComp;
	}

	p->setNode(p->MLNode(maxDepth)); 
    }
}
//-------------------------------------------------------------------------


void LinElemInt:: interpolateMLSolution(Vector<Num>& u, int maxDepth, 
					int maxDepthM1) const
{
    int i, n, node, node0;
    PATCH* p;

    Vector<Num> uPrev(u.high());
    FORALL(u,i) uPrev[i] = u[i];
    u.resize(noOfNodes.Top());
    FORALL(u,i) u[i] = 0.0;


    if (maxDepth > maxDepthM1) 		// a new level was created
    {
	mesh->resetPtIter();

	while ((p=mesh->ptIterAll())) 
	  if (!p->marked()) 					// 'old' node
	  {
	      node  = p->getNode();
	      node0 = p->MLNode(maxDepthM1);
	      FORNCOMP(n)  u[node++] = uPrev[node0++]; 
	  }
    }
    else
    {
	mesh->resetPtIter();

	while ((p=mesh->ptIterAll())) 
	  if (!p->marked()) 					// 'old' node
	  {
	      node = p->getNode();
	      FORNCOMP(n) { u[node] = uPrev[node]; ++node; } 
	  }
    }

    interpolateSolution(u);				// handle new nodes
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
  

void LQuadElemInt:: setHighOrderNodes(int* edNode)
{
    PATCH* patch = 0;

    mesh->resetEdgIter();
    *edNode += 1;

    while ((patch=mesh->edgIterAll())) 
    { 
	patch->setNode(*edNode);  
	*edNode += nComp; 
    }

    *edNode -= 1;
    noOfNodes.push(*edNode);
}
//-------------------------------------------------------------------------
  
  
void HQuadElemInt:: setHighOrderNodes(int* edNode)
{
    PATCH* patch = 0;

    mesh->resetEdgIter();
    *edNode += 1;

    while ((patch=mesh->edgIterAll())) 
    {
	if (patch->isDirichlet()) patch->setNode(0);
	else  			  
	{
	    patch->setNode(*edNode);
	    *edNode += nComp; 
	}
    }

    *edNode -= 1;
    noOfNodes.push(*edNode);
}
//-------------------------------------------------------------------------


void EFTetraInt:: setHighOrderNodes(int* edNode)
{
    int i;
    Vector<PATCH*> neighb(4);
    PATCH* patch = 0;

    mesh->resetEdgIter();
    *edNode += 1;

    while ((patch=mesh->edgIterAll())) 
    {
	if (patch->isDirichlet()) patch->setNode(0);
	else  			  
	{
	    patch->setNode(*edNode); 
	    *edNode += nComp; 
	}
    }

    mesh->resetElemIter();			// first reset faceNodes
    while ((patch=mesh->elemIterAll())) 
    {
	for (i=1; i<=4; ++i)
	{
	    int& node = patch->getFaceNode(i);
	    node = 0;
	}
    }

    mesh->resetElemIter();
    while ((patch=mesh->elemIterAll())) 
    {
	patch->getNeighbours(neighb);
	const Vector<PATCH*>& faces = patch->getFaces();

	FORALL(faces,i)
	{
	    if (faces[i])
	    {
		int& node = patch->getFaceNode(i);

		if (!faces[i]->isDirichlet())
		{
		    node = *edNode; 
		    *edNode += nComp;
		}
	    }
	    else if (neighb[i])
	    {
		int& node = patch->getFaceNode(i);

		if (node == 0)
		{
		    int& neighbNode = neighb[i]->getMyFaceNode(patch);

		    if (neighbNode)
		    {
			cout << "\n*** EFTetraInt:: setHighOrderNodes: "
			     << "faceNode in Neighbour already set\n";
			cout.flush(); abort();
		    }

		    node       = *edNode;
		    neighbNode = *edNode;
		    *edNode += nComp;
		}
	    }
	}
    }

    *edNode -= 1;
    noOfNodes.push(*edNode);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// --  	solution vector to quad. nodal basis: uN = P * uH


void LQuadElemInt:: solToNB(Vector<Num>& u)
{
    int node, n, n1, n2;
    PATCH* edg = 0;

    mesh->resetEdgIter();
    while ((edg=mesh->edgIterAll())) 
    {
	if ((node = edg->getNode()))
	{
	    const Vector<PT*>& p = edg->getPoints();
	    n1 = p[1]->getNode();
	    n2 = p[2]->getNode();
	    FORNCOMP(n) u[node++] += 0.5*(u[n1++] + u[n2++]);
	}
    }
}
//-------------------------------------------------------------------------


// --  	solution vector to quad. hierarchical basis: uH = P**(-1) * uN


void LQuadElemInt:: solToHB(Vector<Num>& u)
{
    int node, n, n1, n2;
    PATCH* edg = 0;

    mesh->resetEdgIter();
    while ((edg=mesh->edgIterAll())) 
    {
	if ((node = edg->getNode()))
	{
	    const Vector<PT*>& p = edg->getPoints();
	    n1 = p[1]->getNode();
	    n2 = p[2]->getNode();
	    FORNCOMP(n) u[node++] -= 0.5*(u[n1++] + u[n2++]);
	}
    }
} 
//-------------------------------------------------------------------------


// --  	right-hand-side to quad. nodal basis: bN = P**(-T) * bH


void LQuadElemInt:: rhsToNB(Vector<Num>& b)
{
    int n, node, nn, n1;
    PATCH* edg = 0;

    mesh->resetEdgIter();
    while ((edg=mesh->edgIterAll())) 
    {
	if ((node = edg->getNode()))
	{
	    const Vector<PT*>& p = edg->getPoints();

	    if (!p[1]->isDirichlet())
	    {
		n1 = p[1]->getNode();
		nn = node;
		FORNCOMP(n) b[n1++] -= 0.5 * b[nn++];

	    }
	    if (!p[2]->isDirichlet())
	    {
		n1 = p[2]->getNode();
		nn = node;
		FORNCOMP(n) b[n1++] -= 0.5 * b[nn++];

	    }
	}
    }
} 
//-------------------------------------------------------------------------


// --  	right-hand-side to quad. hierarchical basis:   bH = P**T * bN


void LQuadElemInt:: rhsToHB(Vector<Num>& b)
{
    int n, node, nn, n1;
    PATCH* edg = 0;

    mesh->resetEdgIter();
    while ((edg=mesh->edgIterAll())) 
    {
	if ((node = edg->getNode()))
	{
	    const Vector<PT*>& p = edg->getPoints();

	    if (!p[1]->isDirichlet())
	    {
		n1 = p[1]->getNode();
		nn = node;
		FORNCOMP(n) b[n1++] += 0.5 * b[nn++];

	    }
	    if (!p[2]->isDirichlet())
	    {
		n1 = p[2]->getNode();
		nn = node;
		FORNCOMP(n) b[n1++] += 0.5 * b[nn++];

	    }
	}
    }
} 
//-------------------------------------------------------------------------


void LQuadElemInt:: setHBGeneration(Generation& gen)
{
    int node;
    PATCH* pt=0, *edg=0;
    
    mesh->resetPtIter();
    mesh->resetEdgIter();

    while ((pt=mesh->ptIterAll())) 
    {
	PointSon* son = new PointSon;
	node = pt->getNode();
	son->father = node;
	gen.insertSon(node, son);
    }

    mesh->resetEdgIter();
    while ((edg=mesh->edgIterAll())) 
    {
	if ((node = edg->getNode()))
	{
	    EdgeSon* eSon = new EdgeSon;
	    
	    const Vector<PT*>& p = edg->getPoints();
	    eSon->father1 = p[1]->getNode();
	    eSon->father2 = p[2]->getNode();
	    gen.insertSon(node, eSon);
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LinElemInt:: updateMGFamilyTree()
{
    PATCH* ed, *pm;
    int    son;

    LinEdgeGeneration* newGen = new LinEdgeGeneration(noOfNodes.Prev()+1, 
						      noOfNodes.Top());
    familyTree->setGeneration(newGen);

    mesh->resetEdgIter();
    while ((ed=mesh->edgIterAllOfHistory())) 
    {
	if (ed->refined())	
	{
	    pm = ed->getMidPoint();

	    if (pm->marked())		// new node here
	    {
		son = pm->getNode();
		const Vector<PT*>& p = ed->getPoints();
		newGen->setFathers(p[1]->getNode(), p[2]->getNode(), son);
	    }
	}
    }
}
//-------------------------------------------------------------------------


void LinElemInt:: updateMLFamilyTree(int maxDepth, int maxDepthM1)
{
    int  depth, nodeNo, father1, father2;
    PATCH* ed, *pm, *pt;


    for (depth=1; depth<=maxDepthM1; ++depth) 
      familyTree->extendGeneration(depth, noOfNodes[depth]);


    if (maxDepth > maxDepthM1)			// propagate 'old' nodes
    {
	Generation* gen;

	if (nComp == 1) gen = new Generation(noOfNodes[maxDepth]);
	else  		gen = new MCGeneration(noOfNodes[maxDepth], nComp);

	familyTree->setGeneration(gen);

	mesh->resetPtIter();
	while ((pt=mesh->ptIterAll()))  
	{
	    if (!pt->marked())				// 'old' node
	    {
		nodeNo  = pt->MLNode(maxDepth);
		father1 = pt->MLNode(maxDepth-1);
		familyTree->newPointSon(nodeNo, father1, maxDepth);
	    }
	}
    }


    // --   care for the new nodes (loop over edges to get the fathers):


    mesh->resetEdgIter();
    while ((ed=mesh->edgIterAllOfHistory()))  
    {
	if (ed->refined())	
	{
	    pm = ed->getMidPoint();

	    if (pm->marked())			// a new node
	    {
		depth  = pm->Depth();
		nodeNo = pm->MLNode(depth);

		const Vector<PT*>& p = ed->getPoints();
		father1 = p[1]->MLNode(depth-1);
		father2 = p[2]->MLNode(depth-1);

		familyTree->newEdgeSon(nodeNo, father1, father2, depth);


		// --  	set the remaining levels:

		for (depth=pm->Depth()+1; depth<=maxDepth; ++depth) 
		{					
		    nodeNo  = pm->MLNode(depth);
		    father1 = pm->MLNode(depth-1);

		    familyTree->newPointSon(nodeNo, father1, depth);
		}
	    }
	}
    }
}