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