/*
$Id: elements2.cc,v 1.5 1996/11/19 09:51:31 bzferdma Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "elements1.h"
#include "elements2.h"
#include "numerics.h"
#include "triang.h"
#include "triangA.h"
#include "materials.h"
#include "integ.h"
//-------------------------------------------------------------------------
static double TriangleSf1(double x, double y) { return 1-x-y; }
static double TriangleSf2(double x, double /*y*/) { return x; }
static double TriangleSf3(double /*x*/, double y) { return y; }
static double TriangleSf1dx(double /*x*/, double /*y*/) { return -1.; }
static double TriangleSf1dy(double /*x*/, double /*y*/) { return -1.; }
static double TriangleSf2dx(double /*x*/, double /*y*/) { return 1.; }
static double TriangleSf2dy(double /*x*/, double /*y*/) { return 0.; }
static double TriangleSf3dx(double /*x*/, double /*y*/) { return 0.; }
static double TriangleSf3dy(double /*x*/, double /*y*/) { return 1.; }
static double HQuadTriangleSf1(double x, double y) { return 1-x-y; }
static double HQuadTriangleSf2(double x, double /*y*/) { return x; }
static double HQuadTriangleSf3(double /*x*/, double y) { return y; }
static double HQuadTriangleSf1dx(double /*x*/, double /*y*/) { return -1.; }
static double HQuadTriangleSf1dy(double /*x*/, double /*y*/) { return -1.; }
static double HQuadTriangleSf2dx(double /*x*/, double /*y*/) { return 1.; }
static double HQuadTriangleSf2dy(double /*x*/, double /*y*/) { return 0.; }
static double HQuadTriangleSf3dx(double /*x*/, double /*y*/) { return 0.; }
static double HQuadTriangleSf3dy(double /*x*/, double /*y*/) { return 1.; }
static double HQuadTriangleSf4 (double x, double y) { return 4.*x*y; }
static double HQuadTriangleSf4dx(double /*x*/, double y) { return 4.*y; }
static double HQuadTriangleSf4dy(double x, double /*y*/) { return 4.*x; }
static double HQuadTriangleSf5 (double x, double y) { return 4.*y*(1-x-y); }
static double HQuadTriangleSf5dx(double /*x*/, double y) { return -4.*y; }
static double HQuadTriangleSf5dy(double x, double y) { return 4.*(1-x-2*y); }
static double HQuadTriangleSf6 (double x, double y) { return 4.*x*(1-x-y); }
static double HQuadTriangleSf6dx(double x, double y) { return 4.*(1-2*x-y); }
static double HQuadTriangleSf6dy(double x, double /*y*/) { return -4.*x; }
static double LQuadTriangleSf1(double x, double y)
{ return 2.*(1.-x-y)*(0.5-x-y); }
static double LQuadTriangleSf1dx(double x, double y) { return 4.*(x+y)-3.; }
static double LQuadTriangleSf1dy(double x, double y) { return 4.*(x+y)-3.; }
static double LQuadTriangleSf2(double x, double /*y*/) { return x*(2.*x-1.); }
static double LQuadTriangleSf2dx(double x, double /*y*/) { return 4.*x-1.; }
static double LQuadTriangleSf2dy(double /*x*/, double /*y*/) { return 0.; }
static double LQuadTriangleSf3(double /*x*/, double y) { return y*(2.*y-1.); }
static double LQuadTriangleSf3dx(double /*x*/, double /*y*/) { return 0.; }
static double LQuadTriangleSf3dy(double /*x*/, double y) { return 4.*y-1.; }
static double LQuadTriangleSf4(double x, double y) { return 4.*x*y; }
static double LQuadTriangleSf4dx(double /*x*/, double y) { return 4.*y; }
static double LQuadTriangleSf4dy(double x, double /*y*/) { return 4.*x; }
static double LQuadTriangleSf5(double x, double y) { return 4.*y*(1.-x-y); }
static double LQuadTriangleSf5dx(double /*x*/, double y) { return -4.*y; }
static double LQuadTriangleSf5dy(double x, double y) { return 4.*(1-x-2.*y); }
static double LQuadTriangleSf6(double x, double y) { return 4.*x*(1.-x-y); }
static double LQuadTriangleSf6dx(double x, double y) { return 4.*(1-2.*x-y); }
static double LQuadTriangleSf6dy(double x, double /*y*/) { return -4.*x; }
//----------------------------- Triangles ------------------------------
Triangle:: Triangle(Material* material0) : StdElement()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void Triangle:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,1); // GaussTriOrder1;
IFMass = new GaussForm(dim,3); // GaussTriOrder3;
IFLumpedMass = new LumpedForm(dim,1,"LumpedTriOrder1");
InfoIntegFormulas("Triangle");
edgeElement = new Line(material);
faceElement = edgeElement;
tmp_nodes = new Vector<int>(faceElement->NoOfNodes()+1);
tmp_Ell = new Matrix<Real>(dim,dim);
tmp_dN = new Matrix<Real>(nnodes,dim);
tmp_u = new Vector<Real>(dim);
tmp_x = new Vector<Real>(dim);
tmp_xBnd = new Vector<Real>(dim+1);
sf[1] = TriangleSf1;
sf[2] = TriangleSf2;
sf[3] = TriangleSf3;
dsf[1][1] = TriangleSf1dx; dsf[1][2] = TriangleSf1dy;
dsf[2][1] = TriangleSf2dx; dsf[2][2] = TriangleSf2dy;
dsf[3][1] = TriangleSf3dx; dsf[3][2] = TriangleSf3dy;
}
//-------------------------------------------------------------------------
HQuadTriangleFast:: HQuadTriangleFast(Material* material0) : HQuadTriangle()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void HQuadTriangleFast:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,2, "GaussTriOrder2OnEdges");
IFMass = new GaussForm(dim,5);
IFLumpedMass = new LumpedForm(dim,2,"LumpedTriOrder2");
InfoIntegFormulas("HQuadTriangleFast");
edgeElement = new HQuadLine(material);
faceElement = edgeElement;
tmp_nodes = new Vector<int>(faceElement->NoOfNodes()+1);
tmp_Ell = new Matrix<Real>(dim,dim);
tmp_dN = new Matrix<Real>(nnodes,dim);
tmp_u = new Vector<Real>(dim);
tmp_x = new Vector<Real>(dim);
tmp_xBnd = new Vector<Real>(dim+1);
tmp_matE = new Vector<Real>(IF->NoOfIPoints());
HQsf[1] = HQuadTriangleSf1;
HQsf[2] = HQuadTriangleSf2;
HQsf[3] = HQuadTriangleSf3;
HQsf[4] = HQuadTriangleSf4;
HQsf[5] = HQuadTriangleSf5;
HQsf[6] = HQuadTriangleSf6;
dHQsf[1][1] = HQuadTriangleSf1dx; dHQsf[1][2] = HQuadTriangleSf1dy;
dHQsf[2][1] = HQuadTriangleSf2dx; dHQsf[2][2] = HQuadTriangleSf2dy;
dHQsf[3][1] = HQuadTriangleSf3dx; dHQsf[3][2] = HQuadTriangleSf3dy;
dHQsf[4][1] = HQuadTriangleSf4dx; dHQsf[4][2] = HQuadTriangleSf4dy;
dHQsf[5][1] = HQuadTriangleSf5dx; dHQsf[5][2] = HQuadTriangleSf5dy;
dHQsf[6][1] = HQuadTriangleSf6dx; dHQsf[6][2] = HQuadTriangleSf6dy;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// the basic matrices for constant coefficients (S1, S2, S3, M, B ...)
// are chosen due to Schwarz's conventions:
void Triangle:: compConstantMatrices()
{
int i,k,ip, noIP;
const int nbnodes = NoOfBaseNodes();
S1 = new Matrix<double>(nbnodes,nbnodes);
S2 = new Matrix<double>(nbnodes,nbnodes);
S3 = new Matrix<double>(nbnodes,nbnodes);
C1 = new Matrix<double>(nbnodes,nbnodes);
C2 = new Matrix<double>(nbnodes,nbnodes);
M = new Matrix<double>(nbnodes,nbnodes);
B = new Vector<double>(nbnodes);
for (i=1; i<=nbnodes; ++i)
{
for (k=1; k<=nbnodes; ++k)
(*S1)(i,k) = (*S2)(i,k) = (*S3)(i,k) =
(*C1)(i,k) = (*C2)(i,k) = (*M)(i,k) = 0.0;
(*B)[i] = 0.0;
}
// partial stiffness matrices:
noIP = IF->NoOfIPoints();
for (ip=1; ip<=noIP; ++ip)
for (i=1; i<=nbnodes; ++i)
for (k=1; k<=i; ++k)
{
(*S1)(i,k) += IF->Weight(ip) * (*dNU)(ip,i,1)*(*dNU)(ip,k,1);
(*S2)(i,k) += IF->Weight(ip) *((*dNU)(ip,i,1)*(*dNU)(ip,k,2) +
(*dNU)(ip,i,2)*(*dNU)(ip,k,1));
(*S3)(i,k) += IF->Weight(ip) * (*dNU)(ip,i,2)*(*dNU)(ip,k,2);
}
// mass matrix and source vector:
noIP = IFMass->NoOfIPoints();
for (ip=1; ip<=IFMass->NoOfIPoints(); ++ip)
for (i=1; i<=nbnodes; ++i)
{
(*B)[i] += IFMass->Weight(ip) * (*NMass)(ip,i);
for (k=1; k<=i; ++k) (*M)(i,k) += IFMass->Weight(ip) *
(*NMass)(ip,i) * (*NMass)(ip,k);
}
for (i=1; i<=nbnodes; ++i) // upper triangles
for (k=1; k<i; ++k)
{
(*S1)(k,i) = (*S1)(i,k);
(*S2)(k,i) = (*S2)(i,k);
(*S3)(k,i) = (*S3)(i,k);
(*M)(k,i) = (*M)(i,k);
}
// partial convection matrices:
noIP = IF->NoOfIPoints();
for (ip=1; ip<=noIP; ++ip)
for (i=1; i<=nbnodes; ++i)
for (k=1; k<=nbnodes; ++k)
{
(*C1)(i,k) += IF->Weight(ip) * (*N)(ip,i)*(*dNU)(ip,k,1);
(*C2)(i,k) += IF->Weight(ip) * (*N)(ip,i)*(*dNU)(ip,k,2);
}
}
//-------------------------------------------------------------------------
Bool HQuadTriangleFast:: assembleEllip(const PATCH& patch, Matrix<Real>& A,
const Jacobian& Jac,
const Matrix<Bool>* pattern) const
{
int i, k, ip, noIP;
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
Real fact, Jabs, sum;
Real sumA, sumB, d, dx, dy, n1, n2, p;
if (!material->EllipticTerm(patch.Class())) return False;
if (pattern == 0) pattern = symPattern;
if (material->constEllipticTerm(patch.Class())) // ident. to HQuadTriangle
{
assembleConstEllip(patch, A, Jac, pattern);
}
else
{
// we have to use the integration formulae GaussTriOrder2OnEdges
Matrix<Real>& dN = *tmp_dN;
Vector<Real>& u = *tmp_u;
Vector<Real>& x = *tmp_x;
Vector<Real>& matE = *tmp_matE;
if (material->isAnisotropic())
{
cout << "\n*** HQuadTriangleFast:: assembleConstEllip: "
"anisotropic material not yet implemented \n";
cout.flush(); abort();
}
else if (!material->isIsotropic())
{
cout << "\n*** HQuadTriangleFast:: assembleConstEllip: no materialtype specified \n";
cout.flush(); abort();
}
Jabs = Abs(Jac.detJ);
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
matE[ip] = material->E(patch.Class(), &x);
fact = Jabs * IF->Weight(ip) * matE[ip];
for (i=1; i<=nnodes; ++i)
{
if ((*pattern)(i,i))
{
sum = 0.0;
for (k=1; k<=dim; ++k) sum += dN(i,k) * dN(i,k);
A(i,i) += sum * fact;
}
}
}
// -- Line integration for normal components:
const Matrix<Real>& JJ = Jac.J;
const Vector<EDG*>& eds = patch.getEdges();
d = 2.0/(3.0*Jac.detJ); // ( JJ(1,1)*JJ(2,2) - JJ(1,2)*JJ(2,1) );
if ( !(eds[1]->isDirichlet()) )
{
dx = JJ(2,1) - JJ(1,1);
dy = JJ(2,2) - JJ(1,2);
p = matE[1]; // value on edge 1
n1 = dy*p;
n2 = -dx*p;
sumA = n1*JJ(2,2)-n2*JJ(2,1);
sumB = n2*JJ(1,1)-n1*JJ(1,2);
A(4,1) += (-d)*(sumA+sumB);
A(4,2) += d*sumA;
A(4,3) += d*sumB;
}
if ( !(eds[2]->isDirichlet()) )
{
dx = -JJ(2,1);
dy = -JJ(2,2);
p = matE[2]; // value on edge 2
n1 = dy*p;
n2 = -dx*p;
sumA = n1*JJ(2,2)-n2*JJ(2,1);
sumB = n2*JJ(1,1)-n1*JJ(1,2);
A(5,1) += (-d)*(sumA+sumB);
A(5,2) += d*sumA;
A(5,3) += d*sumB;
}
if ( !(eds[3]->isDirichlet()) )
{
dx = JJ(1,1);
dy = JJ(1,2);
p = matE[3]; // value on edge 3
n1 = dy*p;
n2 = -dx*p;
sumA = n1*JJ(2,2)-n2*JJ(2,1);
sumB = n2*JJ(1,1)-n1*JJ(1,2);
A(6,1) += (-d)*(sumA+sumB);
A(6,2) += d*sumA;
A(6,3) += d*sumB;
}
}
return True;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
HQuadTriangle:: HQuadTriangle(Material* material0) : Triangle()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void HQuadTriangle:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,3);
IFMass = new GaussForm(dim,5);
IFLumpedMass = new LumpedForm(dim,2,"LumpedTriOrder2");
InfoIntegFormulas("HQuadTriangle");
edgeElement = new HQuadLine(material);
faceElement = edgeElement;
tmp_nodes = new Vector<int>(faceElement->NoOfNodes()+1);
tmp_Ell = new Matrix<Real>(dim,dim);
tmp_dN = new Matrix<Real>(nnodes,dim);
tmp_u = new Vector<Real>(dim);
tmp_x = new Vector<Real>(dim);
tmp_xBnd = new Vector<Real>(dim+1);
HQsf[1] = HQuadTriangleSf1;
HQsf[2] = HQuadTriangleSf2;
HQsf[3] = HQuadTriangleSf3;
HQsf[4] = HQuadTriangleSf4;
HQsf[5] = HQuadTriangleSf5;
HQsf[6] = HQuadTriangleSf6;
dHQsf[1][1] = HQuadTriangleSf1dx; dHQsf[1][2] = HQuadTriangleSf1dy;
dHQsf[2][1] = HQuadTriangleSf2dx; dHQsf[2][2] = HQuadTriangleSf2dy;
dHQsf[3][1] = HQuadTriangleSf3dx; dHQsf[3][2] = HQuadTriangleSf3dy;
dHQsf[4][1] = HQuadTriangleSf4dx; dHQsf[4][2] = HQuadTriangleSf4dy;
dHQsf[5][1] = HQuadTriangleSf5dx; dHQsf[5][2] = HQuadTriangleSf5dy;
dHQsf[6][1] = HQuadTriangleSf6dx; dHQsf[6][2] = HQuadTriangleSf6dy;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
LQuadTriangle:: LQuadTriangle(Material* material0) : HQuadTriangle()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void LQuadTriangle:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,3);
IFMass = new GaussForm(dim,5);
IFLumpedMass = new LumpedForm(dim,2,"LumpedTriOrder2");
InfoIntegFormulas("LQuadTriangle");
edgeElement = new LQuadLine(material);
faceElement = edgeElement;
tmp_nodes = new Vector<int>(faceElement->NoOfNodes()+1);
tmp_Ell = new Matrix<Real>(dim,dim);
tmp_dN = new Matrix<Real>(nnodes,dim);
tmp_u = new Vector<Real>(dim);
tmp_x = new Vector<Real>(dim);
tmp_xBnd = new Vector<Real>(dim+1);
LQsf[1] = LQuadTriangleSf1;
LQsf[2] = LQuadTriangleSf2;
LQsf[3] = LQuadTriangleSf3;
LQsf[4] = LQuadTriangleSf4;
LQsf[5] = LQuadTriangleSf5;
LQsf[6] = LQuadTriangleSf6;
dLQsf[1][1] = LQuadTriangleSf1dx; dLQsf[1][2] = LQuadTriangleSf1dy;
dLQsf[2][1] = LQuadTriangleSf2dx; dLQsf[2][2] = LQuadTriangleSf2dy;
dLQsf[3][1] = LQuadTriangleSf3dx; dLQsf[3][2] = LQuadTriangleSf3dy;
dLQsf[4][1] = LQuadTriangleSf4dx; dLQsf[4][2] = LQuadTriangleSf4dy;
dLQsf[5][1] = LQuadTriangleSf5dx; dLQsf[5][2] = LQuadTriangleSf5dy;
dLQsf[6][1] = LQuadTriangleSf6dx; dLQsf[6][2] = LQuadTriangleSf6dy;
}
//-------------------------------------------------------------------------
void Triangle:: getLocalFaceNodes(int faceNo, Vector<int>& nodes) const
{
if (faceNo == 1) { nodes[1] = 2; nodes[2] = 3; }
else if (faceNo == 2) { nodes[1] = 3; nodes[2] = 1; }
else if (faceNo == 3) { nodes[1] = 1; nodes[2] = 2; }
else
{
cout << "\n*** Triangle:: getLocalFaceNodes: wrong face number: "
<< faceNo << "\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void HQuadTriangle:: getLocalFaceNodes(int faceNo, Vector<int>& nodes) const
{
if (faceNo == 1) { nodes[1] = 2; nodes[2] = 3; nodes[3] = 4; }
else if (faceNo == 2) { nodes[1] = 3; nodes[2] = 1; nodes[3] = 5; }
else if (faceNo == 3) { nodes[1] = 1; nodes[2] = 2; nodes[3] = 6; }
else
{
cout << "\n*** HQuadTriangle:: getLocalFaceNodes: wrong face number: "
<< faceNo << "\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void HQuadTriangle:: setPatterns()
{
int i, k;
const int nnodes = NoOfNodes();
StdElement::setPatterns();
// -- set the assembly pattern for DLY:
DLYPattern = new Matrix<int>(nnodes,nnodes);
FORALL_ROWS(*DLYPattern,i)
FORALL_COLUMNS(*symPattern,k) (*DLYPattern)(i,k) = False;
for (i=4; i<=nnodes; ++i)
{
(*DLYPattern)(i,i) = True;
for (k=1; k<=3; ++k) (*DLYPattern)(i,k) = True;
}
}
syntax highlighted by Code2HTML, v. 0.9.1