/*
 $Id: triang2tr.cc,v 1.2 1996/10/04 15:07:56 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "triang2tr.h"

#include "utils.h"
#include "numerics.h"

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


const PATCH* MESH2:: findPatch(const Vector<Real>& x, Vector<Real>& xUnit, 
			       const PATCH* /*newPatch*/) const
{
    TR2* tr;
    for (tr=trList[0]->first; tr; tr=tr->next)
    {
	if (tr->inPatch(x,xUnit)) return tr->findPatch(x,xUnit);
    }
    patchNotFound();
    return 0;
}
//-------------------------------------------------------------------------


const PATCH* TR2:: findPatch(const Vector<Real>& x, Vector<Real>& xUnit) const
{
    int k, nTri=0;
    TR2* tr;

    if (refined()) 
    {
	nTri = NoOfSons();

	tr = firstSon;
	for (k=1; k<nTri; ++k) 
	{
	    if (tr->inPatch(x,xUnit))  return tr->findPatch(x,xUnit);
	    tr = tr->next;
	}
	return tr->findPatch(x,xUnit);

	// if (tr->inPatch(x,xUnit))  return tr->findPatch(x,xUnit);
	// patchNotFound();
    }
    else  
    {
	unitCoordinates(x,xUnit);
	return this;
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Bool TR2:: inPatch(const Vector<Real>& x, Vector<Real>& xUnit) const
{
    static const Real tiny = -1000.0*machPrec(Real(0.0));
    static const Real approx1 = 1.0 - tiny;

    unitCoordinates(x,xUnit);

    if (xUnit[1] < tiny)  return False;
    if (xUnit[2] < tiny)  return False;
    if (xUnit[1] + xUnit[2] > approx1)  return False;

    return True;
}
//-------------------------------------------------------------------------

/*
Bool TR2:: inPatch(const Vector<Real>& x, Vector<Real>& xUnit) const
{
    static const Real Tiny = 100.0*machPrec(Real(0.0));
    Real det1, det2, det3, tiny;

    PT2* pt[1+3];
    getPoints(pt);

    det1 = det3P(x,pt[1],pt[2]);
    det2 = det3P(x,pt[2],pt[3]);
    det3 = det3P(x,pt[3],pt[1]);

    tiny = -Tiny*det1*det1;

    if (det1 >= tiny && det2 >= tiny && det3 >= tiny)
    {
	unitCoordinates(x,xUnit);
	return True;
    }
    return False;
}
*/
//-------------------------------------------------------------------------

Real TR2:: det3P(const Vector<Real>& x, const PT2* pt2, const PT2* pt3) const
{
    Real j11, j21, j12, j22;
 
    j11 = pt2->x - x[1];	
    j21 = pt3->x - x[1];
    j12 = pt2->y - x[2];
    j22 = pt3->y - x[2];

    return j11*j22 - j12*j21;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


MESH2Trans:: MESH2Trans(const char* inFileName, MESH2* prevMesh0)

	: MESH2(inFileName,False), 
	  prevMesh(prevMesh0), trTrAlloc(ElementsInBlock)
{
    readTriangulation(fileName);
    if (prevMesh) setInitialPartners(prevMesh);
}
//-------------------------------------------------------------------------


void MESH2Trans:: Refine()
{
    MESH2::Refine();
    if (prevMesh) setPartners();
}
//-------------------------------------------------------------------------


const PATCH* MESH2Trans:: findPatch(const Vector<Real>& x, Vector<Real>& xUnit,
				    const PATCH* newPatch) const
{
    return newPatch->Partner()->findPatch(x,xUnit);
}
//-------------------------------------------------------------------------