/* $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(faceElement->NoOfNodes()+1); tmp_Ell = new Matrix(dim,dim); tmp_dN = new Matrix(nnodes,dim); tmp_u = new Vector(dim); tmp_x = new Vector(dim); tmp_xBnd = new Vector(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(nbnodes,nbnodes); S2 = new Matrix(nbnodes,nbnodes); S3 = new Matrix(nbnodes,nbnodes); S4 = new Matrix(nbnodes,nbnodes); S5 = new Matrix(nbnodes,nbnodes); S6 = new Matrix(nbnodes,nbnodes); C1 = new Matrix(nbnodes,nbnodes); C2 = new Matrix(nbnodes,nbnodes); C3 = new Matrix(nbnodes,nbnodes); M = new Matrix(nbnodes,nbnodes); B = new Vector(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; kNoOfIPoints(); 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(faceElement->NoOfNodes()+1); tmp_Ell = new Matrix(dim,dim); tmp_dN = new Matrix(nnodes,dim); tmp_u = new Vector(dim); tmp_x = new Vector(dim); tmp_xBnd = new Vector(dim+1); tmp_matE = new Vector(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(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(faceElement->NoOfNodes()+1); tmp_Ell = new Matrix(dim,dim); tmp_dN = new Matrix(nnodes,dim); tmp_u = new Vector(dim); tmp_x = new Vector(dim); tmp_xBnd = new Vector(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& 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& 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& A, const Jacobian& Jac, const Matrix* 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& dN = *tmp_dN; Vector& u = *tmp_u; Vector& x = *tmp_x; Vector& matE = *tmp_matE; Jabs = Abs(Jac.detJ); const Matrix& 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& 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(faceElement->NoOfNodes()+1); tmp_Ell = new Matrix(dim,dim); tmp_dN = new Matrix(nnodes,dim); tmp_u = new Vector(dim); tmp_x = new Vector(dim); tmp_xBnd = new Vector(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& A, const Jacobian& Jac, const Matrix* 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 dN(nnodes,dim); Vector u(dim), x(dim); Vector matE(IF->NoOfIPoints()); Jabs = Abs(Jac.detJ); const Matrix& 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& 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& 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(); } }