////////////////////////////////////////////////////////////////////////////////
//
//  File:        hgrid.cxx
//  Created by:  Jerzy Majewski - jmajewsk@meil.pw.edu.pl
//  Modified by: 
//    2005-04-19 JM: mbIsOutside flag for HGrdTri was not initialized;
//                   CheckBoundIntegr was creating flat or negative triangles
//
// Copyright notice:
//    Copyright (C) 2000-2003 Jerzy Majewski 
// 
//    This is free software. You can redistribute it and/or
//    modify it under the terms of the GNU Lesser General Public License
//    as published by the Free Software Foundation.
//
//    This software is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY. See the GNU Lesser General Public License 
//    for more details.
//
////////////////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include <time.h>
#include "hgrid.h"
#include "EarThwormDef.h"



MGFloat Angle( const Vect2D &nd, const Vect2D &nd1, const Vect2D &nd2)
{
	static Vect2D	v1, v2;
	static MGFloat	dvect, dsin, dcos;

	//v1 = nd1 - nd;
	//v2 = nd2 - nd;
	v1 = (nd1 - nd).module();
	v2 = (nd2 - nd).module();

	dsin = v1.X()*v2.Y() - v2.X()*v1.Y();
	dcos = v1.X()*v2.X() + v1.Y()*v2.Y();
	if ( fabs( dsin) < ZERO && fabs( dcos) < ZERO)
		return M_PI_2;
	dvect = atan2( dsin, dcos);
	return dvect;
}


bool CheckCrossing( const Vect2D& v1, const Vect2D& v2, const Vect2D& v3, const Vect2D& v4)
{
	Vect2D	vv, vv1, vv2;
	MGFloat	t1, t2, t;
	MGFloat	h1, h2;

	
	vv = v2 - v1;
	vv1 = v3 - v1;
	vv2 = v4 - v1;
	
	ASSERT( fabs( vv.module() ) > ZERO );

	t1 = vv * vv1 / vv.module();
	t2 = vv * vv2 / vv.module();
	
	h1 = (vv.X()*vv1.Y() - vv.Y()*vv1.X()) / vv.module();
	h2 = (vv.X()*vv2.Y() - vv.Y()*vv2.X()) / vv.module();

	if ( fabs( h2 - h1) < ZERO)
		return false;

	t = t1 - (t2 - t1)/(h2 - h1) * h1;
	
	if ( t > 0.0 && t < vv.module() && h1 * h2 < ZERO )
		return true;
	else
		return false;
}


bool AreNeigbours( const IterGCell& ic1, const IterGCell& ic2)
{
	char	sbuf[256];
	bool	b1, b2;

//	TRACE2( "itrcl = %d %d", ic1, *ic1);
//	(*ic1)->DumpTri();
//	TRACE2( "itrcl = %d %d", ic2, *ic2);
//	(*ic2)->DumpTri();

	b1 = b2 = false;
	for ( int i=0; i<NUM_TRI; ++i)
	{
//		if ( (*ic1)->Cell(i) != NULL)
			if ( (*ic1)->Cell(i) == ic2 )
				b1 = true;

//		if ( (*ic2)->Cell(i) != NULL)
			if ( (*ic2)->Cell(i) == ic1 )
				b2 = true;
	}

	sprintf( sbuf, " b1 = %d b2 = %d", static_cast<int>(b1), static_cast<int>(b2) );
	TRACE( sbuf);
	if ( !b1 || !b2)
		THROW_INTERNAL( "Not neighbours found");

	return true;
}

//////////////////////////////////////////////////////////////////////
/// class HGrdTri
//////////////////////////////////////////////////////////////////////
HGrdTri::HGrdTri() : mbIsOutside( false)	
{
	for ( MGInt i=0; i<NUM_TRI; i++)
	{
		mlstNod[i] = IterGPnt();
		mlstCell[i] = IterGCell();
	}
}

bool HGrdTri::IsVisible( const IterGCell& icl, const Vect2D& vct)
{
	static Vect2D	vfac, v1, v2;
	static MGFloat	d1, d2;
	static char		sbuf[1024];

	if ( Cell(0) == icl )
	{
		v1 = *( (*Node(0)) ) - *( (*Node(2)) );
		v2 = *( (*Node(1)) ) - *( (*Node(2)) );
		vfac = vct - *( (*Node(2)) );
	}
	else if ( Cell(1) == icl )
	{
		v1 = *( (*Node(1)) ) - *( (*Node(0)) );
		v2 = *( (*Node(2)) ) - *( (*Node(0)) );
		vfac = vct - *( (*Node(0)) );
	}
	else if ( Cell(2) == icl )
	{
		v1 = *( (*Node(2)) ) - *( (*Node(1)) );
		v2 = *( (*Node(0)) ) - *( (*Node(1)) );
		vfac = vct - *( (*Node(1)) );
	}
	else
	{
		ASSERT( 0);
	}
	d1 = v1.X()*vfac.Y() - v1.Y()*vfac.X();
	d2 = v2.X()*vfac.Y() - v2.Y()*vfac.X();

	d1 /= vfac.module();
	d2 /= vfac.module();
	
    if ( (d2 > ZERO && d1 < -ZERO) )
        return false;
    else
        return true;
}

bool HGrdTri::IsVisibleDump( const IterGCell& icl, const Vect2D& vct)
{
//	static Vect2D	vfac, v1, v2;
//	static MGFloat	d1, d2;
//	static char		sbuf[1024];
	Vect2D	vfac, v1, v2;
	MGFloat	d1, d2;

	THROW_INTERNAL("Should not be used !!!");

	if ( Cell(0) == icl )
	{
		v1 = *( (*Node(0)) ) - *( (*Node(2)) );
		v2 = *( (*Node(1)) ) - *( (*Node(2)) );
		vfac = vct - *( (*Node(2)) );
	}
	else if ( Cell(1) == icl )
	{
		v1 = *( (*Node(1)) ) - *( (*Node(0)) );
		v2 = *( (*Node(2)) ) - *( (*Node(0)) );
		vfac = vct - *( (*Node(0)) );
	}
	else if ( Cell(2) == icl )
	{
		v1 = *( (*Node(2)) ) - *( (*Node(1)) );
		v2 = *( (*Node(0)) ) - *( (*Node(1)) );
		vfac = vct - *( (*Node(1)) );
	}
	else
	{
		ASSERT( 0);
	}
	d1 = v1.X()*vfac.Y() - v1.Y()*vfac.X();
	d2 = v2.X()*vfac.Y() - v2.Y()*vfac.X();
	
    if ( d1 * d2 > ZERO )
	{
		TRACE2( "v1 = (%lg %lg)", v1.X(), v1.Y() );
		TRACE2( "v2 = (%lg %lg)", v2.X(), v2.Y() );
		TRACE2( "vf = (%lg %lg)", vfac.X(), vfac.Y() );
		TRACE2( "d1 = %lg d1 = %lg", d1, d2 );

        return false;
	}
    else
        return true;
}



void HGrdTri::SetNeighbour( const IterGCell& itrcl)
{
	static HGrdTri	*ptr;
	
	ASSERT( itrcl != NULL);
	ptr = (*itrcl);
	ASSERT( ptr);
	
	if ( ( ptr->Node(1) == Node(0) && ptr->Node(0) == Node(1) ) ||
		 ( ptr->Node(0) == Node(0) && ptr->Node(2) == Node(1) ) ||
		 ( ptr->Node(2) == Node(0) && ptr->Node(1) == Node(1) ) )
	{
		rCell(0) = itrcl;
	}
	else
	if ( ( ptr->Node(1) == Node(1) && ptr->Node(0) == Node(2) ) ||
		 ( ptr->Node(0) == Node(1) && ptr->Node(2) == Node(2) ) ||
		 ( ptr->Node(2) == Node(1) && ptr->Node(1) == Node(2) ) )
	{
		rCell(1) = itrcl;
	}
	else
	if ( ( ptr->Node(1) == Node(2) && ptr->Node(0) == Node(0) ) ||
		 ( ptr->Node(0) == Node(2) && ptr->Node(2) == Node(0) ) ||
		 ( ptr->Node(2) == Node(2) && ptr->Node(1) == Node(0) ) )
	{
		rCell(2) = itrcl;
	}
}


void HGrdTri::NullifyThis( HGrdTri *pcl)
{
	for ( MGInt i=0; i<NUM_TRI; i++)
	{
		if ( mlstCell[i] != IterGCell())
			if ( (*mlstCell[i]) == pcl)
				mlstCell[i] = IterGCell();
	}
}


void HGrdTri::NullifyThis( HFroSeg *pseg)
{
	if ( ( pseg->PntLf() == Node(0) && pseg->PntLf() == Node(1) ) ||
		 ( pseg->PntRt() == Node(0) && pseg->PntRt() == Node(1) ) )
	{
		rCell(0) = IterGCell();
	}
	else
	if ( ( pseg->PntLf() == Node(1) && pseg->PntLf() == Node(2) ) ||
		 ( pseg->PntRt() == Node(1) && pseg->PntRt() == Node(2) ) )
	{
		rCell(1) = IterGCell();
	}
	else
	if ( ( pseg->PntLf() == Node(2) && pseg->PntLf() == Node(0) ) ||
		 ( pseg->PntRt() == Node(2) && pseg->PntRt() == Node(0) ) )
	{
		rCell(2) = IterGCell();
	}
}


void HGrdTri::InvalidateNeighb()
{
	for ( MGInt i=0; i<NUM_TRI; i++)
	{
		if ( mlstCell[i] != IterGCell())
			(*mlstCell[i])->NullifyThis( this);
	}
}



bool HGrdTri::IsInside( const Vect2D& vct)
{
// ::TODO:: new and faster algorithm should be introduced
	MGFloat	alf;

	alf  = ::Angle( vct, *(*Node(0)), *(*Node(1)) );
	alf += ::Angle( vct, *(*Node(1)), *(*Node(2)) );
	alf += ::Angle( vct, *(*Node(2)), *(*Node(0)) );

	if ( fabs(alf) < M_PI )
		return false;
	else 
		return true;
}


bool HGrdTri::SetCircCenter()
{
	static MGFloat	x1, y1, x2, y2, x3, y3;
	static MGFloat	xr, yr, d;

	x1 = (*Node(0))->X();
	y1 = (*Node(0))->Y();
	x2 = (*Node(1))->X();
	y2 = (*Node(1))->Y();
	x3 = (*Node(2))->X();
	y3 = (*Node(2))->Y();

	d = y3*(x2 - x1) + y2*(x1 - x3) + y1*(x3 - x2);
	if ( fabs( d) < ZERO)
	{
		DumpTri();
		TRACE1( "d = %lg", d);
		return true;
		// THROW_INTERNAL( "Problem inside SetCircCenter() !!!");
	}

	xr = x1*x1*(y3-y2) + x2*x2*(y1-y3) + x3*x3*(y2-y1) + 
		 y1*y1*(y3-y2) + y2*y2*(y1-y3) + y3*y3*(y2-y1);
	xr *= -0.5/d;

	yr = x1*x1*(x3-x2) + x2*x2*(x1-x3) + x3*x3*(x2-x1) + 
		 y1*y1*(x3-x2) + y2*y2*(x1-x3) + y3*y3*(x2-x1);
	yr *= 0.5/d;

	mCircCenter =  Vect2D( xr, yr);

	return false;
}




Vect2D HGrdTri::Center()
{
	return ( *(*Node(0)) + *(*Node(1)) + *(*Node(2)) )/3.0;
}



IterGCell HGrdTri::NextCell( const Vect2D& vct)
{
	static Vect2D	v1, v2;
	
	if ( Cell(0) != IterGCell())
	{
		v1 = *(*Node(1)) - *(*Node(0));
		v1 = Vect2D( -v1.Y(), v1.X() );
		v2 = ( *(*Node(1)) + *(*Node(0)) )/2.0;
		v2 = vct - v2;
		if ( v1 * v2 < 0.0 )
			return Cell(0);
	}

	if ( Cell(1) != IterGCell())
	{
		v1 = *(*Node(2)) - *(*Node(1));
		v1 = Vect2D( -v1.Y(), v1.X() );
		v2 = ( *(*Node(2)) + *(*Node(1)) )/2.0;
		v2 = vct - v2;
		if ( v1 * v2 < 0.0 )
			return Cell(1);
	}

	if ( Cell(2) != IterGCell())
	{
		v1 = *(*Node(0)) - *(*Node(2));
		v1 = Vect2D( -v1.Y(), v1.X() );
		v2 = ( *(*Node(0)) + *(*Node(2)) )/2.0;
		v2 = vct - v2;
		if ( v1 * v2 < 0.0 )
			return Cell(2);
	}

	return IterGCell();
}

IterGCell  HGrdTri::NextCell( HFroSeg *pseg, const IterGCell& iclprv)
{
	IterGCell	itrnb, itr1, itr2;
	Vect2D		v1, v2, v3, v4;

	v1 = *(*(pseg->PntLf()));
	v2 = *(*(pseg->PntRt()));

	if ( iclprv == IterGCell())
	{

		if ( Node(0) == pseg->PntLf() )
		{
			v3 = *(*Node(1));
			v4 = *(*Node(2));
			itrnb = Cell(1);
			if ( Node(1) == pseg->PntRt() || Node(2) == pseg->PntRt() )
				return IterGCell();
		}
		else if ( Node(1) == pseg->PntLf() )
		{
			v3 = *(*Node(2));
			v4 = *(*Node(0));
			itrnb = Cell(2);
			if ( Node(2) == pseg->PntRt() || Node(0) == pseg->PntRt() )
				return IterGCell();
		}
		else if ( Node(2) == pseg->PntLf() )
		{
			v3 = *(*Node(0));
			v4 = *(*Node(1));
			itrnb = Cell(0);
			if ( Node(0) == pseg->PntRt() || Node(1) == pseg->PntRt() )
				return IterGCell();
		}
		else
		{
			THROW_INTERNAL("NextCell - seg: node not found");
		}

		if ( ::CheckCrossing( v1, v2, v3, v4 ) == true)
		{
			return itrnb;
		}
		else
			return IterGCell();
	}
	else
	{
		int	k;
		for ( int i=0; i<NUM_TRI; ++i)
		{
			if ( iclprv != Cell(i) )
			{
				if ( i == NUM_TRI-1)
					k = 0;
				else
					k = i+1;

				v3 = *(*Node(i));
				v4 = *(*Node(k));

				if ( ::CheckCrossing( v1, v2, v3, v4 ) == true)
					return Cell(i);
			}
		}

		return IterGCell();
	}
}


bool HGrdTri::IsInsideCirc( const Vect2D& vct)
{
	static MGFloat	R2, r2;
	static Vect2D	vtmp, v0;

	v0 = CircCenter();
	vtmp = v0 - *( (*Node(0)) );
	R2 = vtmp * vtmp;

	vtmp = v0 - vct;
	r2 = vtmp * vtmp;

	if ( r2 < R2) 
		return true;
	
	return false;
}




//////////////////////////////////////////////////////////////////////
// class HFront
//////////////////////////////////////////////////////////////////////
HFront::~HFront()
{
	iterator	i;
	for ( i= begin(); i != end(); i++)
		if ( *i != NULL) DELETE_PTR((*i));
}


MGFloat HFront::Angle( const Vect2D& vct)
{
	iterator	itr;
	MGFloat		alf;
	IterGPnt	ipnt1, ipnt2;
	
	alf = 0.0;
	for ( itr = begin(); itr != end(); itr++)
	{
		ipnt1 = (*itr)->PntLf();
		ipnt2 = (*itr)->PntRt();
		alf += ::Angle( vct, *(*ipnt1), *(*ipnt2) );
	}
	
	return alf;
}




//////////////////////////////////////////////////////////////////////
// class HGrid
//////////////////////////////////////////////////////////////////////

HGrid::~HGrid()			
{
	CollGPnt::iterator	itrpnt;
	CollGCell::iterator	itrcell;
	
	for ( itrpnt = mcolPnt.begin(); itrpnt != mcolPnt.end(); itrpnt++)
		if ( (*itrpnt) != NULL) DELETE_PTR((*itrpnt));

	for ( itrcell = mcolCell.begin(); itrcell != mcolCell.end(); itrcell++)
		if ( (*itrcell) != NULL) DELETE_PTR((*itrcell));
}


void HGrid::Init( const std::vector<Vect2D>& tabp, const std::vector<MGInt>& tabn )
{
	MGInt	i, j, nprev;
	IterFro	ifro;

	HGrdPnt		*ppnt;
	IterGPnt	ip0, ipp, ipa;
	HFroSeg		*pfro;
	Vect2D		v0, v1, v2;

	std::map< Vect2D, IterGPnt>				mapNod;
	std::map< Vect2D, IterGPnt>::iterator	imap;

//	char	sbuf[512];
	double	d;

	nprev = 0;
	for ( i=0; i<(MGInt)tabn.size(); ++i)
	{
		v1 = tabp[nprev];
		v2 = tabp[nprev+tabn[i]-1];

		d = (v2-v1).module();
//		sprintf( sbuf, "front %d closure(%lg %lg) (%lg %lg)  = %lg ", i, v1.X(), v1.Y(), v2.X(), v2.Y(), d);
//		TRACE( sbuf);

		if ( (v2-v1).module() < ZERO)
		{
			ifro = mcolFro.insert( mcolFro.end(), HFront() );

			imap = mapNod.find( tabp[nprev]);
			if ( imap != mapNod.end() )
			{
				ip0 = ipp = ipa = (*imap).second;
			}
			else
			{
				ppnt = MGNEW HGrdPnt( tabp[nprev]);
				ppnt->rIndex() = nprev;
				ip0 = ipp = ipa = InsertPoint( ppnt);
				mapNod.insert( std::map< Vect2D, IterGPnt>::value_type( *ppnt, ipa) );
			}

			v0 = *(*ip0); 
		

			for ( j=1; j<tabn[i]; ++j)
			{
				v1 = *(*ipp);
				v2 = tabp[nprev+j];

				if ( (v2 - v1).module() > ZERO)
				{
					if ( j != tabn[i]-1 || (tabp[nprev+j] - v0 ).module() > ZERO)
					{
						imap = mapNod.find( tabp[nprev+j]);
						if ( imap != mapNod.end() )
						{
							ipa = (*imap).second;
							ppnt = *ipa;
						}
						else
						{
							ppnt = MGNEW HGrdPnt( tabp[nprev+j]);
							ppnt->rIndex() = nprev+j;
							ipa = InsertPoint( ppnt);
							mapNod.insert( std::map< Vect2D, IterGPnt>::value_type( *ppnt, ipa) );
						}


						pfro = MGNEW HFroSeg( ipp, ipa);
						(*ifro).insert( (*ifro).end(), pfro);
						ipp = ipa;
					}
				}
			}

			v1 = *(*ipp);
			v2 = *(*ip0);

			if ( (v2 - v1).module() > ZERO)
			{
				pfro = MGNEW HFroSeg( ipp, ip0);
				(*ifro).insert( (*ifro).end(), pfro);
			}
		}

		nprev += tabn[i];
	}



	IterFro		itrfro;
	IterFSeg	itrsg;
	IterGPnt	ip1, ip2;

	for ( itrfro = mcolFro.begin(); itrfro != mcolFro.end(); ++itrfro)
	{
//		TRACE1( "Front size = %d\n", (*itrfro).size() );
		itrsg = (*itrfro).begin();
		ip0 = (*itrsg)->PntLf();
		ip1 = (*itrsg)->PntRt();

		for ( ++itrsg; itrsg != (*itrfro).end(); ++itrsg)
		{
			if ( (*itrsg)->PntLf() != ip1 )
			{
				TRACE( "Front not consistent !!!\n");
			}
			ip1 = (*itrsg)->PntRt();

			v1 = *(*(*itrsg)->PntLf());
			v2 = *(*(*itrsg)->PntRt());

			if ( (v2 - v1).module() < ZERO)
				TRACE1( "seg length = %24.16lg\n", (v2 - v1).module() );
		}

		if ( ip0 != ip1 )
		{
			TRACE( "Front not consistent (closure problem) !!!\n");
		}
	}

//	ASSERT(0);

#ifdef _DEBUG

	FILE	*f = fopen( "front.plt", "wt");

	int isize = 0;
	for ( isize = 0, itrfro = mcolFro.begin(); itrfro != mcolFro.end(); ++itrfro, ++isize)
	{
		fprintf( f, "VARIABLES = \"X\", \"Y\"\n" );
		fprintf( f, "ZONE I=%d F=POINT\n", (*itrfro).size()+1);

		for ( itrsg = (*itrfro).begin(); itrsg != (*itrfro).end(); ++itrsg)
		{
			v1 = *(*(*itrsg)->PntLf());
			v2 = *(*(*itrsg)->PntRt());
			fprintf( f, "%lg %lg\n", v1.X(), v1.Y() );
		}
		fprintf( f, "%lg %lg\n", v2.X(), v2.Y() );
	}

	fclose( f);

#endif // _DEBUG

}



HFroSeg* HGrid::NewFace( MGInt i, IterGCell icl)
{
	HFroSeg*	psg;
	
	THROW_ALLOC( psg = MGNEW HFroSeg() );
	switch ( i)
	{
		case  0:
			psg->rPntLf() = (*icl)->Node(0);
			psg->rPntRt() = (*icl)->Node(1);
			if ( icl != IterGCell())
			{
				psg->rCellUp() = icl;
				psg->rCellLo() = (*icl)->Cell(0);
			}
			break;

		case  1:
			psg->rPntLf() = (*icl)->Node(1);
			psg->rPntRt() = (*icl)->Node(2);
			if ( icl != IterGCell())
			{
				psg->rCellUp() = icl;
				psg->rCellLo() = (*icl)->Cell(1);
			}
			break;
			
		case  2:
			psg->rPntLf() = (*icl)->Node(2);
			psg->rPntRt() = (*icl)->Node(0);
			if ( icl != IterGCell())
			{
				psg->rCellUp() = icl;
				psg->rCellLo() = (*icl)->Cell(2);
			}
			break;
			
		default:
			if ( psg) DELETE_PTR(psg);
			ASSERT( 0);
			return NULL;
	};
	return psg;
}


bool HGrid::CheckNeighb( IterGCell icl, CollFSeg& lstsg, 
						 const Vect2D& vct, const IterGCell& ipvcl)
{
	HGrdTri	*pthis;
	
	pthis = (HGrdTri*)( (*icl) );

	if ( pthis->IsInsideCirc( vct ) == true )
	{
		HGrdTri		*ptri;
		HFroSeg		*pseg;
		IterGCell	itri;
		IterGPnt	ipnt;
		bool		bVis;

		for ( int i=0; i<NUM_TRI; i++)
		{		
			if ( (*icl)->Cell(i) != ipvcl || ipvcl == IterGCell())
			{
				pseg = NewFace( i, icl);	// this allocate memory for pseg !!!
				
				itri = (*icl)->Cell(i);
				bVis = false;
				if ( itri != IterGCell())
					if ( (*itri)->IsOutside() == true)
						itri = IterGCell();
				
				if ( itri != IterGCell())
				{
					ptri = (HGrdTri*)( (*itri) );
					ASSERT( ptri);
	          	  	bVis = ptri->IsVisible( icl, vct);
				}
				
				if ( itri != IterGCell() && bVis)
				{
            	    if ( CheckNeighb( itri, lstsg, vct, icl ) == true )
                	{
						DELETE_PTR(pseg);
                    }
	                else
    	            {
						lstsg.insert( lstsg.end(), pseg );
						pseg->rCellUp() = IterGCell();
					}
	            }
				else
				{
					lstsg.insert( lstsg.end(), pseg );
					pseg->rCellUp() = IterGCell();
				}
			}
		}
					
	    ptri = (HGrdTri*)( (*icl) );
    	ptri->InvalidateNeighb();

        mcolCell.erase( icl);
		DELETE_PTR(ptri);
		return true;
	}

	return false;
}


MGInt HGrid::InsertPointIntoMesh( IterGPnt pntitr)
{
	Vect2D		vct;
	IterGCell	itrcl, itrcl2, itrcl0, itrclout;
	HGrdTri		*ptri;
	CollFSeg	*plstsg;
	IterFSeg	itrsg, itrsg2;
	HFroSeg		*pseg;

	
	static int 	num = 0;
	++num;

	vct = *(*pntitr);
	
//	sprintf( sbuf, "POINT No = %d; x=%14.8lg, y=%14.8lg", num , vct.X(), vct.Y() );
//	TRACE1( "%s", sbuf);
	
	itrcl = mcolCell.begin();

	do
	{
		itrcl0 = (*itrcl)->NextCell( vct);
		if ( itrcl0 == IterGCell())
			break;
		itrcl = itrcl0;
	}
	while ( true);
	
	
	THROW_ALLOC( plstsg = MGNEW CollFSeg );
	
	// next function creates segments bounding Delaunay cavity (stored in plstsg)
	// removes cavity triangles from mcolCell;
	// ALL ITERATORS TO THOSE CELLS ARE THEN INVALID !!!
	// iterators to those cells are set to NULL
	CheckNeighb( itrcl, *plstsg, vct, IterGCell());

	// sorting segments stored in plstsg 
	itrsg = plstsg->end();
	itrsg--;
	do
	{
		for ( itrsg2 = plstsg->begin(); itrsg2 != plstsg->end(); itrsg2++)
		{
			if ( (*itrsg)->PntLf() == (*itrsg2)->PntRt() )
			{
				pseg = (*itrsg2);
				plstsg->erase( itrsg2);
				itrsg = plstsg->insert( itrsg, pseg );
				break;
			}
		}
	}
	while ( itrsg != plstsg->begin() );

	// creating new triangles and connections between triangles in Delaunay cavity	
	itrcl0 = itrcl2 = IterGCell();
	for ( itrsg = plstsg->begin(); itrsg != plstsg->end(); itrsg++)
	{
		THROW_ALLOC( ptri = MGNEW HGrdTri );
		
		itrclout = (*itrsg)->CellLo();
		ptri->rNode(0) = (*itrsg)->PntLf();
		ptri->rNode(1) = (*itrsg)->PntRt();
		ptri->rNode(2) = pntitr;
		ptri->rCell(0) = itrclout;
			
		if ( ptri->SetCircCenter() )
		{
			FILE *f = fopen( "cavity.plt", "wt");
			ExportTECTmp( f);
			fclose( f);
			TRACE1( "num = %d", num);
			TRACE2( "new point = %lg %lg", vct.X(), vct.Y() );
			TRACE1( "no of segs bounding cavity = %d", plstsg->size() );

			FILE *ff = fopen( "cavity_front.plt", "wt");
			fprintf( ff, "VARIABLES = \"X\", \"Y\"\n");
			fprintf( ff, "ZONE I=%d F=POINT\n", plstsg->size()+1 );

			fprintf( ff, "%lg %lg\n", (*(*plstsg->begin())->PntLf())->X(), (*(*plstsg->begin())->PntLf())->Y() );
			for ( itrsg = plstsg->begin(); itrsg != plstsg->end(); itrsg++)
				fprintf( ff, "%lg %lg\n", (*(*itrsg)->PntRt())->X(), (*(*itrsg)->PntRt())->Y() );

			fclose( ff);

			THROW_INTERNAL("Flat triangle !!!");
		}
		

		itrcl = InsertCell( ptri);

		if ( itrclout != IterGCell())
			(*itrclout)->SetNeighbour( itrcl);

		if ( itrcl0 == IterGCell())
			itrcl0 = itrcl;
			
		if ( itrcl2 != IterGCell())
		{
			(*itrcl)->rCell(2)  = itrcl2;
			(*itrcl2)->rCell(1) = itrcl;
		}
		itrcl2 = itrcl;
	}
	if ( itrcl2 != IterGCell() && itrcl0 != IterGCell())
	{
		(*itrcl0)->rCell(2) = itrcl2;
		(*itrcl2)->rCell(1) = itrcl0;
	}	
		
	// removing all segments stored in plstsg
	for ( itrsg = plstsg->begin(); itrsg != plstsg->end(); itrsg++)
		if ( (*itrsg) != NULL ) DELETE_PTR((*itrsg));
		
	if ( plstsg) DELETE_PTR(plstsg);


	return num;
}

bool HGrid::PointExists( const Vect2D& vct)
{
//	Leaf<IterFacGPnt>	*pndlf;
//	IterFacGPnt			ipn;
//
//	pndlf = mPntQTree.ClosestItem( vct );
//	if ( pndlf != NULL)
//	{	
//		ipn = pndlf->Data();
//
//		if ( fabs( (*ipn)->X() - vct.X() ) < ZERO &&
//			 fabs( (*ipn)->Y() - vct.Y() ) < ZERO   )
//		{
//			 return true;
//		}
//	}
//
	return false;
}



// creates basic trinagulation (two triangles) and inserts all boundary points
void HGrid::InitTriangles()
{
	CollGPnt::iterator	itr;
	Vect2D		vmin, vmax, vct;
	bool		bFirst = true;
	HGrdPnt		*pnd1, *pnd2, *pnd3, *pnd4;
	HGrdTri		*ptri1, *ptri2;
	IterGPnt	ind1, ind2, ind3, ind4;
	IterGCell	itri1, itri2;
	
	ind1 = ind2 = ind3 = ind4 = IterGPnt();
	itri1 = itri2 = IterGCell();


	// finding limits
	for ( itr = mcolPnt.begin(); itr != mcolPnt.end(); itr++)
	{
		vct = *(*itr);
		if ( bFirst)
		{
			vmin = vmax = vct;
			bFirst = false;
		}
		else
		{
			if ( vct.X() > vmax.X() ) vmax.rX() = vct.X();
			if ( vct.Y() > vmax.Y() ) vmax.rY() = vct.Y();
			if ( vct.X() < vmin.X() ) vmin.rX() = vct.X();
			if ( vct.Y() < vmin.Y() ) vmin.rY() = vct.Y();
		}
	}
	
	vct = (vmax - vmin)/1.5;
	vmax += vct;
	vmin -= vct;
	
	mBox = HRect( vmin.X(), vmin.Y(), vmax.X(), vmax.Y() ); 
	
	// creating starting triangulation containing two cells and four points
	THROW_ALLOC( pnd1 = MGNEW HGrdPnt( vmin) );
	THROW_ALLOC( pnd2 = MGNEW HGrdPnt( vmax.X(), vmin.Y()) );
	THROW_ALLOC( pnd3 = MGNEW HGrdPnt( vmax) );
	THROW_ALLOC( pnd4 = MGNEW HGrdPnt( vmin.X(), vmax.Y()) );
	
	THROW_ALLOC( ptri1 = MGNEW HGrdTri() );
	THROW_ALLOC( ptri2 = MGNEW HGrdTri() );
		
	mind1 = ind1 = InsertPoint( pnd1);
	mind2 = ind2 = InsertPoint( pnd2);
	mind3 = ind3 = InsertPoint( pnd3);
	mind4 = ind4 = InsertPoint( pnd4);

	itri1 = InsertCell( ptri1);
	itri2 = InsertCell( ptri2);

	ptri1->rNode(0) = ind1;
	ptri1->rNode(1) = ind2;
	ptri1->rNode(2) = ind3;
	
	ptri2->rNode(0) = ind3;
	ptri2->rNode(1) = ind4;
	ptri2->rNode(2) = ind1;
	
	ptri1->rCell(2) = itri2;
	ptri2->rCell(2) = itri1;
	ptri1->SetCircCenter();
	ptri2->SetCircCenter();


	// inserting frontal points into mesh
	IterFro		itrfro;
	IterFSeg	itrsg;

	std::map<HGrdPnt*,int>			mapNod;
	std::map<HGrdPnt*,int>::iterator	imap;

	for ( itrfro = mcolFro.begin(); itrfro != mcolFro.end(); itrfro++)
	{
		for ( itrsg = (*itrfro).begin(); itrsg != (*itrfro).end(); itrsg++)
		{
			itr = (*itrsg)->rPntRt();

			if ( ( imap = mapNod.find( *itr) ) == mapNod.end() )
			{
				InsertPointIntoMesh( itr);
				mapNod.insert( std::map<HGrdPnt*,int>::value_type( *itr, 0) );
			}

		}
	}			

	
#ifdef _DEBUG

	FILE	*ftmp = fopen( "initial.plt", "wt");
	ExportTECTmp( ftmp);
	fclose( ftmp);
	
#endif // _DEBUG
}


bool HGrid::IsOutside( const Vect2D& vct)
{
	IterFro		itrfro;
	IterFSeg	itrsg;
	Vect2D		v1, v2;
	double		x;

//// winding algorithm
//	MGFloat		alf;
//	for ( itrfro = mcolFro.begin(); itrfro != mcolFro.end(); itrfro++)
//	{
//		alf += (*itrfro).Angle( vct);
//	}
//	if ( fabs(alf) < M_PI )
//		return true;
//	else 
//		return false;


// ray casting algorithm
	MGInt	cross = 0;
	for ( itrfro = mcolFro.begin(); itrfro != mcolFro.end(); itrfro++)
		for ( itrsg = (*itrfro).begin(); itrsg != (*itrfro).end(); itrsg++)
		{
			v1 = *(*(*itrsg)->PntLf());
			v2 = *(*(*itrsg)->PntRt());

			if ( ( v1.Y() > vct.Y() && v2.Y() <= vct.Y() ) ||
				 ( v2.Y() > vct.Y() && v1.Y() <= vct.Y() ) )
			{

//				x = ( v1.X()*v2.Y() - v1.Y()*v2.X() ) / ( v2.Y() - v1.Y() );
				x = ( v2.X() - v1.X() ) * ( vct.Y() - v1.Y() ) / ( v2.Y() - v1.Y() ) + v1.X();

				if ( x > vct.X() ) 
					++cross;
			}

		}

	if ( (cross % 2) == 1 )
		return false;
	else 
		return true;
}


void HGrid::FlagOuterTris()
{
	IterFro		itrfro;
	IterFSeg	itrsg;
	IterGCell	itr, itrnb, itrcl;;
	Vect2D		vout, vcnt, vc1, vc2;

	CollGCell	colCell;


//	// flaging all triangles lying outside domain using N^2 algo
//	for ( itr = mcolCell.begin(); itr != mcolCell.end(); itr++)
//	{
//		if( IsOutside( (*itr)->Center() ) )
//			(*itr)->rIsOutside() = true;
//		else
//			(*itr)->rIsOutside() = false;
//	}

	Vect2D	v1, v2, vct;
	MGInt	cross = 0;
	HGrdTri	*ptri;
	MGFloat	x, y1, y2;

	std::multimap<MGFloat, HGrdTri*>				mapCell;
	std::multimap<MGFloat, HGrdTri*>::iterator	imap, ifirst, ilast;

	for ( itr = mcolCell.begin(); itr != mcolCell.end(); itr++)
	{
		vct = (*itr)->Center();
		(*itr)->rCross() = 0;
		mapCell.insert( std::multimap<MGFloat, HGrdTri*>::value_type( vct.Y(), (*itr) ) );
	}


	for ( itrfro = mcolFro.begin(); itrfro != mcolFro.end(); itrfro++)
		for ( itrsg = (*itrfro).begin(); itrsg != (*itrfro).end(); itrsg++)
		{
			v1 = *(*(*itrsg)->PntLf());
			v2 = *(*(*itrsg)->PntRt());
			if ( v1.Y() > v2.Y() )
			{
				y1 = v2.Y();
				y2 = v1.Y();
			}
			else
			{
				y1 = v1.Y();
				y2 = v2.Y();
			}

			ifirst = mapCell.lower_bound( y1 );
			ilast  = mapCell.upper_bound( y2 );

			for ( imap = ifirst; imap != ilast; ++imap)
			{
				ptri = (*imap).second;
				vct = ptri->Center();

				if ( ( v1.Y() > vct.Y() && v2.Y() <= vct.Y() ) ||
					 ( v2.Y() > vct.Y() && v1.Y() <= vct.Y() ) )
				{
					x = ( v2.X() - v1.X() ) * ( vct.Y() - v1.Y() ) / ( v2.Y() - v1.Y() ) + v1.X();

					if ( x > vct.X() ) 
						++(ptri->rCross());
				}
			}
		}
	
	for ( itr = mcolCell.begin(); itr != mcolCell.end(); itr++)
	{
		if ( ((*itr)->rCross() % 2) == 1 )
			(*itr)->rIsOutside() = false;
		else
			(*itr)->rIsOutside() = true;
	}


}


void HGrid::RemoveOuterTris()
{
	IterGCell	itr, itr2;
	HGrdTri		*ptri;

	itr2 = IterGCell();
	for ( itr = mcolCell.begin(); itr != mcolCell.end(); itr++)
	{
		if ( itr2 != IterGCell())
		{
			ptri = *itr2;
    		(*itr2)->InvalidateNeighb();
			mcolCell.erase( itr2);
			DELETE_PTR(ptri);
			itr2 = IterGCell();
		}
		if ( (*itr)->IsOutside() )
			itr2 = itr;
	}
	if ( itr2 != IterGCell())
	{
		ptri = *itr2;
    	(*itr2)->InvalidateNeighb();
		mcolCell.erase( itr2);
		DELETE_PTR(ptri);
		itr2 = IterGCell();
	}

	if ( *mind1) DELETE_PTR(*mind1);
	if ( *mind2) DELETE_PTR(*mind2);
	if ( *mind3) DELETE_PTR(*mind3);
	if ( *mind4) DELETE_PTR(*mind4);
	mcolPnt.erase( mind1);
	mcolPnt.erase( mind2);
	mcolPnt.erase( mind3);
	mcolPnt.erase( mind4);
	mind1 = IterGPnt();
	mind2 = IterGPnt();
	mind3 = IterGPnt();
	mind4 = IterGPnt();
}



bool HGrid::CheckSwapTriangles( HGrdTri *ptri1, HGrdTri *ptri2)
{
	HGrdTri tri1, tri2;

	tri1 = *ptri1;
	tri2 = *ptri2;

	SwapTriangles( &tri1, &tri2, false);

	if ( !tri1.Check() || !tri2.Check() )
		return false;

	if ( tri1.Area() < 0 || tri2.Area() < 0 )
		return false;

	return true;
}

void HGrid::SwapTriangles( HGrdTri *ptri1, HGrdTri *ptri2, bool bgo)
{
	MGInt		ifc1, ifc2;
	IterGPnt	ip1, ip2, ip3, ip4;
	IterGCell	ic1, ic2, ic3, ic4;

	IterGCell itri1, itri2;
	
//	TRACE( "--- swapping !!!");
	
	if ( ptri2->Node(1) == ptri1->Node(0) && ptri2->Node(0) == ptri1->Node(1) )
	{
		ip1 = ptri2->Node(1);
		ip2 = ptri2->Node(2);
		ip3 = ptri2->Node(0);
		ip4 = ptri1->Node(2);
		ifc1 = 0;
		ifc2 = 0;
		ic1 = ptri2->Cell(1);
		ic2 = ptri2->Cell(2);
		ic3 = ptri1->Cell(1);
		ic4 = ptri1->Cell(2);
		itri1 = ptri2->Cell(0);
		itri2 = ptri1->Cell(0);
	}
	else
	if ( ptri2->Node(0) == ptri1->Node(0) && ptri2->Node(2) == ptri1->Node(1) )
	{
		ip1 = ptri2->Node(0);
		ip2 = ptri2->Node(1);
		ip3 = ptri2->Node(2);
		ip4 = ptri1->Node(2);
		ifc1 = 0;
		ifc2 = 2;
		ic1 = ptri2->Cell(0);
		ic2 = ptri2->Cell(1);
		ic3 = ptri1->Cell(1);
		ic4 = ptri1->Cell(2);
		itri1 = ptri2->Cell(2);
		itri2 = ptri1->Cell(0);
	}
	else
	if ( ptri2->Node(2) == ptri1->Node(0) && ptri2->Node(1) == ptri1->Node(1) )
	{
		ip1 = ptri2->Node(2);
		ip2 = ptri2->Node(0);
		ip3 = ptri2->Node(1);
		ip4 = ptri1->Node(2);
		ifc1 = 0;
		ifc2 = 1;
		ic1 = ptri2->Cell(2);
		ic2 = ptri2->Cell(0);
		ic3 = ptri1->Cell(1);
		ic4 = ptri1->Cell(2);
		itri1 = ptri2->Cell(1);
		itri2 = ptri1->Cell(0);
	}
	else
	
	if ( ptri2->Node(1) == ptri1->Node(2) && ptri2->Node(0) == ptri1->Node(0) )
	{
		ip1 = ptri2->Node(1);
		ip2 = ptri2->Node(2);
		ip3 = ptri2->Node(0);
		ip4 = ptri1->Node(1);
		ifc1 = 2;
		ifc2 = 0;
		ic1 = ptri2->Cell(1);
		ic2 = ptri2->Cell(2);
		ic3 = ptri1->Cell(0);
		ic4 = ptri1->Cell(1);
		itri1 = ptri2->Cell(0);
		itri2 = ptri1->Cell(2);
	}
	else
	if ( ptri2->Node(0) == ptri1->Node(2) && ptri2->Node(2) == ptri1->Node(0) )
	{
		ip1 = ptri2->Node(0);
		ip2 = ptri2->Node(1);
		ip3 = ptri2->Node(2);
		ip4 = ptri1->Node(1);
		ifc1 = 2;
		ifc2 = 2;
		ic1 = ptri2->Cell(0);
		ic2 = ptri2->Cell(1);
		ic3 = ptri1->Cell(0);
		ic4 = ptri1->Cell(1);
		itri1 = ptri2->Cell(2);
		itri2 = ptri1->Cell(2);
	}
	else
	if ( ptri2->Node(2) == ptri1->Node(2) && ptri2->Node(1) == ptri1->Node(0) )
	{
		ip1 = ptri2->Node(2);
		ip2 = ptri2->Node(0);
		ip3 = ptri2->Node(1);
		ip4 = ptri1->Node(1);
		ifc1 = 2;
		ifc2 = 1;
		ic1 = ptri2->Cell(2);
		ic2 = ptri2->Cell(0);
		ic3 = ptri1->Cell(0);
		ic4 = ptri1->Cell(1);
		itri1 = ptri2->Cell(1);
		itri2 = ptri1->Cell(2);
	}
	else
	
	if ( ptri2->Node(1) == ptri1->Node(1) && ptri2->Node(0) == ptri1->Node(2) )
	{
		ip1 = ptri2->Node(1);
		ip2 = ptri2->Node(2);
		ip3 = ptri2->Node(0);
		ip4 = ptri1->Node(0);
		ifc1 = 1;
		ifc2 = 0;
		ic1 = ptri2->Cell(1);
		ic2 = ptri2->Cell(2);
		ic3 = ptri1->Cell(2);
		ic4 = ptri1->Cell(0);
		itri1 = ptri2->Cell(0);
		itri2 = ptri1->Cell(1);
	}
	else
	if ( ptri2->Node(0) == ptri1->Node(1) && ptri2->Node(2) == ptri1->Node(2) )
	{
		ip1 = ptri2->Node(0);
		ip2 = ptri2->Node(1);
		ip3 = ptri2->Node(2);
		ip4 = ptri1->Node(0);
		ifc1 = 1;
		ifc2 = 2;
		ic1 = ptri2->Cell(0);
		ic2 = ptri2->Cell(1);
		ic3 = ptri1->Cell(2);
		ic4 = ptri1->Cell(0);
		itri1 = ptri2->Cell(2);
		itri2 = ptri1->Cell(1);
	}
	else
	if ( ptri2->Node(2) == ptri1->Node(1) && ptri2->Node(1) == ptri1->Node(2) )
	{
		ip1 = ptri2->Node(2);
		ip2 = ptri2->Node(0);
		ip3 = ptri2->Node(1);
		ip4 = ptri1->Node(0);
		ifc1 = 1;
		ifc2 = 1;
		ic1 = ptri2->Cell(2);
		ic2 = ptri2->Cell(0);
		ic3 = ptri1->Cell(2);
		ic4 = ptri1->Cell(0);
		itri1 = ptri2->Cell(1);
		itri2 = ptri1->Cell(1);
	}

	ASSERT( itri1 != NULL && itri2 != NULL);
	
	ptri1->rNode(0) = ip2;	
	ptri1->rNode(1) = ip4;	
	ptri1->rNode(2) = ip1;	
	
	ptri1->rCell(0) = itri2;	
	ptri1->rCell(1) = ic4;	
	ptri1->rCell(2) = ic1;	

	
	ptri2->rNode(0) = ip4;	
	ptri2->rNode(1) = ip2;	
	ptri2->rNode(2) = ip3;	
	
	ptri2->rCell(0) = itri1;	
	ptri2->rCell(1) = ic2;	
	ptri2->rCell(2) = ic3;
	
	if ( bgo)
	{
		if ( ic1 != IterGCell() )
		{
			if ( (*ic1)->Cell(0) == itri2)
				(*ic1)->rCell(0) = itri1;
			else if ( (*ic1)->Cell(1) == itri2)
				(*ic1)->rCell(1) = itri1;
			else if ( (*ic1)->Cell(2) == itri2)
				(*ic1)->rCell(2) = itri1;
			else
				ASSERT(0);
		}

		if ( ic3 != IterGCell())
		{
			if ( (*ic3)->Cell(0) == itri1)
				(*ic3)->rCell(0) = itri2;
			else if ( (*ic3)->Cell(1) == itri1)
				(*ic3)->rCell(1) = itri2;
			else if ( (*ic3)->Cell(2) == itri1)
				(*ic3)->rCell(2) = itri2;
			else
				ASSERT(0);
		}
	}

}


void HGrid::CheckBoundIntegr()
{
//	char		sbuf[256];
	IterFro		itrfro;
	IterFSeg	itrsg;
	Vect2D		v1, v2, v3, v4;
//	CollGCell	colCell;
	CollGCell	colPath;
	CollFSeg	colSeg;
	IterFSeg	iseg;
	IterGCell	itr, itrnb, itrcl;
	Vect2D		vcnt;
	bool		bReady, bFound;
	MGInt		k, i=0;

	std::list<IterGCell>				lstCell;
	std::list<IterGCell>::iterator	itrtmp;

	std::multimap< HGrdPnt*, IterGCell >						mapNod;
	std::pair< std::multimap< HGrdPnt*, IterGCell >::iterator, 
	      std::multimap< HGrdPnt*, IterGCell >::iterator>	range[2];	// no of points per froseg

	std::multimap< HGrdPnt*, IterGCell >::iterator			itrmap1, itrmap2;


	// getting number of front segments
	int no = 0;
	for ( itrfro = mcolFro.begin(); itrfro != mcolFro.end(); itrfro++)
		for ( itrsg = (*itrfro).begin(); itrsg != (*itrfro).end(); itrsg++)
			++no;


	for ( itr = mcolCell.begin(); itr != mcolCell.end(); itr++)
		for ( k=0; k<NUM_TRI; ++k)
			mapNod.insert( std::multimap<HGrdPnt*, IterGCell>::value_type( *(*itr)->Node(k), itr) );

	int iii=0;
	int	nrec = 0;
	for ( itrfro = mcolFro.begin(); itrfro != mcolFro.end(); itrfro++)
		for ( itrsg = (*itrfro).begin(); itrsg != (*itrfro).end(); itrsg++)
		{
			++iii;
			//if ( iii - 1000*(iii/1000) == 0)
			//	printf( "triang. iter = %d / %d  nrec = %d\n", iii, no, nrec);

			range[0] = mapNod.equal_range( *(*itrsg)->PntLf() );
			range[1] = mapNod.equal_range( *(*itrsg)->PntRt() );
			
			bReady = false;

			for ( itrmap1=range[0].first; itrmap1!=range[0].second; itrmap1++)
			{
				itrcl = (*itrmap1).second;
				lstCell.insert( lstCell.begin(), itrcl );
				for ( itrmap2=range[1].first; itrmap2!=range[1].second; itrmap2++)
				{
					if ( *itrcl == *(*itrmap2).second )
					{
						bReady = true;
						goto out_of_loops;
					}
				}
			}
			out_of_loops:

	
			// the segment is missing - must be recovered
			if ( ! bReady)
			{
				++nrec;
				//printf( "triang. iter = %d / %d  nrec = %d\n", iii, no, nrec);

				v1 = *(*((*itrsg)->PntLf()));
				v2 = *(*((*itrsg)->PntRt()));


//				TRACE( "\n" );
//				sprintf( sbuf, "seg (%lg, %lg) (%lg, %lg) - lstCell.size = %d", v1.X(), v1.Y(), v2.X(), v2.Y(), lstCell.size() );
//				TRACE( sbuf);

				bFound = false;
				for ( itrtmp = lstCell.begin(); itrtmp != lstCell.end(); ++itrtmp)
				{
					itr = *itrtmp;
					if ( ( itrnb = (*itr)->NextCell( *itrsg, IterGCell()) ) != IterGCell())
					{
						bFound = true;
						break;
					}

				}

				if ( bFound)
				{
					IterGCell	itr1, itr2,	itro1=IterGCell(), itro2=IterGCell();
					std::stack<IterGCell>	stackCell;

//					CheckGrid();

					itr1 = itr;
					itr2 = itrnb;

//					TRACE2( "v1 = %lg %lg", v1.X(), v1.Y() );
//					TRACE2( "v2 = %lg %lg", v2.X(), v2.Y() );

					int niter = 0;

					do 
					{
						++niter;

//						TRACE( "--- loop start");
//						TRACE2( "itrcl = %d / %d", itr1, *itr1);
//						(*itr1)->DumpTri();
//						TRACE2( "itrcl = %d / %d", itr2, *itr2);
//						(*itr2)->DumpTri();

						// check if swap is possible if not then try with another triangles
						while ( /*! (*itr2)->IsVisible( itr1, v1) ||*/ 
							! CheckSwapTriangles( *itr1, *itr2 )
							|| ( itr1 == itro1 && itr2 == itro2 ) ) // avoid swaping the same triangles again
						{
							stackCell.push( itr1);
							itr = itr2;
							itr2 = (*itr)->NextCell( *itrsg, itr1);
							itr1 = itr;
							ASSERT( itr2 != NULL);
						}

						itro1 = itr1;
						itro2 = itr2;


						// pre swapping - remove itr1 and itr2 from map
						for ( k=0; k< NUM_TRI; ++k)
						{
							range[0] = mapNod.equal_range( *(*itr1)->Node(k) );
							for ( itrmap1=range[0].first; itrmap1!=range[0].second; itrmap1++)
							{
								if ( (*itrmap1).second == itr1 )
								{
									mapNod.erase( itrmap1);
									break;
								}
							}

							range[0] = mapNod.equal_range( *(*itr2)->Node(k) );
							for ( itrmap1=range[0].first; itrmap1!=range[0].second; itrmap1++)
							{
								if ( (*itrmap1).second == itr2 )
								{
									mapNod.erase( itrmap1);
									break;
								}
							}
						}

						if ( ! CheckSwapTriangles( *itr1, *itr2 ) )
							THROW_INTERNAL( "Recovery: Swap not possible !!!");
				
						SwapTriangles( *itr1, *itr2 );


						// post swapping - insert modified itr1 and itr2 into map
						for ( k=0; k< NUM_TRI; ++k)
						{
							mapNod.insert(std:: multimap<HGrdPnt*, IterGCell>::value_type( *(*itr1)->Node(k), itr1) );
							mapNod.insert( std::multimap<HGrdPnt*, IterGCell>::value_type( *(*itr2)->Node(k), itr2) );
						}


						if ( stackCell.empty() )
						{
							itrcl = IterGCell();
							itrnb = IterGCell();

							if ( (itr = (*itr1)->NextCell( (*itrsg), IterGCell() )) != IterGCell())
							{
								itrcl = itr1;
								itrnb = itr;
							}

							if ( (itr = (*itr2)->NextCell( (*itrsg), IterGCell() )) != IterGCell())
							{
								itrcl = itr2;
								itrnb = itr;
							}

							itr1 = itrcl;
							itr2 = itrnb;
						}
						else
						{
							itr1 = stackCell.top();
							stackCell.pop();

							if ( stackCell.empty() )
								itr2 = (*itr1)->NextCell( *itrsg, IterGCell());
							else
								itr2 = (*itr1)->NextCell( *itrsg, stackCell.top());
						}


					}
					while ( itr1 != IterGCell() && itr2 != IterGCell() );

					(*itrsg)->mbtmp = true;
				}
				else
				{
					char	sbuf[1024];
					sprintf( sbuf, "recovery problem with seg (%lg, %lg) (%lg, %lg)", v1.X(), v1.Y(), v2.X(), v2.Y() );
					TRACE( sbuf);
					TRACE1( "recovwry lstCell.size = %d", lstCell.size() );
					TRACE1( "%d edges already recovered", nrec);

					FILE *ff = fopen( "lstcell.plt", "wt");
					for ( itrtmp = lstCell.begin(); itrtmp != lstCell.end(); ++itrtmp)
					{
						(*(*itrtmp))->DumpTEC( ff);
					}
					fclose( ff);

					//printf( "%s\n", sbuf);
					//printf( "recovery lstCell.size = %d\n", lstCell.size() );

					(*itrsg)->mbtmp = true;

					FILE *f = fopen( "recovery.plt", "wt");
					ExportTECTmp( f);
					fclose( f);

					THROW_INTERNAL( "Recovery problem !!!");
				}
			}

			lstCell.erase( lstCell.begin(), lstCell.end() );
		}

		//printf( "%d edges recovered\n", nrec);
}

void HGrid::Generate()
{
	clock_t		start, tstart, finish;
	double		duration;

	try
	{
		tstart = clock();

		start = clock();

		InitTriangles();

		finish = clock();
		duration = (double)(finish - start) / CLOCKS_PER_SEC;
		//printf( "InitTriangles() - %2.1f seconds\n\n", duration );

		//FILE *f = fopen( "ini.plt", "wt");
		//ExportTECTmp( f);
		//fclose( f);

		TRACE( "CheckBoundIntegr()");
		//printf( "CheckBoundIntegr()\n");

		start = clock();

		CheckBoundIntegr();

		finish = clock();
		duration = (double)(finish - start) / CLOCKS_PER_SEC;
		//printf( "CheckBoundIntegr() - %2.1f seconds\n\n", duration );

		TRACE( "FlagOuterTris()");
		//printf( "FlagOuterTris()\n");
		FlagOuterTris();

		TRACE( "RemoveOuterTris()");
		//printf( "RemoveOuterTris()\n");
		RemoveOuterTris();

		finish = clock();
		duration = (double)(finish - tstart) / CLOCKS_PER_SEC;
		//printf( "Generation total - %2.1f seconds\n\n", duration );

//		FILE *f = fopen( "param.plt", "wt");
//		ExportTECTmp( f);
//		fclose( f);
	}
	catch ( Except *pe)
	{
		ASSERT( pe);
		TRACE_EXCEPTION( *pe);
		TRACE_TO_STDERR( *pe);
		DELETE_PTR(pe);
	}
}



void HGrid::ExportTECTmp( FILE *f)
{
	IterGPnt	pntitr;
	Vect2D		vct;

	IterGCell	cellitr, itrcl;
	MGInt		id1, id2, id3, id4;
	MGInt		ltri, ltmp;
	std::map<void*,MGInt,std::less<void*> >	tmpMap;

	TRACE( "ExportTEC");


	fprintf( f, "TITLE = \"surface\"\n");
	fprintf( f, "VARIABLES = \"X\", \"Y\"\n");
	fprintf( f, "ZONE T=\"TRIANGLES\", ");
	fprintf( f, "N=%2ld, ", mcolPnt.size() );
	fprintf( f, "E=%2ld, F=FEPOINT, ET=QUADRILATERAL C=BLACK\n ", ltri = (MGInt)mcolCell.size() );
	

	ltmp = 0;
	for ( pntitr = mcolPnt.begin(); pntitr != mcolPnt.end(); pntitr++)
	{
		ltmp++;
		vct = *(*pntitr);
		fprintf( f, "%20.10lg %20.10lg\n", vct.X(), vct.Y() );
		tmpMap.insert( std::map<void*,MGInt,std::less<void*> >::value_type( (void*)(*pntitr), ltmp ));
	}
	
	TRACE( "after Pnt");


	ltri = 0;
	for ( cellitr = mcolCell.begin(); cellitr != mcolCell.end(); cellitr++)
	{	
		ltri++;

		pntitr = (*cellitr)->Node(0);
		id1 = tmpMap[ (*pntitr)];

		pntitr = (*cellitr)->Node(1);
		id2 = tmpMap[ (*pntitr)];

		pntitr = (*cellitr)->Node(2);
		id4 = id3 = tmpMap[ (*pntitr)];
		
		fprintf( f, "%9ld %9ld %9ld %9ld\n", id1, id2, id3, id4 );
	}
	

	TRACE( "after Cells");
}


void HGrid::CheckGrid()
{
//	char		sbuf[256];
	bool		bFound;
	IterGCell	itr, itrnb;

//	TRACE("CheckGrid()\n");
	for ( itr = mcolCell.begin(); itr != mcolCell.end(); ++itr)
	{
//		sprintf( sbuf, "cel = %d   | %d %d %d", itr, (*itr)->Cell(0), (*itr)->Cell(1), (*itr)->Cell(2) );
//		TRACE( sbuf);
		for ( int i=0; i<NUM_TRI; ++i)
		{
			bFound = true;
			itrnb = (*itr)->Cell(i);
			if ( itrnb != IterGCell())
			{
				ASSERT( *itrnb );
				bFound = false;
				for ( int k=0; k<NUM_TRI; ++k)
				{
					if ( (*itrnb)->Cell(k) == itr )
						bFound = true;
				}
			}

			if ( ! bFound)
			{
				THROW_INTERNAL( "Cells pointers are not consistent !");
			}
		}
	}
}



void HGrid::WriteFrontTEC( const char name[])
{
}