/*
 $Id: elementsB.cc,v 1.4 1996/11/19 09:53:55 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"

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


// --	assembling functions for multi-component elements; these ones are
//      practically useless, but may serve as specimen for problem-dependent
//	ones


Bool StdElement:: MCEllip(const PATCH& patch, Matrix<Real>& A, 
			  const Jacobian& Jac, const Matrix<Bool>* pattern, 
			  int nComp, const Matrix<int>& Node) const
{
    int  i, j, k, ip, comp1, comp2, node1, node2; 
    Real fact, Jabs, sum;
    const int nbnodes = NoOfBaseNodes();
    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);
	const int noIP = IF->NoOfIPoints();

	if (material->isAnisotropic())
	{
	    cout << "\n*** StdElement:: MCEllip: "
	    "anisotropic material not yet implemented \n";
	    cout.flush(); abort();
	}
	else if (!material->isIsotropic())
	{
	    cout << "\n*** StdElement:: MCEllip: no materialtype specified \n";
	    cout.flush(); abort();
	}
	for (ip=1; ip<=noIP; ++ip)
	{
	    transformDerivDN(dN, ip, Jac);

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

	    for (comp1=1; comp1<=nComp; ++comp1)
	    for (comp2=1; comp2<=nComp; ++comp2)
	    {
		if (comp1 != comp2) continue;

		fact  = Jabs * IF->Weight(ip) * material->E(patch.Class(),&x); 
	    
		for (i=1; i<=nbnodes; ++i)
		{
		    node1 = Node(i,comp1);

		    for (j=1; j<=nbnodes; ++j)
		    {
			node2 = Node(j,comp2);
			if ((*pattern)(node1,node2))
			{
			    sum = 0.0;
			    for (k=1; k<=dim; ++k) sum += dN(i,k)*dN(j,k);

			    A(node1,node2) += sum * fact; 
			}
		    }
		} 
	    }
	}
    }
    return True;
}
//-------------------------------------------------------------------------


Bool StdElement:: MCConvec(const PATCH& patch, Matrix<Real>& A, 
			   const Jacobian& Jac, const Matrix<Bool>* pattern, 
			   int nComp, const Matrix<int>& Node) const
{
    int  i, j, k, ip, comp1, comp2, node1, node2; 
    Real fact, Jabs;
    Real sum;
    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
    {
	const int noIP = IF->NoOfIPoints();
	const int nbnodes = NoOfBaseNodes();

	Matrix<Real>& dN = *tmp_dN;
	Vector<Real>& u  = *tmp_u;
	Vector<Real>& x  = *tmp_x;
	Vector<Real>& c  = u;		// Use same memory !!!

	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 (comp1=1; comp1<=nComp; ++comp1)
	    for (comp2=1; comp2<=nComp; ++comp2)
	    {
		if (comp1 != comp2) continue;

		for (k=1; k<=dim; ++k) c[k] = material->C(patch.Class(),
								 k, &x); 
	    
		for (i=1; i<=nbnodes; ++i)
		for (j=1; j<=nbnodes; ++j)
		{
		    node1 = Node(i,comp1);
		    node2 = Node(j,comp2);
		    
		    if ((*pattern)(node1,node2))
		    {
			sum = 0.0;
			for (k=1; k<=dim; ++k) 
			  sum += (*N)(ip,i) * c[k] * dN(j,k);
			
			A(node1,node2) += sum * fact; 
		    }
		}
	    } 
	}
    }
    return True;
}
//-------------------------------------------------------------------------


Bool StdElement:: MCSource(const PATCH& patch, Vector<Real>& b, 
			   const Jacobian& Jac, const Matrix<Bool>* pattern, 
			   int nComp, const Matrix<int>& Node, Real time) const
{
    int i, ip, comp, node; 
    Real fact;
    const int nbnodes = NoOfBaseNodes();

    if (!material->SourceTerm(patch.Class())) return False;
    if (pattern == 0) pattern = symPattern;

    Real Jabs = Abs(Jac.detJ);


    if (material->constSourceTerm(patch.Class()))
    {
	for (comp=1; comp<=nComp; ++comp)
	{
	    fact = material->S(patch.Class(), 0, time) * Jabs;

	    for (i=1; i<=nbnodes; ++i) 
	    {
		node = Node(i,comp);
		if ((*pattern)(node,node))  b[node] += 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

	    for (comp=1; comp<=nComp; ++comp)
	    {
		fact = Jabs * IFMass->Weight(ip) * material->S(patch.Class(),
							       &x, time); 
		for (i=1; i<=nbnodes; ++i) 
		{
		    node = Node(i,comp);
		    if ((*pattern)(node,node)) b[node] += fact * (*NMass)(ip,i);
		}
	    }
	}	
    }
    return True;
}
//-------------------------------------------------------------------------


Bool StdElement:: MCMass(const PATCH& patch, Matrix<Real>& A, 
			 const Jacobian& Jac, const Matrix<Bool>* pattern,
			 int nComp, const Matrix<int>& Node ) const
{
    int i, j, ip, comp1, comp2, node1, node2; 
    Real fact;
    const int nbnodes = NoOfBaseNodes();

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

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

    
    if (material->constMassTerm(patch.Class()))
    {
	for (comp1=1; comp1<=nComp; ++comp1)
	for (comp2=1; comp2<=nComp; ++comp2)
	{
	    if (comp1 != comp2) continue;
	
	    fact = material->M(patch.Class()) * Jabs;

	    for (i=1; i<=nbnodes; ++i)
	    for (j=1; j<=nbnodes; ++j)
	    {
		node1 = Node(i,comp1);
		node2 = Node(j,comp2);
		
		if ((*pattern)(node1,node2)) A(node1,node2) += 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


	    for (comp1=1; comp1<=nComp; ++comp1)
	    for (comp2=1; comp2<=nComp; ++comp2)
	    {
		if (comp1 != comp2) continue;
		
		fact = Jabs * IFMass->Weight(ip) * material->M(patch.Class(),&x);
							   
		for (i=1; i<=nbnodes; ++i)
		for (j=1; j<=nbnodes; ++j)
		{
		    node1 = Node(i,comp1);
		    node2 = Node(j,comp2);
		    
		    if ((*pattern)(node1,node2))
		      A(node1,node2) += fact * (*NMass)(ip,i) * (*NMass)(ip,j);
		}
	    }
	}
    }
    return True;
}
//-------------------------------------------------------------------------


Bool StdElement:: MCLumpedMass(const PATCH& patch, Matrix<Real>& A, 
			       const Jacobian& Jac, const Matrix<Bool>* pattern,
			       int nComp, const Matrix<int>& Node) const
{
    int i, comp1, comp2, node1, node2; 
    Real fact;
    const int nbnodes = NoOfBaseNodes();

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

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

    
    if (material->constMassTerm(patch.Class()))
    {
	for (i=1; i<=nbnodes; ++i)
	{
	    for (comp1=1; comp1<=nComp; ++comp1)
	    for (comp2=1; comp2<=nComp; ++comp2)
	    {
		if (comp1 != comp2) continue;
		
		fact = material->M(patch.Class()) * Jabs;
       
		node1 = Node(i,comp1);
		node2 = Node(i,comp2);

		if ((*pattern)(node1,node2))  
		   A(node1,node1) += fact * IFLumpedMass->Weight(i);
	    }
	}
    }
    else
    {
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_x;
	
	for (i=1; i<=nbnodes; ++i)
	{
	    IFLumpedMass->Coordinates(i,u);
	    patch.realCoordinates(u,x);  
	    
	    for (comp1=1; comp1<=nComp; ++comp1)
	    for (comp2=1; comp2<=nComp; ++comp2)
	    {
		if (comp1 != comp2) continue;

		fact = Jabs  * material->M(patch.Class(),&x);
		
		node1 = Node(i,comp1);
		node2 = Node(i,comp2);
		
		if ((*pattern)(node1,node2))
		  A(node1,node2) += fact * IFLumpedMass->Weight(i);
	    }
	}
    }    
    return True;
}
//-------------------------------------------------------------------------


void StdElement:: MCL2Norm(const PATCH& /*patch*/, Matrix<Real>& A, 
			   const Jacobian& Jac, const Matrix<Bool>* pattern,
			   int nComp, const Matrix<int>& Node) const
{
    int i, j, comp1, comp2, node1, node2; 
    const int nbnodes = NoOfBaseNodes();

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


    for (comp1=1; comp1<=nComp; ++comp1)
    for (comp2=1; comp2<=nComp; ++comp2)
    {
	if (comp1 != comp2) continue;
	
	for (i=1; i<=nbnodes; ++i)
	for (j=1; j<=nbnodes; ++j)
	{
	    node1 = Node(i,comp1);
	    node2 = Node(j,comp2);
		    
	    if ((*pattern)(node1,node2))  A(node1,node2) += Jabs * (*M)(i,j);
	}
    }
}
//-------------------------------------------------------------------------


void StdElement:: MCLNorm(const PATCH& /*patch*/, Vector<Real>& b, 
			  const Jacobian& Jac, const Matrix<Bool>* pattern,
			  int nComp, const Matrix<int>& Node) const
{
    int i, comp1, node1; 
    const int nbnodes = NoOfBaseNodes();

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

    for (comp1=1; comp1<=nComp; ++comp1)
    {
	for (i=1; i<=nbnodes; ++i) 
	{
	    node1 = Node(i,comp1);
	    if ((*pattern)(node1,node1))  b[node1] += Jabs * (*B)[i];
	}
    }
}
//-------------------------------------------------------------------------


Bool StdElement:: MCNeumannBCs(const PATCH& patch, Vector<Real>& b, 
			       const Matrix<Bool>* pattern, int nComp, 
			       const Matrix<int>& Node, Real time) const
{
    int k, comp, 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

	for (comp=1; comp<=nComp; ++comp)
	{
	    if (material->isNeumann(face[k]->isNeumann(), face[k]->Class(),
				    comp))
	    {
		getLocalFaceNodes(k, nodes);
		count += faceElement->MCNeumannBCOnBoundary(*(face[k]), nodes, 
							    b, pattern, comp,
							    Node, time);
	    }
	}
    }
    return count;
}
//-------------------------------------------------------------------------


Bool StdElement:: MCCauchyBCs(const PATCH& patch, Matrix<Real>& A, 
			      Vector<Real>& b, const Matrix<Bool>* pattern, 
			      int nComp, const Matrix<int>& Node, 
			      Real time) const
{
    int  k, comp, 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

	for (comp=1; comp<=nComp; ++comp)
	{
	    if (material->isCauchy(face[k]->isCauchy(), face[k]->Class(), comp))
	    {
		++count;
		getLocalFaceNodes(k, nodes);
		faceElement->MCCauchyBCOnBoundary(*(face[k]), nodes, A, b, 
						  pattern, comp, Node, time);
	    }
	}
    }	
    return count;
}
//-------------------------------------------------------------------------


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


Bool StdElement:: MCNeumannBCOnBoundary(const PATCH& tr, Vector<int>& nodes, 
					Vector<Real>& b,
					const Matrix<Bool>* pattern,
					int comp, const Matrix<int>& Node,
					Real time) const
{
    int  i, ip, node; 
    Real fact;
    const int nbnodes = NoOfBaseNodes();

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

    Real Jabs = tr.volume();				// length of edge
    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<=nbnodes; ++i) 
	{
	    node = Node(nodes[i],comp);
	    if ((*pattern)(node,node)) b[node] += fact * (*B)[i];
	}
    }
    else
    {
	const int nnodes = NoOfNodes();
	const int noIP = IFMass->NoOfIPoints();
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_xBnd;

	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 = Node(nodes[i],comp);
		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:: MCCauchyBCOnBoundary(const PATCH& tr, Vector<int>& nodes, 
				       Matrix<Real>& A, Vector<Real>& b,
				       const Matrix<Bool>* pattern,
				       int comp, const Matrix<int>& Node,
				       Real time)  const
{
    int  i, j, ip, node1, node2; 
    Real fact;
    const int nnodes = NoOfNodes();

    MCNeumannBCOnBoundary(tr, nodes, b, pattern, comp, Node, 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)
	{
	    node1 = Node(nodes[i],comp);

	    for (j=1; j<=nnodes; ++j)
	    {
		node2 = Node(nodes[j],comp);
		
		if ((*pattern)(node1,node2)) A(node1,node2) += fact * (*M)(i,j); 
	    }
	}
    }
    else
    {
	const int noIP = IFMass->NoOfIPoints();
	Vector<Real>& u = *tmp_u;
	Vector<Real>& x = *tmp_xBnd;

	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->Cauchy(tr.Class(),
								&x, time); 
	    for (i=1; i<=nnodes; ++i)
	    {
		node1 = Node(nodes[i],comp);

		for (j=1; j<=nnodes; ++j)
		{
		    node2 = Node(nodes[j],comp);

		    if ((*pattern)(node1,node2)) 
		      A(node1,node2) += fact * (*NMass)(ip,i) * (*NMass)(ip,j);
		}
	    }
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void StdElement:: MCConstEllip1(const PATCH& patch, Matrix<Real>& A, 
				const Jacobian& Jac, 
				const Matrix<Bool>* pattern,
				int nComp, const Matrix<int>& Node) const
{
    int i, j, comp1, comp2, node1, node2; 
    Real fact, Jabs;
    const int nbnodes = NoOfBaseNodes();

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

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

    for (comp1=1; comp1<=nComp; ++comp1)
    for (comp2=1; comp2<=nComp; ++comp2)
    {
	if (comp1 != comp2) continue;

	for (i=1; i<=nbnodes; ++i)
	{
	    node1 = Node(i,comp1);
	    for (j=1; j<=nbnodes; ++j)
	    {
		node2 = Node(j,comp2);
		if ((*pattern)(node1,node2)) A(node1,node2) += fact * (*S1)(i,j);
	    }
	}
    }
}
//-------------------------------------------------------------------------


void StdElement:: MCConstConvec1(const PATCH& patch, Matrix<Real>& A, 
				 const Jacobian& Jac,
				 const Matrix<Bool>* pattern,
				 int nComp, const Matrix<int>& Node) const
{
    int  i, j, comp1, comp2, node1, node2; 
    Real a, c, Jabs;
    const int nbnodes = NoOfBaseNodes();

    Jabs = Abs(Jac.detJ);

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

    for (comp1=1; comp1<=nComp; ++comp1)
    for (comp2=1; comp2<=nComp; ++comp2)
    {
	if (comp1 != comp2) continue;

	for (i=1; i<=nbnodes; ++i)
	{
	    node1 = Node(i,comp1);
	    for (j=1; j<=nbnodes; ++j)
	    {
		node2 = Node(j,comp2);

		if ((*pattern)(node1,node2))  A(node1,node2) += a * (*C1)(i,j);
	    }
	}
    }
}
//-------------------------------------------------------------------------


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

    if (pattern == 0) pattern = symPattern;

    for (k=1; k<=2; k++)
    {
	for (comp=1; comp<=nComp; ++comp)
	{
	    if (material->isNeumann(pt[k]->isNeumann(), pt[k]->Class(), comp))
	    {
		node = Node(pt[k]->getNode(), comp);

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


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

    if (pattern == 0) pattern = symPattern;

    for (k=1; k<=2; k++)
    {
	for (comp=1; comp<=nComp; ++comp)
	{
	    if (material->isCauchy(pt[k]->isCauchy(), pt[k]->Class(), comp))
	    {
		node = Node(pt[k]->getNode(), comp);

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


void StdElement:: MCConstEllip2(const PATCH& patch, Matrix<Real>& A, 
				const Jacobian& Jac, 
				const Matrix<Bool>* pattern,
				int nComp, const Matrix<int>& Node) const
{
    int i, j, comp1, comp2, node1, node2; 
    Real a1, a2, a3, fact, Jabs;
    const int nbnodes = NoOfBaseNodes();

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

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

    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());
    a1 *= fact;
    a2 *= fact;
    a3 *= fact;


    for (comp1=1; comp1<=nComp; ++comp1)
    for (comp2=1; comp2<=nComp; ++comp2)
    {
	if (comp1 != comp2) continue;

	for (i=1; i<=nbnodes; ++i)
	{
	    node1 = Node(i,comp1);
	    for (j=1; j<=nbnodes; ++j)
	    {
		node2 = Node(j,comp2);
		if ((*pattern)(node1,node2))
		 A(node1,node2) += a1*(*S1)(i,j) + a2*(*S2)(i,j) + a3*(*S3)(i,j);
	    }
	}
    }
}
//-------------------------------------------------------------------------


void StdElement:: MCConstConvec2(const PATCH& patch, Matrix<Real>& A, 
				 const Jacobian& Jac,
				 const Matrix<Bool>* pattern,
				 int nComp, const Matrix<int>& Node) const
{
    int  i, j, k, comp1, comp2, node1, node2; 

    const int nbnodes = NoOfBaseNodes();
    const int dim = SpaceDim();

    Real a1, a2, Jabs;
    Vector<Real>& c = *tmp_u;

    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 (comp1=1; comp1<=nComp; ++comp1)
    for (comp2=1; comp2<=nComp; ++comp2)
    {
	if (comp1 != comp2) continue;

	for (i=1; i<=nbnodes; ++i)
	{
	    node1 = Node(i,comp1);
	    for (j=1; j<=nbnodes; ++j)
	    {
		node2 = Node(j,comp2);

		if ((*pattern)(node1,node2))  
		  A(node1,node2) += a1*(*C1)(i,j) + a2*(*C2)(i,j);
	    }
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void StdElement:: MCConstEllip3(const PATCH& patch, Matrix<Real>& A, 
				const Jacobian& Jac, 
				const Matrix<Bool>* pattern,
				int nComp, const Matrix<int>& Node) const
{
    int i, j, comp1, comp2, node1, node2; 
    Real a[7], fact, Jabs;
    const int nbnodes = NoOfBaseNodes();

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

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

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

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

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


    for (comp1=1; comp1<=nComp; ++comp1)
    for (comp2=1; comp2<=nComp; ++comp2)
    {
	if (comp1 != comp2) continue;

	for (i=1; i<= nbnodes; ++i)
	{
	    node1 = Node(i,comp1);
	    for (j=1; j<=nbnodes; ++j)
	    {
		node2 = Node(j,comp2);
		if ((*pattern)(node1,node2))
		   A(node1,node2) += 
		    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:: MCConstConvec3(const PATCH& patch, Matrix<Real>& A, 
				 const Jacobian& Jac,
				 const Matrix<Bool>* pattern,
				 int nComp, const Matrix<int>& Node) const
{
    int  i, j, k, comp1, comp2, node1, node2; 
    Real a1, a2, a3, Jabs;
    Vector<Real> c(SpaceDim());
    const int nbnodes = NoOfBaseNodes();
    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 (comp1=1; comp1<=nComp; ++comp1)
    for (comp2=1; comp2<=nComp; ++comp2)
    {
	if (comp1 != comp2) continue;

	for (i=1; i<=nbnodes; ++i)
	{
	    node1 = Node(i,comp1);
	    for (j=1; j<=nbnodes; ++j)
	    {
		node2 = Node(j,comp2);

		if ((*pattern)(node1,node2))  
		  A(node1,node2) += 
		    a1*(*C1)(i,j) + a2*(*C2)(i,j) + a3*(*C3)(i,j);
	    }
	}
    }
}