*** empty log message ***

Numeric/Iso.cpp

@@ -13,6 +13,25 @@ static void affect(double *xi, double *yi, double *zi, int i,
   zi[i] = zp[j];
 }
 
+double InterpolateIso(double *X, double *Y, double *Z,
+                      double *Val, double V, int I1, int I2,
+                      double *XI, double *YI, double *ZI)
+{
+  if(Val[I1] == Val[I2]) {
+    *XI = X[I1];
+    *YI = Y[I1];
+    *ZI = Z[I1];
+    return 0;
+  }
+  else {
+    double coef = (V - Val[I1]) / (Val[I2] - Val[I1]);
+    *XI = coef * (X[I2] - X[I1]) + X[I1];
+    *YI = coef * (Y[I2] - Y[I1]) + Y[I1];
+    *ZI = coef * (Z[I2] - Z[I1]) + Z[I1];
+    return coef;
+  }
+}
+
 // Compute an iso-point in a line
 
 int IsoLine(double *X, double *Y, double *Z, double *Val, double V,

Numeric/Iso.h

@@ -6,6 +6,10 @@
 #ifndef _ISO_H_
 #define _ISO_H_
 
+double InterpolateIso(double *X, double *Y, double *Z, 
+                      double *Val, double V, int I1, int I2, 
+                      double *XI, double *YI ,double *ZI);
+
 int IsoLine(double *X, double *Y, double *Z, double *Val, double V,
             double *Xp, double *Yp, double *Zp);
 

Numeric/Numeric.cpp

@@ -337,25 +337,6 @@ float char2float(char c)
   else return f;
 }
 
-double InterpolateIso(double *X, double *Y, double *Z,
-                      double *Val, double V, int I1, int I2,
-                      double *XI, double *YI, double *ZI)
-{
-  if(Val[I1] == Val[I2]) {
-    *XI = X[I1];
-    *YI = Y[I1];
-    *ZI = Z[I1];
-    return 0;
-  }
-  else {
-    double coef = (V - Val[I1]) / (Val[I2] - Val[I1]);
-    *XI = coef * (X[I2] - X[I1]) + X[I1];
-    *YI = coef * (Y[I2] - Y[I1]) + Y[I1];
-    *ZI = coef * (Z[I2] - Z[I1]) + Z[I1];
-    return coef;
-  }
-}
-
 void gradSimplex(double *x, double *y, double *z, double *v, double *grad)
 {
   // p = p1 * (1-u-v-w) + p2 u + p3 v + p4 w

Numeric/Numeric.h

@@ -65,9 +65,6 @@ float char2float(char c);
 void eigenvalue(double mat[3][3], double re[3]);
 void FindCubicRoots(const double coeff[4], double re[3], double im[3]);
 void eigsort(double d[3]);
-double InterpolateIso(double *X, double *Y, double *Z, 
-                      double *Val, double V, int I1, int I2, 
-                      double *XI, double *YI ,double *ZI);
 void gradSimplex(double *x, double *y, double *z, double *v, double *grad);
 double ComputeVonMises(double *val);
 double ComputeScalarRep(int numComp, double *val);

Plugin/Levelset.cpp

@@ -6,6 +6,7 @@
 #include "Levelset.h"
 #include "MakeSimplex.h"
 #include "Numeric.h"
+#include "Iso.h"
 #include "adaptiveData.h"
 #include "GmshDefines.h"
 

Solver/elasticitySolver.cpp

@@ -60,7 +60,7 @@ void elasticitySolver::setMesh(const std::string &meshFileName)
 
 void elasticitySolver::readInputFile(const std::string &fn)
 {
-  FILE *f = fopen (fn.c_str(),"r");
+  FILE *f = fopen(fn.c_str(), "r");
   char what[256];
   while(!feof(f)){
     fscanf(f, "%s", what);
@@ -142,7 +142,7 @@ void elasticitySolver::solve()
   for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
        it != nodalDisplacements.end(); ++it){
     MVertex *v = pModel->getMeshVertexByTag(it->first.first);
-    pAssembler-> fixVertex(v , it->first.second, _tag, it->second);    
+    pAssembler->fixVertex(v , it->first.second, _tag, it->second);    
     printf("-- Fixing node %3d comp %1d to %8.5f\n", 
            it->first.first, it->first.second, it->second);
   }
@@ -172,11 +172,11 @@ void elasticitySolver::solve()
     int physical = it->first;
     std::vector<MVertex *> v;     
     pModel->getMeshVerticesForPhysicalGroup(_dim, physical, v);
-    printf("Physical %d, dim: %d, nb vert: %d\n", physical, _dim, v.size());
+    printf("Physical %d, dim: %d, nb vert: %d\n", physical, _dim, (int)v.size());
     for (unsigned int i = 0; i < v.size(); i++){  
-      pAssembler->numberVertex (v[i], 0, _tag);  
-      pAssembler->numberVertex (v[i], 1, _tag);  
-      pAssembler->numberVertex (v[i], 2, _tag);  
+      pAssembler->numberVertex(v[i], 0, _tag);  
+      pAssembler->numberVertex(v[i], 1, _tag);  
+      pAssembler->numberVertex(v[i], 2, _tag);  
     }
   }
   
@@ -189,9 +189,9 @@ void elasticitySolver::solve()
     SVector3 f = it->second;
     std::vector<GEntity*> ent = groups[0][iVertex];
     for (unsigned int i = 0; i < ent.size(); i++){      
-      pAssembler-> assemble(ent[i]->mesh_vertices[0] , 0, _tag, f.x());
-      pAssembler-> assemble(ent[i]->mesh_vertices[0] , 1, _tag, f.y());
-      pAssembler-> assemble(ent[i]->mesh_vertices[0] , 2, _tag, f.z());
+      pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
+      pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
+      pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
       printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n", 
              ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
     }
@@ -236,7 +236,7 @@ void elasticitySolver::solve()
   }
   
   // assembling the stifness matrix
-  for (std::map<int, std::pair<double, double> > :: iterator it = elasticConstants.begin();
+  for (std::map<int, std::pair<double, double> >::iterator it = elasticConstants.begin();
        it != elasticConstants.end() ; it++){
     int physical = it->first;
     double E = it->second.first;

Solver/femTerm.h

@@ -22,6 +22,8 @@ class femTerm {
  protected:
   GModel *_gm;
  public:
+  femTerm(GModel *gm) : _gm(gm) {}
+  virtual ~femTerm (){}
   // return the number of columns of the element matrix
   virtual int sizeOfC(MElement*) const = 0;
   // return the number of rows of the element matrix
@@ -29,27 +31,28 @@ class femTerm {
   // in a given element, return the dof associated to a given row (column)
   // of the local element matrix
   virtual Dof getLocalDofR(MElement *e, int iRow) const = 0;
-  // default behavior : symmetric
+  // default behavior: symmetric
   virtual Dof getLocalDofC(MElement *e, int iCol) const
   { 
     return getLocalDofR(e, iCol); 
   }
- public:
-  femTerm(GModel *gm) : _gm(gm) {}
-  virtual ~femTerm (){}
+  // compute the elementary matrix
   virtual void elementMatrix(MElement *e, fullMatrix<dataMat> &m) const = 0;
   virtual void elementVector(MElement *e, fullVector<dataVec> &m) const 
   {
     m.scale(0.0);
   }
-  void addToMatrix(dofManager<dataVec,dataMat> &dm, GEntity *ge) const
+  // add the contribution from all the elements in the entity ge to
+  // the dof manager
+  void addToMatrix(dofManager<dataVec, dataMat> &dm, GEntity *ge) const
   {
     for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
       MElement *e = ge->getMeshElement(i);
       addToMatrix(dm, e);
     }
   }
-  void addToMatrix(dofManager<dataVec,dataMat> &dm, MElement *e) const
+  // add the contribution from a single element to the dof manager
+  void addToMatrix(dofManager<dataVec, dataMat> &dm, MElement *e) const
   {
     const int nbR = sizeOfR(e);
     const int nbC = sizeOfC(e);
@@ -57,7 +60,7 @@ class femTerm {
     elementMatrix(e, localMatrix);
     addToMatrix(dm, localMatrix, e);
   }
-  void addToMatrix(dofManager<dataVec,dataMat> &dm, 
+  void addToMatrix(dofManager<dataVec, dataMat> &dm, 
                    fullMatrix<dataMat> &localMatrix, 
                    MElement *e) const
   {
@@ -67,7 +70,7 @@ class femTerm {
       Dof R = getLocalDofR(e, j);
       for (int k = 0; k < nbC; k++){
         Dof C = getLocalDofC(e, k);
-        dm.assemble(R,C, localMatrix(j, k));
+        dm.assemble(R, C, localMatrix(j, k));
       }
     }
   }
@@ -83,7 +86,7 @@ class femTerm {
   }
   void neumannNodalBC(int physical, int dim, int comp,int field,
                       const simpleFunction<dataVec> &fct,
-                      dofManager<dataVec,dataMat> &dm)
+                      dofManager<dataVec, dataMat> &dm)
   {
     std::map<int, std::vector<GEntity*> > groups[4];
     GModel *m = _gm;
@@ -120,12 +123,12 @@ class femTerm {
   void addToRightHandSide(dofManager<dataVec,dataMat> &dm, GEntity *ge) const 
   {
     for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
-      MElement *e = ge->getMeshElement (i);
+      MElement *e = ge->getMeshElement(i);
       int nbR = sizeOfR(e);
-      fullVector<dataVec> V (nbR);
-      elementVector (e, V);
+      fullVector<dataVec> V(nbR);
+      elementVector(e, V);
       // assembly
-      for (int j=0;j<nbR;j++)dm.assemble(getLocalDofR(e,j),V(j));
+      for (int j = 0; j < nbR; j++) dm.assemble(getLocalDofR(e, j), V(j));
     }
   }
 };

Solver/helmholtzTerm.h

@@ -23,6 +23,10 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
   const int _iFieldR;
   int _iFieldC ;
  public:
+  helmoltzTerm(GModel *gm, int iFieldR, int iFieldC, simpleFunction<scalar> *k,
+               simpleFunction<scalar> *a) 
+    : femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR), 
+      _iFieldC(iFieldC) {}
   // one dof per vertex (nodal fem)
   virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
   virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
@@ -34,11 +38,6 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
   {
     return Dof(e->getVertex(iRow)->getNum(), _iFieldC);
   }
-  public:
-  helmoltzTerm(GModel *gm, int iFieldR, int iFieldC, 
-               simpleFunction<scalar> *k, 
-               simpleFunction<scalar> *a) : 
-    femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR), _iFieldC(iFieldC){}
   virtual void elementMatrix(MElement *e, fullMatrix<scalar> &m) const
   {
     // compute integration rule