diff --git a/Numeric/Makefile b/Numeric/Makefile
index a1135dba541f9dc0ace51b8b1455322f31541f9d..6e6bd9c8fa2bb2ae60314b5e06cde1eca892e7a1 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -84,7 +84,8 @@ gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
+ Numeric.h NumericEmbedded.h
gmshElasticity${OBJEXT}: gmshElasticity.cpp gmshElasticity.h \
gmshTermOfFormulation.h GmshMatrix.h ../Common/GmshMessage.h \
../Common/Gmsh.h ../Common/GmshMessage.h ../Geo/GModel.h \
@@ -98,7 +99,8 @@ gmshElasticity${OBJEXT}: gmshElasticity.cpp gmshElasticity.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
+ Numeric.h NumericEmbedded.h
EigSolve${OBJEXT}: EigSolve.cpp
FunctionSpace${OBJEXT}: FunctionSpace.cpp FunctionSpace.h GmshMatrix.h \
../Common/GmshMessage.h ../Common/GmshDefines.h
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index c5a253350533ad35ab93031e1f552abee0e00b82..43f439ceb32448b9a26e653960c7778793d7ebae 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -4,21 +4,20 @@
// bugs and problems to <gmsh@geuz.org>.
#include "gmshElasticity.h"
-
-extern double inv3x3 (double a[3][3], double b[3][3]);
+#include "Numeric.h"
void gmshElasticityTerm::elementMatrix(MElement *e, Double_Matrix & m) const
{
int nbNodes = e->getNumVertices();
- int integrationOrder = 2*(e->getPolynomialOrder()-1);
+ int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
int npts;
IntPt *GP;
- double jac [3][3];
- double invjac [3][3];
- double Grads[256][3],grads[100][3];
+ double jac[3][3];
+ double invjac[3][3];
+ double Grads[256][3], grads[100][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
- m.set_all(0.0);
+ m.set_all(0.);
double FACT = _E / (1 + _nu);
double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
@@ -32,43 +31,46 @@ void gmshElasticityTerm::elementMatrix(MElement *e, Double_Matrix & m) const
{ 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };
- Double_Matrix H(6,6);
- Double_Matrix B (6,3*nbNodes);
- Double_Matrix BTH (3*nbNodes,6);
- Double_Matrix BT (3*nbNodes,6);
+ Double_Matrix H(6, 6);
+ Double_Matrix B(6, 3 * nbNodes);
+ Double_Matrix BTH(3 * nbNodes, 6);
+ Double_Matrix BT(3 * nbNodes, 6);
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
H(i, j) = C[i][j];
for (int i = 0; i < npts; i++){
- const double u = GP[i].pt[0];
- const double v = GP[i].pt[1];
- const double w = GP[i].pt[2];
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
const double weight = GP[i].weight;
const double detJ = e->getJacobian(u, v, w, jac);
e->getGradShapeFunctions(u, v, w, grads);
inv3x3(jac, invjac);
- B.set_all(0.0);
- BT.set_all(0.0);
+ B.set_all(0.);
+ BT.set_all(0.);
for (int j = 0; j < nbNodes; j++){
- Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] + invjac[0][2] * grads[j][2];
- Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] + invjac[1][2] * grads[j][2];
- Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] + invjac[2][2] * grads[j][2];
+ Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ invjac[0][2] * grads[j][2];
+ Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][2] * grads[j][2];
+ Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ invjac[2][2] * grads[j][2];
- BT(j,0) = B(0,j) = Grads[j][0];
- BT(j,3) = B(3,j) = Grads[j][1];
- BT(j,4) = B(4,j) = Grads[j][2];
+ BT(j, 0) = B(0, j) = Grads[j][0];
+ BT(j, 3) = B(3, j) = Grads[j][1];
+ BT(j, 4) = B(4, j) = Grads[j][2];
- BT(j+nbNodes,1) = B(1,j+nbNodes) = Grads[j][1];
- BT(j+nbNodes,3) = B(3,j+nbNodes) = Grads[j][0];
- BT(j+nbNodes,5) = B(5,j+nbNodes) = Grads[j][2];
+ BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
+ BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
+ BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
- BT(j+2*nbNodes,2) = B(2,j+2*nbNodes) = Grads[j][2];
- BT(j+2*nbNodes,4) = B(4,j+2*nbNodes) = Grads[j][0];
- BT(j+2*nbNodes,5) = B(5,j+2*nbNodes) = Grads[j][1];
+ BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
+ BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
+ BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
}
- BTH.set_all(0.0);
+ BTH.set_all(0.);
BTH.blas_dgemm(BT, H);
- m.blas_dgemm(BTH, B, weight * detJ, 1.0);
+ m.blas_dgemm(BTH, B, weight * detJ, 1.);
}
}
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index 1410d342654a555d378c8e8c901e0ecb4ade3318..176cab13aad6b77dcd05599ef11c726c414fb61e 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -4,49 +4,49 @@
// bugs and problems to <gmsh@geuz.org>.
#include "gmshLaplace.h"
-
-extern double inv3x3 (double a[3][3], double b[3][3]);
+#include "Numeric.h"
void gmshLaplaceTerm::elementMatrix(MElement *e, Double_Matrix & m) const
{
int nbNodes = e->getNumVertices();
- int integrationOrder = 2*(e->getPolynomialOrder()-1);
+ int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
int npts;
IntPt *GP;
- double jac [3][3];
- double invjac [3][3];
- double Grads[256][3],grads[256][3];
+ double jac[3][3];
+ double invjac[3][3];
+ double Grads[256][3], grads[256][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
- m.set_all(0.0);
+ m.set_all(0.);
- for (int i=0;i<npts;i++){
- const double u = GP[i].pt[0];
- const double v = GP[i].pt[1];
- const double w = GP[i].pt[2];
+ for (int i = 0; i < npts; i++){
+ const double u = GP[i].pt[0];
+ const double v = GP[i].pt[1];
+ const double w = GP[i].pt[2];
const double weight = GP[i].weight;
- const double detJ = e->getJacobian (u,v,w,jac);
- SPoint3 p; e->pnt(u,v,w,p);
- const double _diff = (*_diffusivity)(p.x(),p.y(),p.z());
- inv3x3 ( jac, invjac) ;
- e->getGradShapeFunctions (u,v,w,grads);
- for (int j=0;j<nbNodes;j++){
- Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] + invjac[0][2] * grads[j][2];
- Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] + invjac[1][2] * grads[j][2];
- Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] + invjac[2][2] * grads[j][2];
+ const double detJ = e->getJacobian(u, v, w, jac);
+ SPoint3 p; e->pnt(u, v, w, p);
+ const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
+ inv3x3(jac, invjac) ;
+ e->getGradShapeFunctions(u, v, w, grads);
+ for (int j = 0; j < nbNodes; j++){
+ Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ invjac[0][2] * grads[j][2];
+ Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][2] * grads[j][2];
+ Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ invjac[2][2] * grads[j][2];
}
- for (int j=0;j<nbNodes;j++){
- for (int k=0;k<=j;k++){
- m(j,k) += (Grads[j][0]*Grads[k][0]+
- Grads[j][1]*Grads[k][1]+
- Grads[j][2]*Grads[k][2]) * weight * detJ * _diff*0.5;
+ for (int j = 0; j < nbNodes; j++){
+ for (int k = 0; k <= j; k++){
+ m(j, k) += (Grads[j][0] * Grads[k][0] +
+ Grads[j][1] * Grads[k][1] +
+ Grads[j][2] * Grads[k][2]) * weight * detJ * _diff * 0.5;
}
}
}
- for (int j=0;j<nbNodes;j++)
- for (int k=0;k<j;k++)
- m(k,j) = m(j,k);
-
+ for (int j = 0; j < nbNodes; j++)
+ for (int k = 0; k < j; k++)
+ m(k, j) = m(j, k);
}
-
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index f43a722bd0c6ebc4aa13294c968de5e2bf0e2ccb..f184c6725b8904ab602daed12fabf2b65eaa51d3 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -32,5 +32,4 @@ class gmshLaplaceTerm : public gmshNodalFemTerm {
virtual void elementMatrix(MElement *e, Double_Matrix &m) const;
};
-
#endif