diff --git a/CMakeLists.txt b/CMakeLists.txt
index 89bfd3a0ce46bfcef0bc49303f1458e4c48fe9b8..716c8f82b8696de8fe6cdab3e3ad7e1e352f88ac 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -28,6 +28,7 @@ option(ENABLE_FOURIER_MODEL "Enable Fourier geometrical models" OFF)
option(ENABLE_GMM "Enable GMM linear algebra solvers" ON)
option(ENABLE_GRAPHICS "Compile-in OpenGL graphics even if there is no GUI" OFF)
option(ENABLE_KBIPACK "Enable Kbipack for homology solver" ON)
+option(ENABLE_LUA "Enable the Programming Language Lua " ON)
option(ENABLE_MATHEX "Enable MathEx expression parser" ON)
option(ENABLE_MED "Enable MED mesh and post-processing file formats" ON)
option(ENABLE_MESH "Build the mesh module" ON)
@@ -75,7 +76,7 @@ set(GMSH_API
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h
Mesh/meshGFaceDelaunayInsertion.h
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
- Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemFull.h
+ Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h Solver/linearSystemFull.h
Solver/elasticitySolver.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h
Graphics/drawContext.h
@@ -547,6 +548,19 @@ if(ENABLE_TAUCS)
endif(TAUCS_LIB)
endif(ENABLE_TAUCS)
+if(ENABLE_LUA)
+ find_library(LUA_LIB lua PATH_SUFFIXES lib)
+ if(LUA_LIB)
+ find_path(LUA_INC "lua.h" PATH_SUFFIXES src include)
+ if(LUA_INC)
+ set(HAVE_LUA TRUE)
+ list(APPEND CONFIG_OPTIONS "Lua")
+ list(APPEND EXTERNAL_LIBRARIES ${LUA_LIB})
+ list(APPEND EXTERNAL_INCLUDES ${LUA_INC})
+ endif(LUA_INC)
+ endif(LUA_LIB)
+endif(ENABLE_LUA)
+
if(ENABLE_PETSC)
set(ENV_PETSC_DIR $ENV{PETSC_DIR})
set(ENV_PETSC_ARCH $ENV{PETSC_ARCH})
@@ -989,13 +1003,16 @@ message("Run 'ccmake ${CMAKE_SOURCE_DIR}' to fine-tune the configuration.")
message("")
mark_as_advanced(BISON FLEX GMP_LIB GMSH_EXTRA_VERSION HDF5_LIB MAKEINFO
- MED_LIB OCC_INC SZ_LIB TAUCS_LIB TEXI2PDF)
+ MED_LIB OCC_INC SZ_LIB TAUCS_LIB LUA_LIB TEXI2PDF)
-add_executable(gmshdg EXCLUDE_FROM_ALL Solver/dgMainTest.cpp ${GMSH_SRC})
+add_executable(gmshdg EXCLUDE_FROM_ALL Solver/dgMainTest.cpp ${GMSH_SRC})
+add_executable(gmshlua EXCLUDE_FROM_ALL Solver/dgMainLua.cpp ${GMSH_SRC})
target_link_libraries(gmshdg ${LINK_LIBRARIES})
+target_link_libraries(gmshlua ${LINK_LIBRARIES})
add_executable(gmshdgsw EXCLUDE_FROM_ALL Solver/dgMainShallowWater2d.cpp ${GMSH_SRC})
target_link_libraries(gmshdgsw ${LINK_LIBRARIES})
add_executable(gmshdgwave EXCLUDE_FROM_ALL Solver/dgMainWaveEquation.cpp ${GMSH_SRC})
target_link_libraries(gmshdgwave ${LINK_LIBRARIES})
+
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index 0659077523f181eba21f808f2620593700476ae2..245960ab20f22ce30901e2f7b1f02284db09591f 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -21,6 +21,7 @@
#cmakedefine HAVE_LIBJPEG
#cmakedefine HAVE_LIBPNG
#cmakedefine HAVE_LIBZ
+#cmakedefine HAVE_LUA
#cmakedefine HAVE_MATHEX
#cmakedefine HAVE_MED
#cmakedefine HAVE_MESH
diff --git a/Common/luna.h b/Common/luna.h
new file mode 100755
index 0000000000000000000000000000000000000000..67d066d637cafc0a5b0d2fda4122559f2cbeda79
--- /dev/null
+++ b/Common/luna.h
@@ -0,0 +1,131 @@
+/*
+Copyright (c) 2005 Leonardo Palozzi
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+/*
+ * Lenny Palozzi - lenny.palozzi@home.com
+ */
+#ifndef _LUNA_H_
+#define _LUNA_H_
+extern "C" {
+#include <lua.h>
+}
+extern "C" {
+#include "lua.h"
+#include "lauxlib.h"
+}
+
+template <typename T> class Luna {
+ typedef struct { T *pT; } userdataType;
+public:
+ typedef int (T::*mfp)(lua_State *L);
+ typedef struct { const char *name; mfp mfunc; } RegType;
+
+ static void Register(lua_State *L) {
+ lua_newtable(L);
+ int methods = lua_gettop(L);
+
+ luaL_newmetatable(L, T::className);
+ int metatable = lua_gettop(L);
+
+ // store method table in globals so that
+ // scripts can add functions written in Lua.
+ lua_pushstring(L, T::className);
+ lua_pushvalue(L, methods);
+ lua_settable(L, LUA_GLOBALSINDEX);
+
+ lua_pushliteral(L, "__metatable");
+ lua_pushvalue(L, methods);
+ lua_settable(L, metatable); // hide metatable from Lua getmetatable()
+
+ lua_pushliteral(L, "__index");
+ lua_pushvalue(L, methods);
+ lua_settable(L, metatable);
+
+ lua_pushliteral(L, "__tostring");
+ lua_pushcfunction(L, tostring_T);
+ lua_settable(L, metatable);
+
+ lua_pushliteral(L, "__gc");
+ lua_pushcfunction(L, gc_T);
+ lua_settable(L, metatable);
+
+ lua_newtable(L); // mt for method table
+ int mt = lua_gettop(L);
+ lua_pushliteral(L, "__call");
+ lua_pushcfunction(L, new_T);
+ lua_pushliteral(L, "new");
+ lua_pushvalue(L, -2); // dup new_T function
+ lua_settable(L, methods); // add new_T to method table
+ lua_settable(L, mt); // mt.__call = new_T
+ lua_setmetatable(L, methods);
+
+ // fill method table with methods from class T
+ for (RegType *l = T::methods; l->name; l++) {
+ /* edited by Snaily: shouldn't it be const RegType *l ... ? */
+ lua_pushstring(L, l->name);
+ lua_pushlightuserdata(L, (void*)l);
+ lua_pushcclosure(L, thunk, 1);
+ lua_settable(L, methods);
+ }
+
+ lua_pop(L, 2); // drop metatable and method table
+ }
+
+ // get userdata from Lua stack and return pointer to T object
+ static T *check(lua_State *L, int narg) {
+ userdataType *ud =
+ static_cast<userdataType*>(luaL_checkudata(L, narg, T::className));
+ if(!ud) luaL_typerror(L, narg, T::className);
+ return ud->pT; // pointer to T object
+ }
+
+private:
+ Luna(); // hide default constructor
+
+ // member function dispatcher
+ static int thunk(lua_State *L) {
+ // stack has userdata, followed by method args
+ T *obj = check(L, 1); // get 'self', or if you prefer, 'this'
+ lua_remove(L, 1); // remove self so member function args start at index 1
+ // get member function from upvalue
+ RegType *l = static_cast<RegType*>(lua_touserdata(L, lua_upvalueindex(1)));
+ return (obj->*(l->mfunc))(L); // call member function
+ }
+
+ // create a new T object and
+ // push onto the Lua stack a userdata containing a pointer to T object
+ static int new_T(lua_State *L) {
+ lua_remove(L, 1); // use classname:new(), instead of classname.new()
+ T *obj = new T(L); // call constructor for T objects
+ userdataType *ud =
+ static_cast<userdataType*>(lua_newuserdata(L, sizeof(userdataType)));
+ ud->pT = obj; // store pointer to object in userdata
+ luaL_getmetatable(L, T::className); // lookup metatable in Lua registry
+ lua_setmetatable(L, -2);
+ return 1; // userdata containing pointer to T object
+ }
+
+ // garbage collection metamethod
+ static int gc_T(lua_State *L) {
+ userdataType *ud = static_cast<userdataType*>(lua_touserdata(L, 1));
+ T *obj = ud->pT;
+ delete obj; // call destructor for T objects
+ return 0;
+ }
+
+ static int tostring_T (lua_State *L) {
+ char buff[32];
+ userdataType *ud = static_cast<userdataType*>(lua_touserdata(L, 1));
+ T *obj = ud->pT;
+ sprintf(buff, "%p", obj);
+ lua_pushfstring(L, "%s (%s)", T::className, buff);
+ return 1;
+ }
+};
+#endif
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 89e0c22ee78eb24687ec7723fc55c3ab0e93b3f0..d119667e11b5395ee6e04288b44fcc7b83377bc1 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -25,6 +25,7 @@ set(SRC
findLinks.cpp
SOrientedBoundingBox.cpp
GeomMeshMatcher.cpp
+ LuaBindings_Geo.cpp
MVertex.cpp
MFace.cpp
MElement.cpp
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index da343bbda8c210c1ba4a11cb5f1b057339f10bb2..de9f1ce6509b96ebbeb8d41ecde90650588433b4 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -33,6 +33,7 @@
#include "Context.h"
#include "discreteFace.h"
#include "eigenSolver.h"
+#include "multiscaleLaplace.h"
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
{
@@ -934,6 +935,20 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
return;
}
+ // TEST --------------------------------------------------
+ if(step == ITERV){
+ std::vector<MElement*> _elements;
+ std::list<GFace*>::const_iterator it = _compound.begin();
+ for( ; it != _compound.end(); ++it){
+ for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ _elements.push_back(t);
+ }
+ }
+ // multiscaleLaplace(_elements,ordered,coords);
+ }
+ // TEST --------------------------------------------------
+
for(unsigned int i = 0; i < ordered.size(); i++){
MVertex *v = ordered[i];
const double theta = 2 * M_PI * coords[i];
@@ -1217,6 +1232,7 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->curvatureMax(pv);
}
+ // emi fait qqch ici ...
return 0.;
}
@@ -1248,10 +1264,10 @@ GPoint GFaceCompound::point(double par1, double par2) const
gp.setNoSuccess();
return gp;
}
- if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
- SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
- return lt->gf->point(pv.x(),pv.y());
- }
+ // if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
+ // SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
+ // return lt->gf->point(pv.x(),pv.y());
+ // }
const bool LINEARMESH = true; //false
diff --git a/Geo/GModel.h b/Geo/GModel.h
index c528a3135b18c8436dc62e5725685e0b77d8e775..719d9b6245ae8b300ea08ebb389fb8d5c405d61f 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -120,6 +120,18 @@ class GModel
// return the current model, and sets the current model index if
// index >= 0
static GModel *current(int index=-1);
+
+ // sets a model to current
+ static int setCurrent(GModel *m){
+ for (int i=0;i<list.size();i++){
+ if (list[i] ==m){
+ _current = i;
+ return i;
+ }
+ }
+ }
+
+
// find a model by name
static GModel *findByName(std::string name);
diff --git a/Geo/LuaBindings_Geo.cpp b/Geo/LuaBindings_Geo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9cd89477dd7211be48c338955579094f57956f9
--- /dev/null
+++ b/Geo/LuaBindings_Geo.cpp
@@ -0,0 +1,51 @@
+#include "LuaBindings_Geo.h"
+#include "OpenFile.h"
+#include "CreateFile.h"
+// define the name of the object that will be used
+// in Lua
+
+const char LuaGModel::className[] = "GModel";
+
+// Define the methods we will expose to Lua
+#define _method(class, name) {#name, &class::name}
+
+Luna<LuaGModel>::RegType LuaGModel::methods[] = {
+ _method(LuaGModel, mesh),
+ _method(LuaGModel, load),
+ _method(LuaGModel, save),
+ {0,0}
+};
+
+LuaGModel::LuaGModel(lua_State *L){
+ _gm = new GModel;
+}
+
+LuaGModel::~LuaGModel(){
+ delete _gm;
+}
+
+void LuaGModel::Register (lua_State *L){
+ Luna<LuaGModel>::Register(L);
+}
+
+int LuaGModel::load(lua_State *L){
+ GModel *temp = GModel::current();
+ GModel::setCurrent(_gm);
+ MergeFile(luaL_checkstring(L, 1));
+ GModel::setCurrent(temp);
+ return 1;
+}
+
+int LuaGModel::mesh(lua_State *L){
+ _gm -> mesh((int)luaL_checknumber(L, 1));
+ return 1;
+}
+
+int LuaGModel::save(lua_State *L){
+ GModel *temp = GModel::current();
+ GModel::setCurrent(_gm);
+ int guess = GuessFileFormatFromFileName(std::string(lua_tostring(L, 1)));
+ CreateOutputFile (std::string(lua_tostring(L, 1)), guess);
+ GModel::setCurrent(temp);
+ return 1;
+}
diff --git a/Geo/LuaBindings_Geo.h b/Geo/LuaBindings_Geo.h
new file mode 100644
index 0000000000000000000000000000000000000000..be6dda96f04fdd5b0f79320222cf9859f8a75c70
--- /dev/null
+++ b/Geo/LuaBindings_Geo.h
@@ -0,0 +1,39 @@
+#ifndef _LUA_BINDINGS_GEO_
+#define _LUA_BINDINGS_GEO_
+#include "GmshConfig.h"
+
+#if defined(HAVE_LUA)
+// include lua stuff
+extern "C" {
+#include "lua.h"
+#include "lauxlib.h"
+#include "lualib.h"
+}
+// Use luna for bindings
+#include "luna.h"
+// The header file for the GModel
+#include "GModel.h"
+
+// we do not use the GModel::current()
+// this should not be a class anyway ...
+// this should be removed !!!!!! and be put
+// inside GModel itself
+class LuaGModel{
+ GModel *_gm;
+public:
+ static const char className[];
+ static Luna<LuaGModel>::RegType methods[];
+ static void Register(lua_State *L);
+ LuaGModel(lua_State *L);
+ ~LuaGModel();
+ // Methods we will bind
+ int load (lua_State *L); // load a file
+ int mesh (lua_State *L); // make the mesh
+ int save (lua_State *L); // save a file
+ GModel * getGModel () const {return _gm;}
+private:
+};
+
+#endif // HAVE LUA
+
+#endif // _LUA_BINDINGS_GEO_
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 52c4c867092e34164e88e144cac849ad5d6910c9..0312dd4b40440f8fa33c684cb537c248deda1126 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -277,4 +277,58 @@ bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S)
return false;
}
+#if defined(HAVE_LUA)
+template<>
+int fullMatrix<double>::gemm(lua_State *L)
+{
+#if defined(HAVE_BLAS)
+ int n = lua_gettop(L);
+ if (n < 2)throw;
+
+ fullMatrix<double> *a = Luna<fullMatrix<double> >::check(L,1);
+ fullMatrix<double> *b = Luna<fullMatrix<double> >::check(L,2);
+
+ if(!a || !b)throw;
+ // printf("%d %d\n",a->size1(),b->size1());
+ double alpha=1, beta=1;
+ if (n > 2) alpha = luaL_checknumber(L, 3);
+ if (n > 3) beta = luaL_checknumber(L, 4);
+ int M = size1(), N = size2(), K = a->size2();
+ int LDA = a->size1(), LDB = b->size1(), LDC = size1();
+ F77NAME(dgemm)("N", "N", &M, &N, &K, &alpha, a->_data, &LDA, b->_data, &LDB,
+ &beta, _data, &LDC);
+ return 0;
+#else
+ throw;
+ return 0;
+#endif
+}
+
+// define the name of the object that will be used in Lua
+template <>
+const char fullMatrix<double> ::className[] = "fullMatrix";
+template <>
+const char fullMatrix<float> ::className[] = "fullMatrixFloat";
+
+// Define the methods we will expose to Lua
+#define _method(class, name) {#name, &class::name}
+template <>
+Luna<fullMatrix<double> >::RegType fullMatrix<double>::methods[] = {
+ _method(fullMatrix<double> , size1),
+ _method(fullMatrix<double> , size2),
+ _method(fullMatrix<double> , get),
+ _method(fullMatrix<double> , set),
+ _method(fullMatrix<double> , gemm),
+ {0,0}
+};
+
+// this function has to be called in the main in order to register
+// the names defined above
+template<>
+void fullMatrix<double>::Register (lua_State *L){
+ Luna<fullMatrix<double> >::Register(L);
+}
+#endif
+
+
#endif
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index b4561cb25cd847ba7c556af82d7c3de6627e57c2..18d0efef5faedf1ff334c6e59802924f13d38157 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -11,6 +11,17 @@
#include "GmshConfig.h"
#include "GmshMessage.h"
+#if defined(HAVE_LUA)
+// include lua stuff
+extern "C" {
+#include "lua.h"
+#include "lauxlib.h"
+#include "lualib.h"
+}
+// Use luna for c++ class member functions bindings
+#include "luna.h"
+#endif // HAVE LUA
+
template <class scalar> class fullMatrix;
template <class scalar>
@@ -79,6 +90,51 @@ class fullMatrix
int _r, _c;
scalar *_data;
public:
+#if defined(HAVE_LUA)
+ static const char className[];
+ static Luna<fullMatrix<double> >::RegType methods[];
+ static void Register(lua_State *L);
+
+ fullMatrix(lua_State *L){
+ int n = lua_gettop(L);
+ if (n == 1){
+ fullMatrix<scalar> * _ud = (fullMatrix<scalar> *) lua_touserdata(L, 1);
+ _c = _ud->_c;
+ _r = _ud->_r;
+ _own_data = false;
+ _data = _ud->_data;
+ return;
+ }
+ if (n != 2)throw;
+ _r = luaL_checkint(L, 1);
+ _c = luaL_checkint(L, 2);
+ _data = new scalar[_r * _c];
+ _own_data = true;
+ scale(0.);
+ }
+ int size1(lua_State *L){
+ lua_pushinteger(L, _r);
+ return 1;
+ }
+ int size2(lua_State *L){
+ lua_pushinteger(L, _c);
+ return 1;
+ }
+ int get(lua_State *L){
+ int r = luaL_checkint(L, 1);
+ int c = luaL_checkint(L, 2);
+ lua_pushnumber(L, (*this)(r,c));
+ return 1;
+ }
+ int set(lua_State *L){
+ int r = luaL_checkint(L, 1);
+ int c = luaL_checkint(L, 2);
+ scalar val = luaL_checknumber(L, 3);
+ (*this)(r,c) = val;
+ return 0;
+ }
+ int gemm(lua_State *L);
+#endif // HAVE LUA
fullMatrix(fullMatrix<scalar> &original, int c_start, int c){
_c = c;
_r = original._r;
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index b052c9abd604c1c3b6069bcd16860b2938b5d0a5..60898af6e31f95eb5e2e0deb4f1e2c8e654ea68a 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -15,6 +15,7 @@ set(SRC
dgAlgorithm.cpp
dgConservationLaw.cpp
dgConservationLawAdvection.cpp
+ dgSystemOfEquations.cpp
dgConservationLawShallowWater2d.cpp
dgConservationLawWaveEquation.cpp
function.cpp
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 5fd24ac4c00ae45a863b77e438077072f70bfd92..9a1345a569fafef732ca4ddba58bd2fa2927adba 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -121,11 +121,11 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfFaces &group,
- const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &solutionRight,
- fullMatrix<double> &residual // residual !! at faces nodes
- )
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &solutionRight,
+ fullMatrix<double> &residual // residual !! at faces nodes
+ )
{
// A) global operations before entering the loop over the faces
// A1 ) copy locally some references from the group for the sake of lisibility
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 7e23dfe336bcd01a5c2e0328725491b8ebb31800..3482bbb4da76bcead48db0b2351c33230f46f59e 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -84,6 +84,7 @@ public:
};
class dgGroupOfFaces {
+ // normals always point outside left to right
// only used if this is a group of boundary faces
std::string _boundaryTag;
const dgGroupOfElements &_groupLeft,&_groupRight;
@@ -149,6 +150,7 @@ public:
//keep this outside the Algorithm because this is the only place where data overlap
void mapToInterface(int nFields, const fullMatrix<double> &vLeft, const fullMatrix<double> &vRight, fullMatrix<double> &v);
void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
+ const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
};
diff --git a/Solver/dgMainTest.cpp b/Solver/dgMainTest.cpp
index 48fdacd8d8ac92fe28aff3456a53f87d8fb8d6fc..bb5c8a25a95ccc877e10076bb5a0fe7af83105be 100644
--- a/Solver/dgMainTest.cpp
+++ b/Solver/dgMainTest.cpp
@@ -57,9 +57,6 @@ int main(int argc, char **argv){
algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
}
-
-
-
fullMatrix<double> advectionSpeed(3,1);
advectionSpeed(0,0)=0.15;
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4faccae24e068cf76a2d749249e5551596abd1fe
--- /dev/null
+++ b/Solver/dgSystemOfEquations.cpp
@@ -0,0 +1,142 @@
+#include "dgSystemOfEquations.h"
+#include "LuaBindings_Geo.h"
+
+#if defined(HAVE_LUA)
+// define the name of the object that will be used in Lua
+const char dgSystemOfEquations::className[] = "dgSystemOfEquations";
+// Define the methods we will expose to Lua
+#define _method(class, name) {#name, &class::name}
+Luna<dgSystemOfEquations>::RegType dgSystemOfEquations::methods[] = {
+ _method(dgSystemOfEquations, getRHS),
+ _method(dgSystemOfEquations, computeRHS),
+ _method(dgSystemOfEquations, getSolution),
+ _method(dgSystemOfEquations, exportSolution),
+ _method(dgSystemOfEquations, setConservationLaw),
+ _method(dgSystemOfEquations, setOrder),
+ _method(dgSystemOfEquations, setup),
+ _method(dgSystemOfEquations, addBoundaryCondition),
+ {0,0}
+};
+
+// this function has to be called in the main in order to register
+// the names defined above
+void dgSystemOfEquations::Register (lua_State *L){
+ Luna<dgSystemOfEquations>::Register(L);
+}
+
+// system of equations is build using a mesh
+dgSystemOfEquations::dgSystemOfEquations(lua_State *L){
+ // get the number of arguments
+ int n = lua_gettop(L);
+ if (n < 1)throw;
+ LuaGModel *obj =Luna<LuaGModel>::check(L, 1);
+ if (!obj)throw;
+ _gm = obj->getGModel();
+ _dimension = _gm->getNumRegions() ? 3 : 2;
+}
+
+// set the conservation law as a string (for now)
+int dgSystemOfEquations::setConservationLaw(lua_State *L){
+ _claw = 0;
+ _cLawName = std::string (luaL_checkstring(L, 1));
+ if (_cLawName == "WaveEquation")
+ _claw = dgNewConservationLawWaveEquation();
+ else if (_cLawName == "ShallowWater2d")
+ _claw = dgNewConservationLawShallowWater2d();
+ if (!_claw)throw;
+ return 0;
+}
+
+// set the order of interpolation
+int dgSystemOfEquations::setOrder(lua_State *L){
+ _order = luaL_checkint(L, 1);
+}
+
+// add a boundary Condition
+int dgSystemOfEquations::addBoundaryCondition (lua_State *L){
+ std::string physicalName(luaL_checkstring(L, 1));
+ std::string bcName(luaL_checkstring(L, 2));
+ if (_cLawName == "WaveEquation"){
+ if (bcName == "Wall"){
+ _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionWaveEquationWall());
+ }
+ else throw;
+ }
+ else if (_cLawName == "ShallowWater2d"){
+ if (bcName == "Wall"){
+ _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionShallowWater2dWall());
+ }
+ else throw;
+ }
+ else throw;
+}
+
+
+// ok, we can setup the groups and create solution vectors
+int dgSystemOfEquations::setup(lua_State *L){
+ if (!_claw) throw;
+ _algo->buildGroups(_gm,
+ _dimension,
+ _order,
+ _elementGroups,
+ _faceGroups,
+ _boundaryGroups);
+ //for now, we suppose there is only one group of elements
+ int nbNodes = _elementGroups[0]->getNbNodes();
+ int nbElements = _elementGroups[0]->getNbElements();
+ // allocate solution and RHS
+ int nbFields = _claw->nbFields();
+ _solution = new fullMatrix<double> (nbNodes,nbFields*nbElements);
+ _rightHandSide = new fullMatrix<double> (nbNodes,nbFields*nbElements);
+}
+
+
+int dgSystemOfEquations::getSolution(lua_State *L){
+ lua_pushlightuserdata (L, _solution);
+ return 1;
+}
+
+int dgSystemOfEquations::computeRHS(lua_State *L){
+ _algo->residual(*_claw, _elementGroups, _faceGroups, _boundaryGroups, *_solution, *_rightHandSide);
+ lua_pushlightuserdata (L, _solution);
+ return 1;
+}
+
+int dgSystemOfEquations::getRHS(lua_State *L){
+ lua_pushlightuserdata (L, _rightHandSide);
+ return 1;
+}
+
+int dgSystemOfEquations::exportSolution(lua_State *L){
+}
+#endif // HAVE_LUA
+
+dgSystemOfEquations::dgSystemOfEquations(GModel *obj, dgConservationLaw *claw,int order){
+ setup(obj,claw,order);
+}
+
+dgSystemOfEquations::~dgSystemOfEquations(){
+ delete _solution;
+ delete _rightHandSide;
+}
+
+void dgSystemOfEquations::setup(GModel *obj, dgConservationLaw *claw, int order){
+ _gm = obj;
+ _claw = claw;
+ _algo = new dgAlgorithm;
+ _order= order;
+ _dimension = _gm->getNumRegions() ? 3 : 2;
+ _algo->buildGroups(_gm,
+ _dimension,
+ _order,
+ _elementGroups,
+ _faceGroups,
+ _boundaryGroups);
+ //for now, we suppose there is only one group of elements
+ int nbNodes = _elementGroups[0]->getNbNodes();
+ int nbElements = _elementGroups[0]->getNbElements();
+ // allocate solution and RHS
+ int nbFields = _claw->nbFields();
+ _solution = new fullMatrix<double> (nbNodes,nbFields*nbElements);
+ _rightHandSide = new fullMatrix<double> (nbNodes,nbFields*nbElements);
+}
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
new file mode 100644
index 0000000000000000000000000000000000000000..879a6998b28c93e63064b221b6ebba6c681bdf52
--- /dev/null
+++ b/Solver/dgSystemOfEquations.h
@@ -0,0 +1,70 @@
+#ifndef _DG_SYSTEM_OF_EQUATIONS_
+#define _DG_SYSTEM_OF_EQUATIONS_
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "dgAlgorithm.h"
+#include "dgGroupOfElements.h"
+#include "dgAlgorithm.h"
+#include "dgConservationLaw.h"
+#include "Gmsh.h"
+
+#if defined(HAVE_LUA)
+// include lua stuff
+extern "C" {
+#include "lua.h"
+#include "lauxlib.h"
+#include "lualib.h"
+}
+// Use luna for c++ class member functions bindings
+#include "luna.h"
+#endif // HAVE LUA
+
+class dgSystemOfEquations {
+ // the mesh and the model
+ GModel *_gm;
+ // the algorithm that computes DG operators
+ dgAlgorithm *_algo;
+ // the conservation law
+ dgConservationLaw *_claw;
+ std::string _cLawName;
+ // polynomial order (should be more general)
+ int _order;
+ // dimension of the problem
+ int _dimension;
+ // the solution
+ fullMatrix<double> *_solution;
+ // the right hand side (a temporary vector)
+ fullMatrix<double> *_rightHandSide;
+ // groups of elements (volume terms)
+ std::vector<dgGroupOfElements*> _elementGroups;
+ // groups of faces (interface terms)
+ std::vector<dgGroupOfFaces*> _faceGroups;
+ // groups of faces (boundary conditions)
+ std::vector<dgGroupOfFaces*> _boundaryGroups;
+ // sets up everything
+ void setup(GModel *gm, dgConservationLaw *claw, int order);
+ dgSystemOfEquations(const dgSystemOfEquations &) {}
+public:
+ // lua stuff
+#if defined(HAVE_LUA)
+ static const char className[];
+ static Luna<dgSystemOfEquations>::RegType methods[];
+ static void Register(lua_State *L);
+ int computeRHS (lua_State *L); // get the Right Hand
+ int getRHS (lua_State *L); // get the Right Hand
+ int getSolution (lua_State *L); // get the Solution
+ int exportSolution (lua_State *L); // export the solution
+ int setOrder (lua_State *L); // set the polynomial order
+ int setConservationLaw (lua_State *L); // set the conservationLaw
+ int addBoundaryCondition (lua_State *L); // add a boundary condition : "physical name", "type", [options]
+ int setup (lua_State *L); // setup the groups and allocate
+ dgSystemOfEquations(lua_State *L);
+#endif // HAVE LUA
+ dgSystemOfEquations(GModel *gm, dgConservationLaw *claw, int order);
+ ~dgSystemOfEquations();
+
+private:
+};
+
+#endif // _DG_SYSTEM_OF_EQUATIONS_
+
diff --git a/Solver/function.h b/Solver/function.h
index e3be565594ee3b23669d5e5e670faf9dee4894bb..df56c4a4542915f1455079abb0b3b27752e09610 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -8,7 +8,7 @@
class dataCacheMap;
class MElement;
-// those classes manage complex function dependcies and keep their values in cache so that they are not recomputed when it is not necessary. To do this, we use three classes : function, dataCache and dataCacheMap. The workflow is :
+// those classes manage complex function dependencies and keep their values in cache so that they are not recomputed when it is not necessary. To do this, we use three classes : function, dataCache and dataCacheMap. The workflow is :
//
// a) before parsing the input file : a few general function that are always available are created (xyz for example).
//
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 0f5c09f5795a7ccd5705658a0c73e11fc950df2e..7f3e0ed9afc952f1ce329dccf7c56ed1834125fa 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -1,3 +1,4 @@
+#include "Context.h"
#include "multiscaleLaplace.h"
#include "GmshConfig.h"
#include "GmshDefines.h"
@@ -9,50 +10,291 @@
#include "laplaceTerm.h"
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
+#include "MTriangle.h"
+#include "robustPredicates.h"
+
+/*
+//-- returns 0 if no intersection occurs --\\
+//-- returns 1
+p1 (1-u) + p2 u = q1 (1-t) + q2 t
+
+we search the intersection into segment q1-q2
+
+robustPredicates
+
+if ( isLeftOf (
+
+*/
+static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
+ SPoint2 &q1, SPoint2 &q2,
+ double x[2]){
+ double P1[2] = {p1.x(),p1.y()};
+ double P2[2] = {p2.x(),p2.y()};
+ double Q1[2] = {q1.x(),q1.y()};
+ double Q2[2] = {q2.x(),q2.y()};
+
+ double PQ1 = robustPredicates::orient2d(P1,P2,Q1);
+ double PQ2 = robustPredicates::orient2d(P1,P2,Q2);
+
+ double QP1 = robustPredicates::orient2d(Q1,Q2,P1);
+ double QP2 = robustPredicates::orient2d(Q1,Q2,P2);
+
+}
+
+
+static int intersection_segments (SPoint2 &p1, SPoint2 &p2,
+ SPoint2 &q1, SPoint2 &q2,
+ double x[2]){
+ double A[2][2];
+ A[0][0] = p2.x()-p1.x();
+ A[0][1] = q1.x()-q2.x();
+ A[1][0] = p2.y()-p1.y();
+ A[1][1] = q1.y()-q2.y();
+ double b[2] = {q1.x()-p1.x(),q1.y()-p1.y()};
+ sys2x2(A,b,x);
+
+ return (x[0] >= 0.0 && x[0] <= 1. &&
+ x[1] >= 0.0 && x[1] <= 1.);
+
+}
+
+static void recur_compute_centers_ (double R, double a1, double a2,
+ multiscaleLaplaceLevel * root ){
+
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
+ multiscaleLaplaceLevel* zero = 0;
+
+ SPoint2 PL (R*cos(a1),R*sin(a1));
+ SPoint2 PR (R*cos(a2),R*sin(a2));
+ centers.push_back(std::make_pair(PL,zero));
+ centers.push_back(std::make_pair(PR,zero));
+ for (int i=0;i<root->childeren.size();i++){
+ centers.push_back(std::make_pair(root->childeren[i]->center,root->childeren[i]));
+ multiscaleLaplaceLevel* m = root->childeren[i];
+ m->radius = 0.0;
+ for (std::map<MVertex*,SPoint2>::iterator it = m->coordinates.begin();
+ it != m->coordinates.end() ; ++it){
+ const SPoint2 &p = it->second;
+ m->radius = std::max(m->radius,sqrt ((m->center.x() - p.x())*(m->center.x() - p.x())+
+ (m->center.y() - p.y())*(m->center.y() - p.y())));
+ }
+ }
+ std::sort(centers.begin(),centers.end());
+
+ for (int i=1;i<centers.size()-1;i++){
+ multiscaleLaplaceLevel* m1 = centers[i-1].second;
+ multiscaleLaplaceLevel* m2 = centers[i].second;
+ multiscaleLaplaceLevel* m3 = centers[i+1].second;
+ if (m2){
+ a1 = atan2 (m2->center.y() - centers[i-1].first.y(),m2->center.x() - centers[i-1].first.x());
+ a2 = atan2 (m2->center.y() - centers[i+1].first.y(),m2->center.x() - centers[i+1].first.x());
+ recur_compute_centers_ (m2->radius, a1, a2, m2);
+ }
+ }
+}
+
+static void recur_cut_edges_ (multiscaleLaplaceLevel * root,
+ std::multimap<MEdge,MElement*,Less_Edge> &e2e,
+ std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
+ std::set<MVertex*> &cutVertices){
+ std::set<MEdge,Less_Edge> allEdges;
+ for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.begin();
+ it != e2e.end() ; ++it){
+ allEdges.insert(it->first);
+ }
+
+ const double EPS = 0.001;
+
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
+ for (int i=0;i<centers.size()-1;i++){
+ SPoint2 p1 = centers[i].first;
+ SPoint2 p2 = centers[i+1].first;
+ for (std::set <MEdge,Less_Edge>::iterator it = allEdges.begin();
+ it != allEdges.end() ; ++it){
+ if(e2e.count(*it) == 2 && cutEdges.find(*it) == cutEdges.end()){
+ std::map<MVertex *, SPoint2>::iterator it0 = root->coordinates.find(it->getVertex(0));
+ std::map<MVertex *, SPoint2>::iterator it1 = root->coordinates.find(it->getVertex(1));
+ if (it0 != root->coordinates.end() && it1 != root->coordinates.end()){
+ SPoint2 q1 = root->coordinates[it->getVertex(0)];
+ SPoint2 q2 = root->coordinates[it->getVertex(1)];
+ double x[2];
+ int inters = intersection_segments (p1,p2,q1,q2,x);
+ if (inters && x[1] > EPS && x[1] < 1.-EPS){
+ // printf("%g %g -- %g %g -- %g %g -- %g %g\n",p1.x(),p1.y(),p2.x(),p2.y(),q1.x(),q1.y(),q2.x(),q2.y());
+ MVertex *newv = new MVertex ((1.-x[1])*it->getVertex(0)->x() + x[1]*it->getVertex(1)->x(),
+ (1.-x[1])*it->getVertex(0)->y() + x[1]*it->getVertex(1)->y(),
+ (1.-x[1])*it->getVertex(0)->z() + x[1]*it->getVertex(1)->z());
+ cutEdges[*it] = newv;
+ root->coordinates[newv] = q1*(1.-x[1]) + q2*x[1] ;
+ }
+ else if (inters && x[1] <= EPS)cutVertices.insert(it->getVertex(0));
+ else if (inters && x[1] >= 1.-EPS)cutVertices.insert(it->getVertex(1));
+ }
+ }
+ }
+ }
+ for (int i=0;i<centers.size();i++){
+ multiscaleLaplaceLevel* m2 = centers[i].second;
+ if (m2){
+ recur_cut_edges_ (m2,e2e,cutEdges,cutVertices);
+ }
+ }
+}
+
+static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
+ std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
+ std::set<MVertex*> &cutVertices,
+ std::set<MEdge,Less_Edge> &theCut,
+ std::vector<MElement*> &_all){
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
+ std::vector<MElement*> newElements;
+ for (int i=0;i<root->elements.size();i++){
+ MVertex *c[3] = {0,0,0};
+ for (int j=0;j<3;j++){
+ MEdge ed = root->elements[i]->getEdge(j);
+ // if (cutVertices.find (ed.getVertex(0)) != cutVertices.end() &&
+ // cutVertices.find (ed.getVertex(1)) != cutVertices.end() )theCut.insert(ed);
+
+ std::map<MEdge,MVertex*,Less_Edge> :: iterator it = cutEdges.find(ed);
+ if (it != cutEdges.end()){
+ c[j] = it->second;
+ }
+ }
+ if (c[0] && c[1]){
+ newElements.push_back(new MTriangle (c[0],root->elements[i]->getVertex(1),c[1]));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],c[1]));
+ theCut.insert(MEdge(c[0],c[1]));
+ // FIXME should be done !!!!!!
+ //delete root->elements[i];
+ }
+ else if (c[0] && c[2]){
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],c[2]));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[0],root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],c[2]));
+ theCut.insert(MEdge(c[0],c[2]));
+ //delete root->elements[i];
+ }
+ else if (c[1] && c[2]){
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[2],c[1]));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),root->elements[i]->getVertex(1),c[2]));
+ newElements.push_back(new MTriangle (c[2],root->elements[i]->getVertex(1),c[1]));
+ theCut.insert(MEdge(c[1],c[2]));
+ //delete root->elements[i];
+ }
+ else if (c[0]){
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],root->elements[i]->getVertex(1)));
+ if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end())
+ theCut.insert(MEdge(c[0],root->elements[i]->getVertex(0)));
+ else if (cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end())
+ theCut.insert(MEdge(c[0],root->elements[i]->getVertex(1)));
+ else
+ theCut.insert(MEdge(c[0],root->elements[i]->getVertex(2)));
+ //delete root->elements[i];
+ }
+ else if (c[1]){
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[1],root->elements[i]->getVertex(0)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[1],root->elements[i]->getVertex(2)));
+ if (cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end())
+ theCut.insert(MEdge(c[1],root->elements[i]->getVertex(1)));
+ else if (cutVertices.find (root->elements[i]->getVertex(2)) != cutVertices.end())
+ theCut.insert(MEdge(c[1],root->elements[i]->getVertex(2)));
+ else
+ theCut.insert(MEdge(c[1],root->elements[i]->getVertex(0)));
+ //delete root->elements[i];
+ }
+ else if (c[2]){
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[2],root->elements[i]->getVertex(1)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[2],root->elements[i]->getVertex(2)));
+ if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end())
+ theCut.insert(MEdge(c[2],root->elements[i]->getVertex(0)));
+ else if (cutVertices.find (root->elements[i]->getVertex(2)) != cutVertices.end())
+ theCut.insert(MEdge(c[2],root->elements[i]->getVertex(2)));
+ else
+ theCut.insert(MEdge(c[2],root->elements[i]->getVertex(1)));
+ //delete root->elements[i];
+ }
+ else newElements.push_back(root->elements[i]);
+ }
+ root->elements = newElements;
+ _all.insert(_all.end(),newElements.begin(),newElements.end());
+ for (int i=0;i<centers.size();i++){
+ multiscaleLaplaceLevel* m2 = centers[i].second;
+ if (m2){
+ recur_cut_elements_ (m2,cutEdges,cutVertices,theCut,_all);
+ }
+ }
+}
+
+static void recur_split_ (MElement *e,
+ std::multimap<MEdge,MElement*,Less_Edge> &e2e,
+ std::set<MElement*> &group,
+ std::set<MEdge,Less_Edge> &theCut){
+ if (group.find(e) != group.end())return;
+ group.insert(e);
+ for (int i=0;i<e->getNumEdges();i++){
+ MEdge ed = e->getEdge(i);
+ if (theCut.find(ed) == theCut.end()){
+ for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.lower_bound(ed);
+ it != e2e.upper_bound(ed) ; ++it){
+ if (it->second != e)recur_split_ (it->second,e2e,group,theCut);
+ }
+ }
+ }
+}
+
// starting form a list of elements, returns
// lists of lists that are all simply connected
-static void recur_connect (MVertex *v,
- std::multimap<MVertex*,MElement*> &v2e,
+static void recur_connect (const MEdge &e,
+ std::multimap<MEdge,MElement*,Less_Edge> &e2e,
std::set<MElement*> &group,
- std::set<MVertex*> &touched){
- if (touched.find(v) != touched.end())return;
- touched.insert(v);
- for (std::multimap <MVertex*,MElement*>::iterator it = v2e.lower_bound(v);
- it != v2e.upper_bound(v) ; ++it){
+ std::set<MEdge,Less_Edge> &touched){
+ if (touched.find(e) != touched.end())return;
+ touched.insert(e);
+ for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.lower_bound(e);
+ it != e2e.upper_bound(e) ; ++it){
group.insert(it->second);
- for (int i=0;i<it->second->getNumVertices();++i){
- recur_connect (it->second->getVertex(i),v2e,group,touched);
+ for (int i=0;i<it->second->getNumEdges();++i){
+ recur_connect (it->second->getEdge(i),e2e,group,touched);
}
}
}
+
static void connectedRegions (std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions)
{
- // build vertex 2 elements
- std::multimap<MVertex*,MElement*> v2e;
+ std::multimap<MEdge,MElement*,Less_Edge> e2e;
for (int i=0;i<elements.size();++i){
- for (int j=0;j<elements[i]->getNumVertices();j++){
- v2e.insert(std::make_pair(elements[i]->getVertex(j),elements[i]));
+ for (int j=0;j<elements[i]->getNumEdges();j++){
+ e2e.insert(std::make_pair(elements[i]->getEdge(j),elements[i]));
}
}
- while (!v2e.empty()){
+ while (!e2e.empty()){
std::set<MElement*> group;
- std::set<MVertex*> touched;
- recur_connect (v2e.begin()->first,v2e,group,touched);
+ std::set<MEdge,Less_Edge> touched;
+ recur_connect (e2e.begin()->first,e2e,group,touched);
std::vector<MElement*> temp;
temp.insert(temp.begin(), group.begin(), group.end());
regions.push_back(temp);
- for ( std::set<MVertex*>::iterator it = touched.begin() ; it != touched.end();++it)
- v2e.erase(*it);
+ for ( std::set<MEdge>::iterator it = touched.begin() ; it != touched.end();++it)
+ e2e.erase(*it);
}
}
-
static void printLevel ( const char* fn,
std::vector<MElement *> &elements,
- std::map<MVertex*,SPoint2> &coordinates,
- double version, bool param){
+ std::map<MVertex*,SPoint2> *coordinates,
+ double version){
+
+
+ std::set<MVertex*> vs;
+ for (int i=0;i<elements.size();i++){
+ for (int j=0;j<elements[i]->getNumVertices();j++){
+ vs.insert(elements[i]->getVertex(j));
+ }
+ }
bool binary = false;
FILE *fp = fopen (fn, "w");
@@ -60,14 +302,15 @@ static void printLevel ( const char* fn,
fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
fprintf(fp, "$EndMeshFormat\n");
- fprintf(fp,"$Nodes\n%d\n",coordinates.size());
- std::map<MVertex*,SPoint2> :: iterator it = coordinates.begin();
+ fprintf(fp,"$Nodes\n%d\n",vs.size());
+ std::set<MVertex*> :: iterator it = vs.begin();
int index = 1;
- for (; it != coordinates.end() ; ++it){
- it->first->setIndex(index++);
- if (param) fprintf(fp,"%d %g %g 0\n",it->first->getIndex(),it->second[0],it->second[1]);
- else fprintf(fp,"%d %g %g %g\n",it->first->getIndex(),
- it->first->x(),it->first->y(),it->first->z());
+ for (; it != vs.end() ; ++it){
+ (*it)->setIndex(index++);
+ SPoint2 p = (coordinates) ? (*coordinates)[*it] : SPoint2(0,0);
+ if (coordinates) fprintf(fp,"%d %g %g 0\n",(*it)->getIndex(),p.x(),p.y());
+ else fprintf(fp,"%d %g %g %g\n",(*it)->getIndex(),
+ (*it)->x(),(*it)->y(),(*it)->z());
}
fprintf(fp,"$EndNodes\n",elements.size());
@@ -80,11 +323,11 @@ static void printLevel ( const char* fn,
fclose(fp);
}
-
multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
std::vector<MVertex*> &boundaryNodes,
std::vector<double> &linearAbscissa) {
+ if (!CTX::instance()->mesh.smoothInternalEdges)return;
#if defined(HAVE_TAUCS)
_lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_GMM)
@@ -95,17 +338,88 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
_lsys = new linearSystemFull<double>;
#endif
- multiscaleLaplaceLevel *level = new multiscaleLaplaceLevel;
- level->elements.insert(level->elements.begin(),elements.begin(),elements.end());
+ root = new multiscaleLaplaceLevel;
+ root->elements = elements;
for(unsigned int i = 0; i < boundaryNodes.size(); i++){
MVertex *v = boundaryNodes[i];
const double theta = 2 * M_PI * linearAbscissa[i];
- level->coordinates[v] = SPoint2(cos(theta),sin(theta));
+ root->coordinates[v] = SPoint2(cos(theta),sin(theta));
+ }
+ root->recur = 0;
+ root->region = 0;
+ parametrize(*root);
+
+ std::multimap<MEdge,MElement*,Less_Edge> e2e;
+ for (int i=0;i<elements.size();++i){
+ for (int j=0;j<elements[i]->getNumEdges();j++){
+ e2e.insert(std::make_pair(elements[i]->getEdge(j),elements[i]));
+ }
+ }
+ std::map<MEdge,MVertex*,Less_Edge> cutEdges;
+ std::set<MVertex*> cutVertices;
+ recur_compute_centers_ (1.0,0, M_PI, root);
+ recur_cut_edges_ (root,e2e,cutEdges,cutVertices);
+
+ if (1){
+ // DEBUG : EXPORT -----------------
+ std::map<MEdge,MVertex*,Less_Edge>::iterator it = cutEdges.begin();
+ FILE *f = fopen("points.pos","w");
+ fprintf(f,"View\"\"{\n");
+ for ( ; it != cutEdges.end();++it){
+ fprintf(f,"SP(%g,%g,%g){1.0};\n",it->second->x(),it->second->y(),it->second->z());
+ }
+ fprintf(f,"};\n");
+ fclose(f);
+ // ------END DEBUG
}
- levels.push_back(level);
- level->recur = 0;
- level->region = 0;
- parametrize(*level);
+ std::set<MEdge,Less_Edge> theCut;
+ std::vector<MElement*> _all;
+ recur_cut_elements_ (root,cutEdges,cutVertices,theCut,_all);
+ if (1){
+ // DEBUG : EXPORT -----------------
+ std::set<MEdge,Less_Edge>::iterator it = theCut.begin();
+ FILE *f = fopen("edges.pos","w");
+ fprintf(f,"View\"\"{\n");
+ for ( ; it != theCut.end();++it){
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g){1.0,1.0};\n",it->getVertex(0)->x(),it->getVertex(0)->y(),it->getVertex(0)->z(),
+ it->getVertex(1)->x(),it->getVertex(1)->y(),it->getVertex(1)->z());
+ }
+ fprintf(f,"};\n");
+ fclose(f);
+ // ------END DEBUG
+ }
+ e2e.clear();
+ for (int i=0;i<_all.size();++i){
+ for (int j=0;j<_all[i]->getNumEdges();j++){
+ e2e.insert(std::make_pair(_all[i]->getEdge(j),_all[i]));
+ }
+ }
+ // std::set<MElement*> leftSet;
+ // recur_split_ (_all[0],e2e,leftSet,theCut);
+
+ std::vector<MElement*> left,right;
+ /*
+ for (int i=0;i<_all.size();i++){
+ if(leftSet.find(_all[i]) == leftSet.end())right.push_back(_all[i]);
+ else left.push_back(_all[i]);
+ }
+ */
+ cut (left,right);
+
+ printLevel ("multiscale_left.msh",left,0,1.0);
+ printLevel ("multiscale_right.msh",right,0,1.0);
+ printLevel ("multiscale_all.msh",_all,0,1.0);
+ // FIXME !!!!!
+ throw;
+}
+
+static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
+ SBoundingBox3d local;
+ for(unsigned int j = 0; j<e->getNumVertices(); ++j){
+ SPoint2 p = solution[e->getVertex(j)];
+ local += SPoint3(p.x(),p.y(),0.0);
+ }
+ return local.max().distance(local.min());
}
void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
@@ -161,44 +475,59 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
double global_size = bbox.max().distance(bbox.min());
// check elements that are too small
std::vector<MElement*> tooSmall, goodSize;
- std::set<MVertex*> goodSizev;
+
for(unsigned int i = 0; i < level.elements.size(); ++i){
MElement *e = level.elements[i];
std::vector<SPoint2> localCoord;
- SBoundingBox3d local;
- for(unsigned int j = 0; j<e->getNumVertices(); ++j){
- SPoint2 p = solution[e->getVertex(j)];
- local += SPoint3(p.x(),p.y(),0.0);
- }
- double local_size = local.max().distance(local.min());
+ double local_size = localSize(e,solution);
if (local_size < 1.e-6 * global_size){
tooSmall.push_back(e);
}
else{
goodSize.push_back(e);
- for(unsigned int j = 0; j<e->getNumVertices(); ++j){
- goodSizev.insert(e->getVertex(j));
- level.coordinates[e->getVertex(j)] = solution[e->getVertex(j)];
- }
}
}
// only keep the right elements and nodes
level.elements = goodSize;
+ std::vector<std::vector<MElement*> > regions, regions_;
+ connectedRegions (tooSmall,regions_);
+
+ // remove small regions
+
+ for (int i=0;i< regions_.size() ; i++){
+ bool region_has_really_small_elements = false;
+ for (int k=0; k<regions_[i].size() ; k++){
+ MElement *e = regions_[i][k];
+ double local_size = localSize(e,solution);
+ if (local_size < 1.e-8 * global_size){
+ region_has_really_small_elements = true;
+ }
+ }
+ if(region_has_really_small_elements) regions.push_back(regions_[i]);
+ else level.elements.insert(level.elements.begin(), regions_[i].begin(), regions_[i].end() );
+ }
+
+ std::set<MVertex*> goodSizev;
+ for(unsigned int i = 0; i < level.elements.size(); ++i){
+ MElement *e = level.elements[i];
+ for(unsigned int j = 0; j<e->getNumVertices(); ++j){
+ goodSizev.insert(e->getVertex(j));
+ level.coordinates[e->getVertex(j)] = solution[e->getVertex(j)];
+ }
+ }
+
//END PARAMETRIZE ---------------------------------
// DEBUG
char name[245];
- sprintf(name,"multiscale_level%d_region%d_param.msh",level.recur, level.region);
- printLevel (name,level.elements,level.coordinates,1.0,true);
sprintf(name,"multiscale_level%d_region%d_real.msh",level.recur, level.region);
- printLevel (name,level.elements,level.coordinates,1.0,false);
+ printLevel (name,level.elements,0,1.0);
+ sprintf(name,"multiscale_level%d_region%d_param.msh",level.recur, level.region);
+ printLevel (name,level.elements,&level.coordinates,1.0);
// END DEBUG
- std::vector<std::vector<MElement*> > regions;
- connectedRegions (tooSmall,regions);
-
Msg::Info("%d connected regions\n",regions.size());
-
+
for (int i=0;i< regions.size() ; i++){
// check nodes that appear both on too small and good size
// and set it to the next level
@@ -233,7 +562,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
// recursively continue if tooSmall is not empty
if (!tooSmall.empty()){
Msg::Info("Multiscale Laplace, level %d region %d, %d too small",level.recur,level.region,tooSmall.size());
- levels.push_back(nextLevel);
+ level.childeren.push_back(nextLevel);
parametrize (*nextLevel);
}
else {
@@ -242,3 +571,86 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
}
}
}
+
+static void recur_cut_ (double R, double a1, double a2,
+ multiscaleLaplaceLevel * root,
+ std::vector<MElement *> &left,
+ std::vector<MElement *> &right ){
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
+ centers.clear();
+ multiscaleLaplaceLevel* zero = 0;
+
+ SPoint2 PL (R*cos(a1),R*sin(a1));
+ SPoint2 PR (R*cos(a2),R*sin(a2));
+ centers.push_back(std::make_pair(PL,zero));
+ centers.push_back(std::make_pair(PR,zero));
+ for (int i=0;i<root->childeren.size();i++){
+ centers.push_back(std::make_pair(root->childeren[i]->center,root->childeren[i]));
+ multiscaleLaplaceLevel* m = root->childeren[i];
+ m->radius = 0.0;
+ for (std::map<MVertex*,SPoint2>::iterator it = m->coordinates.begin();
+ it != m->coordinates.end() ; ++it){
+ const SPoint2 &p = it->second;
+ m->radius = std::max(m->radius,sqrt ((m->center.x() - p.x())*(m->center.x() - p.x())+
+ (m->center.y() - p.y())*(m->center.y() - p.y())));
+ }
+ }
+ std::sort(centers.begin(),centers.end());
+ /*
+ printf("%d %d ",root->recur,root->region);
+ for (int j=0;j<centers.size();j++){
+ printf("(%g %g) ",centers[j].first.x(),centers[j].first.y());
+ }
+ printf("\n");
+ */
+ double d = sqrt((PL.x()-PR.x())*(PL.x()-PR.x())+
+ (PL.y()-PR.y())*(PL.y()-PR.y()));
+ SPoint2 farLeft (0.5*(PL.x()+PR.x()) - (PR.y()-PL.y())/d ,
+ 0.5*(PL.y()+PR.y()) + (PR.x()-PL.x())/d );
+
+ for (int i=0;i<root->elements.size();i++){
+ SPoint2 pp (0,0);
+ for (int j=0; j<root->elements[i]->getNumVertices(); j++){
+ pp += root->coordinates[root->elements[i]->getVertex(j)];
+ }
+ pp *= 1./(double)root->elements[i]->getNumVertices();
+ int nbIntersect = 0;
+ for (int j=0;j<centers.size()-1;j++){
+ double x[2];
+ nbIntersect += intersection_segments (centers[j].first,centers[j+1].first,pp,farLeft,x);
+ }
+ if (root->recur != 0){
+ if (nbIntersect %2 == 0)
+ left.push_back(root->elements[i]);
+ else
+ right.push_back(root->elements[i]);
+ }
+ else{
+ if (nbIntersect %2 != 0)
+ left.push_back(root->elements[i]);
+ else
+ right.push_back(root->elements[i]);
+ }
+ }
+
+ for (int i=1;i<centers.size()-1;i++){
+ multiscaleLaplaceLevel* m1 = centers[i-1].second;
+ multiscaleLaplaceLevel* m2 = centers[i].second;
+ multiscaleLaplaceLevel* m3 = centers[i+1].second;
+ if (m2){
+ /*center of the local system is always 0,0
+ its relative position to its parent is center
+ only 2 angles have to be computed for in and out*/
+ a1 = atan2 (m2->center.y() - centers[i-1].first.y(),m2->center.x() - centers[i-1].first.x());
+ a2 = atan2 (m2->center.y() - centers[i+1].first.y(),m2->center.x() - centers[i+1].first.x());
+ recur_cut_ (m2->radius, a1, a2, m2, left, right);
+ }
+ }
+}
+
+
+void multiscaleLaplace::cut (std::vector<MElement *> &left,
+ std::vector<MElement *> &right)
+{
+ recur_cut_ (1.0,0, M_PI, root,left,right);
+}
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index 1661c692610938d89b5aa80fc5bda72bb47fdde7..90f662d261c658a7b957a3e244f39c8a5080165b 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -4,7 +4,7 @@
#include <vector>
#include <map>
#include "SPoint2.h"
-#include "linearSystemGMM.h"
+#include "linearSystem.h"
class MElement;
class MVertex;
@@ -12,20 +12,23 @@ class MVertex;
struct multiscaleLaplaceLevel {
SPoint2 center;
double scale;
+ double radius;
int recur,region;
+ std::vector<multiscaleLaplaceLevel*> childeren;
std::vector<MElement *> elements;
std::map<MVertex*,SPoint2> coordinates;
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > cut;
};
class multiscaleLaplace{
linearSystem<double> *_lsys;
- std::vector<multiscaleLaplaceLevel*> levels;
- void parametrize (multiscaleLaplaceLevel &);
-
+ multiscaleLaplaceLevel* root;
+ void parametrize (multiscaleLaplaceLevel &);
public:
multiscaleLaplace (std::vector<MElement *> &elements,
std::vector<MVertex*> &boundaryNodes,
std::vector<double> &linearAbscissa) ;
-
+ void cut (std::vector<MElement *> &left,
+ std::vector<MElement *> &right);
};
#endif
diff --git a/benchmarks/extrude/Cylinder.geo b/benchmarks/extrude/Cylinder.geo
index 55323ba11083895f75a1ac6e8b404a2c48d2d633..3a1f9aeed0f263458475a256350abb30a4558ae6 100644
--- a/benchmarks/extrude/Cylinder.geo
+++ b/benchmarks/extrude/Cylinder.geo
@@ -1,175 +1,16 @@
-// Cylinder //
-//----------//
-
-Mesh.Algorithm = 2;
-
-N = 0.5;
-
-bl = 0.025;
-NBL = 10*N;
-NBL_PROG = 1.1;
-N_CYL = 100*N;
-CYL_PROG = 0.99;
-Delta_0 = 0.5*3.1415*(0.5+bl)*(CYL_PROG-1)/(CYL_PROG*(CYL_PROG^N_CYL-1));
-Delta_N = Delta_0*CYL_PROG^N_CYL;
-h_square = 0.75;
-l_rectangle = 10.;
-m_box = 0.02/N;
-N_LARGEBOX = 11*N;
-L_LARGEBOX = 100;
-
-// CYLINDER //
-
-Point(11) = { 0, 0,0,1};
-Point(12) = { 0, 0.5,0,1};
-Point(13) = { 0.5, 0,0,1};
-Point(14) = { 0,-0.5,0,1};
-Point(15) = {-0.5, 0,0,1};
-
-Circle(21) = {12,11,13};
-Circle(22) = {13,11,14};
-Circle(23) = {14,11,15};
-Circle(24) = {15,11,12};
-
-// 15 , 24 //
-
-
-// BOUNDARY LAYER //
-
-Point(16) = { 0, 0.5+bl,0,Delta_0};
-Point(17) = { 0.5+bl, 0,0,Delta_N};
-Point(18) = { 0,-0.5-bl,0,Delta_0};
-Point(19) = {-0.5-bl, 0,0,Delta_N};
-
-Line(25) = {12,16};
-Line(26) = {13,17};
-Line(27) = {14,18};
-Line(28) = {15,19};
-
-Circle(29) = {16,11,17};
-Circle(210)= {17,11,18};
-Circle(211)= {18,11,19};
-Circle(212)= {19,11,16};
-
-Line Loop(213) = {-25,21,26,-29};
-Line Loop(214) = {-26,22,27,-210};
-Line Loop(215) = {-27,23,28,-211};
-Line Loop(216) = {-28,24,25,-212};
-
-Ruled Surface(31) = {213};
-Ruled Surface(32) = {214};
-Ruled Surface(33) = {215};
-Ruled Surface(34) = {216};
-
-Transfinite Line{29,-210,211,-212,21,-22,23,-24} = Floor(N_CYL) Using Progression CYL_PROG;
-Transfinite Line{25,26,27,28} = Floor(NBL) Using Progression NBL_PROG;
-
-Transfinite Surface{31} = {16,12,13,17};
-Transfinite Surface{32} = {17,13,14,18};
-Transfinite Surface{33} = {18,14,15,19};
-Transfinite Surface{34} = {19,15,12,16};
-
-Recombine Surface{31,32,33,34};
-
-// 19 , 216 , 34 //
-
-
-// SQUARE //
-
-Point(111)={-h_square,-h_square,0,m_box*3};
-Point(112)={-h_square, h_square,0,m_box*3};
-Point(113)={ h_square, h_square,0,m_box};
-Point(114)={ h_square,-h_square,0,m_box};
-
-Line(217) = {111,112};
-Line(218) = {112,113};
-Line(219) = {113,114};
-Line(220) = {114,111};
-
-Line Loop(221) = {-217,-218,-219,-220};
-Line Loop(222) = {-29,-210,-211,-212};
-
-Plane Surface(35) = {221,222};
-
-// 114 , 222 , 35 //
-
-
-// RECTANGLE //
-
-Point(115) = {h_square+l_rectangle,-h_square,0,m_box};
-Point(116) = {h_square+l_rectangle, h_square,0,m_box};
-
-Line(223) = {114,115};
-Line(224) = {115,116};
-Line(225) = {116,113};
-
-Line Loop(226) = {223,224,225,219};
-
-Ruled Surface(36) = {226};
-
-Transfinite Line{224,219} = Floor(2*h_square/m_box) Using Progression 1;
-Transfinite Line{223,225} = Floor(l_rectangle/m_box) Using Progression 1;
-
-Transfinite Surface{36} = {114,115,116,113};
-Recombine Surface{36};
-
-// 116 , 226 , 36 //
-
-
-
-// EXTERIOR //
-
-Point(117) = {-L_LARGEBOX,-L_LARGEBOX,0,2*L_LARGEBOX/N_LARGEBOX};
-Point(118) = { L_LARGEBOX,-L_LARGEBOX,0,2*L_LARGEBOX/N_LARGEBOX};
-Point(119) = { L_LARGEBOX, L_LARGEBOX,0,2*L_LARGEBOX/N_LARGEBOX};
-Point(120) = {-L_LARGEBOX, L_LARGEBOX,0,2*L_LARGEBOX/N_LARGEBOX};
-
-Line(227) = {117,118};
-Line(228) = {118,119};
-Line(229) = {119,120};
-Line(230) = {120,117};
-
-Line(231) = {117,111};
-Line(232) = {118,115};
-Line(233) = {119,116};
-Line(234) = {120,112};
-
-
-Line Loop(235) = {-231,227,232,-223,220};
-Line Loop(236) = {-232,228,233,-224};
-Line Loop(237) = {-233,229,234,218,-225};
-Line Loop(238) = {-234,230,231,217};
-
-Plane Surface(37) = {235};
-Plane Surface(38) = {236};
-Plane Surface(39) = {237};
-Plane Surface(310) = {238};
-
-Recombine Surface{37,38,39,310};
-
-// 120 , 238 , 310 //
-
-
-// EXTRUSION //
-
-out1[] = Extrude{0,0,0.1}{Surface{31}; Layers{4,1} ; Recombine;};
-out2[] = Extrude{0,0,0.1}{Surface{32}; Layers{4,1} ; Recombine;};
-out3[] = Extrude{0,0,0.1}{Surface{33}; Layers{4,1} ; Recombine;};
-out4[] = Extrude{0,0,0.1}{Surface{34}; Layers{4,1} ; Recombine;};
-out5[] = Extrude{0,0,0.1}{Surface{35}; Layers{4,1} ; Recombine;};
-out6[] = Extrude{0,0,0.1}{Surface{36}; Layers{4,1} ; Recombine;};
-out7[] = Extrude{0,0,0.1}{Surface{37}; Layers{4,1} ; Recombine;};
-out8[] = Extrude{0,0,0.1}{Surface{38}; Layers{4,1} ; Recombine;};
-out9[] = Extrude{0,0,0.1}{Surface{39}; Layers{4,1} ; Recombine;};
-out10[] = Extrude{0,0,0.1}{Surface{310}; Layers{4,1} ; Recombine;};
-
-Physical Volume(1) = {out1[1],out2[1],out3[1],out4[1],out5[1],out6[1],out7[1],out8[1],out9[1],out10[1]};
-Physical Surface(2) = {31,32,33,34,35,36,37,38,39,310}; // Z_DOWN //
-Physical Surface(3) = {out1[0],out2[0],out3[0],out4[0],out5[0],out6[0],out7[0],out8[0],out9[0],out10[0]}; // Z_UP //
-Physical Surface(4) = {476}; // Y_DOWN //
-Physical Surface(5) = {525}; // Y_UP //
-Physical Surface(6) = {551}; // INLET //
-Physical Surface(7) = {502}; // OUTLET //
-Physical Surface(8) = {389,323,345,367}; // CYLINDER //
-
+Point(1) = {0, 0, 0};
+Point(2) = {1, 0, 0};
+Point(3) = {0, 1, 0};
+Point(4) = {0, -1, 0};
+Point(5) = {-1, 0, 0};
+Circle(1) = {3, 1, 2};
+Circle(2) = {2, 1, 4};
+Circle(3) = {4, 1, 5};
+Circle(4) = {5, 1, 3};
+Extrude {0, 0, 13} {
+ Line{4, 1, 2, 3};
+}
+
+
+Compound Surface(10000) = {20, 16, 12, 8};
diff --git a/benchmarks/testsuite/testsuite.sh b/benchmarks/testsuite/testsuite.sh
index bb63797faf1f1db081d55cfb5c593ce866ac3855..22411f4087419b006b31b1aa113b2a6cca96701c 100755
--- a/benchmarks/testsuite/testsuite.sh
+++ b/benchmarks/testsuite/testsuite.sh
@@ -26,3 +26,5 @@ ${GMSH} Cone_1.brep -clscale 1 -2 -v 1 > /dev/null
${GMSH} Cylinder_1.brep -clscale 1 -2 -v 1 > /dev/null
${GMSH} Torus_1.brep -clscale 300 -2 -v 1 > /dev/null
${GMSH} A319.geo -clscale 3 -2 -v 1 > /dev/null
+${GMSH} capot.geo -clscale .2 -2 -saveall -v 1 > /dev/null
+${GMSH} bouda_.geo -clscale .2 -2 -saveall -v 1 > /dev/null
diff --git a/utils/api_demos/Makefile b/utils/api_demos/Makefile
index e0475823274e65331bb4cf33b7d2f59ad3294916..929e1bc7a0de5c7b4668283a460053d81a947908 100644
--- a/utils/api_demos/Makefile
+++ b/utils/api_demos/Makefile
@@ -6,10 +6,18 @@ GLUT = -framework OpenGL -framework GLUT -framework Cocoa -framework Application
OCC = -L/usr/local/opencascade/lib -L/Users/remacle/SOURCES/opencascade/lib -lTKSTEP -lTKSTEP209 -lTKSTEPAttr -lTKSTEPBase -lTKIGES -lTKXSBase -lTKOffset -lTKFeat -lTKFillet -lTKBool -lTKShHealing -lTKMesh -lTKHLR -lTKBO -lTKPrim -lTKTopAlgo -lTKGeomAlgo -lTKBRep -lTKGeomBase -lTKG3d -lTKG2d -lTKAdvTools -lTKMath -lTKernel
TAUCS = -L/Users/remacle/SOURCES/gmsh/contrib/taucs/lib -L/home/gaetan/develop/gmsh/contrib/taucs/lib/linux -ltaucs
+mainLua: mainLua.cpp
+ g++ ${FLAGS} -I../../Geo/ -I../../lib/Common -I../../Common -DHAVE_LUA -o mainLua ${INC} mainLua.cpp\
+ ${GMSH} -llapack -lblas -lm ${TAUCS} ${OCC} -llua
+
mainElasticity: mainElasticity.cpp
g++ ${FLAGS} -o mainElasticity ${INC} mainElasticity.cpp\
${GMSH} -llapack -lblas -lm -lmetis ${TAUCS}
+mainFDPlate: mainFDPlate.cpp
+ g++ ${FLAGS} -o mainFDPlate ${INC} mainFDPlate.cpp\
+ ${GMSH} -llapack -lblas -lm ${TAUCS} ${OCC}
+
mainCartesian: mainCartesian.cpp
g++ ${FLAGS} -o mainCartesian ${INC} mainCartesian.cpp\
${GMSH} -llapack -lblas -lm ${OCC} ${GLUT}