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

#include "elementsA.h"

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

#include "cmdpars.h"
extern CmdPars Cmd;

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


StdElement:: StdElement() : Element()
{ 
    material = 0;

    IF 		= 0; 
    IFMass 	= 0; 
    IFLumpedMass= 0; 

    N  		= 0;
    dNU 	= 0; 
    NMass 	= 0; 

    edgeElement = 0;
    faceElement = 0;

    M = 0;
    B = 0;
    S1 = S2 = S3 = S4 = S5 = S6 = 0;
    C1 = C2 = C3 = 0;

    tmp_nodes	= 0;		// defined in basicInit(), cf. elements[123].cc
    tmp_Ell	= 0;		//                 "                 "
    tmp_dN	= 0;		//                 "                 "
    tmp_u	= 0;		//                 "                 "
    tmp_x	= 0;		//                 "                 "
    tmp_xBnd	= 0;		//                 "                 "
}
//-------------------------------------------------------------------------


StdElement:: ~StdElement() 
{ 
    delete IF; 
    delete IFMass; 
    delete IFLumpedMass; 

    delete N; 
    delete dNU; 
    delete NMass; 

    if (edgeElement != faceElement)  delete edgeElement;
    delete faceElement;

    delete M;
    delete B;
    delete S1; delete S2; delete S3; delete S4; delete S5; delete S6;
    delete C1; delete C2; delete C3;

    delete tmp_nodes;
    delete tmp_Ell;
    delete tmp_dN;
    delete tmp_u;
    delete tmp_x;
    delete tmp_xBnd;
}
//-------------------------------------------------------------------------

void StdElement:: InfoIntegFormulas(const char* name) const
{
    if (Cmd.isTrue("infoInteg"))
    {    
	cout << "\n  " << name << ": Integration Formulas";
	if (IF)  	      cout << "\n\tIF:"; IF->Info();
	if (IFMass)       cout << "\n\tIFMass:"; IFMass->Info();
	if (IFLumpedMass) cout << "\n\tIFLumpedMass:"; IFLumpedMass->Info();
    }
}
//-------------------------------------------------------------------------

void StdElement:: getLocalFaceNodes(int /*k*/, Vector<int>& /*nodes*/) const
{
    notImplemented("getLocalFaceNodes");
}
//-------------------------------------------------------------------------


void StdElement:: assembleConstEllip(const PATCH& /*patch*/, Matrix<Real>& /*A*/, 
				     const Jacobian& /*Jac*/, 
				     const Matrix<Bool>* /*pattern*/) const
{
    notImplemented("assembleConstEllip");
}
				     
void StdElement:: assembleConstConvec (const PATCH& /*patch*/, Matrix<Real>& /*A*/, 
				       const Jacobian& /*Jac*/, 
				       const Matrix<Bool>* /*pattern*/) const
{
    notImplemented("assembleConstConvec");
}
//-------------------------------------------------------------------------


    // --  compute shape function values at integration points:


void StdElement:: setShapeFunctionsAtIPs(int /*nComp*/)
{
    int k, m, ip, noIP;
    const int nbnodes = NoOfBaseNodes();
    const int dim = SpaceDim();

    Vector<Real> x(dim);

    N     = new Matrix<Real>(IF->NoOfIPoints(),nbnodes);
    dNU   = new Array3<Real>(IF->NoOfIPoints(),nbnodes,dim);
    NMass = new Matrix<Real>(IFMass->NoOfIPoints(),nbnodes);

    noIP = IF->NoOfIPoints();
    for (ip=1; ip<=noIP; ++ip)
    {
	IF->Coordinates(ip,x);

        for (k=1; k<=nbnodes; ++k)
	{
	     (*N)(ip,k) = SF(k,x);
	     for (m=1; m<=dim; ++m) (*dNU)(ip,k,m) = dSF(k,m,x);
	 }
    }

    noIP = IFMass->NoOfIPoints();
    for (ip=1; ip<=noIP; ++ip)
    {
	IFMass->Coordinates(ip,x);

        for (k=1; k<=nbnodes; ++k)  (*NMass)(ip,k) = SF(k,x);
    }
}
//-------------------------------------------------------------------------


Num StdElement:: valueAt(const Vector<Real>& unitCoord, const Vector<Num>& sol) 
									const
{
    int k;
    const int nnodes  = NoOfNodes();

    Num sum = 0.0;

    checkUnitCoord(unitCoord);

    for (k=1; k<=nnodes; ++k) sum += SF(k,unitCoord) * sol[k];
    return sum;
}
//-------------------------------------------------------------------------


void StdElement:: valueAt(const Vector<Real>& unitCoord, const Vector<Num>& sol,
			  Vector<Num>& uInt, int baseNode, int nComp) const
{
    int k, n, node, comp;
    const int nbnodes = NoOfBaseNodes();
    Real SFunct;

    checkUnitCoord(unitCoord);

    node = baseNode;
    for (comp=1; comp<=nComp; ++comp) uInt[node++] = 0.0;
  
    n = 1;
    for (k=1; k<=nbnodes; ++k) 
    {
	SFunct = SF(k,unitCoord);
	node = baseNode; 
	for (comp=1; comp<=nComp; ++comp)  uInt[node++] += SFunct * sol[n++];
    }
}
//-------------------------------------------------------------------------


void StdElement:: checkUnitCoord(const Vector<Real>& unitCoord) const
{
    int k;
    Real csum = 0.0;
    static const Real tiny = -100.0*machPrec(Real(0.0));
    static const Real approx1 = 1.0 - tiny;

    FORALL(unitCoord,k)
    { 
	csum += unitCoord[k]; 
	if (unitCoord[k] < tiny)
	{
	    cout << "\n** Element::valueAt: wrong unit-coordinate "
	         << unitCoord[k] << "\n";
	    unitCoord.print();
 	}
    }
    if (csum > approx1)
    {
	cout << "\n** Element::valueAt: wrong unit-coordinates: "
	     << "sum " << csum << " > 1.0\n";
	unitCoord.print();
    }
}
//-------------------------------------------------------------------------


void StdElement:: gradientAt(int ip, const Jacobian& Jac, 
			     const Vector<Num>& sol, Vector<Num>& grad)  const
{
    int i, k;
    const int nnodes = NoOfNodes();
    const int dim = SpaceDim();

    static Matrix<Real> dN(nnodes,dim);

    transformDerivDN(dN, ip, Jac);
    
    for (i=1; i<=dim; ++i)  grad[i] = 0.0;

    for (k=1; k<nnodes; ++k)
        for (i=1; i<=dim; ++i)  grad[i] += sol[k] * dN(k,i);
}
//-------------------------------------------------------------------------


// --   transform derivatives of unit element at integration point ip
//	to "real" element by chain rule:


void StdElement:: transformDerivDN(Matrix<Real>& dN, int ip, 
				   const Jacobian& Jac) const
{
    int i, k, m;
    const int nbnodes = NoOfBaseNodes();
    const int dim = SpaceDim();
    Real sum;

    for (k=1; k<=nbnodes; ++k)
	for (i=1; i<=dim; ++i) 
	{
	    sum = 0.0;
	    for (m=1; m<=dim; ++m) sum += (Jac.Jinv)(i,m)*(*dNU)(ip,k,m);
	    dN(k,i) = sum;
	}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


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

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

    if (pattern == 0) pattern = symPattern;


    if (material->constEllipticTerm(patch.Class())) 
    {
	assembleConstEllip(patch, A, Jac, pattern);
    }
    else
    {
	Matrix<Real>& dN = *tmp_dN;
	Vector<Real>& u  = *tmp_u;
	Vector<Real>& x  = *tmp_x;

	Jabs = Abs(Jac.detJ);

	if (material->isAnisotropic())
	{
	    Matrix<Real>& Ell = *tmp_Ell;

	    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

		for (k=1; k<=dim; ++k) 
		for (int l=1; l<=dim; ++l) 
		    Ell(k,l) = material->E(patch.Class(),&x,k,l);

		fact  = Jabs * IF->Weight(ip);

		for (i=1; i<=nnodes; ++i)
		for (j=1; j<=nnodes; ++j)
		{
		    if ((*pattern)(i,j))
		    {
			sum = 0.0;
			for (k=1; k<=dim; ++k) 
			for (int l=1; l<=dim; ++l) 
			    sum += Ell(k,l) * dN(i,l) * dN(j,k);
			A(i,j) += sum * fact; 
		    }
		} 
	    }
        }
	else if (material->isIsotropic())
	{
	    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

		fact  = Jabs * IF->Weight(ip) * material->E(patch.Class(),&x); 
	    
		for (i=1; i<=nnodes; ++i)
		for (j=1; j<=nnodes; ++j)
		{
		    if ((*pattern)(i,j))
		    {
			sum = 0.0;
			for (k=1; k<=dim; ++k)  sum += dN(i,k) * dN(j,k);
			A(i,j) += sum * fact; 
		    }
		} 
	    }
	}
	else 
	{
	    cout << "\n*** StdElement:: assembleEllip: no materialtype specified\n";
	    cout.flush(); abort();
	}
    }

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


Bool StdElement:: assembleSource(const PATCH& patch, Vector<Real>& b, 
				 const Jacobian& Jac, 
				 const Matrix<Bool>* pattern, 
				 Real time) const
{
    int i, ip; 
    Real fact, factS;
    const int nnodes = NoOfNodes();

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

    if (pattern == 0) pattern = symPattern;

    Real Jabs = Abs(Jac.detJ);



    if(material->constSourceTerm(patch.Class()))
    {
	fact = material->S(patch.Class(), 0, time) * Jabs;
	for (i=1; i<=nnodes; ++i) 
	  if ((*pattern)(i,i))  b[i] += fact * (*B)[i];
    }
    else
    {
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_x;
	const int noIP = IFMass->NoOfIPoints();

	for (ip=1; ip<=noIP; ++ip)
	{
            IFMass->Coordinates(ip,u);
	    patch.realCoordinates(u,x);  	// location x in real space

            factS = material->S(patch.Class(), &x, time);
	    fact  = Jabs * IFMass->Weight(ip) * factS; 

	    for (i=1; i<=nnodes; ++i) 
	      if ((*pattern)(i,i)) b[i] += fact * (*NMass)(ip,i);
	}
    }

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


Bool StdElement:: assembleConvec(const PATCH& patch, Matrix<Real>& A, 
				 const Jacobian& Jac, 
				 const Matrix<Bool>* pattern) const
{
    int  i, j, k, ip; 
    Real fact, Jabs;
    Real sum;
    
    const int nnodes = NoOfNodes();
    const int dim = SpaceDim();

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

    if (pattern == 0) pattern = asymPattern;


    if (material->constConvectionTerm(patch.Class())) 
    {
	assembleConstConvec(patch, A, Jac, pattern);
    }
    else 
    {
	Matrix<Real>& dN = *tmp_dN;
	Vector<Real>& u  = *tmp_u;
	Vector<Real>& x  = *tmp_x;
	Vector<Real>& c  = u;		// Use same memory !!!

	const int noIP = IF->NoOfIPoints();
	Jabs = Abs(Jac.detJ);

	for (ip=1; ip<=noIP; ++ip)
	{
	    transformDerivDN(dN, ip, Jac);	// location xy in real space:

	    IF->Coordinates(ip,u);
	    patch.realCoordinates(u,x);  	// location x in real space

	    fact  = Jabs * IF->Weight(ip); 
	    for (k=1; k<=dim; ++k) c[k] = material->C(patch.Class(),k,&x); 
	    
	    for (i=1; i<=nnodes; ++i)
	    for (j=1; j<=nnodes; ++j)
	    {
		if ((*pattern)(i,j))
		{
		    sum = 0.0;
		    for (k=1; k<=dim; ++k) sum += (*N)(ip,i)*c[k]*dN(j,k);
		    A(i,j) += sum * fact; 
		}
	    } 
	}
    }

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


Bool StdElement:: assembleMass(const PATCH& patch, Matrix<Real>& A, 
			       const Jacobian& Jac,
			       const Matrix<Bool>* pattern) const
{
    int i, j, ip; 
    Real fact;
    const int nnodes = NoOfNodes();
    const int dim = SpaceDim();

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

    if (pattern == 0) pattern = symPattern;

    Real Jabs = Abs(Jac.detJ);

    if (material->constMassTerm(patch.Class()))
    {
	fact = material->M(patch.Class()) * Jabs;

	for (i=1; i<=nnodes; ++i)
	for (j=1; j<=nnodes; ++j)
	  if ((*pattern)(i,j)) 
	    {
	      A(i,j) +=  fact * (*M)(i,j); 
	    }
    }
    else
    {
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_x;
	const int noIP = IFMass->NoOfIPoints();

	for (ip=1; ip<=noIP; ++ip)
	{
            IFMass->Coordinates(ip,u);
	    patch.realCoordinates(u,x);  	// location x in real space

	    fact = Jabs * IFMass->Weight(ip) * material->M(patch.Class(),&x);
		   
	    for (i=1; i<=nnodes; ++i)
	    for (j=1; j<=nnodes; ++j)
	    {
		if ((*pattern)(i,j))  
		  A(i,j) += fact * (*NMass)(ip,i) * (*NMass)(ip,j);
	    }
	}
    }

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


Bool StdElement:: assembleLumpedMass(const PATCH& patch, Matrix<Real>& A, 
				     const Jacobian& Jac, 
				     const Matrix<Bool>* pattern) const
{
    int i; 
    Real fact;
    const int nnodes = NoOfNodes();

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

    if (pattern == 0) pattern = symPattern;

    Real Jabs = Abs(Jac.detJ);

    if (material->constMassTerm(patch.Class()))
    {
	fact = material->M(patch.Class()) * Jabs;

	for (i=1; i<=nnodes; ++i)
	  if ((*pattern)(i,i))  A(i,i) += fact * IFLumpedMass->Weight(i);
    }
    else
    {
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_x;

	for (i=1; i<=nnodes; ++i)
	{
	    if ((*pattern)(i,i)) 
	    {
                IFLumpedMass->Coordinates(i,u);
		patch.realCoordinates(u,x);  	// location x in real space

		fact = Jabs  * material->M(patch.Class(),&x);
		A(i,i) += fact * IFLumpedMass->Weight(i);
	    }
	}
    }
    
    return True;
}
//-------------------------------------------------------------------------


Bool StdElement:: assembleP(const PATCH& patch, Matrix<Real>& A, 
				const Jacobian& Jac,
				const Matrix<Bool>* pattern) const
{
    int i, j, ip; 
    Real fact;
    const int nnodes = NoOfNodes();

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

    if (pattern == 0) pattern = symPattern;

    Real Jabs = Abs(Jac.detJ);

    if (material->constPTerm(patch.Class()))
    {
	fact = material->P(patch.Class()) * Jabs;

	for (i=1; i<=nnodes; ++i)
	for (j=1; j<=nnodes; ++j)
       	  if ((*pattern)(i,j))  A(i,j) +=  fact * (*M)(i,j); 
    }
    else
    {
	const int noIP = IFMass->NoOfIPoints();
	Vector<Real>& u  = *tmp_u;
	Vector<Real>& x  = *tmp_x;

	for (ip=1; ip<=noIP; ++ip)
	{
            IFMass->Coordinates(ip,u);
	    patch.realCoordinates(u,x);  	// location x in real space

	    fact = Jabs * IFMass->Weight(ip) * material->P(patch.Class(),&x);
							   
	    for (i=1; i<=nnodes; ++i)
	    for (j=1; j<=nnodes; ++j)
	    {
		if ((*pattern)(i,j))  
		  A(i,j) += fact * (*NMass)(ip,i) * (*NMass)(ip,j);
	    }
	}
    }

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


Bool StdElement:: assembleLumpedP(const PATCH& patch, Matrix<Real>& A, 
				     const Jacobian& Jac, 
				     const Matrix<Bool>* pattern) const
{
    int i; 
    Real fact;
    const int nnodes = NoOfNodes();

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

    if (pattern == 0) pattern = symPattern;

    Real Jabs = Abs(Jac.detJ);

    if (material->constPTerm(patch.Class()))
    {
	fact = material->P(patch.Class()) * Jabs;

	for (i=1; i<=nnodes; ++i)
	  if ((*pattern)(i,i))  A(i,i) += fact * IFLumpedMass->Weight(i);
    }
    else
    {
	Vector<Real>& u  = *tmp_u;
	Vector<Real>& x  = *tmp_x;

	for (i=1; i<=nnodes; ++i)
	{
	    if ((*pattern)(i,i)) 
	    {
                IFLumpedMass->Coordinates(i,u);
		patch.realCoordinates(u,x);  	// location x in real space

		fact = Jabs  * material->P(patch.Class(),&x);
		A(i,i) += fact * IFLumpedMass->Weight(i);
	    }
	}
    }
    
    return True;
}
//-------------------------------------------------------------------------


Bool StdElement:: assembleNeumannBCs(const PATCH& patch, Vector<Real>& b,
				     const Matrix<Bool>* pattern,Real time) const
{
    int k, count = 0;
    const Vector<PATCH*>& face = patch.getFaces();
    Vector<int>&         nodes = *tmp_nodes;

    if (pattern == 0) pattern = symPattern;

    FORALL(face,k)
    {
	if (face[k] == 0) continue;			// interior face

	if (material->isNeumann(face[k]->isNeumann(), face[k]->Class()))
	{
	    getLocalFaceNodes(k, nodes);
	    count += faceElement->assNeumannBCOnBoundary(*(face[k]), nodes, 
							 b, pattern, time);
	}
    }
    return count;
}
//-------------------------------------------------------------------------
Bool StdElement:: assembleInnerBCs(const PATCH& patch, Vector<Real>& b,
				     const Matrix<Bool>* pattern,Real time) const
{
    int k, count = 0;
    Vector<PATCH*> face(NoOfFaces()), neighb(NoOfFaces());    

    Vector<int>&         nodes = *tmp_nodes;

    if (pattern == 0) pattern = symPattern;

    if (patch.Class() != 1) return 1;

    patch.getAllFaces(face);
    patch.getNeighbours(neighb);

    FORALL(face,k)
    {
	if ( neighb[k] == 0 ) continue;
	if ( (neighb[k]->Class()) == 2)
	{
	    getLocalFaceNodes(k, nodes);
	    count += faceElement->assBConInnerBoundary(*(face[k]), nodes, 
							 b, pattern, time);
	}
    }

    patch.returnFaces(face);
    return count;
}
//-------------------------------------------------------------------------


Bool StdElement:: assembleCauchyBCs(const PATCH& patch, Matrix<Real>& A, 
				    Vector<Real>& b,
				    const Matrix<Bool>* pattern, Real time) const
{
    int  k, count = 0;
    const Vector<PATCH*>& face = patch.getFaces();
    Vector<int>&         nodes = *tmp_nodes;

    if (pattern == 0) pattern = symPattern;

    FORALL(face,k)
    {
	if (face[k] == 0) continue;			// interior face

	if (material->isCauchy(face[k]->isCauchy(), face[k]->Class()))
	{
	    ++count;
	    getLocalFaceNodes(k, nodes);
	    faceElement->assCauchyBCOnBoundary(*(face[k]), nodes, A, b, 
					       pattern, time);
	}
    }

    return count;
}
//-------------------------------------------------------------------------


// --   this routine may be called by a 'hyperelement':: assembleNeumannBCs,
//	where this is a face element


Bool StdElement:: assNeumannBCOnBoundary(const PATCH& tr, Vector<int>& nodes, 
					 Vector<Real>& b, 
					 const Matrix<Bool>* pattern, 
					 Real time) const
{
    int  i, ip, node; 
    Real fact;
    static Real tiny = machMin(Real(0.0));
    const int nnodes = NoOfNodes();

    Real Jabs = tr.volume();		// length of edge, area of triangle
    if (SpaceDim() == 2) Jabs *= 2.0;   // 2D:integration points normed on 0.5

    if (material->constNeumannTerm(tr.Class()))
    {
	fact = material->Neumann(tr.Class(), 0, time);
	if (Abs(fact) <= tiny) return 0;    // a zero term need not be assembled

	fact *= Jabs;
	for (i=1; i<=nnodes; ++i) 
	{
	    node = nodes[i];
	    if ((*pattern)(node,node))  b[node] += fact * (*B)[i];
	}
    }
    else
    {
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_xBnd;
	const int noIP = IFMass->NoOfIPoints();

	for (ip=1; ip<=noIP; ++ip)
	{
            IFMass->Coordinates(ip,u);

	    tr.realCoordinates(u,x);  	// location x in real space

	    fact = Jabs * IFMass->Weight(ip) * material->Neumann(tr.Class(),
								 &x, time); 
	    for (i=1; i<=nnodes; ++i) 
	    {
		node = nodes[i];
		if ((*pattern)(node,node))  b[node] += fact * (*NMass)(ip,i);
	    }
	}
    }

    return 1;
}
//-------------------------------------------------------------------------

Bool StdElement:: assBConInnerBoundary(const PATCH& tr, Vector<int>& nodes, 
					 Vector<Real>& b, 
					 const Matrix<Bool>* pattern, 
					 Real time) const
{
    int  i, ip, node; 
    Real fact;
    const int nnodes = NoOfNodes();

    static Real tiny = machMin(Real(0.0));

    Real Jabs = tr.volume();		// length of edge, area of triangle
    if (SpaceDim() == 2) Jabs *= 2.0;   // 2D:integration points normed on 0.5

    Vector<Real>& u = *tmp_u;
    Vector<Real>& x = *tmp_xBnd;

    const int noIP = IFMass->NoOfIPoints();
    for (ip=1; ip<=noIP; ++ip)
    {
	IFMass->Coordinates(ip,u);

	tr.realCoordinates(u,x);  	// location x in real space

	fact = Jabs * IFMass->Weight(ip) * material->Inner(tr.Class(), &x, time); 
	for (i=1; i<=nnodes; ++i) 
	{
	    node = nodes[i];
	    if ((*pattern)(node,node))  b[node] += fact * (*NMass)(ip,i);
	}
    }

    return 1;
}
//-------------------------------------------------------------------------

// --   this routine may be called by a 'hyperelement':: assembleCauchyBCs,
//	where this is a face element


void StdElement:: assCauchyBCOnBoundary(const PATCH& tr, Vector<int>& nodes, 
					Matrix<Real>& A, Vector<Real>& b,
					const Matrix<Bool>* pattern, 
					Real time)  const
{
    int  i, j, ip, n1, n2;
    Real fact;
    int nnodes = NoOfNodes();

    assNeumannBCOnBoundary(tr, nodes, b, pattern, time);

    Real Jabs = tr.volume();
    if (SpaceDim() == 2) Jabs *= 2.0;   // 2D:integration points normed on 0.5

    if (material->constCauchyTerm(tr.Class()))
    {
	fact = material->Cauchy(tr.Class(), 0, time) * Jabs;

	for (i=1; i<=nnodes; ++i)
	{
	    n1 = nodes[i]; 
	    for (j=1; j<=nnodes; ++j)  
	    {
		n2 = nodes[j];
		if ((*pattern)(n1,n2))  A(n1,n2) += fact * (*M)(i,j); 
	    }
	}
    }
    else
    {
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_xBnd;

	for (ip=1; ip<=IFMass->NoOfIPoints(); ++ip)
	{
            IFMass->Coordinates(ip,u);
	    tr.realCoordinates(u,x);  	// location x in real space

	    fact = Jabs * IFMass->Weight(ip) * material->Cauchy(tr.Class(),
								&x, time); 
	    for (i=1; i<=nnodes; ++i)
	    {
		n1 = nodes[i]; 
		for (j=1; j<=nnodes; ++j)  
		{
		    n2 = nodes[j]; 
		    if ((*pattern)(n1,n2)) 
		      {
			A(n1,n2) += fact * (*NMass)(ip,i) * (*NMass)(ip,j);
		      }
		}
	    }
	}
    }
}
//-------------------------------------------------------------------------


void StdElement:: assembleL2Norm(const PATCH& /*patch*/, Matrix<Real>& A, 
				 const Jacobian& Jac, 
				 const Matrix<Bool>* pattern) const
{
    int i, j; 
    int nnodes = NoOfNodes();

    if (pattern == 0) pattern = symPattern;
    Real Jabs = Abs(Jac.detJ);

    for (i=1; i<=nnodes; ++i)
        for (j=1; j<=i; ++j) if ((*pattern)(i,j))  A(i,j) += Jabs * (*M)(i,j); 
}
//-------------------------------------------------------------------------


void StdElement:: assembleLNorm(const PATCH& /*patch*/, Vector<Real>& b, 
				const Jacobian& Jac, 
				const Matrix<Bool>* pattern) const
{
    int i; 
    int nnodes = NoOfNodes();

    if (pattern == 0) pattern = symPattern;
    Real Jabs = Abs(Jac.detJ);

    for (i=1; i<=nnodes; ++i) 
	if ((*pattern)(i,i))  b[i] += Jabs * (*B)[i];
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// -- 			   1D - routines:


void StdElement:: assembleConstEllip1(const PATCH& patch, Matrix<Real>& A, 
				      const Jacobian& Jac, 
				      const Matrix<Bool>* pattern) const
{
    int i, j; 
    Real fact, Jabs;
    int nnodes = NoOfNodes();

    Jabs = Abs(Jac.detJ);			// length of line

    fact = material->E(patch.Class()) / Jabs;

    for (i=1; i<=nnodes; ++i)
    for (j=1; j<=nnodes; ++j)
    {
	if ((*pattern)(i,j))  A(i,j) += fact*(*S1)(i,j);
    }
}
//-------------------------------------------------------------------------


void StdElement:: assembleConstConvec1(const PATCH& patch, Matrix<Real>& A, 
				       const Jacobian& Jac,
				       const Matrix<Bool>* pattern) const
{
    int  i, j; 
    Real a, c, Jabs;
    int nnodes = NoOfNodes();

    Jabs = Abs(Jac.detJ);

    c = material->C(patch.Class(),1) / Jabs; 
    a = Jabs*c;

    for (i=1; i<=nnodes; ++i)
    for (j=1; j<=nnodes; ++j)
    {
	if ((*pattern)(i,j))  A(i,j) += a * (*C1)(i,j);
    }
}
//-------------------------------------------------------------------------


Bool StdElement:: assembleNeumannBCs1(const PATCH& patch, Vector<Real>& b, 
				      const Matrix<Bool>* pattern, 
				      Real time) const
{
    int k, count=0;
    const Vector<PATCH*>& pt = patch.getFaces();

    if (pattern == 0) pattern = symPattern;

    for (k=1; k<=2; k++)
    {
	if (material->isNeumann(pt[k]->isNeumann(), pt[k]->Class()))
	{
	    if ((*pattern)(k,k))
	    {
		b[k] += material->Neumann(pt[k]->Class(), 0, time);
		++count;
	    }
	}
    }
    return count;
}
//-------------------------------------------------------------------------


Bool StdElement:: assembleCauchyBCs1(const PATCH& patch, Matrix<Real>& A, 
				     Vector<Real>& b, 
				     const Matrix<Bool>* pattern,
				     Real time) const
{
    int k, count=0;
    const Vector<PATCH*>& pt = patch.getFaces();

    if (pattern == 0) pattern = symPattern;

    for (k=1; k<=2; k++)
    {
	if (material->isCauchy(pt[k]->isCauchy(), pt[k]->Class()))
	{
	    if ((*pattern)(k,k))
	    {
		b[k] += material->Neumann(pt[k]->Class(), 0, time);
		A(k,k) += material->Cauchy(pt[k]->Class(), 0, time);
		++count;
	    }
	}
    }
    return count;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// -- 			   2D - routines:


void StdElement:: assembleConstEllip2(const PATCH& patch, Matrix<Real>& A, 
				      const Jacobian& Jac,
				      const Matrix<Bool>* pattern) const
{
    int i, j; 
    Real a1, a2, a3, fact, Jabs;
    const int nnodes = NoOfNodes();

    Jabs = Abs(Jac.detJ);

    const Matrix<Real>& J = Jac.Jinv;

    if ( material->isAnisotropic() )
    {
	// Use of *S2 and definition of a2 only valid for symmetric tensor E!

	a1 = 0.0;
	a2 = 0.0;
	a3 = 0.0;

	for (int k=1; k<=2; ++k) 
	for (int l=1; l<=2; ++l) 
	{
	    Real  Ekl = material->E(patch.Class(),0,k,l);
	    a1 += Ekl * J(l,1)*J(k,1);
	    a2 += Ekl * J(l,1)*J(k,2); // Only valid for symmetric tensor E!
	    a3 += Ekl * J(l,2)*J(k,2);
	}
	
	fact = Jabs;
    }
    else  if ( material->isIsotropic() )
    {
	a1 = (J(1,1)*J(1,1) + J(2,1)*J(2,1));
	a2 = (J(1,1)*J(1,2) + J(2,1)*J(2,2));
	a3 = (J(1,2)*J(1,2) + J(2,2)*J(2,2));
    
	fact = Jabs * material->E(patch.Class());
    }
    else
    {
	cout << "\n*** StdElement:: assembleConstEllip2: no materialtype specified\n ";
	cout.flush(); abort();
    }


    a1 *= fact;
    a2 *= fact;
    a3 *= fact;
    
    for (i=1; i<=nnodes; ++i)
    for (j=1; j<=nnodes; ++j)
    {
	if ((*pattern)(i,j))
	  A(i,j) += a1*(*S1)(i,j) + a2*(*S2)(i,j) + a3*(*S3)(i,j);
    }
}
//-------------------------------------------------------------------------


void StdElement:: assembleConstConvec2(const PATCH& patch, Matrix<Real>& A, 
				       const Jacobian& Jac, 
				       const Matrix<Bool>* pattern) const
{
    int  i, j, k; 
    Real a1, a2, Jabs;
    Vector<Real>& c = *tmp_u;
    const int nnodes = NoOfNodes();
    const int dim = SpaceDim();

    Jabs = Abs(Jac.detJ);

    const Matrix<Real>& J = Jac.Jinv;

    for (k=1; k<=dim; ++k) c[k] = Jabs * material->C(patch.Class(),k); 
    
    a1 = J(1,1)*c[1] + J(2,1)*c[2];
    a2 = J(1,2)*c[1] + J(2,2)*c[2];
    
    for (i=1; i<=nnodes; ++i)
    for (j=1; j<=nnodes; ++j)
    {
	if ((*pattern)(i,j))  A(i,j) += a1*(*C1)(i,j) + a2*(*C2)(i,j);
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// -- 			   3D - routines:


void StdElement:: assembleConstEllip3(const PATCH& patch, Matrix<Real>& A, 
				      const Jacobian& Jac,
				      const Matrix<Bool>* pattern) const
{
    int i, j; 
    Real a[7], fact, Jabs;
    const int nnodes = NoOfNodes();


    Jabs = Abs(Jac.detJ);

    const Matrix<Real>& J = Jac.Jinv;

    for (i=1; i<=6; ++i) a[i] = 0.0;

    if (material->isAnisotropic())
    {
	// Use of *S2,*S3,*S5 and definition of a2,a3,a5
	// only valid for symmetric tensor E!

	for (int k=1; k<=3; ++k) 
	for (int l=1; l<=3; ++l) 
	{
	    Real    Ekl = material->E(patch.Class(),0,k,l);
	    a[1] += Ekl * J(l,1)*J(k,1);
	    a[2] += Ekl * J(l,1)*J(k,2); // Only valid for symmetric tensor E!
	    a[3] += Ekl * J(l,1)*J(k,3); // Only valid for symmetric tensor E!
	    a[4] += Ekl * J(l,2)*J(k,2);
	    a[5] += Ekl * J(l,2)*J(k,3); // Only valid for symmetric tensor E!
	    a[6] += Ekl * J(l,3)*J(k,3);
	}
	
	fact = Jabs;
    }
    else if (material->isIsotropic())
    {
	for (j=1; j<=3; ++j) 
	{
	    a[1] += J(j,1)*J(j,1);
	    a[2] += J(j,1)*J(j,2);
	    a[3] += J(j,1)*J(j,3);
	    a[4] += J(j,2)*J(j,2);
	    a[5] += J(j,2)*J(j,3);
	    a[6] += J(j,3)*J(j,3);
	}
    
	fact = Jabs * material->E(patch.Class());
    }
    else
    {
	cout << "\n*** StdElement:: assembleConstEllip3: no materialtype specified\n ";
	cout.flush(); abort();
    }


    for (i=1; i<=6; ++i) a[i] *= fact;
    
    for (i=1; i<=nnodes; ++i)
    for (j=1; j<=nnodes; ++j)
    {
	if ((*pattern)(i,j))
	  A(i,j) += a[1]*(*S1)(i,j) + a[2]*(*S2)(i,j) + a[3]*(*S3)(i,j) + 
	            a[4]*(*S4)(i,j) + a[5]*(*S5)(i,j) + a[6]*(*S6)(i,j);
    }
}
//-------------------------------------------------------------------------


void StdElement:: assembleConstConvec3(const PATCH& patch, Matrix<Real>& A, 
				       const Jacobian& Jac,
				       const Matrix<Bool>* pattern) const
{
    int i, j, k; 
    Real a1, a2, a3, Jabs;
    Vector<Real>& c = *tmp_u;
    const int nnodes = NoOfNodes();
    const int dim = SpaceDim();

    Jabs = Abs(Jac.detJ);

    const Matrix<Real>& J = Jac.Jinv;

    for (k=1; k<=dim; ++k) c[k] = Jabs * material->C(patch.Class(),k); 

    a1 = J(1,1)*c[1] + J(2,1)*c[2] + J(3,1)*c[3];
    a2 = J(1,2)*c[1] + J(2,2)*c[2] + J(3,2)*c[3];
    a3 = J(1,3)*c[1] + J(2,3)*c[2] + J(3,3)*c[3];
    
    for (i=1; i<=nnodes; ++i)
    for (j=1; j<=nnodes; ++j)
    {
	if ((*pattern)(i,j))
	  A(i,j) += a1*(*C1)(i,j) + a2*(*C2)(i,j) + a3*(*C3)(i,j);
    }
}