lua bindings for elasticity

Common/LuaBindings.cpp

@@ -40,6 +40,10 @@
 #include "linearSystemCSR.h"
 #endif
 
+#if defined(HAVE_POST)
+#include "PView.h"
+#endif
+
 extern "C" {
 #include "lua.h"
 #include "lualib.h"
@@ -412,6 +416,9 @@ binding::binding()
   linearSystemCSRGmm<double>::registerBindings(this);
   elasticitySolverRegisterBindings(this); 
 #endif
+#if defined(HAVE_POST)
+  PView::registerBindings(this);
+#endif
 }
 
 binding *binding::_instance=NULL;

Common/LuaBindings.h

@@ -17,6 +17,7 @@
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "Bindings.h"
+#include "SVector3.h"
 
 #if defined(HAVE_LUA)
 
@@ -108,6 +109,13 @@ class luaStack<int>{
   static std::string getName(){ return "int"; }
 };
 
+template<>
+class luaStack<bool>{
+ public:
+  static int get(lua_State *L, int ia){ return lua_toboolean(L, ia); }
+  static void push(lua_State *L, bool i){ lua_pushboolean(L, i); }
+  static std::string getName(){ return "bool"; }
+};
 template<>
 class luaStack<unsigned int>{
  public:
@@ -178,6 +186,7 @@ class luaStack<std::list<lType > >{
   }
 };
 
+
 template<class type>
 class luaStack<std::vector<type> &>{
  public:
@@ -245,6 +254,36 @@ class luaStack<double>{
   static std::string getName(){ return "double"; }
 };
 
+template <>
+class luaStack<SVector3>{
+ public:
+  static SVector3 get(lua_State *L, int ia)
+  {
+    double v[3];
+    size_t size = lua_objlen(L, ia);
+    if (size!=3)
+      luaL_typerror(L, ia, "SVector3");
+    for(size_t i = 0; i< size; i++){
+      lua_rawgeti(L, ia, i + 1);
+      v[i] = luaStack<double>::get(L, -1);
+      lua_pop(L, 1);
+    }
+    return SVector3(v[0],v[1],v[2]);
+  }
+  static void push(lua_State *L, const SVector3& v)
+  {
+    lua_createtable(L, 3, 0);
+    for(size_t i = 0; i < 3; i++){
+      luaStack<double>::push(L, v[i]);
+      lua_rawseti(L, 2, i + 1);
+    }
+  }
+  static std::string getName()
+  {
+    return "3-uple of double";
+  }
+};
+
 template<>
 class luaStack<std::string>{
   public:

Numeric/simpleFunction.h

@@ -24,6 +24,29 @@ class simpleFunction {
   { dfdx = dfdy = dfdz = 0.0; }*/
 };
 
+#include "GmshConfig.h"
+#ifdef HAVE_LUA
+#include "LuaBindings.h"
+template <class scalar>
+class simpleFunctionLua:public simpleFunction<scalar> {
+  lua_State *_L;
+  std::string _luaFunctionName;
+  public:
+  scalar operator () (double x, double y, double z) const {
+    lua_getfield(_L, LUA_GLOBALSINDEX, _luaFunctionName.c_str());
+    luaStack<double>::push(_L, x);
+    luaStack<double>::push(_L, y);
+    luaStack<double>::push(_L, z);
+    lua_call(_L, 3, 1);
+    return luaStack<scalar>::get(_L,-1);
+  }
+  simpleFunctionLua (lua_State *L, const std::string luaFunctionName, scalar s):simpleFunction<scalar>(s) {
+    _L = L;
+    _luaFunctionName = luaFunctionName;
+  }
+};
+#endif
+
 template <class scalar>
 class simpleFunctionOnElement : public simpleFunction<scalar>
 {

Post/PView.cpp

@@ -311,3 +311,12 @@ PView *PView::getViewByNum(int num, int timeStep, int partition)
   return 0;
 }
 
+#include "Bindings.h"
+void PView::registerBindings(binding *b) {
+  classBinding *cb = b->addClass<PView>("PView");
+  cb->setDescription("A post-processing view");
+  methodBinding *cm;
+  cm = cb->addMethod("write",&PView::write);
+  cm->setArgNames("fileName","format","append",NULL);
+  cm->setDescription("write data to a file. Format can be : 0 for ascii pos file, 1 for binary pos file, 2 for parsed pos file, 3 for STL, 4 for TXT, 5 for MSH and 6 for MED files. 'append' option is only supported for pos format.");
+}

Post/PView.h

@@ -19,6 +19,7 @@ class GModel;
 class GMSH_PostPlugin;
 class ConnectionManager;
 
+class binding;
 // A post-processing view.
 class PView{
  private:
@@ -124,6 +125,7 @@ class PView{
 
   // smoothed normals
   smooth_normals *normals;
+  static void registerBindings(binding *b);
 };
 
 // this is the maximum number of nodes of elements we actually *draw*

Solver/elasticitySolver.cpp

@@ -24,6 +24,7 @@
 #include "PViewData.h"
 #endif
 
+/********************* deprecated api ? ****************************************/
 static void printLinearSystem(linearSystemCSRTaucs<double> *lsys)
 {
   int *startIndex;
@@ -65,13 +66,50 @@ void elasticitySolver::setMesh(const std::string &meshFileName)
   LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1);
 }
 
+void elasticitySolver::solve()
+{
+  linearSystemFull<double> *lsys = new linearSystemFull<double>;
+  /*
+#if defined(HAVE_TAUCS)
+  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+  linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#else
+  linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+  lsys->setNoisy(2);
+#endif
+  */
+  assemble(lsys);
+  lsys->systemSolve();
+  printf("-- done solving!\n");
+  GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+
+//   printLinearSystem(lsys);
+
+  double energ=0;
+  for (unsigned int i = 0; i < elasticFields.size(); i++)
+  {
+    SolverField<SVector3> Field(pAssembler, LagSpace);
+    IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
+  }
+  printf("elastic energy=%f\n",energ);
+}
+
 void elasticitySolver::readInputFile(const std::string &fn)
 {
   FILE *f = fopen(fn.c_str(), "r");
   char what[256];
   while(!feof(f)){
     if(fscanf(f, "%s", what) != 1) return;
-    if (!strcmp(what, "ElasticDomain")){
+    if(what[0]=='#'){
+      char *line=NULL;
+      size_t l = 0;
+      int r = getline(&line,&l,f);
+      free(line);
+    }
+    else if (!strcmp(what, "ElasticDomain")){
       elasticField field;
       int physical;
       if(fscanf(f, "%d %lf %lf", &physical, &field._E, &field._nu) != 3) return;
@@ -79,7 +117,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       field.g = new groupOfElements (_dim, physical);
       elasticFields.push_back(field);
     }
-    if (!strcmp(what, "LagrangeMultipliers")){
+    else if (!strcmp(what, "LagrangeMultipliers")){
       LagrangeMultiplierField field;
       int physical;
       double d1, d2, d3, val;
@@ -106,7 +144,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       if(fscanf(f, "%d %d %lf", &node, &comp, &val) != 3) return;
       dirichletBC diri;
       diri.g = new groupOfElements (0, node);
-      diri._f= simpleFunction<double>(val);
+      diri._f= new simpleFunction<double>(val);
       diri._comp=comp;
       diri._tag=node;
       diri.onWhat=BoundaryCondition::ON_VERTEX;
@@ -118,7 +156,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       if(fscanf(f, "%d %d %lf", &edge, &comp, &val) != 3) return;
       dirichletBC diri;
       diri.g = new groupOfElements (1, edge);
-      diri._f= simpleFunction<double>(val);
+      diri._f= new simpleFunction<double>(val);
       diri._comp=comp;
       diri._tag=edge;
       diri.onWhat=BoundaryCondition::ON_EDGE;
@@ -130,7 +168,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       if(fscanf(f, "%d %d %lf", &face, &comp, &val) != 3) return;
       dirichletBC diri;
       diri.g = new groupOfElements (2, face);
-      diri._f= simpleFunction<double>(val);
+      diri._f= new simpleFunction<double>(val);
       diri._comp=comp;
       diri._tag=face;
       diri.onWhat=BoundaryCondition::ON_FACE;
@@ -142,7 +180,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
       neumannBC neu;
       neu.g = new groupOfElements (0, node);
-      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._f= new simpleFunction<SVector3>(SVector3(val1, val2, val3));
       neu._tag=node;
       neu.onWhat=BoundaryCondition::ON_VERTEX;
       allNeumann.push_back(neu);
@@ -153,7 +191,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       if(fscanf(f, "%d %lf %lf %lf", &edge, &val1, &val2, &val3) != 4) return;
       neumannBC neu;
       neu.g = new groupOfElements (1, edge);
-      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._f= new simpleFunction<SVector3>(SVector3(val1, val2, val3));
       neu._tag=edge;
       neu.onWhat=BoundaryCondition::ON_EDGE;
       allNeumann.push_back(neu);
@@ -164,7 +202,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       if(fscanf(f, "%d %lf %lf %lf", &face, &val1, &val2, &val3) != 4) return;
       neumannBC neu;
       neu.g = new groupOfElements (2, face);
-      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._f= new simpleFunction<SVector3>(SVector3(val1, val2, val3));
       neu._tag=face;
       neu.onWhat=BoundaryCondition::ON_FACE;
       allNeumann.push_back(neu);
@@ -175,7 +213,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
       if(fscanf(f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3) != 4) return;
       neumannBC neu;
       neu.g = new groupOfElements (3, volume);
-      neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+      neu._f= new simpleFunction<SVector3>(SVector3(val1, val2, val3));
       neu._tag=volume;
       neu.onWhat=BoundaryCondition::ON_VOLUME;
       allNeumann.push_back(neu);
@@ -200,26 +238,96 @@ void elasticitySolver::readInputFile(const std::string &fn)
       pModel = pModel->buildCutGModel(&ls);
     }
     else {
-      Msg::Error("Invalid input : %s", what);
+      Msg::Error("Invalid input : '%s'", what);
 //      return;
     }
   }
   fclose(f);
 }
 
-void elasticitySolver::solve()
+/********************* end deprecated api ****************************************/
+
+void elasticitySolver::addDirichletBC (int dim, int entityId, int component, double value) {
+  dirichletBC diri;
+  diri.g = new groupOfElements (dim, entityId);
+  diri._f= new simpleFunction<double>(value);
+  diri._comp=component;
+  diri._tag=entityId;
+  switch (dim) {
+    case 0 : diri.onWhat=BoundaryCondition::ON_VERTEX; break;
+    case 1 : diri.onWhat=BoundaryCondition::ON_EDGE; break;
+    case 2 : diri.onWhat=BoundaryCondition::ON_FACE; break;
+    default : return;
+  }
+  allDirichlet.push_back(diri);
+}
+
+void elasticitySolver::addDirichletBCLua (int dim, int entityId, int component, std::string luaFunctionName, lua_State *L) {
+  dirichletBC diri;
+  diri.g = new groupOfElements (dim, entityId);
+  diri._f= new simpleFunctionLua<double>(L, luaFunctionName,0.);
+  diri._comp=component;
+  diri._tag=entityId;
+  switch (dim) {
+    case 0 : diri.onWhat=BoundaryCondition::ON_VERTEX; break;
+    case 1 : diri.onWhat=BoundaryCondition::ON_EDGE; break;
+    case 2 : diri.onWhat=BoundaryCondition::ON_FACE; break;
+    default : return;
+  }
+  allDirichlet.push_back(diri);
+}
+
+void elasticitySolver::addNeumannBC (int dim, int entityId, const std::vector<double> value) {
+  if(value.size()!=3) return;
+  neumannBC neu;
+  neu.g = new groupOfElements (dim, entityId);
+  neu._f= new simpleFunction<SVector3>(SVector3(value[0], value[1], value[2]));
+  neu._tag=entityId;
+  switch (dim) {
+    case 0 : neu.onWhat=BoundaryCondition::ON_VERTEX; break;
+    case 1 : neu.onWhat=BoundaryCondition::ON_EDGE; break;
+    case 2 : neu.onWhat=BoundaryCondition::ON_FACE; break;
+    default : return;
+  }
+  allNeumann.push_back(neu);
+}
+
+void elasticitySolver::addNeumannBCLua (int dim, int entityId, std::string luaFunctionName, lua_State *L) {
+  neumannBC neu;
+  neu.g = new groupOfElements (dim, entityId);
+  neu._f= new simpleFunctionLua<SVector3>(L, luaFunctionName, SVector3(0,0,0));
+  neu._tag=entityId;
+  switch (dim) {
+    case 0 : neu.onWhat=BoundaryCondition::ON_VERTEX; break;
+    case 1 : neu.onWhat=BoundaryCondition::ON_EDGE; break;
+    case 2 : neu.onWhat=BoundaryCondition::ON_FACE; break;
+    default : return;
+  }
+  allNeumann.push_back(neu);
+}
+
+void elasticitySolver::addElasticDomain (int physical, double e, double nu){
+  elasticField field;
+  field._tag = _tag;
+  field._E = e;
+  field._nu = nu;
+  field.g = new groupOfElements (_dim, physical);
+  elasticFields.push_back(field);
+}
+
+elasticitySolver::elasticitySolver(GModel *model, int tag)
+{
+  pModel = model;
+  _dim = pModel->getNumRegions() ? 3 : 2;
+  _tag = tag;
+  pAssembler = NULL;
+  if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
+  if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+  LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1);
+}
+
+void elasticitySolver::assemble(linearSystem<double> *lsys)
 {
-  linearSystemFull<double> *lsys = new linearSystemFull<double>;
-  /*
-#if defined(HAVE_TAUCS)
-  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
-#elif defined(HAVE_PETSC)
-  linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
-#else
-  linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
-  lsys->setNoisy(2);
-#endif
-  */
   if (pAssembler) delete pAssembler;
   pAssembler = new dofManager<double>(lsys);
 
@@ -231,7 +339,7 @@ void elasticitySolver::solve()
   for (unsigned int i = 0; i < allDirichlet.size(); i++)
   {
     FilterDofComponent filter(allDirichlet[i]._comp);
-    FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
+    FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,*allDirichlet[i]._f,filter);
   }
   // LagrangeMultipliers
   for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i)
@@ -255,7 +363,7 @@ void elasticitySolver::solve()
   std::cout <<  "Neumann BC"<< std::endl;
   for (unsigned int i = 0; i < allNeumann.size(); i++)
   {
-    LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
+    LoadTerm<SVector3> Lterm(*LagSpace,*allNeumann[i]._f);
     Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
   }
   // Assemble cross term, laplace term and rhs term for LM
@@ -281,32 +389,18 @@ void elasticitySolver::solve()
     Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
   }
 
-  for (int i=0;i<pAssembler->sizeOfR();i++){
+  /*for (int i=0;i<pAssembler->sizeOfR();i++){
     for (int j=0;j<pAssembler->sizeOfR();j++){
       double d;lsys->getFromMatrix(i,j,d);
       printf("%12.5E ",d);
     }
     double d;lsys->getFromRightHandSide(i,d);
     printf(" |  %12.5E\n",d);
-  }
+  }*/
 
   printf("nDofs=%d\n",pAssembler->sizeOfR());
   printf("nFixed=%d\n",pAssembler->sizeOfF());
   printf("-- done assembling!\n");
-  lsys->systemSolve();
-  printf("-- done solving!\n");
-
-//   printLinearSystem(lsys);
-
-  double energ=0;
-  for (unsigned int i = 0; i < elasticFields.size(); i++)
-  {
-    SolverField<SVector3> Field(pAssembler, LagSpace);
-    IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
-    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
-    Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
-  }
-  printf("elastic energy=%f\n",energ);
 }
 
 
@@ -352,7 +446,7 @@ static double vonMises(dofManager<double> *a, MElement *e,
   return ComputeVonMises(s);
 }
 
-PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
+PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
 {
   std::cout <<  "build Displacement View"<< std::endl;
   std::set<MVertex*> v;
@@ -396,7 +490,7 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
   return pv;
 }
 
-PView* elasticitySolver::buildLagrangeMultiplierView (const std::string &postFileName)
+PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFileName)
 {
   std::cout <<  "build Lagrange Multiplier View"<< std::endl;
   if(!LagrangeMultiplierSpace) return new PView();
@@ -424,7 +518,7 @@ PView* elasticitySolver::buildLagrangeMultiplierView (const std::string &postFil
   return pv;
 }
 
-PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
+PView *elasticitySolver::buildElasticEnergyView(const std::string postFileName)
 {
   std::cout <<  "build Elastic Energy View"<< std::endl;
   std::map<int, std::vector<double> > data;
@@ -454,7 +548,7 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
   return pv;
 }
 
-PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
+PView *elasticitySolver::buildVonMisesView(const std::string postFileName)
 {std::cout <<  "build elastic view"<< std::endl;
   std::map<int, std::vector<double> > data;
   GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
@@ -506,7 +600,6 @@ PView* elasticitySolver::buildVonMisesView(const std::string &postFileName)
 }
 #endif
 
-
 #include "Bindings.h"
 void elasticitySolverRegisterBindings(binding *b) 
 {
@@ -519,12 +612,52 @@ void elasticitySolverRegisterBindings(binding *b)
   cm->setDescription ("reads an input file");
   cm->setArgNames("fileName",NULL);
 
-  cm = cb->addMethod("solve", &elasticitySolver::solve);
-  cm->setDescription ("solves the problem");
+  cm = cb->addMethod("assemble", &elasticitySolver::assemble);
+  cm->setDescription ("assembles the problem");
+  cm->setArgNames ("linearSystem",NULL);
 
-  cm = cb->setConstructor<elasticitySolver,int>();
-  cm->setDescription ("A new elasticitySolver. The parameter is the unknowns tag");
-  cm->setArgNames("tag",NULL);
+  cm = cb->addMethod("addDirichletBC", &elasticitySolver::addDirichletBC);
+  cm->setDescription ("add a Dirichlet (displacement) boundary condition on a given entity");
+  cm->setArgNames ("dim", "entityId", "component", "value", NULL);
+
+  cm = cb->addMethod("addDirichletBCLua", &elasticitySolver::addDirichletBCLua);
+  cm->setDescription ("add a Dirichlet (displacement) boundary condition on a given entity. The lua function takes x,y,z as arguments and return the displacement value.");
+  cm->setArgNames ("dim", "entityId", "component", "luaFunctionName", NULL);
 
+  cm = cb->addMethod("addNeumannBC", &elasticitySolver::addNeumannBC);
+  cm->setDescription ("add a Neumann (force) boundary condition on a given entity. Size of value has to be 3.");
+  cm->setArgNames ("dim", "entityId", "value", NULL);
+
+  cm = cb->addMethod("addNeumannBCLua", &elasticitySolver::addNeumannBCLua);
+  cm->setDescription ("add a Neumann (force) boundary condition on a given entity. The lua function takes x,y,z as arguments and return a vector with the 3 components of the force.");
+  cm->setArgNames ("dim", "entityId", "luaFunctionName", NULL);
+
+  cm = cb->addMethod("addElasticDomain", &elasticitySolver::addElasticDomain);
+  cm->setDescription ("add an elastic region with parameters E and nu");
+  cm->setArgNames ("tag", "E", "nu", NULL);
+
+#if defined HAVE_POST
+
+  cm = cb->addMethod ("buildVonMisesView", &elasticitySolver::buildVonMisesView);
+  cm->setDescription ("create a new view.");
+  cm->setArgNames ("fileName", NULL);
+
+  cm = cb->addMethod ("buildDisplacementView", &elasticitySolver::buildDisplacementView);
+  cm->setDescription ("create a new view."); 
+  cm->setArgNames ("fileName", NULL);
+
+  cm = cb->addMethod ("buildLagrangeMultiplierView", &elasticitySolver::buildLagrangeMultiplierView);
+  cm->setDescription ("create a new view.");
+  cm->setArgNames ("fileName", NULL);
+
+  cm = cb->addMethod ("buildElasticEnergyView", &elasticitySolver::buildElasticEnergyView);
+  cm->setDescription ("create a new view.");
+  cm->setArgNames ("fileName", NULL);
+
+#endif
+
+  cm = cb->setConstructor<elasticitySolver,GModel*,int>();
+  cm->setDescription ("A new elasticitySolver. The parameter is the unknowns tag");
+  cm->setArgNames("model","tag",NULL);
 }
 

Solver/elasticitySolver.h

@@ -45,14 +45,15 @@ struct BoundaryCondition
 struct dirichletBC : public BoundaryCondition
 {
   int _comp; // component
-  simpleFunction<double> _f;
-  dirichletBC () :BoundaryCondition(),_comp(0),_f(0){}
+  simpleFunction<double> *_f;
+  dirichletBC ():BoundaryCondition(),_comp(0),_f(0){}
+  dirichletBC (int dim, int entityId, int component, double value);
 };
 
 struct neumannBC  : public BoundaryCondition
 {
-  simpleFunction<SVector3> _f;
-  neumannBC () : BoundaryCondition(),_f(SVector3(0,0,0)){}
+  simpleFunction<SVector3> *_f;
+  neumannBC () : BoundaryCondition(),_f(NULL){}
 };
 
 // an elastic solver ...
@@ -75,20 +76,29 @@ class elasticitySolver
   
  public:
   elasticitySolver(int tag) : _tag(tag),LagSpace(0),pAssembler(0),LagrangeMultiplierSpace(0) {}
+
+  elasticitySolver(GModel *model, int tag);
+  void addDirichletBC (int dim, int entityId, int component, double value);
+  void addDirichletBCLua (int dim, int entityId, int component, std::string luaFunctionName, lua_State *L);
+  void addNeumannBC (int dim, int entityId, const std::vector<double> value);
+  void addNeumannBCLua (int dim, int entityId, std::string luaFunctionName, lua_State *L);
+  void addElasticDomain (int tag, double e, double nu);
+
   virtual ~elasticitySolver()
   {
     if (LagSpace) delete LagSpace;
     if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
     if (pAssembler) delete pAssembler;
   }
+  void assemble (linearSystem<double> *lsys);
   void readInputFile(const std::string &meshFileName);
   void read(const std::string s) {readInputFile(s.c_str());}
   virtual void setMesh(const std::string &meshFileName);
   void solve();
-  virtual PView *buildDisplacementView(const std::string &postFileName);
-  virtual PView *buildLagrangeMultiplierView(const std::string &posFileName);
-  virtual PView *buildElasticEnergyView(const std::string &postFileName);
-  virtual PView *buildVonMisesView(const std::string &postFileName);
+  virtual PView *buildDisplacementView(const std::string postFileName);
+  virtual PView *buildLagrangeMultiplierView(const std::string posFileName);
+  virtual PView *buildElasticEnergyView(const std::string postFileName);
+  virtual PView *buildVonMisesView(const std::string postFileName);
   // std::pair<PView *, PView*> buildErrorEstimateView
   //   (const std::string &errorFileName, double, int);
   // std::pair<PView *, PView*> buildErrorEstimateView

benchmarks/elasticity/conge.geo

+
+Mesh.CharacteristicLengthFromCurvature = 0;
+Mesh.CharacteristicLengthFromPoints = 1;
+Mesh.CharacteristicLengthExtendFromBoundary = 1;
+
+unit = 1.0e-02 ;
+
+e1 =  4.5 * unit ;
+e2 =  6.0 * unit / 2.0 ;
+e3 =  5.0 * unit / 2.0 ;
+h1 =  5.0 * unit ;
+h2 = 10.0 * unit ;
+h3 =  5.0 * unit ;
+//h4 =  1.0 * unit ;
+h4 =  2.0 * unit ;
+h5 =  4.5 * unit ;
+R1 =  1.0 * unit ;
+R2 =  1.5 * unit ;
+//R2 =  1.5 * unit ;
+r  =  2 * unit ;
+ccos = (-h5*R1+e2* (h5*h5+e2*e2-R1*R1)^0.5) / (h5*h5+e2*e2) ;
+ssin = ( 1.0 - ccos*ccos )^0.5 ;
+
+eps = 0.01 * unit;
+
+Lc1 = 0.001 ;
+Lc2 = 0.001 ;
+
+Point(1) = { -e1-e2, 0.0  , 0.0 , Lc1};
+Point(2) = { -e1-e2, h1   , 0.0 , Lc1};
+Point(3) = { -e3-r , h1   , 0.0 , Lc2};
+Point(4) = { -e3-r , h1+r , 0.0 , Lc2};
+Point(5) = { -e3   , h1+r , 0.0 , Lc2};
+Point(6) = { -e3   , h1+h2-eps, 0.0 , Lc1};
+Point(7) = {  e3   , h1+h2, 0.0 , Lc1};
+Point(8) = {  e3   , h1+r , 0.0 , Lc2};
+Point(9) = {  e3+r , h1+r , 0.0 , Lc2};
+Point(10)= {  e3+r , h1   , 0.0 , Lc2};
+Point(11)= {  e1+e2, h1   , 0.0 , Lc1};
+Point(12)= {  e1+e2, 0.0  , 0.0 , Lc1};
+Point(13)= {  e2   , 0.0  , 0.0 , Lc1};
+
+Point(14)= {  R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
+Point(15)= {  0.0       , h5        , 0.0 , Lc2};
+Point(16)= { -R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
+Point(17)= { -e2        , 0.0       , 0.0 , Lc1};
+
+Point(18)= { -R2  , h1+h3   , 0.0 , Lc2};
+Point(19)= { -R2  , h1+h3+h4, 0.0 , Lc2};
+Point(20)= {  0.0 , h1+h3+h4, 0.0 , Lc2};
+Point(21)= {  R2  , h1+h3+h4, 0.0 , Lc2};
+Point(22)= {  R2  , h1+h3   , 0.0 , Lc2};
+Point(23)= {  0.0 , h1+h3   , 0.0 , Lc2};
+
+Point(24)= {  0 , h1+h3+h4+R2, 0.0 , Lc2};
+Point(25)= {  0 , h1+h3-R2,    0.0 , Lc2};
+
+Line(1)  = {1 ,17};    /* ux=uy=0 */
+Line(2)  = {17,16};
+Circle(3) = {14,15,16};
+Line(4)  = {14,13};
+Line(5)  = {13,12};    /* ux=uy=0 */
+Line(6)  = {12,11};
+Line(7)  = {11,10};
+Circle(8) = { 8, 9,10};
+Line(9)  = { 8, 7};
+Line(10) = { 7, 6};    /* T=10000 N */
+Line(11) = { 6, 5};
+Circle(12) = { 3, 4, 5};
+Line(13) = { 3, 2};
+Line(14) = { 2, 1};
+
+Line(15) = {18,19};
+Circle(16) = {21,20,24};
+Circle(17) = {24,20,19};
+Circle(18) = {18,23,25};
+Circle(19) = {25,23,22};
+Line(20) = {21,22};
+
+Line Loop(21) = {17,-15,18,19,-20,16};
+//Plane Surface(22) = {21};
+Line Loop(23) = {11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10};
+Plane Surface(24) = {23,21};
+
+//Physical Line(25) = {9,1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,18,19,20,10};
+//Physical Surface(26) = {22,24};
+Physical Line(9) = {11};
+Physical Line(10) = {10};
+Physical Line(8) = {1, 5};
+Physical Surface(7) = {24};
+Physical Point(20) = {2};
+Physical Point(21) = {11};
+// Physical Line(25) = {11, 12, 13, 14};
+

benchmarks/elasticity/conge.lua

+function neumannCondition (x,y,z)
+  return {-100e6,0,0}
+end
+
+m = GModel()
+m:load("conge.geo")
+m:load("conge.msh")
+e = elasticitySolver(m,1)
+e:addElasticDomain (7, 200e9, 0.3)
+e:addDirichletBC (1,8,0,0)
+e:addDirichletBC (1,8,1,0)
+e:addDirichletBC (1,8,2,0)
+e:addNeumannBCLua (1,9,"neumannCondition")
+sys = linearSystemCSRTaucs()
+e:assemble(sys)
+sys:systemSolve()
+view = e:buildVonMisesView("vonMises")
+view = e:buildDisplacementView("displacement")
+view = e:buildLagrangeMultiplierView("lagrangeMultiplier")
+view = e:buildElasticEnergyView("elasticEnergy")
+--view:write("test.pos",2,false)