*** empty log message ***

Solver/convexCombinationTerm.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _CONVEX_COMBINATION_TERM_H_
+#define _CONVEX_COMBINATION_TERM_H_
+
+#include "femTerm.h"
+#include "simpleFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "SElement.h"
+#include "fullMatrix.h"
+
+class convexCombinationTerm : public femTerm<double, double> {
+ protected:
+  const simpleFunction<double> *_diffusivity;
+  const int _iField;
+ public:
+  convexCombinationTerm(GModel *gm, int iField, simpleFunction<double> *diffusivity)
+    : femTerm<double, double>(gm), _iField(iField), _diffusivity(diffusivity) {}
+  virtual int sizeOfR(SElement *se) const
+  {
+    return se->getMeshElement()->getNumVertices();
+  }
+  virtual int sizeOfC(SElement *se) const 
+  { 
+    return se->getMeshElement()->getNumVertices();
+  }
+  void getLocalDofR(SElement *se, int iRow) const
+  {
+    return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iField);
+  }
+  virtual void elementMatrix(SElement *e, fullMatrix<double> &m) const
+  {
+    m.set_all(0.);
+    int nbNodes = e->getMeshElement()->getNumVertices();
+    //double x = 0.3*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x());
+    //double y = 0.3*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y());
+    //double z = 0.3*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z());
+    const double _diff = 1.0; 
+    for (int j = 0; j < nbNodes; j++){
+      for (int k = 0; k < nbNodes; k++){
+        m(j,k) = -1.*_diff;
+      }
+      m(j,j) = (nbNodes - 1) * _diff;
+    }
+  }
+};
+
+#endif

Solver/crossConfTerm.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _CROSS_CONF_TERM_H_
+#define _CROSS_CONF_TERM_H_
+
+#include "femTerm.h"
+#include "simpleFunction.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "SElement.h"
+#include "fullMatrix.h"
+#include "Numeric.h"
+
+class crossConfTerm : public femTerm<double, double> {
+ protected:
+  const simpleFunction<double> *_diffusivity;
+  const int _iFieldR;
+  int _iFieldC;
+ public:
+  crossConfTerm(GModel *gm, int iFieldR, int iFieldC, 
+                simpleFunction<double> *diffusivity)
+    : femTerm<double, double>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC), 
+      _diffusivity(diffusivity) {}
+  virtual int sizeOfR(SElement *se) const 
+  {
+    return se->getMeshElement()->getNumVertices(); 
+  }
+  virtual int sizeOfC(SElement *se) const
+  {
+    return se->getMeshElement()->getNumVertices(); 
+  }
+  Dof getLocalDofR(SElement *se, int iRow) const
+  {
+    return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldR);
+  }
+  void getLocalDofC(SElement *se, int iRow) const
+  {
+    return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldC);
+  }
+  virtual void elementMatrix(SElement *e, fullMatrix<double> &m) const
+  {
+    MElement *e = se->getMeshElement();
+    int nbNodes = e->getNumVertices();
+    int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
+    int npts;
+    IntPt *GP;
+    double jac[3][3];
+    double invjac[3][3];
+    SVector3 Grads [256];
+    double grads[256][3];
+    e->getIntegrationPoints(integrationOrder, &npts, &GP);
+  
+    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];
+      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] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] + 
+                            invjac[0][2] * grads[j][2],
+                            invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] + 
+                            invjac[1][2] * grads[j][2],
+                            invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+                            invjac[2][2] * grads[j][2]);
+      }
+      SVector3 N (jac[2][0], jac[2][1], jac[2][2]);
+      for (int j = 0; j < nbNodes; j++){
+        for (int k = 0; k <= j; k++){
+          m(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ * _diff;
+        }
+      }
+    }
+    for (int j = 0; j < nbNodes; j++)
+      for (int k = 0; k < j; k++)
+        m(k, j) = -1.*m(j, k);    
+  }
+};
+
+#endif

Solver/helmholtzTerm.h

@@ -25,8 +25,8 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
  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) {}
+    : femTerm<scalar, scalar>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
+      _k(k), _a(a) {}
   // one dof per vertex (nodal fem)
   virtual int sizeOfR(SElement *se) const 
   {