/*
$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);
}
}
}
}
syntax highlighted by Code2HTML, v. 0.9.1