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

#include "triangA.h"

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

#include "cmdpars.h"
extern CmdPars Cmd;


//-------------------------------------------------------------------------
//------------------------   0D objects  ---------------------------------


ostream& operator << (ostream& os, const PT& p) { p.print(os); return os; }

Bool operator == (const PT& p1, const PT& p2) 
{ 
    return p1.equal(p2);
}

Bool PT:: equal(const PT& /*pt*/) const 
{ 
    notImplemented("equal");  return 0;
}
//-------------------------------------------------------------------------

void PT:: print(ostream& /*os*/) const  {   notImplemented("print"); }


//-------------------------------------------------------------------------
//------------------------   1D objects  ---------------------------------


void EDG:: notImplemented(const char* s) const
{
    cout << "\n*** class EDG: " << s << " not implemented\n";  cout.flush(); abort();
}

ostream& operator << (ostream& os, const EDG& ed) { ed.print(os); return os; }

Bool operator == (const EDG& t1, const EDG& t2)
{ 
    return t1.equal(t2);
}

Bool EDG:: equal(const EDG& /*ed*/)  const {  notImplemented("equal"); return 0; }

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

void EDG:: print(ostream& /*os*/) const
{
    notImplemented("print");
}
//-------------------------------------------------------------------------

/*
  void EDG1:: getInnerNodes(Vector<int>& nodes) const { node[1] = node; }
  void EDG2:: getInnerNodes(Vector<int>& nodes) const { node[1] = node; }
  void EDG3:: getInnerNodes(Vector<int>& nodes) const { node[1] = node; }
  
  void EDG1:: getPointNodes(Vector<int>& nodes) const
  {
  nodes[1] = p1->getNode();
  nodes[2] = p2->getNode();
  }
  void EDG2:: getPointNodes(Vector<int>& nodes) const
  {
  nodes[1] = p1->getNode();
  nodes[2] = p2->getNode();
  }
  void EDG3:: getPointNodes(Vector<int>& nodes) const
  {
  nodes[1] = p1->getNode();
  nodes[2] = p2->getNode();
  }
  */
//-------------------------------------------------------------------------

Real EDG:: volume() const
{
    int dim = spaceDim();
    Vector<Real> cP1(dim), cP2(dim);
    
    getCoordinates(cP1,cP2);
    
    if (dim==1) return Abs(cP2[1]-cP1[1]);
    else
    {
	Real l = 0.0;
	for (int i=1; i<=dim; ++i) l += sqr(cP2[i]-cP1[i]);
	return sqrt(l);
    }
}

Real EDG:: lengthSqr() const
{
    int dim = spaceDim();
    Vector<Real> cP1(dim), cP2(dim);
    
    getCoordinates(cP1,cP2);
    
    Real l = 0.0;
    for (int i=1; i<=dim; ++i) l += sqr(cP2[i]-cP1[i]);
    return l;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

void EDG:: centerOfGravity(Vector<Real>& xG)  const
{
    int dim = spaceDim();
    Vector<Real> cP1(dim), cP2(dim);

    getCoordinates(cP1,cP2);
    for (int i=1; i<=dim; ++i)  xG[i] = 0.5 * (cP1[i] + cP2[i]);
}

//-------------------------------------------------------------------------
//------------------------   2D objects  ---------------------------------


Bool operator == (const TR& t1, const TR& t2)
{ 
    return t1.equal(t2);
}

Bool TR:: equal(const TR& /*lr*/)  const {  notImplemented("equal"); return 0; }

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

void TR:: notImplemented(const char* s) const
{
    cout << "\n*** class TR: " << s << " not implemented\n";  cout.flush(); abort();
}
//-------------------------------------------------------------------------

ostream& operator << (ostream& os, const TR& t) { t.print(os); return os; }

void TR:: print(ostream& os) const 
{
    const Vector<PT*>& p = getPoints();
    os << "Triangle points";
    for (int i=1; i<=3; ++i) os << " " << p[i]->getNode();
}

void TR:: centerOfGravity(Vector<Real>& xG) const
{
    int dim = spaceDim();
    Vector<Real> cP1(dim), cP2(dim), cP3(dim);
    const Real OneThird = 1./3.;

    getCoordinates(cP1,cP2,cP3);
    for (int i=1; i<=dim; ++i) xG[i] = OneThird * (cP1[i]+cP2[i]+cP3[i]);
}
//-------------------------------------------------------------------------
/*
  void TR2:: getPointNodes(Vector<int>& nodes) const
  {
  nodes[1] = p1->getNode();
  nodes[2] = p2->getNode();
  nodes[3] = p3->getNode();
  }
  void TR3:: getPointNodes(Vector<int>& nodes) const
  {
  nodes[1] = p1->getNode();
  nodes[2] = p2->getNode();
  nodes[3] = p3->getNode();
  }
  
  void TR2:: getEdgeNodes(Vector<int>& nodes) const
  {
  nodes[4] = e1->getNode();	// e1 is opposite to p1
  nodes[5] = e2->getNode();
  nodes[6] = e3->getNode();
  }
  void TR3:: getEdgeNodes(Vector<int>& nodes) const
  {
  nodes[4] = e1->getNode();	// e1 is opposite to p1
  nodes[5] = e2->getNode();
  nodes[6] = e3->getNode();
  }
  */

//-------------------------------------------------------------------------
//------------------------   3D objects  ---------------------------------


ostream& operator << (ostream& os, const TET& t) { t.print(os); return os; }

Bool operator == (const TET& t1, const TET& t2)  { return t1.equal(t2); }

Bool TET:: equal(const TET& /*t*/)  const { notImplemented("equal"); return 0; }

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

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

void TET:: centerOfGravity(Vector<Real>& xG) const
{
    int dim = spaceDim();
    Vector<Real> cP1(dim), cP2(dim), cP3(dim), cP4(dim);

    getCoordinates(cP1,cP2,cP3,cP4);
    for (int i=1; i<=dim; ++i) xG[i] = 0.25 * (cP1[i]+cP2[i]+cP3[i]+cP4[i]);
}

void TET:: print(ostream& os) const 
{
    const Vector<PT*>& p = getPoints();

    os << "Tetra points";
    for (int i=1; i<=4; ++i) os << " " << p[i]->getNode();
}

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