/*
$Id: elements3.cc,v 1.5 1996/11/19 09:52:45 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 "elements3.h"
#include "utils.h"
#include "numerics.h"
#include "triang.h"
#include "materials.h"
#include "integ.h"
//-------------------------------------------------------------------------
// -- linear tetrahedron:
static double TetraSf1(double x, double y, double z) { return 1-x-y-z; }
static double TetraSf2(double x, double /*y*/, double /*z*/) { return x; }
static double TetraSf3(double /*x*/, double y, double /*z*/) { return y; }
static double TetraSf4(double /*x*/, double /*y*/, double z) { return z; }
static double TetraSf1dx(double /*x*/, double /*y*/, double /*z*/) { return -1.; }
static double TetraSf1dy(double /*x*/, double /*y*/, double /*z*/) { return -1.; }
static double TetraSf1dz(double /*x*/, double /*y*/, double /*z*/) { return -1.; }
static double TetraSf2dx(double /*x*/, double /*y*/, double /*z*/) { return 1.; }
static double TetraSf2dy(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double TetraSf2dz(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double TetraSf3dx(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double TetraSf3dy(double /*x*/, double /*y*/, double /*z*/) { return 1.; }
static double TetraSf3dz(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double TetraSf4dx(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double TetraSf4dy(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double TetraSf4dz(double /*x*/, double /*y*/, double /*z*/) { return 1.; }
// -- quadratic Lagrange-tetrahedron:
static double LQuadTetraSf1(double x, double y, double z)
{ return 2.*(1-x-y-z)*(0.5-x-y-z); }
static double LQuadTetraSf2(double x, double /*y*/, double /*z*/) { return x*(2.*x-1); }
static double LQuadTetraSf3(double /*x*/, double y, double /*z*/) { return y*(2.*y-1); }
static double LQuadTetraSf4(double /*x*/, double /*y*/, double z) { return z*(2.*z-1); }
static double LQuadTetraSf1dx(double x, double y, double z)
{ return 4.*(x+y+z)-3.; }
static double LQuadTetraSf1dy(double x, double y, double z)
{ return 4.*(x+y+z)-3.; }
static double LQuadTetraSf1dz(double x, double y, double z)
{ return 4.*(x+y+z)-3.; }
static double LQuadTetraSf2dx(double x, double /*y*/, double /*z*/) { return 4.*x-1.; }
static double LQuadTetraSf2dy(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double LQuadTetraSf2dz(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double LQuadTetraSf3dx(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double LQuadTetraSf3dy(double /*x*/, double y, double /*z*/) { return 4.*y-1.; }
static double LQuadTetraSf3dz(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double LQuadTetraSf4dx(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double LQuadTetraSf4dy(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double LQuadTetraSf4dz(double /*x*/, double /*y*/, double z) { return 4.*z-1.; }
static double LQuadTetraSf5(double x, double y, double z)
{ return 4.0*y*(1.0-x-y-z); }
static double LQuadTetraSf6(double x, double y, double /*z*/) { return 4.0*x*y; }
static double LQuadTetraSf7(double x, double y, double z)
{ return 4.0*x*(1.0-x-y-z); }
static double LQuadTetraSf8(double /*x*/, double y, double z) { return 4.0*y*z; }
static double LQuadTetraSf9(double x, double /*y*/, double z) { return 4.0*x*z; }
static double LQuadTetraSf10(double x, double y, double z)
{ return 4.0*z*(1.0-x-y-z); }
static double LQuadTetraSf5dx (double /*x*/, double y, double /*z*/) { return -4.0*y; }
static double LQuadTetraSf6dx (double /*x*/, double y, double /*z*/) { return 4.0*y; }
static double LQuadTetraSf7dx (double x, double y, double z)
{ return 4.0*(1-2.0*x-y-z); }
static double LQuadTetraSf8dx (double /*x*/, double /*y*/, double /*z*/) { return 0.0; }
static double LQuadTetraSf9dx (double /*x*/, double /*y*/, double z) { return 4.0*z; }
static double LQuadTetraSf10dx(double /*x*/, double /*y*/, double z) { return -4.0*z; }
static double LQuadTetraSf5dy(double x, double y, double z)
{ return 4.0*(1-x-2.0*y-z); }
static double LQuadTetraSf6dy(double x, double /*y*/, double /*z*/) { return 4.0*x; }
static double LQuadTetraSf7dy(double x, double /*y*/, double /*z*/) { return -4.0*x; }
static double LQuadTetraSf8dy(double /*x*/, double /*y*/, double z) { return 4.0*z; }
static double LQuadTetraSf9dy(double /*x*/, double /*y*/, double /*z*/) { return 0.0; }
static double LQuadTetraSf10dy(double /*x*/, double /*y*/, double z) { return -4.0*z; }
static double LQuadTetraSf5dz(double /*x*/, double y, double /*z*/) { return -4.0*y; }
static double LQuadTetraSf6dz(double /*x*/, double /*y*/, double /*z*/) { return 0.0; }
static double LQuadTetraSf7dz(double x, double /*y*/, double /*z*/) { return -4.0*x; }
static double LQuadTetraSf8dz(double /*x*/, double y, double /*z*/) { return 4.0*y; }
static double LQuadTetraSf9dz(double x, double /*y*/, double /*z*/) { return 4.0*x; }
static double LQuadTetraSf10dz(double x, double y, double z)
{ return 4.0*(1-x-y-2.0*z); }
//-------------------------------------------------------------------------
static double Csf1(double x, double y, double z) { return 1-x-y-z; }
static double Csf2(double x, double /*y*/, double /*z*/) { return x; }
static double Csf3(double /*x*/, double y, double /*z*/) { return y; }
static double Csf4(double /*x*/, double /*y*/, double z) { return z; }
static double Csf1dx(double /*x*/, double /*y*/, double /*z*/) { return -1.; }
static double Csf1dy(double /*x*/, double /*y*/, double /*z*/) { return -1.; }
static double Csf1dz(double /*x*/, double /*y*/, double /*z*/) { return -1.; }
static double Csf2dx(double /*x*/, double /*y*/, double /*z*/) { return 1.; }
static double Csf2dy(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double Csf2dz(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double Csf3dx(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double Csf3dy(double /*x*/, double /*y*/, double /*z*/) { return 1.; }
static double Csf3dz(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double Csf4dx(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double Csf4dy(double /*x*/, double /*y*/, double /*z*/) { return 0.; }
static double Csf4dz(double /*x*/, double /*y*/, double /*z*/) { return 1.; }
static double Csf5(double x, double y, double z) { return 4.0*y*(1.0-x-y-z); }
static double Csf6(double x, double y, double /*z*/) { return 4.0*x*y; }
static double Csf7(double x, double y, double z) { return 4.0*x*(1.0-x-y-z); }
static double Csf8(double /*x*/, double y, double z) { return 4.0*y*z; }
static double Csf9(double x, double /*y*/, double z) { return 4.0*x*z; }
static double Csf10(double x, double y, double z) { return 4.0*z*(1.0-x-y-z); }
static double Csf5dx(double /*x*/, double y, double /*z*/) { return -4.0*y; }
static double Csf6dx(double /*x*/, double y, double /*z*/) { return 4.0*y; }
static double Csf7dx(double x, double y, double z) { return 4.0*(1-2.0*x-y-z);}
static double Csf8dx(double /*x*/, double /*y*/, double /*z*/) { return 0.0; }
static double Csf9dx(double /*x*/, double /*y*/, double z) { return 4.0*z; }
static double Csf10dx(double /*x*/, double /*y*/, double z) { return -4.0*z; }
static double Csf5dy(double x, double y, double z) { return 4.0*(1-x-2.0*y-z);}
static double Csf6dy(double x, double /*y*/, double /*z*/) { return 4.0*x; }
static double Csf7dy(double x, double /*y*/, double /*z*/) { return -4.0*x; }
static double Csf8dy(double /*x*/, double /*y*/, double z) { return 4.0*z; }
static double Csf9dy(double /*x*/, double /*y*/, double /*z*/) { return 0.0; }
static double Csf10dy(double /*x*/, double /*y*/, double z) { return -4.0*z; }
static double Csf5dz(double /*x*/, double y, double /*z*/) { return -4.0*y; }
static double Csf6dz(double /*x*/, double /*y*/, double /*z*/) { return 0.0; }
static double Csf7dz(double x, double /*y*/, double /*z*/) { return -4.0*x; }
static double Csf8dz(double /*x*/, double y, double /*z*/) { return 4.0*y; }
static double Csf9dz(double x, double /*y*/, double /*z*/) { return 4.0*x; }
static double Csf10dz(double x, double y, double z) { return 4.0*(1-x-y-2.0*z);}
static double Csf11(double x, double y, double z) { return 27.0*x*y*z; }
static double Csf12(double x, double y, double z) { return 27.0*(1-x-y-z)*y*z; }
static double Csf13(double x, double y, double z) { return 27.0*(1-x-y-z)*x*z; }
static double Csf14(double x, double y, double z) { return 27.0*(1-x-y-z)*x*y; }
static double Csf11dx(double /*x*/, double y, double z) { return 27.0*y*z; }
static double Csf12dx(double /*x*/, double y, double z) { return -27.0*y*z; }
static double Csf13dx(double x, double y, double z) { return 27.0*(1.0-2.0*x-y-z)*z; }
static double Csf14dx(double x, double y, double z) { return 27.0*(1.0-2.0*x-y-z)*y; }
static double Csf11dy(double x, double /*y*/, double z) { return 27.0*x*z; }
static double Csf12dy(double x, double y, double z) { return 27.0*(1.0-x-2.0*y-z)*z; }
static double Csf13dy(double x, double /*y*/, double z) { return -27.0*x*z; }
static double Csf14dy(double x, double y, double z) { return 27.0*(1.0-x-2.0*y-z)*x; }
static double Csf11dz(double x, double y, double /*z*/) { return 27.0*x*y; }
static double Csf12dz(double x, double y, double z) { return 27.0*(1.0-x-y-2.0*z)*y; }
static double Csf13dz(double x, double y, double z) { return 27.0*(1.0-x-y-2.0*z)*x; }
static double Csf14dz(double x, double y, double /*z*/) { return -27.0*x*y; }
//----------------------------- Tetrahedra ------------------------------
Tetra:: Tetra(Material* material0) : StdElement()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void Tetra:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,1); // GaussTetraOrder1;
IFMass = new GaussForm(dim,3); // GaussTetraOrder3;
IFLumpedMass = new LumpedForm(dim,1,"LumpedTetraOrder1");
InfoIntegFormulas("Tetra");
edgeElement = new Line(material);
faceElement = new Triangle(material);
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] = TetraSf1;
sf[2] = TetraSf2;
sf[3] = TetraSf3;
sf[4] = TetraSf4;
dsf[1][1] = TetraSf1dx; dsf[1][2] = TetraSf1dy; dsf[1][3] = TetraSf1dz;
dsf[2][1] = TetraSf2dx; dsf[2][2] = TetraSf2dy; dsf[2][3] = TetraSf2dz;
dsf[3][1] = TetraSf3dx; dsf[3][2] = TetraSf3dy; dsf[3][3] = TetraSf3dz;
dsf[4][1] = TetraSf4dx; dsf[4][2] = TetraSf4dy; dsf[4][3] = TetraSf4dz;
}
//-------------------------------------------------------------------------
// the basic matrices for constant coefficients (S1, S2, S3, M, B ...)
// are chosen due to Schwarz's conventions:
void Tetra:: 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);
S4 = new Matrix<double>(nbnodes,nbnodes);
S5 = new Matrix<double>(nbnodes,nbnodes);
S6 = new Matrix<double>(nbnodes,nbnodes);
C1 = new Matrix<double>(nbnodes,nbnodes);
C2 = new Matrix<double>(nbnodes,nbnodes);
C3 = 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) = (*S4)(i,k) = (*S5)(i,k) =
(*S6)(i,k) =
(*C1)(i,k) = (*C2)(i,k) = (*C3)(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,1)*(*dNU)(ip,k,3) +
(*dNU)(ip,i,3)*(*dNU)(ip,k,1));
(*S4)(i,k) += IF->Weight(ip) * (*dNU)(ip,i,2)*(*dNU)(ip,k,2);
(*S5)(i,k) += IF->Weight(ip) *((*dNU)(ip,i,2)*(*dNU)(ip,k,3) +
(*dNU)(ip,i,3)*(*dNU)(ip,k,2));
(*S6)(i,k) += IF->Weight(ip) * (*dNU)(ip,i,3)*(*dNU)(ip,k,3);
}
// 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);
(*S4)(k,i) = (*S4)(i,k);
(*S5)(k,i) = (*S5)(i,k);
(*S6)(k,i) = (*S6)(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);
(*C3)(i,k) += IF->Weight(ip) * (*N)(ip,i)*(*dNU)(ip,k,3);
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
HQuadTetra:: HQuadTetra(Material* material0) : Tetra()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void HQuadTetra:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,3);
IFMass = new GaussForm(dim,5);
IFLumpedMass = new LumpedForm(dim,2,"LumpedTetraOrder2");
InfoIntegFormulas("HQuadTetra");
edgeElement = new HQuadLine(material);
faceElement = new HQuadTriangle(material);
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] = TetraSf1;
HQsf[2] = TetraSf2;
HQsf[3] = TetraSf3;
HQsf[4] = TetraSf4;
HQsf[5] = LQuadTetraSf5;
HQsf[6] = LQuadTetraSf6;
HQsf[7] = LQuadTetraSf7;
HQsf[8] = LQuadTetraSf8;
HQsf[9] = LQuadTetraSf9;
HQsf[10]= LQuadTetraSf10;
dHQsf[1][1] = TetraSf1dx;
dHQsf[1][2] = TetraSf1dy;
dHQsf[1][3] = TetraSf1dz;
dHQsf[2][1] = TetraSf2dx;
dHQsf[2][2] = TetraSf2dy;
dHQsf[2][3] = TetraSf2dz;
dHQsf[3][1] = TetraSf3dx;
dHQsf[3][2] = TetraSf3dy;
dHQsf[3][3] = TetraSf3dz;
dHQsf[4][1] = TetraSf4dx;
dHQsf[4][2] = TetraSf4dy;
dHQsf[4][3] = TetraSf4dz;
dHQsf[5][1] = LQuadTetraSf5dx;
dHQsf[5][2] = LQuadTetraSf5dy;
dHQsf[5][3] = LQuadTetraSf5dz;
dHQsf[6][1] = LQuadTetraSf6dx;
dHQsf[6][2] = LQuadTetraSf6dy;
dHQsf[6][3] = LQuadTetraSf6dz;
dHQsf[7][1] = LQuadTetraSf7dx;
dHQsf[7][2] = LQuadTetraSf7dy;
dHQsf[7][3] = LQuadTetraSf7dz;
dHQsf[8][1] = LQuadTetraSf8dx;
dHQsf[8][2] = LQuadTetraSf8dy;
dHQsf[8][3] = LQuadTetraSf8dz;
dHQsf[9][1] = LQuadTetraSf9dx;
dHQsf[9][2] = LQuadTetraSf9dy;
dHQsf[9][3] = LQuadTetraSf9dz;
dHQsf[10][1]= LQuadTetraSf10dx;
dHQsf[10][2]= LQuadTetraSf10dy;
dHQsf[10][3]= LQuadTetraSf10dz;
}
//-------------------------------------------------------------------------
void HQuadTetra:: 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(*DLYPattern,k) (*DLYPattern)(i,k) = False;
for (i=5; i<=nnodes; ++i)
{
(*DLYPattern)(i,i) = True;
for (k=1; k<=4; ++k) (*DLYPattern)(i,k) = True;
}
}
//-------------------------------------------------------------------------
LQuadTetra:: LQuadTetra(Material* material0) : HQuadTetra()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void LQuadTetra:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,3);
IFMass = new GaussForm(dim,5);
IFLumpedMass = new LumpedForm(dim,2,"LumpedTetraOrder2");
InfoIntegFormulas("LQuadTetra");
edgeElement = new LQuadLine(material);
faceElement = new LQuadTriangle(material);
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] = LQuadTetraSf1;
LQsf[2] = LQuadTetraSf2;
LQsf[3] = LQuadTetraSf3;
LQsf[4] = LQuadTetraSf4;
LQsf[5] = LQuadTetraSf5;
LQsf[6] = LQuadTetraSf6;
dLQsf[1][1] = LQuadTetraSf1dx; dLQsf[1][2] = LQuadTetraSf1dy;
dLQsf[2][1] = LQuadTetraSf2dx; dLQsf[2][2] = LQuadTetraSf2dy;
dLQsf[3][1] = LQuadTetraSf3dx; dLQsf[3][2] = LQuadTetraSf3dy;
dLQsf[4][1] = LQuadTetraSf4dx; dLQsf[4][2] = LQuadTetraSf4dy;
dLQsf[5][1] = LQuadTetraSf5dx; dLQsf[5][2] = LQuadTetraSf5dy;
dLQsf[6][1] = LQuadTetraSf6dx; dLQsf[6][2] = LQuadTetraSf6dy;
LQsf[1] = LQuadTetraSf1;
LQsf[2] = LQuadTetraSf2;
LQsf[3] = LQuadTetraSf3;
LQsf[4] = LQuadTetraSf4;
LQsf[5] = LQuadTetraSf5;
LQsf[6] = LQuadTetraSf6;
LQsf[7] = LQuadTetraSf7;
LQsf[8] = LQuadTetraSf8;
LQsf[9] = LQuadTetraSf9;
LQsf[10]= LQuadTetraSf10;
dLQsf[1][1] = LQuadTetraSf1dx;
dLQsf[1][2] = LQuadTetraSf1dy;
dLQsf[1][3] = LQuadTetraSf1dz;
dLQsf[2][1] = LQuadTetraSf2dx;
dLQsf[2][2] = LQuadTetraSf2dy;
dLQsf[2][3] = LQuadTetraSf2dz;
dLQsf[3][1] = LQuadTetraSf3dx;
dLQsf[3][2] = LQuadTetraSf3dy;
dLQsf[3][3] = LQuadTetraSf3dz;
dLQsf[4][1] = LQuadTetraSf4dx;
dLQsf[4][2] = LQuadTetraSf4dy;
dLQsf[4][3] = LQuadTetraSf4dz;
dLQsf[5][1] = LQuadTetraSf5dx;
dLQsf[5][2] = LQuadTetraSf5dy;
dLQsf[5][3] = LQuadTetraSf5dz;
dLQsf[6][1] = LQuadTetraSf6dx;
dLQsf[6][2] = LQuadTetraSf6dy;
dLQsf[6][3] = LQuadTetraSf6dz;
dLQsf[7][1] = LQuadTetraSf7dx;
dLQsf[7][2] = LQuadTetraSf7dy;
dLQsf[7][3] = LQuadTetraSf7dz;
dLQsf[8][1] = LQuadTetraSf8dx;
dLQsf[8][2] = LQuadTetraSf8dy;
dLQsf[8][3] = LQuadTetraSf8dz;
dLQsf[9][1] = LQuadTetraSf9dx;
dLQsf[9][2] = LQuadTetraSf9dy;
dLQsf[9][3] = LQuadTetraSf9dz;
dLQsf[10][1]= LQuadTetraSf10dx;
dLQsf[10][2]= LQuadTetraSf10dy;
dLQsf[10][3]= LQuadTetraSf10dz;
}
//-------------------------------------------------------------------------
void Tetra:: 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] = 1; nodes[2] = 4; nodes[3] = 3; }
else if (faceNo == 3) { nodes[1] = 1; nodes[2] = 2; nodes[3] = 4; }
else if (faceNo == 4) { nodes[1] = 1; nodes[2] = 3; nodes[3] = 2; }
else
{
cout << "\n*** Tetra:: getLocalFaceNodes: wrong face number: "
<< faceNo << "\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void HQuadTetra:: getLocalFaceNodes(int faceNo, Vector<int>& nodes) const
{
if (faceNo == 1) { nodes[1] = 2; nodes[2] = 3; nodes[3] = 4;
nodes[4] = 8; nodes[5] = 9; nodes[6] = 6;
}
else if (faceNo == 2) { nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
nodes[4] = 8; nodes[5] = 5; nodes[6] = 10;
}
else if (faceNo == 3) { nodes[1] = 1; nodes[2] = 2; nodes[3] = 4;
nodes[4] = 9; nodes[5] = 10;nodes[6] = 7;
}
else if (faceNo == 4) { nodes[1] = 1; nodes[2] = 3; nodes[3] = 2;
nodes[4] = 6; nodes[5] = 7; nodes[6] = 5;
}
else
{
cout << "\n*** HQuadTetra:: getLocalFaceNodes: wrong face number: "
<< faceNo << "\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool HQuadTetraFast:: 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 p, di, sum1, sum2, sum3, sum4;
Real nx, ny, nz;
if (!material->EllipticTerm(patch.Class())) return False;
if (pattern == 0) pattern = symPattern;
if (material->constEllipticTerm(patch.Class()))
{
assembleConstEllip(patch, A, Jac, pattern);
}
else
{
// we have to use the integration formulae GaussTetraOrder2OnPoints
Matrix<Real>& dN = *tmp_dN;
Vector<Real>& u = *tmp_u;
Vector<Real>& x = *tmp_x;
Vector<Real>& matE = *tmp_matE;
Jabs = Abs(Jac.detJ);
const Matrix<Real>& J = Jac.Jinv;
if (material->isAnisotropic())
{
cout << "\n*** HQuadTetraFast:: assembleEllip: "
"anisotropic material not yet implemented \n";
cout.flush(); abort();
}
else if (!material->isIsotropic())
{
cout << "\n*** HQuadTetraFast:: assembleEllip: no materialtype specified \n";
cout.flush(); abort();
}
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;
}
}
}
// constant integration on faces
const Vector<PATCH*>& tri = patch.getFaces();
di = Jabs;
for (k=1; k<=NoOfFaces(); k++)
{
if ( !(tri[k]->isDirichlet()) )
{
switch (k)
{
case 1:
p=matE[1]; // value on face 1
nx = J(1,1) + J(1,2) + J(1,3);
ny = J(2,1) + J(2,2) + J(2,3);
nz = J(3,1) + J(3,2) + J(3,3);
sum1 = -di*p*( nx*nx + ny*ny + nz*nz )/6.0;
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(6,1) += sum1;
A(6,2) += sum2;
A(6,3) += sum3;
A(6,4) += sum4;
A(8,1) += sum1;
A(8,2) += sum2;
A(8,3) += sum3;
A(8,4) += sum4;
A(9,1) += sum1;
A(9,2) += sum2;
A(9,3) += sum3;
A(9,4) += sum4;
break;
case 2:
p=matE[2]; // value on face 2
nx = -J(1,1);
ny = -J(2,1);
nz = -J(3,1);
sum1 = -di*p/6.0*
( nx*(J(1,1)+J(1,2)+J(1,3) )
+ny*(J(2,1)+J(2,2)+J(2,3) )
+nz*(J(3,1)+J(3,2)+J(3,3) )
);
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(5,1) += sum1;
A(5,2) += sum2;
A(5,3) += sum3;
A(5,4) += sum4;
A(8,1) += sum1;
A(8,2) += sum2;
A(8,3) += sum3;
A(8,4) += sum4;
A(10,1) += sum1;
A(10,2) += sum2;
A(10,3) += sum3;
A(10,4) += sum4;
break;
case 3:
p=matE[3]; // value on face 3
nx = -J(1,2);
ny = -J(2,2);
nz = -J(3,2);
sum1 = -di*p/6.0*
( nx*(J(1,1)+J(1,2)+J(1,3) )
+ny*(J(2,1)+J(2,2)+J(2,3) )
+nz*(J(3,1)+J(3,2)+J(3,3) )
);
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(7,1) += sum1;
A(7,2) += sum2;
A(7,3) += sum3;
A(7,4) += sum4;
A(9,1) += sum1;
A(9,2) += sum2;
A(9,3) += sum3;
A(9,4) += sum4;
A(10,1) += sum1;
A(10,2) += sum2;
A(10,3) += sum3;
A(10,4) += sum4;
break;
case 4:
p=matE[4]; // value on face 4
nx = -J(1,3);
ny = -J(2,3);
nz = -J(3,3);
sum1 = -di*p/6.0*
( nx*(J(1,1)+J(1,2)+J(1,3) )
+ny*(J(2,1)+J(2,2)+J(2,3) )
+nz*(J(3,1)+J(3,2)+J(3,3) )
);
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(5,1) += sum1;
A(5,2) += sum2;
A(5,3) += sum3;
A(5,4) += sum4;
A(6,1) += sum1;
A(6,2) += sum2;
A(6,3) += sum3;
A(6,4) += sum4;
A(7,1) += sum1;
A(7,2) += sum2;
A(7,3) += sum3;
A(7,4) += sum4;
break;
}
}
}
}
return True;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
EFTetraFast:: EFTetraFast(Material* material0) : Tetra()
{
material = material0;
basicInit();
setShapeFunctionsAtIPs();
compConstantMatrices();
setPatterns();
}
//-------------------------------------------------------------------------
void EFTetraFast:: basicInit()
{
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
IF = new GaussForm(dim,5);
IFMass = new GaussForm(dim,7);
IFLumpedMass = new LumpedForm(dim,2,"LumpedTetraOrder2");
InfoIntegFormulas("EFTetraFast");
edgeElement = new HQuadLine(material);
faceElement = new HQuadTriangle(material);
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);
Csf[1] = Csf1;
Csf[2] = Csf2;
Csf[3] = Csf3;
Csf[4] = Csf4;
Csf[5] = Csf5;
Csf[6] = Csf6;
Csf[7] = Csf7;
Csf[8] = Csf8;
Csf[9] = Csf9;
Csf[10]= Csf10;
Csf[11]= Csf11;
Csf[12]= Csf12;
Csf[13]= Csf13;
Csf[14]= Csf14;
dCsf[1][1] = Csf1dx; dCsf[1][2] = Csf1dy; dCsf[1][3] = Csf1dz;
dCsf[2][1] = Csf2dx; dCsf[2][2] = Csf2dy; dCsf[2][3] = Csf2dz;
dCsf[3][1] = Csf3dx; dCsf[3][2] = Csf3dy; dCsf[3][3] = Csf3dz;
dCsf[4][1] = Csf4dx; dCsf[4][2] = Csf4dy; dCsf[4][3] = Csf4dz;
dCsf[5][1] = Csf5dx; dCsf[5][2] = Csf5dy; dCsf[5][3] = Csf5dz;
dCsf[6][1] = Csf6dx; dCsf[6][2] = Csf6dy; dCsf[6][3] = Csf6dz;
dCsf[7][1] = Csf7dx; dCsf[7][2] = Csf7dy; dCsf[7][3] = Csf7dz;
dCsf[8][1] = Csf8dx; dCsf[8][2] = Csf8dy; dCsf[8][3] = Csf8dz;
dCsf[9][1] = Csf9dx; dCsf[9][2] = Csf9dy; dCsf[9][3] = Csf9dz;
dCsf[10][1]= Csf10dx; dCsf[10][2]= Csf10dy; dCsf[10][3]= Csf10dz;
dCsf[11][1]= Csf11dx; dCsf[11][2]= Csf11dy; dCsf[11][3]= Csf11dz;
dCsf[12][1]= Csf12dx; dCsf[12][2]= Csf12dy; dCsf[12][3]= Csf12dz;
dCsf[13][1]= Csf13dx; dCsf[13][2]= Csf13dy; dCsf[13][3]= Csf13dz;
dCsf[14][1]= Csf14dx; dCsf[14][2]= Csf14dy; dCsf[14][3]= Csf14dz;
}
//-------------------------------------------------------------------------
Bool EFTetraFast:: assembleEllip(const PATCH& patch, Matrix<Real>& A,
const Jacobian& Jac,
const Matrix<Bool>* pattern) const
{
int i, k, ip, noIP=IF->NoOfIPoints();
const int nnodes = NoOfNodes();
const int dim = SpaceDim();
Real fact, Jabs, sum;
Real p, di, sum1, sum2, sum3, sum4;
Real nx, ny, nz;
if (!material->EllipticTerm(patch.Class())) return False;
if (pattern == 0) pattern = symPattern;
if (material->constEllipticTerm(patch.Class()))
{
assembleConstEllip(patch, A, Jac, pattern);
}
else
{
// we have to use the integration formulae GaussTetraOrder2OnPoints
Matrix<Real> dN(nnodes,dim);
Vector<Real> u(dim), x(dim);
Vector<Real> matE(IF->NoOfIPoints());
Jabs = Abs(Jac.detJ);
const Matrix<Real>& J = Jac.Jinv;
if (material->isAnisotropic())
{
cout << "\n*** EFTetraFast:: assembleEllip: "
"anisotropic material not yet implemented \n";
cout.flush(); abort();
}
else if (!material->isIsotropic())
{
cout << "\n*** EFTetraFast:: assembleEllip: no materialtype specified " << endl;
abort();
}
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;
}
}
}
// quadratic and cubic bumbs: constant integration on faces
const Vector<PATCH*>& tri = patch.getFaces();
di = Jabs;
Real iConst = 9.0/40.0;
for (k=1; k<=NoOfFaces(); k++)
{
if ( !(tri[k]->isDirichlet()) )
{
switch (k)
{
case 1:
// quadratic bumbs
p=matE[1]; // value on face 1
nx = J(1,1) + J(1,2) + J(1,3);
ny = J(2,1) + J(2,2) + J(2,3);
nz = J(3,1) + J(3,2) + J(3,3);
sum1 = -di*p*( nx*nx + ny*ny + nz*nz )/6.0;
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(6,1) += sum1;
A(6,2) += sum2;
A(6,3) += sum3;
A(6,4) += sum4;
A(8,1) += sum1;
A(8,2) += sum2;
A(8,3) += sum3;
A(8,4) += sum4;
A(9,1) += sum1;
A(9,2) += sum2;
A(9,3) += sum3;
A(9,4) += sum4;
// cubic bumb
sum1 = -iConst*di*p*(nx*nx + ny*ny + nz*nz);
sum2 = iConst*di*p*(nx*J(1,1) + ny*J(2,1) + nz*J(3,1));
sum3 = iConst*di*p*(nx*J(1,2) + ny*J(2,2) + nz*J(3,2));
sum4 = iConst*di*p*(nx*J(1,3) + ny*J(2,3) + nz*J(3,3));
A(11,1) += sum1;
A(11,2) += sum2;
A(11,3) += sum3;
A(11,4) += sum4;
break;
case 2:
// quadratic bumbs
p=matE[2]; // value on face 2
nx = -J(1,1);
ny = -J(2,1);
nz = -J(3,1);
sum1 = -di*p/6.0*
( nx*(J(1,1)+J(1,2)+J(1,3))
+ny*(J(2,1)+J(2,2)+J(2,3))
+nz*(J(3,1)+J(3,2)+J(3,3))
);
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(5,1) += sum1;
A(5,2) += sum2;
A(5,3) += sum3;
A(5,4) += sum4;
A(8,1) += sum1;
A(8,2) += sum2;
A(8,3) += sum3;
A(8,4) += sum4;
A(10,1) += sum1;
A(10,2) += sum2;
A(10,3) += sum3;
A(10,4) += sum4;
// cubic bumb
sum1 = -iConst*di*p*
( nx*(J(1,1)+J(1,2)+J(1,3))
+ny*(J(2,1)+J(2,2)+J(2,3))
+nz*(J(3,1)+J(3,2)+J(3,3))
);
sum2 = iConst*di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) );
sum3 = iConst*di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) );
sum4 = iConst*di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) );
A(12,1) += sum1;
A(12,2) += sum2;
A(12,3) += sum3;
A(12,4) += sum4;
break;
case 3:
// quadratic bumbs
p=matE[3]; // value on face 3
nx = -J(1,2);
ny = -J(2,2);
nz = -J(3,2);
sum1 = -di*p/6.0*
( nx*(J(1,1)+J(1,2)+J(1,3) )
+ny*(J(2,1)+J(2,2)+J(2,3) )
+nz*(J(3,1)+J(3,2)+J(3,3) )
);
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(7,1) += sum1;
A(7,2) += sum2;
A(7,3) += sum3;
A(7,4) += sum4;
A(9,1) += sum1;
A(9,2) += sum2;
A(9,3) += sum3;
A(9,4) += sum4;
A(10,1) += sum1;
A(10,2) += sum2;
A(10,3) += sum3;
A(10,4) += sum4;
// cubic bumb
sum1 = -iConst*di*p*
( nx*(J(1,1)+J(1,2)+J(1,3))
+ny*(J(2,1)+J(2,2)+J(2,3))
+nz*(J(3,1)+J(3,2)+J(3,3))
);
sum2 = iConst*di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) );
sum3 = iConst*di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) );
sum4 = iConst*di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) );
A(13,1) += sum1;
A(13,2) += sum2;
A(13,3) += sum3;
A(13,4) += sum4;
break;
case 4:
// quadratic bumbs
p=matE[4]; // value on face 4
nx = -J(1,3);
ny = -J(2,3);
nz = -J(3,3);
sum1 = -di*p/6.0*
( nx*(J(1,1)+J(1,2)+J(1,3) )
+ny*(J(2,1)+J(2,2)+J(2,3) )
+nz*(J(3,1)+J(3,2)+J(3,3) )
);
sum2 = di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) )/6.0;
sum3 = di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) )/6.0;
sum4 = di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) )/6.0;
A(5,1) += sum1;
A(5,2) += sum2;
A(5,3) += sum3;
A(5,4) += sum4;
A(6,1) += sum1;
A(6,2) += sum2;
A(6,3) += sum3;
A(6,4) += sum4;
A(7,1) += sum1;
A(7,2) += sum2;
A(7,3) += sum3;
A(7,4) += sum4;
// cubic bumb
sum1 = -iConst*di*p*
( nx*(J(1,1)+J(1,2)+J(1,3))
+ny*(J(2,1)+J(2,2)+J(2,3))
+nz*(J(3,1)+J(3,2)+J(3,3))
);
sum2 = iConst*di*p*( nx*J(1,1) + ny*J(2,1) + nz*J(3,1) );
sum3 = iConst*di*p*( nx*J(1,2) + ny*J(2,2) + nz*J(3,2) );
sum4 = iConst*di*p*( nx*J(1,3) + ny*J(2,3) + nz*J(3,3) );
A(14,1) += sum1;
A(14,2) += sum2;
A(14,3) += sum3;
A(14,4) += sum4;
break;
}
}
}
}
return True;
}
//-------------------------------------------------------------------------
void EFTetraFast:: getLocalFaceNodes(int faceNo, Vector<int>& nodes) const
{
if (faceNo == 1) { nodes[1] = 2; nodes[2] = 3; nodes[3] = 4;
nodes[4] = 8; nodes[5] = 9; nodes[6] = 6;
nodes[7] = 11;
}
else if (faceNo == 2) { nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
nodes[4] = 8; nodes[5] = 5; nodes[6] = 10;
nodes[7] = 12;
}
else if (faceNo == 3) { nodes[1] = 1; nodes[2] = 2; nodes[3] = 4;
nodes[4] = 9; nodes[5] = 10;nodes[6] = 7;
nodes[7] = 13;
}
else if (faceNo == 4) { nodes[1] = 1; nodes[2] = 3; nodes[3] = 2;
nodes[4] = 6; nodes[5] = 7; nodes[6] = 5;
nodes[7] = 14;
}
else
{
cout << "\n*** CubTetraFast:: getLocalFaceNodes: wrong face number: "
<< faceNo << "\n";
cout.flush(); abort();
}
}
syntax highlighted by Code2HTML, v. 0.9.1