*** empty log message ***

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

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.);
   } 
 }

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);
 } 
-

Numeric/gmshLaplace.h

@@ -32,5 +32,4 @@ class gmshLaplaceTerm : public gmshNodalFemTerm {
   virtual void elementMatrix(MElement *e, Double_Matrix &m) const;
 };
 
-
 #endif