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

#include "elements1.h"

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

//-------------------------- shape functions -----------------------------
	

// --			Linear Line:

static double LineSf1(double x)   { return 1.-x; }
static double LineSf2(double x)   { return  x;  }
static double LineSf1dx(double /*x*/) { return -1.; }
static double LineSf2dx(double /*x*/) { return  1.; }


		// -- quadratic Lagrange element:

static double LQuadLineSf1(double x)  { return (1.-x)*(1.-2.*x); }
//{ double L1 = 1.-x, L2 = x;  return L2*(L2-L1); }

static double LQuadLineSf2(double x)  { return  x*(2.*x-1.); }
//{ double L1 = 1.-x, L2 = x;  return L1*(L1-L2); }

static double LQuadLineSf3(double x)  { return 4.*x*(1.-x); }
//{ double L1 = 1.-x, L2 = x;  return 4.*L1*L2; }

static double LQuadLineSf1dx(double x) { return 4.*x-3.; }
//{ double L1 = 1.-x, L2 = x;  return 3.*L2-L1; }

static double LQuadLineSf2dx(double x) { return 4.*x-1.; }
//{ double L1 = 1.-x, L2 = x;  return 3.*L1-L2; }

static double LQuadLineSf3dx(double x)  { return 4.*(1.-2.*x); }
//{ double L1 = 1.-x, L2 = x;  return 4.*(L1-L2); }

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

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

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

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

    IF 	     	 = new GaussForm(dim,1);              // GaussLineOrder1; 
    IFMass 	 = new GaussForm(dim,3);              // GaussLineOrder3; 
    IFLumpedMass = new LumpedForm(dim,1,"LumpedLineOrder1"); 
    
    InfoIntegFormulas("Line");

    edgeElement = 0;  
    faceElement = 0; 
    

    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] = LineSf1; 
    sf[2] = LineSf2; 
    dsf[1][1] = LineSf1dx;
    dsf[2][1] = LineSf2dx;
}
//-------------------------------------------------------------------------

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


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

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

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

		    // stiffness matrix:
    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);


		    // mass matrix and source vector:

    noIP = IFMass->NoOfIPoints();
    for (ip=1; ip<=noIP; ++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);
	    (*M)(k,i) = (*M)(i,k);
	}


		    //   convection matrix:

    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);
	    }

    //cout << "\n S, C, M, B : \n", 
    //     S->print(cout); C->print(cout);  M->print(cout);  B->print(cout);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

	
HQuadLine:: HQuadLine (Material* material0) : Line()	
{ 
    material = material0;	

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


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

    IF 	     	 = new GaussForm(dim,3);           
    IFMass 	 = new GaussForm(dim,5);           
    IFLumpedMass = new LumpedForm(dim,2,"LumpedLineOrder2"); 
    
    InfoIntegFormulas("HQuadLine");
    
    edgeElement = 0;  
    faceElement = 0;  

    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] = LineSf1; 
    HQsf[2] = LineSf2; 
    HQsf[3] = LQuadLineSf3; 

    dHQsf[1][1] = LineSf1dx;
    dHQsf[2][1] = LineSf2dx;
    dHQsf[3][1] = LQuadLineSf3dx;
}
//-------------------------------------------------------------------------


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

    StdElement::setPatterns();
	

    // -- 	set the assembly pattern for DLY:

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

	
LQuadLine:: LQuadLine(Material* material0) : HQuadLine()	
{ 
    material = material0;	

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


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

    IF 	     	 = new GaussForm(dim,3);           
    IFMass 	 = new GaussForm(dim,5);           
    IFLumpedMass = new LumpedForm(dim,2,"LumpedLineOrder2"); 
    
    InfoIntegFormulas("LQuadLine");
    
    edgeElement = 0;  
    faceElement = 0;  

    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] = LQuadLineSf1; 
    LQsf[2] = LQuadLineSf2; 
    LQsf[3] = LQuadLineSf3; 

    dLQsf[1][1] = LQuadLineSf1dx;
    dLQsf[2][1] = LQuadLineSf2dx;
    dLQsf[3][1] = LQuadLineSf3dx;
}