/*
 $Id: elements2.cc,v 1.5 1996/11/19 09:51:31 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "elements1.h"
#include "elements2.h"

#include "numerics.h"
#include "triang.h"
#include "triangA.h"
#include "materials.h"
#include "integ.h"

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

static double TriangleSf1(double x, double y)   { return 1-x-y; }
static double TriangleSf2(double x, double /*y*/)   { return  x;  }
static double TriangleSf3(double /*x*/, double y)   { return  y;  }
static double TriangleSf1dx(double /*x*/, double /*y*/) { return -1.; }
static double TriangleSf1dy(double /*x*/, double /*y*/) { return -1.; }
static double TriangleSf2dx(double /*x*/, double /*y*/) { return  1.; }
static double TriangleSf2dy(double /*x*/, double /*y*/) { return  0.; }
static double TriangleSf3dx(double /*x*/, double /*y*/) { return  0.; }
static double TriangleSf3dy(double /*x*/, double /*y*/) { return  1.; }

static double HQuadTriangleSf1(double x, double y)   { return 1-x-y; }
static double HQuadTriangleSf2(double x, double /*y*/)   { return  x;  }
static double HQuadTriangleSf3(double /*x*/, double y)   { return  y;  }
static double HQuadTriangleSf1dx(double /*x*/, double /*y*/) { return -1.; }
static double HQuadTriangleSf1dy(double /*x*/, double /*y*/) { return -1.; }
static double HQuadTriangleSf2dx(double /*x*/, double /*y*/) { return  1.; }
static double HQuadTriangleSf2dy(double /*x*/, double /*y*/) { return  0.; }
static double HQuadTriangleSf3dx(double /*x*/, double /*y*/) { return  0.; }
static double HQuadTriangleSf3dy(double /*x*/, double /*y*/) { return  1.; }
static double HQuadTriangleSf4  (double x, double y) { return 4.*x*y; }
static double HQuadTriangleSf4dx(double /*x*/, double y) { return 4.*y; }
static double HQuadTriangleSf4dy(double x, double /*y*/) { return 4.*x; }
static double HQuadTriangleSf5  (double x, double y) { return  4.*y*(1-x-y); }
static double HQuadTriangleSf5dx(double /*x*/, double y) { return -4.*y; }
static double HQuadTriangleSf5dy(double x, double y) { return  4.*(1-x-2*y); }
static double HQuadTriangleSf6  (double x, double y) { return  4.*x*(1-x-y); }
static double HQuadTriangleSf6dx(double x, double y) { return  4.*(1-2*x-y); }
static double HQuadTriangleSf6dy(double x, double /*y*/) { return -4.*x; }

static double LQuadTriangleSf1(double x, double y)   
						{ return 2.*(1.-x-y)*(0.5-x-y); }
static double LQuadTriangleSf1dx(double x, double y) { return 4.*(x+y)-3.; }
static double LQuadTriangleSf1dy(double x, double y) { return 4.*(x+y)-3.; }  
static double LQuadTriangleSf2(double x, double /*y*/)   { return x*(2.*x-1.); }
static double LQuadTriangleSf2dx(double x, double /*y*/) { return 4.*x-1.; }
static double LQuadTriangleSf2dy(double /*x*/, double /*y*/) { return 0.; }
static double LQuadTriangleSf3(double /*x*/, double y)   { return y*(2.*y-1.); }
static double LQuadTriangleSf3dx(double /*x*/, double /*y*/) { return 0.; }
static double LQuadTriangleSf3dy(double /*x*/, double y) { return 4.*y-1.; }
static double LQuadTriangleSf4(double x, double y)   { return 4.*x*y; }
static double LQuadTriangleSf4dx(double /*x*/, double y) { return 4.*y; }
static double LQuadTriangleSf4dy(double x, double /*y*/) { return 4.*x; }
static double LQuadTriangleSf5(double x, double y)   { return 4.*y*(1.-x-y); }
static double LQuadTriangleSf5dx(double /*x*/, double y) { return -4.*y; }
static double LQuadTriangleSf5dy(double x, double y) { return 4.*(1-x-2.*y); }
static double LQuadTriangleSf6(double x, double y)   { return 4.*x*(1.-x-y); }
static double LQuadTriangleSf6dx(double x, double y) { return 4.*(1-2.*x-y); }
static double LQuadTriangleSf6dy(double x, double /*y*/) { return -4.*x; }


//----------------------------- Triangles ------------------------------


Triangle:: Triangle(Material* material0) : StdElement()	     
{ 
    material = material0;	

    basicInit();
    setShapeFunctionsAtIPs();
    compConstantMatrices();
    setPatterns();
} 
//-------------------------------------------------------------------------


void Triangle:: basicInit()
{
    const int nnodes = NoOfNodes();
    const int dim    = SpaceDim();

    IF      = new GaussForm(dim,1);   		// GaussTriOrder1; 
    IFMass  = new GaussForm(dim,3);        	// GaussTriOrder3; 
    IFLumpedMass = new LumpedForm(dim,1,"LumpedTriOrder1"); 	
 
    InfoIntegFormulas("Triangle");

    edgeElement = new Line(material);
    faceElement = edgeElement;

    tmp_nodes	= new Vector<int>(faceElement->NoOfNodes()+1);

    tmp_Ell	= new Matrix<Real>(dim,dim);
    tmp_dN	= new Matrix<Real>(nnodes,dim);
    tmp_u	= new Vector<Real>(dim);
    tmp_x	= new Vector<Real>(dim);
    tmp_xBnd	= new Vector<Real>(dim+1);

    sf[1] = TriangleSf1; 
    sf[2] = TriangleSf2; 
    sf[3] = TriangleSf3;
    dsf[1][1] = TriangleSf1dx;  dsf[1][2] = TriangleSf1dy;
    dsf[2][1] = TriangleSf2dx;  dsf[2][2] = TriangleSf2dy;
    dsf[3][1] = TriangleSf3dx;  dsf[3][2] = TriangleSf3dy;
}
//-------------------------------------------------------------------------

HQuadTriangleFast:: HQuadTriangleFast(Material* material0) :  HQuadTriangle()
{ 
    material = material0;	

    basicInit();
    setShapeFunctionsAtIPs();
    compConstantMatrices();
    setPatterns();
}
//-------------------------------------------------------------------------


void HQuadTriangleFast:: basicInit()
{
    const int nnodes = NoOfNodes();
    const int dim    = SpaceDim();

    IF      = new GaussForm(dim,2, "GaussTriOrder2OnEdges");
    IFMass  = new GaussForm(dim,5);      
    IFLumpedMass = new LumpedForm(dim,2,"LumpedTriOrder2"); 	
 
    InfoIntegFormulas("HQuadTriangleFast");

    edgeElement = new HQuadLine(material);
    faceElement = edgeElement;

    tmp_nodes	= new Vector<int>(faceElement->NoOfNodes()+1);

    tmp_Ell	= new Matrix<Real>(dim,dim);
    tmp_dN	= new Matrix<Real>(nnodes,dim);
    tmp_u	= new Vector<Real>(dim);
    tmp_x	= new Vector<Real>(dim);
    tmp_xBnd	= new Vector<Real>(dim+1);
    tmp_matE    = new Vector<Real>(IF->NoOfIPoints());

    HQsf[1] = HQuadTriangleSf1; 
    HQsf[2] = HQuadTriangleSf2; 
    HQsf[3] = HQuadTriangleSf3;
    HQsf[4] = HQuadTriangleSf4; 
    HQsf[5] = HQuadTriangleSf5; 
    HQsf[6] = HQuadTriangleSf6;
    dHQsf[1][1] = HQuadTriangleSf1dx;  dHQsf[1][2] = HQuadTriangleSf1dy;
    dHQsf[2][1] = HQuadTriangleSf2dx;  dHQsf[2][2] = HQuadTriangleSf2dy;
    dHQsf[3][1] = HQuadTriangleSf3dx;  dHQsf[3][2] = HQuadTriangleSf3dy;
    dHQsf[4][1] = HQuadTriangleSf4dx;  dHQsf[4][2] = HQuadTriangleSf4dy;
    dHQsf[5][1] = HQuadTriangleSf5dx;  dHQsf[5][2] = HQuadTriangleSf5dy;
    dHQsf[6][1] = HQuadTriangleSf6dx;  dHQsf[6][2] = HQuadTriangleSf6dy;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


    // the basic matrices for constant coefficients (S1, S2, S3, M, B ...) 
    // are chosen due to Schwarz's conventions:


void Triangle:: compConstantMatrices()
{
    int i,k,ip, noIP;
    const int nbnodes = NoOfBaseNodes();

    S1 = new Matrix<double>(nbnodes,nbnodes);
    S2 = new Matrix<double>(nbnodes,nbnodes);
    S3 = new Matrix<double>(nbnodes,nbnodes);
    C1 = new Matrix<double>(nbnodes,nbnodes);
    C2 = new Matrix<double>(nbnodes,nbnodes);
    M  = new Matrix<double>(nbnodes,nbnodes);
    B  = new Vector<double>(nbnodes);

    for (i=1; i<=nbnodes; ++i) 
    { 
	for (k=1; k<=nbnodes; ++k) 
		(*S1)(i,k) = (*S2)(i,k) = (*S3)(i,k) =
		(*C1)(i,k) = (*C2)(i,k) = (*M)(i,k) = 0.0; 
	(*B)[i] = 0.0;
    }

		    // partial stiffness matrices:

    noIP = IF->NoOfIPoints();
    for (ip=1; ip<=noIP; ++ip)
        for (i=1; i<=nbnodes; ++i)
	    for (k=1; k<=i; ++k) 
	    {
		(*S1)(i,k) += IF->Weight(ip) * (*dNU)(ip,i,1)*(*dNU)(ip,k,1);
		(*S2)(i,k) += IF->Weight(ip) *((*dNU)(ip,i,1)*(*dNU)(ip,k,2) + 
					       (*dNU)(ip,i,2)*(*dNU)(ip,k,1));
		(*S3)(i,k) += IF->Weight(ip) * (*dNU)(ip,i,2)*(*dNU)(ip,k,2);
	    }

		    // mass matrix and source vector:

    noIP = IFMass->NoOfIPoints();
    for (ip=1; ip<=IFMass->NoOfIPoints(); ++ip)
	for (i=1; i<=nbnodes; ++i) 
	{
	    (*B)[i] += IFMass->Weight(ip) * (*NMass)(ip,i);
	    for (k=1; k<=i; ++k)  (*M)(i,k) += IFMass->Weight(ip) *
	    					(*NMass)(ip,i) * (*NMass)(ip,k);
	}

    for (i=1; i<=nbnodes; ++i) 		// upper triangles
    	for (k=1; k<i; ++k)  
	{
	    (*S1)(k,i) = (*S1)(i,k);
	    (*S2)(k,i) = (*S2)(i,k);
	    (*S3)(k,i) = (*S3)(i,k);
	     (*M)(k,i) =  (*M)(i,k);
	}

		    // partial convection matrices:

    noIP = IF->NoOfIPoints();
    for (ip=1; ip<=noIP; ++ip)
        for (i=1; i<=nbnodes; ++i)
	    for (k=1; k<=nbnodes; ++k) 
	    {
		(*C1)(i,k) += IF->Weight(ip) * (*N)(ip,i)*(*dNU)(ip,k,1);
		(*C2)(i,k) += IF->Weight(ip) * (*N)(ip,i)*(*dNU)(ip,k,2);
	    }
}
//-------------------------------------------------------------------------


Bool HQuadTriangleFast:: assembleEllip(const PATCH& patch, Matrix<Real>& A, 
				      const Jacobian& Jac, 
				      const Matrix<Bool>* pattern) const
{
    int i, k, ip, noIP; 
    const int nnodes = NoOfNodes();
    const int dim = SpaceDim();

    Real fact, Jabs, sum;
    Real sumA, sumB, d, dx, dy, n1, n2, p;


    if (!material->EllipticTerm(patch.Class())) return False;

    if (pattern == 0) pattern = symPattern;


    if (material->constEllipticTerm(patch.Class()))    // ident. to HQuadTriangle
    {
	assembleConstEllip(patch, A, Jac, pattern);
    }
    else
    {
        // we have to use the integration formulae GaussTriOrder2OnEdges

	Matrix<Real>& dN = *tmp_dN;
	Vector<Real>& u  = *tmp_u;
	Vector<Real>& x  = *tmp_x;
        Vector<Real>& matE = *tmp_matE;

	if (material->isAnisotropic())
	{
	    cout << "\n*** HQuadTriangleFast:: assembleConstEllip: "
	    "anisotropic material not yet implemented \n";
	    cout.flush(); abort();
	}
	else if (!material->isIsotropic())
	{
	    cout << "\n*** HQuadTriangleFast:: assembleConstEllip: no materialtype specified \n";
	    cout.flush(); abort();
	}

	Jabs = Abs(Jac.detJ);

	noIP = IF->NoOfIPoints();
        for (ip=1; ip<=noIP; ++ip)
	{
	    transformDerivDN(dN, ip, Jac);

            IF->Coordinates(ip,u);
	    patch.realCoordinates(u,x);  	// location x in real space
	    
            matE[ip] = material->E(patch.Class(), &x);
	    fact  = Jabs * IF->Weight(ip) * matE[ip]; 
	    
	    for (i=1; i<=nnodes; ++i)
	    {
		if ((*pattern)(i,i))
		{
		    sum = 0.0;
		    for (k=1; k<=dim; ++k)  sum += dN(i,k) * dN(i,k);
		    A(i,i) += sum * fact; 
		}
	    } 
	}

	// --      Line integration for normal components:

        const Matrix<Real>& JJ = Jac.J;
        const Vector<EDG*>& eds = patch.getEdges();

        d = 2.0/(3.0*Jac.detJ);		// ( JJ(1,1)*JJ(2,2) - JJ(1,2)*JJ(2,1) );

        if ( !(eds[1]->isDirichlet()) )
	{
            dx = JJ(2,1) - JJ(1,1);
            dy = JJ(2,2) - JJ(1,2);
            p  = matE[1];  // value on edge 1
            n1 =  dy*p;
            n2 = -dx*p;

            sumA = n1*JJ(2,2)-n2*JJ(2,1);
            sumB = n2*JJ(1,1)-n1*JJ(1,2);
            A(4,1) += (-d)*(sumA+sumB);
            A(4,2) += d*sumA;
            A(4,3) += d*sumB;
	}

        if ( !(eds[2]->isDirichlet()) )
	{
            dx  = -JJ(2,1);
            dy  = -JJ(2,2);
            p  = matE[2];   // value on edge 2
            n1 =  dy*p;
            n2 = -dx*p;

            sumA = n1*JJ(2,2)-n2*JJ(2,1);
            sumB = n2*JJ(1,1)-n1*JJ(1,2);
            A(5,1) += (-d)*(sumA+sumB);
            A(5,2) += d*sumA;
            A(5,3) += d*sumB;
	}

        if ( !(eds[3]->isDirichlet()) )
	{
            dx  = JJ(1,1);
            dy  = JJ(1,2);
            p  = matE[3];   // value on edge 3
            n1 =  dy*p;
            n2 = -dx*p;

            sumA = n1*JJ(2,2)-n2*JJ(2,1);
            sumB = n2*JJ(1,1)-n1*JJ(1,2);
            A(6,1) += (-d)*(sumA+sumB);
            A(6,2) += d*sumA;
            A(6,3) += d*sumB;
	}
    }

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


HQuadTriangle:: HQuadTriangle(Material* material0) : Triangle()	
{ 
    material = material0;	

    basicInit();
    setShapeFunctionsAtIPs();
    compConstantMatrices();
    setPatterns();
} 
//-------------------------------------------------------------------------


void HQuadTriangle:: basicInit()
{
    const int nnodes = NoOfNodes();
    const int dim    = SpaceDim();

    IF      = new GaussForm(dim,3);   
    IFMass  = new GaussForm(dim,5);   
    IFLumpedMass = new LumpedForm(dim,2,"LumpedTriOrder2"); 	
 
    InfoIntegFormulas("HQuadTriangle");

    edgeElement = new HQuadLine(material);
    faceElement = edgeElement;

    tmp_nodes	= new Vector<int>(faceElement->NoOfNodes()+1);

    tmp_Ell	= new Matrix<Real>(dim,dim);
    tmp_dN	= new Matrix<Real>(nnodes,dim);
    tmp_u	= new Vector<Real>(dim);
    tmp_x	= new Vector<Real>(dim);
    tmp_xBnd	= new Vector<Real>(dim+1);

    HQsf[1] = HQuadTriangleSf1; 
    HQsf[2] = HQuadTriangleSf2; 
    HQsf[3] = HQuadTriangleSf3;
    HQsf[4] = HQuadTriangleSf4; 
    HQsf[5] = HQuadTriangleSf5; 
    HQsf[6] = HQuadTriangleSf6;
    dHQsf[1][1] = HQuadTriangleSf1dx;  dHQsf[1][2] = HQuadTriangleSf1dy;
    dHQsf[2][1] = HQuadTriangleSf2dx;  dHQsf[2][2] = HQuadTriangleSf2dy;
    dHQsf[3][1] = HQuadTriangleSf3dx;  dHQsf[3][2] = HQuadTriangleSf3dy;
    dHQsf[4][1] = HQuadTriangleSf4dx;  dHQsf[4][2] = HQuadTriangleSf4dy;
    dHQsf[5][1] = HQuadTriangleSf5dx;  dHQsf[5][2] = HQuadTriangleSf5dy;
    dHQsf[6][1] = HQuadTriangleSf6dx;  dHQsf[6][2] = HQuadTriangleSf6dy;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


LQuadTriangle:: LQuadTriangle(Material* material0) : HQuadTriangle()	
{ 
    material = material0;	

    basicInit();
    setShapeFunctionsAtIPs();
    compConstantMatrices();
    setPatterns();
}
//-------------------------------------------------------------------------


void LQuadTriangle:: basicInit()
{
    const int nnodes = NoOfNodes();
    const int dim    = SpaceDim();

    IF      = new GaussForm(dim,3);   
    IFMass  = new GaussForm(dim,5);   
    IFLumpedMass = new LumpedForm(dim,2,"LumpedTriOrder2"); 	
 
    InfoIntegFormulas("LQuadTriangle");

    edgeElement = new LQuadLine(material);
    faceElement = edgeElement;

    tmp_nodes	= new Vector<int>(faceElement->NoOfNodes()+1);

    tmp_Ell	= new Matrix<Real>(dim,dim);
    tmp_dN	= new Matrix<Real>(nnodes,dim);
    tmp_u	= new Vector<Real>(dim);
    tmp_x	= new Vector<Real>(dim);
    tmp_xBnd	= new Vector<Real>(dim+1);

    LQsf[1] = LQuadTriangleSf1; 
    LQsf[2] = LQuadTriangleSf2; 
    LQsf[3] = LQuadTriangleSf3;
    LQsf[4] = LQuadTriangleSf4; 
    LQsf[5] = LQuadTriangleSf5; 
    LQsf[6] = LQuadTriangleSf6;
    dLQsf[1][1] = LQuadTriangleSf1dx;  dLQsf[1][2] = LQuadTriangleSf1dy;
    dLQsf[2][1] = LQuadTriangleSf2dx;  dLQsf[2][2] = LQuadTriangleSf2dy;
    dLQsf[3][1] = LQuadTriangleSf3dx;  dLQsf[3][2] = LQuadTriangleSf3dy;
    dLQsf[4][1] = LQuadTriangleSf4dx;  dLQsf[4][2] = LQuadTriangleSf4dy;
    dLQsf[5][1] = LQuadTriangleSf5dx;  dLQsf[5][2] = LQuadTriangleSf5dy;
    dLQsf[6][1] = LQuadTriangleSf6dx;  dLQsf[6][2] = LQuadTriangleSf6dy;
}
//-------------------------------------------------------------------------


void Triangle:: getLocalFaceNodes(int faceNo, Vector<int>& nodes) const 
{
    if      (faceNo == 1) { nodes[1] = 2;  nodes[2] = 3; }
    else if (faceNo == 2) { nodes[1] = 3;  nodes[2] = 1; }
    else if (faceNo == 3) { nodes[1] = 1;  nodes[2] = 2; }
    else 
    {
	cout << "\n*** Triangle:: getLocalFaceNodes: wrong face number: "
	<< faceNo << "\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------

void HQuadTriangle:: getLocalFaceNodes(int faceNo, Vector<int>& nodes) const
{
    if      (faceNo == 1) { nodes[1] = 2;  nodes[2] = 3;  nodes[3] = 4; }
    else if (faceNo == 2) { nodes[1] = 3;  nodes[2] = 1;  nodes[3] = 5; }
    else if (faceNo == 3) { nodes[1] = 1;  nodes[2] = 2;  nodes[3] = 6; }
    else 
    {
	cout << "\n*** HQuadTriangle:: getLocalFaceNodes: wrong face number: "
	<< faceNo << "\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------

void HQuadTriangle:: setPatterns()
{
    int i, k;
    const int nnodes = NoOfNodes();

    StdElement::setPatterns();

    
    // -- 	set the assembly pattern for DLY:

    DLYPattern  = new Matrix<int>(nnodes,nnodes);
    
    FORALL_ROWS(*DLYPattern,i)
    FORALL_COLUMNS(*symPattern,k)  (*DLYPattern)(i,k) = False;
 
    for (i=4; i<=nnodes; ++i)
    {
	(*DLYPattern)(i,i) = True;
	for (k=1; k<=3; ++k) (*DLYPattern)(i,k) = True;
    }
}