diff --git a/Common/Bindings.h b/Common/Bindings.h new file mode 100644 index 0000000000000000000000000000000000000000..031c3dc40027e0f1d6ef42a50bd19cba07cec82c --- /dev/null +++ b/Common/Bindings.h @@ -0,0 +1,457 @@ +#ifndef _LUNA_SIGNATURE_H_ +#define _LUNA_SIGNATURE_H_ + +extern "C" { +#include "lua.h" +#include "lauxlib.h" +} +#include <vector> + + +class binding { +}; + +class methodBinding { + public: + std::string _luaname; + virtual int call (lua_State *L)=0; + methodBinding(const std::string luaname){ + _luaname=luaname; + } +}; + +class constructorBinding { + public: + virtual int call (lua_State *L)=0; +}; + +template <typename T> class classBinding { + typedef struct { T *pT; bool owned;} userdataType; +public: + static void push(lua_State *L, T* obj, bool owned) { + userdataType *ud = static_cast<userdataType*>(lua_newuserdata(L, sizeof(userdataType))); + ud->pT=obj; + ud->owned=owned; + luaL_getmetatable(L, T::className); // lookup metatable in Lua registry + lua_setmetatable(L, -2); + } + static void Register(lua_State *L) { + lua_newtable(L); + int methods = lua_gettop(L); //-- create methods , i.e. the class itself + + luaL_newmetatable(L, T::className); + int metatable = lua_gettop(L); // -- create metatable of methods + + // 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); // -- global[T::className] = methods + +// 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); // metatable.__index=methods + + 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_pushlightuserdata(L,(void*)T::constructorMethod); + lua_pushcclosure(L, callConstructor, 1); + + 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 + + if(T::parentClassName!=std::string("")){ + lua_pushliteral(L,"__index"); + lua_pushstring(L,T::parentClassName); + lua_gettable(L,LUA_GLOBALSINDEX); + lua_settable(L,mt); // mt.__index = global[T::parentClassName] // this is the inheritance bit + } + lua_setmetatable(L, methods); // setmetatable(methods, mt) + + // fill method table with methods from class T + for (methodBinding **l = T::methods; *l; l++) { + lua_pushstring(L, (*l)->_luaname.c_str()); + lua_pushlightuserdata(L, (void*)*l); + lua_pushcclosure(L, callMethod, 1); + lua_settable(L, methods); //methods.(l->name) = callMethod(l) + } + 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*>(lua_touserdata(L, narg)); + if(!ud) luaL_typerror(L, narg, T::className); + return ud->pT; // pointer to T object*/ + } + +private: + classBinding(); // hide default constructor + + // member function dispatcher + static int callMethod(lua_State *L) { + methodBinding *l = static_cast<methodBinding*>(lua_touserdata(L, lua_upvalueindex(1))); + return l->call(L); // call member function + } + static int callConstructor(lua_State *L) { + constructorBinding *l = static_cast<constructorBinding*>(lua_touserdata(L, lua_upvalueindex(1))); + if(!l){ + printf("this class does not have constructor\n"); + //throw(); + } + return l->call(L); // call member function + } + + // garbage collection metamethod + static int gc_T(lua_State *L) { + userdataType *ud = static_cast<userdataType*>(lua_touserdata(L, 1)); + /*if(ud->owned) //disable gc + delete ud->pT; // 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 *t = static_cast<T*>(lua_touserdata(L, 1)); + sprintf(buff, "%p", ud->pT); + lua_pushfstring(L, "%s (%s)", T::className, buff); + return 1; + } +}; + +template<class type> +class lua_template { + public: + static type get(lua_State *L, int ia){ + printf("error cannot get generic class in lua, only pointers are implemented\n"); + } + static void push(lua_State *L, type obj){ + printf("error cannot push generic class in lua, only pointers are implemented\n"); + } +}; + +template <class type> +class lua_template<type*>{ + public: + static type* get(lua_State *L, int ia){ + type *a=classBinding<type>::check(L,ia); + return a; + + } + static void push(lua_State *L, type *obj){ + classBinding<type>::push(L,obj,false); + } +}; + +template<> +class lua_template<lua_State *>{ + public: + static lua_State *get(lua_State *L, int ia){ + return L; + } + static void push(lua_State *L, int i){ + printf("error cannot push a lua_State in lua\n"); + } +}; + +template<> +class lua_template<int>{ + public: + static int get(lua_State *L, int ia){ + int a= luaL_checkint(L,ia); + return a; + } + static void push(lua_State *L, int i){ + lua_pushinteger(L,i); + } +}; + +template<class type> +class lua_template<std::vector<type > >{ + public: + static std::vector<type> get(lua_State *L, int ia){ + std::vector<type> v; + size_t size=lua_objlen(L,ia); + v.resize(size); + for(size_t i=0;i<size;i++){ + lua_pushinteger(L,i+1); + lua_gettable(L,ia); + v[i]=lua_template<type>::get(L,-1); + lua_pop(L,1); + } + return v; + } + static void push(lua_State *L, std::vector<type> a){ + } +}; + +template<> +class lua_template<double>{ + public: + static double get(lua_State *L, int ia){ + return luaL_checknumber(L,ia); + } + static void push(lua_State *L, double v){ + lua_pushnumber(L,v); + } +}; + +template<> +class lua_template<std::string>{ + public: + static std::string get(lua_State *L, int ia){ + return luaL_checkstring(L,ia); + } + static void push(lua_State *L, std::string s){ + lua_pushstring(L,s.c_str()); + } +}; + +//full : 4 args with return +template <class objectType, class returnType=void, class arg0Type=void, class arg1Type=void, class arg2Type=void, class arg3Type=void> +class methodBindingTemplate:public methodBinding { + typedef returnType (objectType::*callback)(arg0Type,arg1Type,arg2Type,arg3Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + arg1Type a1 = lua_template<arg1Type>::get(L,3); + arg2Type a2 = lua_template<arg2Type>::get(L,4); + arg3Type a3 = lua_template<arg3Type>::get(L,5); + returnType r=(obj->*(_f))(a0,a1,a2,a3); + lua_template<returnType>::push(L,r); + return 1; + } +}; +//3 args with return +template <class objectType, class returnType, class arg0Type, class arg1Type, class arg2Type> +class methodBindingTemplate<objectType,returnType,arg0Type,arg1Type,arg2Type,void>:public methodBinding { + typedef returnType (objectType::*callback)(arg0Type,arg1Type,arg2Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + arg1Type a1 = lua_template<arg1Type>::get(L,3); + arg2Type a2 = lua_template<arg2Type>::get(L,4); + returnType r=(obj->*(_f))(a0,a1,a2); + lua_template<returnType>::push(L,r); + return 1; + } +}; +//2 args with return +template <class objectType, class returnType, class arg0Type, class arg1Type> +class methodBindingTemplate<objectType,returnType,arg0Type,arg1Type,void,void>:public methodBinding { + typedef returnType (objectType::*callback)(arg0Type,arg1Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + arg1Type a1 = lua_template<arg1Type>::get(L,3); + returnType r=(obj->*(_f))(a0,a1); + lua_template<returnType>::push(L,r); + return 1; + } +}; +//1 arg with return +template <class objectType, class returnType, class arg0Type> +class methodBindingTemplate<objectType,returnType,arg0Type,void,void,void>:public methodBinding { + typedef returnType (objectType::*callback)(arg0Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + returnType r=(obj->*(_f))(a0); + lua_template<returnType>::push(L,r); + return 1; + } +}; +//0 arg with return +template <class objectType, class returnType> +class methodBindingTemplate<objectType,returnType,void,void,void,void>:public methodBinding { + typedef returnType (objectType::*callback)(); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + returnType r=(obj->*(_f))(); + lua_template<returnType>::push(L,r); + return 1; + } +}; +//4 args without return +template <class objectType, class arg0Type, class arg1Type, class arg2Type, class arg3Type> +class methodBindingTemplate<objectType,void,arg0Type,arg1Type,arg2Type,arg3Type>:public methodBinding { + typedef void (objectType::*callback)(arg0Type,arg1Type,arg2Type,arg3Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + arg1Type a1 = lua_template<arg1Type>::get(L,3); + arg2Type a2 = lua_template<arg2Type>::get(L,4); + arg3Type a3 = lua_template<arg3Type>::get(L,5); + (obj->*(_f))(a0,a1,a2,a3); + return 0; + } +}; +//3 args without return +template <class objectType, class arg0Type, class arg1Type, class arg2Type> +class methodBindingTemplate<objectType,void,arg0Type,arg1Type,arg2Type,void>:public methodBinding { + typedef void (objectType::*callback)(arg0Type,arg1Type,arg2Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + arg1Type a1 = lua_template<arg1Type>::get(L,3); + arg2Type a2 = lua_template<arg2Type>::get(L,4); + (obj->*(_f))(a0,a1,a2); + return 0; + } +}; +//2 args without return +template <class objectType, class arg0Type, class arg1Type> +class methodBindingTemplate<objectType,void,arg0Type,arg1Type,void,void>:public methodBinding { + typedef void (objectType::*callback)(arg0Type,arg1Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + arg1Type a1 = lua_template<arg1Type>::get(L,3); + (obj->*(_f))(a0,a1); + return 0; + } +}; +//1 arg without return +template <class objectType, class arg0Type> +class methodBindingTemplate<objectType,void,arg0Type,void,void,void>:public methodBinding { + typedef void (objectType::*callback)(arg0Type); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,2); + (obj->*(_f))(a0); + return 0; + } +}; +//0 arg without return +template <class objectType> +class methodBindingTemplate<objectType,void,void,void,void,void>:public methodBinding { + typedef void (objectType::*callback)(); + callback _f; + public: + methodBindingTemplate(const std::string luaname,callback f):methodBinding(luaname){ + _f=f; + } + int call (lua_State *L) { + objectType *obj = classBinding<objectType>::check(L,1); + (obj->*(_f))(); + return 0; + } +}; + +//constructor 4 args +template<class objectType, class arg0Type=void, class arg1Type=void, class arg2Type=void, class arg3Type=void> +class constructorBindingTemplate:public constructorBinding { + int call (lua_State *L){ + lua_remove(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,1); + arg1Type a1 = lua_template<arg1Type>::get(L,2); + arg2Type a2 = lua_template<arg2Type>::get(L,3); + arg3Type a3 = lua_template<arg3Type>::get(L,4); + classBinding<objectType>::push(L, new objectType(a0,a1,a2,a3),true); + return 1; + } +}; +//constructor 3 args +template<class objectType, class arg0Type, class arg1Type, class arg2Type> +class constructorBindingTemplate<objectType,arg0Type,arg1Type,arg2Type,void>:public constructorBinding { + int call (lua_State *L){ + lua_remove(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,1); + arg1Type a1 = lua_template<arg1Type>::get(L,2); + arg2Type a2 = lua_template<arg2Type>::get(L,3); + classBinding<objectType>::push(L, new objectType(a0,a1,a2),true); + return 1; + } +}; +//constructor 2 args +template<class objectType, class arg0Type, class arg1Type> +class constructorBindingTemplate<objectType,arg0Type,arg1Type,void,void>:public constructorBinding { + int call (lua_State *L){ + lua_remove(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,1); + arg1Type a1 = lua_template<arg1Type>::get(L,2); + classBinding<objectType>::push(L,new objectType(a0,a1),true); + return 1; + } +}; +//constructor 1 args +template<class objectType, class arg0Type> +class constructorBindingTemplate<objectType,arg0Type,void,void,void>:public constructorBinding { + int call (lua_State *L){ + lua_remove(L,1); + arg0Type a0 = lua_template<arg0Type>::get(L,1); + classBinding<objectType>::push(L, new objectType(a0),true); + return 1; + } +}; +//constructor 0 args +template<class objectType> +class constructorBindingTemplate<objectType,void,void,void,void>:public constructorBinding { + int call (lua_State *L){ + lua_remove(L,1); + classBinding<objectType>::push(L, new objectType(),true); + return 1; + } +}; +#endif diff --git a/Common/luna.h b/Common/luna.h deleted file mode 100755 index 67d066d637cafc0a5b0d2fda4122559f2cbeda79..0000000000000000000000000000000000000000 --- a/Common/luna.h +++ /dev/null @@ -1,131 +0,0 @@ -/* -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 d119667e11b5395ee6e04288b44fcc7b83377bc1..89e0c22ee78eb24687ec7723fc55c3ab0e93b3f0 100644 --- a/Geo/CMakeLists.txt +++ b/Geo/CMakeLists.txt @@ -25,7 +25,6 @@ set(SRC findLinks.cpp SOrientedBoundingBox.cpp GeomMeshMatcher.cpp - LuaBindings_Geo.cpp MVertex.cpp MFace.cpp MElement.cpp diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp index 6949444b6940c2d06fe3b09f67b7909a58e4ed30..453dae3999666308c496b6715564a7d3d0268136 100644 --- a/Geo/GModel.cpp +++ b/Geo/GModel.cpp @@ -27,6 +27,9 @@ #include "Context.h" #include "OS.h" +#include "OpenFile.h" +#include "CreateFile.h" + #if defined(HAVE_MESH) #include "Field.h" #include "Generator.h" @@ -1360,3 +1363,30 @@ GModel *GModel::buildCutGModel(gLevelset *ls) return cutGM; } +void GModel::load(std::string fileName){ + GModel *temp = GModel::current(); + GModel::setCurrent(this); + MergeFile(fileName.c_str()); + GModel::setCurrent(temp); +} + +void GModel::save(std::string fileName){ + GModel *temp = GModel::current(); + GModel::setCurrent(this); + int guess = GuessFileFormatFromFileName(fileName); + CreateOutputFile (fileName, guess); + GModel::setCurrent(temp); +} + +#ifdef HAVE_LUA +#include "Bindings.h" +const char GModel::className[]="GModel"; +const char GModel::parentClassName[]=""; +methodBinding *GModel::methods[]={ + new methodBindingTemplate<GModel,int,int>("mesh",&GModel::mesh), + new methodBindingTemplate<GModel,void,std::string>("load",&GModel::load), + new methodBindingTemplate<GModel,void,std::string>("save",&GModel::save), + 0 +}; +constructorBinding *GModel::constructorMethod=new constructorBindingTemplate<GModel>(); +#endif diff --git a/Geo/GModel.h b/Geo/GModel.h index b43b3d33d6e492da559f76e6d7566fdf6e34782b..80c4a8597fc651a9fbdae8f730261bdee28a1c04 100644 --- a/Geo/GModel.h +++ b/Geo/GModel.h @@ -19,6 +19,11 @@ #include "SBoundingBox3d.h" #include "discreteFace.h" +#ifdef HAVE_LUA + class methodBinding; + class constructorBinding; +#endif + class Octree; class FM_Internals; class GEO_Internals; @@ -433,6 +438,17 @@ class GModel int readDIFF(const std::string &name); int writeDIFF(const std::string &name, bool binary=false, bool saveAll=false, double scalingFactor=1.0); + + + void save(std::string fileName); + void load(std::string fileName); + +#ifdef HAVE_LUA + static const char className[]; + static const char parentClassName[]; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; +#endif }; #endif diff --git a/Geo/LuaBindings_Geo.cpp b/Geo/LuaBindings_Geo.cpp deleted file mode 100644 index f237b3d0cdef5017b3e4450f72c5d4ffa1eaacc7..0000000000000000000000000000000000000000 --- a/Geo/LuaBindings_Geo.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include "LuaBindings_Geo.h" -#include "OpenFile.h" -#include "CreateFile.h" - -#if defined(HAVE_LUA) - -// 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; -} - -#endif diff --git a/Geo/LuaBindings_Geo.h b/Geo/LuaBindings_Geo.h deleted file mode 100644 index be6dda96f04fdd5b0f79320222cf9859f8a75c70..0000000000000000000000000000000000000000 --- a/Geo/LuaBindings_Geo.h +++ /dev/null @@ -1,39 +0,0 @@ -#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 6bb80183857e1b094861f555b56deaae6d0f9f7a..2b169cb54851a34c6aec2c82a07a7f7eee0fb5dd 100644 --- a/Numeric/fullMatrix.cpp +++ b/Numeric/fullMatrix.cpp @@ -290,56 +290,26 @@ bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S) } #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} +#include "Bindings.h" +template<> +const char fullMatrix<double>::className[]="fullMatrix"; +template<> +const char fullMatrix<float>::className[]="fullMatrixFloat"; +template<> +const char fullMatrix<double>::parentClassName[]=""; +template<> +const char fullMatrix<float>::parentClassName[]=""; +template<> +methodBinding *fullMatrix<double>::methods[]={ + new methodBindingTemplate<const fullMatrix<double>,int>("size1",&fullMatrix<double>::size1), + new methodBindingTemplate<const fullMatrix<double>,int>("size2",&fullMatrix<double>::size2), + new methodBindingTemplate<const fullMatrix<double>,double,int,int>("get",&fullMatrix<double>::get), + new methodBindingTemplate<fullMatrix<double>,void,int,int,double>("set",&fullMatrix<double>::set), + new methodBindingTemplate<fullMatrix<double>,void,const fullMatrix<double>*,const fullMatrix<double> *>("gemm",&fullMatrix<double>::gemm), + 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); -} +template<> +constructorBinding *fullMatrix<double>::constructorMethod = new constructorBindingTemplate<fullMatrix<double>,int,int>(); #endif diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h index 6d3f28b87605db3ca4174b0c189414544a198994..4a3a2c451b4a05317061487d6e15fefe07889a88 100644 --- a/Numeric/fullMatrix.h +++ b/Numeric/fullMatrix.h @@ -11,17 +11,6 @@ #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> @@ -104,6 +93,10 @@ class fullVector printf("\n"); } }; +#if defined(HAVE_LUA) +class methodBinding; +class constructorBinding; +#endif template <class scalar> class fullMatrix @@ -113,51 +106,17 @@ 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; + inline scalar get (int r, int c)const { + return (*this)(r,c); } - 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; + inline void set(int r, int c, scalar v){ + (*this)(r,c)=v; } - - 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); +#if defined(HAVE_LUA) + static const char className[]; + static const char parentClassName[]; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; #endif // HAVE LUA fullMatrix(scalar *original, int r, int c){ _r = r; @@ -260,6 +219,9 @@ class fullMatrix } #endif ; + inline void gemm (const fullMatrix<scalar> *a, const fullMatrix<scalar> *b){ + gemm(*a,*b); + } void gemm(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b, scalar alpha=1., scalar beta=1.) #if !defined(HAVE_BLAS) diff --git a/Solver/TESTCASES/Advection1D.lua b/Solver/TESTCASES/Advection1D.lua index f472e9d558b041590acea174740168442c7b4a75..86a3f399324f325e1e061bf8298cf2670431f056 100644 --- a/Solver/TESTCASES/Advection1D.lua +++ b/Solver/TESTCASES/Advection1D.lua @@ -2,7 +2,7 @@ model = GModel () model:load ('edge.geo') --model:mesh (1) --model:save ('edge_new.msh') -model:load ('edge_new.msh') +model:load ('edge.msh') dg = dgSystemOfEquations (model) dg:setOrder(1) @@ -17,20 +17,21 @@ v:set(2,0,0) nu=fullMatrix(1,1); nu:set(0,0,0) -dg:setConservationLaw('AdvectionDiffusion',createFunction.constant(v),createFunction.constant(nu)) +law = ConservationLawAdvection(FunctionConstant(v):getName(),FunctionConstant(nu):getName()) + +dg:setConservationLaw(law) -- boundary condition outside=fullMatrix(1,1) outside:set(0,0,0.) -dg:addBoundaryCondition('Left','OutsideValues',createFunction.constant(outside)) -dg:addBoundaryCondition('Right','OutsideValues',createFunction.constant(outside)) +bndcondition=law:newOutsideValueBoundary(FunctionConstant(outside):getName()) +law:addBoundaryCondition('Left',bndcondition) +law:addBoundaryCondition('Right',bndcondition) dg:setup() -- initial condition -function initial_condition( _x , _f ) - xyz = fullMatrix(_x) - f = fullMatrix(_f) +function initial_condition( xyz , f ) for i=0,xyz:size1()-1 do x = xyz:get(i,0) y = xyz:get(i,1) @@ -38,14 +39,14 @@ function initial_condition( _x , _f ) f:set (i, 0, math.exp(-100*((x+0.8)^2))) end end -dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ')) +dg:L2Projection(FunctionLua(1,'initial_condition',{'XYZ'}):getName()) dg:exportSolution('output/Adv1D_00000') -- main loop n = 5 -for i=1,50*n do - norm = dg:RK44(0.0325) +for i=1,100*n do + norm = dg:RK44(0.03) if (i % n == 0) then print('iter',i,norm) dg:exportSolution(string.format("output/Adv1D-%05d", i)) diff --git a/Solver/TESTCASES/AdvectionDiffusion.lua b/Solver/TESTCASES/AdvectionDiffusion.lua index 90a85807b86d6b0a7c41c114e30a045a5922c036..36e1f84f688cb5b81131da72af52f8e9d4566dc6 100644 --- a/Solver/TESTCASES/AdvectionDiffusion.lua +++ b/Solver/TESTCASES/AdvectionDiffusion.lua @@ -1,33 +1,41 @@ model = GModel () model:load ('square.geo') model:load ('square.msh') + +model2 = GModel() +model2:load('') + +vmodel = {GModel(),GModel()} +vmodel[1]:load('') +vmodel[2]:load('') + dg = dgSystemOfEquations (model) dg:setOrder(5) - -- conservation law + -- advection speed v=fullMatrix(3,1); v:set(0,0,0.15) v:set(1,0,0.05) v:set(2,0,0) + -- diffusivity nu=fullMatrix(1,1); nu:set(0,0,0.001) -dg:setConservationLaw('AdvectionDiffusion',createFunction.constant(v),createFunction.constant(nu)) +law = ConservationLawAdvection(FunctionConstant(v):getName(),FunctionConstant(nu):getName()) +dg:setConservationLaw(law) -- boundary condition outside=fullMatrix(1,1) outside:set(0,0,0.) -dg:addBoundaryCondition('Border','OutsideValues',createFunction.constant(outside)) +law:addBoundaryCondition('Border',law:newOutsideValueBoundary(FunctionConstant(outside):getName())) dg:setup() -- initial condition -function initial_condition( _x , _f ) - xyz = fullMatrix(_x) - f = fullMatrix(_f) +function initial_condition( xyz , f ) for i=0,xyz:size1()-1 do x = xyz:get(i,0) y = xyz:get(i,1) @@ -35,7 +43,7 @@ function initial_condition( _x , _f ) f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2))) end end -dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ')) +dg:L2Projection(FunctionLua(1,'initial_condition',{'XYZ'}):getName()) dg:exportSolution('output/Advection_00000') diff --git a/Solver/TESTCASES/BackwardFacingStep.lua b/Solver/TESTCASES/BackwardFacingStep.lua index 00bc95c90faa79457141cc7fdd07894e68e56332..7002a931fe069633ad95d1401e8b32259b459570 100644 --- a/Solver/TESTCASES/BackwardFacingStep.lua +++ b/Solver/TESTCASES/BackwardFacingStep.lua @@ -7,9 +7,7 @@ SOUND = V/MACH --[[ Function for initial conditions --]] -function free_stream( x, f ) - FCT = fullMatrix(f) - XYZ = fullMatrix(x) +function free_stream( XYZ, FCT ) for i=0,XYZ:size1()-1 do FCT:set(i,0,RHO) FCT:set(i,1,RHO*V) @@ -30,12 +28,13 @@ myModel:load ('step.msh') print'*** Create a dg solver ***' DG = dgSystemOfEquations (myModel) DG:setOrder(order) -DG:setConservationLaw('PerfectGas2d') -DG:addBoundaryCondition('Walls','Wall') +law=ConservationLawPerfectGas2d() +DG:setConservationLaw(law) +law:addBoundaryCondition('Walls',law:newWallBoundary()) -FS = createFunction.lua(4, 'free_stream', 'XYZ') +FS = FunctionLua(4, 'free_stream', {'XYZ'}):getName() -DG:addBoundaryCondition('LeftRight','FreeStream',FS) +law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS)) DG:setup() @@ -56,7 +55,7 @@ print('DT=',dt) for i=1,1000 do norm = DG:RK44(dt) print('*** ITER ***',i,norm) - if (i % 20 == 0) then + if (i % 10 == 0) then DG:exportSolution(string.format("solution-%03d", i)) end end diff --git a/Solver/TESTCASES/Diffusion.lua b/Solver/TESTCASES/Diffusion.lua index 99f20e84de282191a574fd989ba94d3305c5eec6..4b969dc01e2a50ea87dc33a0363c98c17660fd37 100644 --- a/Solver/TESTCASES/Diffusion.lua +++ b/Solver/TESTCASES/Diffusion.lua @@ -9,19 +9,18 @@ dg:setOrder(5) -- advection speed nu=fullMatrix(1,1); nu:set(0,0,0.01) -dg:setConservationLaw('AdvectionDiffusion','',createFunction.constant(nu)) +law = ConservationLawAdvection('',FunctionConstant(nu):getName()) +dg:setConservationLaw(law) -- boundary condition outside=fullMatrix(1,1) outside:set(0,0,0.) -dg:addBoundaryCondition('Border','0Flux') +law:addBoundaryCondition('Border',law:new0FluxBoundary()) dg:setup() -- initial condition -function initial_condition( _x , _f ) - xyz = fullMatrix(_x) - f = fullMatrix(_f) +function initial_condition( xyz , f ) for i=0,xyz:size1()-1 do x = xyz:get(i,0) y = xyz:get(i,1) @@ -29,7 +28,7 @@ function initial_condition( _x , _f ) f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2))) end end -dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ')) +dg:L2Projection(FunctionLua(1,'initial_condition',{'XYZ'}):getName()) dg:exportSolution('output/Diffusion_00000') diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua index 1d9e3b2ec90d376b6e039fae0b545fa3b4e89199..b5017e97c941d5cf8703f8efab0df1fdf243cea9 100644 --- a/Solver/TESTCASES/Stommel.lua +++ b/Solver/TESTCASES/Stommel.lua @@ -1,23 +1,18 @@ - model = GModel() model:load ('stommel_square.msh') -print('load..') - -order = 3 dg = dgSystemOfEquations (model) +order=1 dg:setOrder (order) - - -dg:setConservationLaw('ShallowWater2d') -dg:addBoundaryCondition('Wall','Wall') +claw = ShallowWater2d() +claw:addBoundaryCondition('Wall',claw:newWallBoundary()) +dg:setConservationLaw(claw) dg:setup() +dg:exportSolution('output/init') -dg:exportSolution('output/solution_00000') - -for i=1,1000 do +for i=1,10000 do norm = dg:RK44(80*(3/(2.*order+1))) - if ( i%10 ==0 ) then + if ( i%100 ==0 ) then print ('iter ', i, norm) end if ( i%100 ==0 ) then diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua index 815eb769e9ebb4af100cb057a8d0772b5ca8f7d1..0ef157ba57622cc354f45d21315b79bb459d95cd 100644 --- a/Solver/TESTCASES/WavePulse.lua +++ b/Solver/TESTCASES/WavePulse.lua @@ -1,9 +1,7 @@ --[[ Function for initial conditions --]] -function initial_condition( _x , _f ) - xyz = fullMatrix(_x) - f = fullMatrix(_f) +function initial_condition( xyz , f ) for i=0,xyz:size1()-1 do X = xyz:get(i,0) - .5 Y = xyz:get(i,1) - .5 @@ -15,24 +13,22 @@ function initial_condition( _x , _f ) end end ---[[ - Example of a lua program driving the DG code ---]] print'*** Loading the mesh and the model ***' myModel = GModel () -myModel:load('box.geo') -myModel:load('box.msh') +--myModel:load('box.geo') +--myModel:load('box.msh') +myModel:load('square.msh') print'*** Create a dg solver ***' DG = dgSystemOfEquations (myModel) DG:setOrder(1) -DG:setConservationLaw('WaveEquation') - -DG:addBoundaryCondition('Border','Wall') +law=ConservationLawWaveEquation(2) +DG:setConservationLaw(law) +law:addBoundaryCondition('Border',law:newWallBoundary()) DG:setup() -initialCondition = createFunction.lua(3,'initial_condition','XYZ') +initialCondition = FunctionLua(3,'initial_condition',{'XYZ'}):getName() print'*** setting the initial solution ***' @@ -45,12 +41,12 @@ DG:exportSolution('output/solution_000') print'*** solve ***' dt = 0.00125; -N = 1; +N = 1000; for i=1,N do norm = DG:RK44(dt) print('*** ITER ***',i,norm) if (i % 10 == 0) then - DG:exportSolution(string.format("solution-%03d", i)) + DG:exportSolution(string.format("output/solution-%03d", i)) end end diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp index e4c7d1cb6d77a306aee59fe0b6c58b0c90e083ab..f1e8cd05def2f64a22f0f373eb95ac6ce3c8e0eb 100644 --- a/Solver/dgConservationLaw.cpp +++ b/Solver/dgConservationLaw.cpp @@ -26,7 +26,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition { if(riemannSolver){ for(int i=0;i<_value.size1(); i++) for(int j=0;j<_value.size2(); j++) - _value(i,j) = (*riemannSolver)(i,j*2); + _value(i,j) = (*riemannSolver)(i,j); } } }; @@ -55,10 +55,23 @@ class dgBoundaryCondition0Flux : public dgBoundaryCondition { } }; -dgBoundaryCondition *dgBoundaryCondition::newOutsideValueCondition(dgConservationLaw &claw,const std::string outsideValueFunctionName) { - return new dgBoundaryConditionOutsideValue(claw,outsideValueFunctionName); +dgBoundaryCondition *dgConservationLaw::newOutsideValueBoundary(const std::string outsideValueFunctionName) { + return new dgBoundaryConditionOutsideValue(*this,outsideValueFunctionName); } -dgBoundaryCondition *dgBoundaryCondition::new0FluxCondition(dgConservationLaw &claw) { - return new dgBoundaryCondition0Flux(claw); +dgBoundaryCondition *dgConservationLaw::new0FluxBoundary() { + return new dgBoundaryCondition0Flux(*this); } +#include "Bindings.h" +const char dgConservationLaw::className[]="ConservationLaw"; +const char dgConservationLaw::parentClassName[]=""; +methodBinding * dgConservationLaw::methods[]={ + new methodBindingTemplate<dgConservationLaw,void,std::string,dgBoundaryCondition*>("addBoundaryCondition",&dgConservationLaw::addBoundaryCondition), + new methodBindingTemplate<dgConservationLaw,dgBoundaryCondition*>("new0FluxBoundary",&dgConservationLaw::new0FluxBoundary), + new methodBindingTemplate<dgConservationLaw,dgBoundaryCondition*,std::string>("newOutsideValueBoundary",&dgConservationLaw::newOutsideValueBoundary), +0}; +constructorBinding * dgConservationLaw::constructorMethod=NULL; +const char dgBoundaryCondition::className[]="BoundaryCondition"; +const char dgBoundaryCondition::parentClassName[]=""; +methodBinding * dgBoundaryCondition::methods[]={0}; +constructorBinding * dgBoundaryCondition::constructorMethod=NULL; diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h index f8c9f840ecfc6842f2fa95211b20dab1dfe46696..f500f9587847e3ecdc98063a755d2f46812cf581 100644 --- a/Solver/dgConservationLaw.h +++ b/Solver/dgConservationLaw.h @@ -9,6 +9,10 @@ #include "fullMatrix.h" class dataCacheDouble; class dataCacheMap; +#ifdef HAVE_LUA +class constructorBinding; +class methodBinding; +#endif class dgConservationLaw; @@ -16,9 +20,12 @@ class dgBoundaryCondition { public: virtual ~dgBoundaryCondition () {} virtual dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const = 0; - //a generic boundary condition using the Riemann solver of the conservation Law - static dgBoundaryCondition *newOutsideValueCondition(dgConservationLaw &claw,const std::string outsideValueFunctionName); - static dgBoundaryCondition *new0FluxCondition(dgConservationLaw &claw); +#if defined(HAVE_LUA) + static const char className[]; + static const char parentClassName[]; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; +#endif }; class dgConservationLaw { @@ -45,21 +52,26 @@ class dgConservationLaw { return it->second; } - inline void addBoundaryCondition(const std::string tag, dgBoundaryCondition * condition) { + void addBoundaryCondition(std::string tag, dgBoundaryCondition * condition) { if(_boundaryConditions.find(tag)!=_boundaryConditions.end()) throw; _boundaryConditions[tag]=condition; } + //a generic boundary condition using the Riemann solver of the conservation Law + dgBoundaryCondition *newOutsideValueBoundary(std::string outsideValueFunctionName); + dgBoundaryCondition *new0FluxBoundary(); + + #ifdef HAVE_LUA + static const char className[]; + static const char parentClassName[]; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; + #endif }; -dgConservationLaw *dgNewConservationLawAdvection(const std::string vname,const std::string nuname); -dgConservationLaw *dgNewConservationLawShallowWater2d(); -dgConservationLaw *dgNewConservationLawWaveEquation(int); dgConservationLaw *dgNewPerfectGasLaw2d(); -dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall(int); -dgBoundaryCondition *dgNewBoundaryConditionShallowWater2dWall(); dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dWall(); dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dFreeStream(std::string&); diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp index ac44af08c3908c9e6ea10516db4e08b081ef2130..a97805fca2510274a6412bce64eff5d821bcc30b 100644 --- a/Solver/dgConservationLawAdvection.cpp +++ b/Solver/dgConservationLawAdvection.cpp @@ -4,98 +4,98 @@ #include "SPoint3.h" #include "MElement.h" #include "function.h" +#include "dgConservationLawAdvection.h" -class dgConservationLawAdvection : public dgConservationLaw { - std::string _vFunctionName,_nuFunctionName; - class advection : public dataCacheDouble { - dataCacheDouble &sol, &v; - public: - advection(std::string vFunctionName, dataCacheMap &cacheMap): - dataCacheDouble(cacheMap.getNbEvaluationPoints(),3), - sol(cacheMap.get("Solution",this)), - v(cacheMap.get(vFunctionName,this)) - {}; - void _eval () { - for(int i=0; i< sol().size1(); i++) { - _value(i,0) = sol(i,0)*v(i,0); - _value(i,1) = sol(i,0)*v(i,1); - _value(i,2) = sol(i,0)*v(i,2); - } - } - }; - class riemann : public dataCacheDouble { - dataCacheDouble &normals, &solLeft, &solRight,&v; - public: - riemann(std::string vFunctionName, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight): - dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),2), - normals(cacheMapLeft.get("Normals", this)), - solLeft(cacheMapLeft.get("Solution", this)), - solRight(cacheMapRight.get("Solution", this)), - v(cacheMapLeft.get(vFunctionName,this)) - {}; - void _eval () { - for(int i=0; i< _value.size1(); i++) { - double un=v(i,0)*normals(0,i)+v(i,1)*normals(1,i)+v(i,2)*normals(2,i); - if(un>0){ - _value(i,0) = -solLeft(i,0)*un; - _value(i,1) = solLeft(i,0)*un; - }else{ - _value(i,0) = -solRight(i,0)*un; - _value(i,1) = solRight(i,0)*un; - } - } +class dgConservationLawAdvection::advection : public dataCacheDouble { + dataCacheDouble &sol, &v; + public: + advection(std::string vFunctionName, dataCacheMap &cacheMap): + dataCacheDouble(cacheMap.getNbEvaluationPoints(),3), + sol(cacheMap.get("Solution",this)), + v(cacheMap.get(vFunctionName,this)) + {}; + void _eval () { + for(int i=0; i< sol().size1(); i++) { + _value(i,0) = sol(i,0)*v(i,0); + _value(i,1) = sol(i,0)*v(i,1); + _value(i,2) = sol(i,0)*v(i,2); } - }; - class diffusion : public dataCacheDouble { - dataCacheDouble &solgrad, ν - public: - diffusion(std::string nuFunctionName, dataCacheMap &cacheMap): - solgrad(cacheMap.get("SolutionGradient",this)), - nu(cacheMap.get(nuFunctionName,this)) - {}; - void _eval () { - if(_value.size1() != solgrad().size1()/3) - _value=fullMatrix<double>(solgrad().size1()/3,3); - for(int i=0; i< solgrad().size1()/3; i++) { - _value(i,0) = -solgrad(i*3,0)*nu(i,0); - _value(i,1) = -solgrad(i*3+1,0)*nu(i,0); - _value(i,2) = -solgrad(i*3+1,0)*nu(i,0); + } +}; +class dgConservationLawAdvection::riemann : public dataCacheDouble { + dataCacheDouble &normals, &solLeft, &solRight,&v; + public: + riemann(std::string vFunctionName, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight): + dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),2), + normals(cacheMapLeft.get("Normals", this)), + solLeft(cacheMapLeft.get("Solution", this)), + solRight(cacheMapRight.get("Solution", this)), + v(cacheMapLeft.get(vFunctionName,this)) + {}; + void _eval () { + for(int i=0; i< _value.size1(); i++) { + double un=v(i,0)*normals(0,i)+v(i,1)*normals(1,i)+v(i,2)*normals(2,i); + if(un>0){ + _value(i,0) = -solLeft(i,0)*un; + _value(i,1) = solLeft(i,0)*un; + }else{ + _value(i,0) = -solRight(i,0)*un; + _value(i,1) = solRight(i,0)*un; } } - }; - public: - dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const { - if( !_vFunctionName.empty()) - return new advection(_vFunctionName,cacheMap); - else - return NULL; - } - dataCacheDouble *newMaximumDiffusivity( dataCacheMap &cacheMap) const { - if( !_nuFunctionName.empty()) - return &cacheMap.get(_nuFunctionName); - else - return NULL; - } - dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { - if( !_vFunctionName.empty()) - return new riemann(_vFunctionName,cacheMapLeft, cacheMapRight); - else - return NULL; - } - dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const { - if( !_nuFunctionName.empty()) - return new diffusion(_nuFunctionName,cacheMap); - else - return NULL; } - dgConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName) - { - _vFunctionName = vFunctionName; - _nuFunctionName = nuFunctionName; - _nbf = 1; +}; +class dgConservationLawAdvection::diffusion : public dataCacheDouble { + dataCacheDouble &solgrad, ν + public: + diffusion(std::string nuFunctionName, dataCacheMap &cacheMap): + solgrad(cacheMap.get("SolutionGradient",this)), + nu(cacheMap.get(nuFunctionName,this)) + {}; + void _eval () { + if(_value.size1() != solgrad().size1()/3) + _value=fullMatrix<double>(solgrad().size1()/3,3); + for(int i=0; i< solgrad().size1()/3; i++) { + _value(i,0) = -solgrad(i*3,0)*nu(i,0); + _value(i,1) = -solgrad(i*3+1,0)*nu(i,0); + _value(i,2) = -solgrad(i*3+2,0)*nu(i,0); + } } }; - -dgConservationLaw *dgNewConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName) { - return new dgConservationLawAdvection(vFunctionName,nuFunctionName); +dataCacheDouble *dgConservationLawAdvection::newConvectiveFlux( dataCacheMap &cacheMap) const { + if( !_vFunctionName.empty()) + return new advection(_vFunctionName,cacheMap); + else + return NULL; +} +dataCacheDouble *dgConservationLawAdvection::newMaximumDiffusivity( dataCacheMap &cacheMap) const { + if( !_nuFunctionName.empty()) + return &cacheMap.get(_nuFunctionName); + else + return NULL; } +dataCacheDouble *dgConservationLawAdvection::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { + if( !_vFunctionName.empty()) + return new riemann(_vFunctionName,cacheMapLeft, cacheMapRight); + else + return NULL; +} +dataCacheDouble *dgConservationLawAdvection::newDiffusiveFlux( dataCacheMap &cacheMap) const { + if( !_nuFunctionName.empty()) + return new diffusion(_nuFunctionName,cacheMap); + else + return NULL; +} +dgConservationLawAdvection::dgConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName) +{ + _vFunctionName = vFunctionName; + _nuFunctionName = nuFunctionName; + _nbf = 1; +} + +#include "Bindings.h" +const char * dgConservationLawAdvection::className = "ConservationLawAdvection"; +const char * dgConservationLawAdvection::parentClassName = "ConservationLaw"; +constructorBinding *dgConservationLawAdvection::constructorMethod = new constructorBindingTemplate<dgConservationLawAdvection,std::string,std::string>(); +methodBinding *dgConservationLawAdvection::methods []={0}; + diff --git a/Solver/dgConservationLawAdvection.h b/Solver/dgConservationLawAdvection.h new file mode 100644 index 0000000000000000000000000000000000000000..69c11578a2bb322280158e1ed91a1326a3ea47d5 --- /dev/null +++ b/Solver/dgConservationLawAdvection.h @@ -0,0 +1,22 @@ +#ifndef _DG_CONSERVATION_LAW_ADVECTION_H +#define _DG_CONSERVATION_LAW_ADVECTION_H +#include "dgConservationLaw.h" +class constructorBinding; +class methodBinding; +class dgConservationLawAdvection : public dgConservationLaw { + std::string _vFunctionName,_nuFunctionName; + class advection; + class riemann; + class diffusion; + public: + dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const; + dataCacheDouble *newMaximumDiffusivity( dataCacheMap &cacheMap) const; + dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const; + dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const; + dgConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName); + static const char * className; + static const char * parentClassName; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; +}; +#endif diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp index e90235a0d4da4bdd6e2b250e58f6ad6571756d1c..4776a51719042691995e5fe9b4e433a5e398055d 100644 --- a/Solver/dgConservationLawPerfectGas.cpp +++ b/Solver/dgConservationLawPerfectGas.cpp @@ -1,4 +1,4 @@ -#include "dgConservationLaw.h" +#include "dgConservationLawPerfectGas.h" #include "function.h" const double GAMMA = 1.4; @@ -127,107 +127,97 @@ static inline void _ROE2D (const double &_GAMMA, FLUX[k] = -lflux; } } +// perfect gas law, GAMMA is the only parameter -class dgPerfectGasLaw2d : public dgConservationLaw { +class dgPerfectGasLaw2d::advection : public dataCacheDouble { + dataCacheDouble / + public: + advection(dataCacheMap &cacheMap): + sol(cacheMap.get("Solution",this)) + {}; + void _eval () { + const int nQP = sol().size1(); + if(_value.size1() != nQP) + _value=fullMatrix<double>(nQP,8); + const double GM1 = GAMMA - 1.0; + for (size_t k = 0 ; k < nQP; k++ ){ + // printf("%d %g %g %g %g\n",k,sol(k,0),sol(k,1),sol(k,2),sol(k,3)); + const double invrho = 1./sol(k,0); - // perfect gas law, GAMMA is the only parameter + const double q12 = sol(k,1)*sol(k,2)*invrho; + const double q11 = sol(k,1)*sol(k,1)*invrho; + const double q22 = sol(k,2)*sol(k,2)*invrho; - class advection : public dataCacheDouble { - dataCacheDouble / - public: - advection(dataCacheMap &cacheMap): - sol(cacheMap.get("Solution",this)) - {}; - void _eval () { - const int nQP = sol().size1(); - if(_value.size1() != nQP) - _value=fullMatrix<double>(nQP,8); - const double GM1 = GAMMA - 1.0; - for (size_t k = 0 ; k < nQP; k++ ){ - // printf("%d %g %g %g %g\n",k,sol(k,0),sol(k,1),sol(k,2),sol(k,3)); - const double invrho = 1./sol(k,0); - - const double q12 = sol(k,1)*sol(k,2)*invrho; - const double q11 = sol(k,1)*sol(k,1)*invrho; - const double q22 = sol(k,2)*sol(k,2)*invrho; - - const double p = GM1*sol(k,3) - 0.5*GM1*(q11+q22); - const double qq = invrho*(sol(k,3)+p); - - _value(k,0) = sol(k,1); - _value(k,1) = q11+p; - _value(k,2) = q12; - _value(k,3) = sol(k,1)*qq; + const double p = GM1*sol(k,3) - 0.5*GM1*(q11+q22); + const double qq = invrho*(sol(k,3)+p); - _value(k,0+4) = sol(k,2); - _value(k,1+4) = q12; - _value(k,2+4) = q22+p; - _value(k,3+4) = sol(k,2)*qq; + _value(k,0) = sol(k,1); + _value(k,1) = q11+p; + _value(k,2) = q12; + _value(k,3) = sol(k,1)*qq; - /* _value(k,8) = 0; - _value(k,9) = 0; - _value(k,10) = 0; - _value(k,11) = 0;*/ - } - } - }; - - class riemann : public dataCacheDouble { - dataCacheDouble &normals, &solL, &solR; - public: - riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight): - normals(cacheMapLeft.get("Normals", this)), - solL(cacheMapLeft.get("Solution", this)), - solR(cacheMapRight.get("Solution", this)) - {}; - void _eval () { - int nQP = solL().size1(); - if(_value.size1() != nQP) - _value = fullMatrix<double>(nQP,8); - - for(int i=0; i< nQP; i++) { - const double solLeft [4] = {solL(i,0),solL(i,1),solL(i,2),solL(i,3)}; - const double solRight[4] = {solR(i,0),solR(i,1),solR(i,2),solR(i,3)}; - double FLUX[4] ; - const double nx = normals(0,i); - const double ny = normals(1,i); - _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX); - - /* - printf("RSOLL %g %g %g %g\n",solLeft[0],solLeft[1],solLeft[2], solLeft[3]); - printf("RSOLR %g %g %g %g\n",solRight[0],solRight[1],solRight[2], solRight[3]); - printf("RN %g %g\n",nx,ny); - printf("RFLUX %g %g %g %g\n",FLUX[0],FLUX[1],FLUX[2],FLUX[3]); - */ - _value(i,0) = FLUX[0]; - _value(i,1) = FLUX[1]; - _value(i,2) = FLUX[2]; - _value(i,3) = FLUX[3]; - _value(i,4) = -_value(i,0); - _value(i,5) = -_value(i,1); - _value(i,6) = -_value(i,2); - _value(i,7) = -_value(i,3); - } + _value(k,0+4) = sol(k,2); + _value(k,1+4) = q12; + _value(k,2+4) = q22+p; + _value(k,3+4) = sol(k,2)*qq; + + /* _value(k,8) = 0; + _value(k,9) = 0; + _value(k,10) = 0; + _value(k,11) = 0;*/ } - }; - public: - dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const { - return new advection(cacheMap); - } - dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { - return new riemann(cacheMapLeft, cacheMapRight); - } - dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const { - return 0; } - dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const { - return 0; - } - dgPerfectGasLaw2d() - { - _nbf = 4; // \rho \rho u \rho v \rho e +}; + +class dgPerfectGasLaw2d::riemann : public dataCacheDouble { + dataCacheDouble &normals, &solL, &solR; + public: + riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight): + normals(cacheMapLeft.get("Normals", this)), + solL(cacheMapLeft.get("Solution", this)), + solR(cacheMapRight.get("Solution", this)) + {}; + void _eval () { + int nQP = solL().size1(); + if(_value.size1() != nQP) + _value = fullMatrix<double>(nQP,8); + + for(int i=0; i< nQP; i++) { + const double solLeft [4] = {solL(i,0),solL(i,1),solL(i,2),solL(i,3)}; + const double solRight[4] = {solR(i,0),solR(i,1),solR(i,2),solR(i,3)}; + double FLUX[4] ; + const double nx = normals(0,i); + const double ny = normals(1,i); + _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX); + + _value(i,0) = FLUX[0]; + _value(i,1) = FLUX[1]; + _value(i,2) = FLUX[2]; + _value(i,3) = FLUX[3]; + _value(i,4) = -_value(i,0); + _value(i,5) = -_value(i,1); + _value(i,6) = -_value(i,2); + _value(i,7) = -_value(i,3); + } } }; +dataCacheDouble *dgPerfectGasLaw2d::newConvectiveFlux( dataCacheMap &cacheMap) const { + return new advection(cacheMap); +} +dataCacheDouble *dgPerfectGasLaw2d::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { + return new riemann(cacheMapLeft, cacheMapRight); +} +dataCacheDouble *dgPerfectGasLaw2d::newDiffusiveFlux( dataCacheMap &cacheMap) const { + return 0; +} +dataCacheDouble *dgPerfectGasLaw2d::newSourceTerm (dataCacheMap &cacheMap) const { + return 0; +} +dgPerfectGasLaw2d::dgPerfectGasLaw2d() +{ + _nbf = 4; // \rho \rho u \rho v \rho e +} + class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition { class term : public dataCacheDouble { @@ -268,7 +258,12 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition { return new term(cacheMapLeft); } }; +dgBoundaryCondition *dgPerfectGasLaw2d::newWallBoundary() const { + return new dgBoundaryConditionPerfectGasLaw2dWall(); +} +#if 0 // JF : I commented this out as I think this equivalent to the generic OutsideValue condition +can you confirm and remove it ? class dgBoundaryConditionPerfectGasLaw2dFreeStream : public dgBoundaryCondition { class term : public dataCacheDouble { dataCacheDouble &sol,&normals,&freeStream; @@ -313,16 +308,13 @@ class dgBoundaryConditionPerfectGasLaw2dFreeStream : public dgBoundaryCondition return new term(cacheMapLeft,_freeStreamName); } }; +#endif +#include "Bindings.h" +const char *dgPerfectGasLaw2d::className = "ConservationLawPerfectGas2d"; +const char *dgPerfectGasLaw2d::parentClassName = "ConservationLaw"; +methodBinding *dgPerfectGasLaw2d::methods[] ={ + new methodBindingTemplate<const dgPerfectGasLaw2d,dgBoundaryCondition*>("newWallBoundary",&dgPerfectGasLaw2d::newWallBoundary), +0}; +constructorBinding *dgPerfectGasLaw2d::constructorMethod=new constructorBindingTemplate<dgPerfectGasLaw2d>(); -dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dWall() { - return new dgBoundaryConditionPerfectGasLaw2dWall(); -} - -dgBoundaryCondition *dgNewBoundaryConditionPerfectGasLaw2dFreeStream(std::string &freeStreamName) { - return new dgBoundaryConditionPerfectGasLaw2dFreeStream(freeStreamName); -} - -dgConservationLaw *dgNewPerfectGasLaw2d() { - return new dgPerfectGasLaw2d(); -} diff --git a/Solver/dgConservationLawPerfectGas.h b/Solver/dgConservationLawPerfectGas.h new file mode 100644 index 0000000000000000000000000000000000000000..434e9e602a51a2adf66ea016798b1d4ff2651539 --- /dev/null +++ b/Solver/dgConservationLawPerfectGas.h @@ -0,0 +1,20 @@ +#ifndef _DG_CONSERVATION_LAW_PERFECT_GAS_H_ +#define _DG_CONSERVATION_LAW_PERFECT_GAS_H_ +#include "dgConservationLaw.h" +class dgPerfectGasLaw2d : public dgConservationLaw { + // perfect gas law, GAMMA is the only parameter + class advection; + class riemann; + public: + dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const; + dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const; + dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const; + dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const; + dgPerfectGasLaw2d(); + static const char * className; + static const char * parentClassName; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; + dgBoundaryCondition *newWallBoundary()const ; +}; +#endif diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp index dbf70dded98eb8b57bdbe9e93d221e1ee5ad5a6e..76ddc2f910b0c58d65417cbb44a546e926e63f83 100644 --- a/Solver/dgConservationLawShallowWater2d.cpp +++ b/Solver/dgConservationLawShallowWater2d.cpp @@ -1,118 +1,100 @@ -#include "dgConservationLaw.h" +#include "dgConservationLawShallowWater2d.h" #include "function.h" static double g = 9.81; static double h = 1000; -class dgConservationLawShallowWater2d : public dgConservationLaw { - class advection : public dataCacheDouble { - dataCacheDouble / - public: - advection(dataCacheMap &cacheMap): - sol(cacheMap.get("Solution",this)) - {}; - void _eval () { - int nQP = sol().size1(); - if(_value.size1() != nQP) - _value=fullMatrix<double>(nQP,9); - for(int i=0; i< nQP; i++) { - double eta = sol(i,0); - double u = sol(i,1); - double v = sol(i,2); - // flux_x - _value(i,0) = h*u; - _value(i,1) = g*eta; - _value(i,2) = 0; - // flux_y - _value(i,3) = h*v; - _value(i,4) = 0; - _value(i,5) = g*eta; - // flux_z - _value(i,6) = 0; - _value(i,7) = 0; - _value(i,8) = 0; - } - } - }; - class source : public dataCacheDouble { - dataCacheDouble &xyz, &solution; - public : - source(dataCacheMap &cacheMap) : - xyz(cacheMap.get("XYZ",this)), - solution(cacheMap.get("Solution",this)) - {} - void _eval () { - int nQP = xyz().size1(); - if(_value.size1() != nQP) - _value = fullMatrix<double>(nQP,3); - for (int i=0; i<nQP; i++) { - double eta = solution(i,0); - double u = solution(i,1); - double v = solution(i,2); - double wind = 0.1*sin(xyz(i,1)/1e6*M_PI)/1000/h; - double f = 1e-4+xyz(i,1)*2e-11; - double gamma = 1e-6; - _value (i,0) = 0; - _value (i,1) = f*v + - gamma*u + wind; - _value (i,2) = -f*u - gamma*v ; - } +class dgConservationLawShallowWater2d::advection: public dataCacheDouble { + dataCacheDouble / + public: + advection(dataCacheMap &cacheMap): + sol(cacheMap.get("Solution",this)) + {}; + void _eval () { + int nQP = sol().size1(); + if(_value.size1() != nQP) + _value=fullMatrix<double>(nQP,9); + for(int i=0; i< nQP; i++) { + double eta = sol(i,0); + double u = sol(i,1); + double v = sol(i,2); + // flux_x + _value(i,0) = h*u; + _value(i,1) = g*eta; + _value(i,2) = 0; + // flux_y + _value(i,3) = h*v; + _value(i,4) = 0; + _value(i,5) = g*eta; + // flux_z + _value(i,6) = 0; + _value(i,7) = 0; + _value(i,8) = 0; } - }; - class riemann : public dataCacheDouble { - dataCacheDouble &normals, &solL, &solR; - public: - riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight): - normals(cacheMapLeft.get("Normals", this)), - solL(cacheMapLeft.get("Solution", this)), - solR(cacheMapRight.get("Solution", this)) - {}; - void _eval () { - int nQP = solL().size1(); - if(_value.size1() != nQP) - _value = fullMatrix<double>(nQP,6); - for(int i=0; i< nQP; i++) { - double nx = normals(0,i); - double ny = normals(1,i); - double unL = nx*solL(i,1) + ny*solL(i,2); - double unR = nx*solR(i,1) + ny*solR(i,2); - double utL = ny*solL(i,1) - nx*solL(i,2); - double utR = ny*solR(i,1) - nx*solR(i,2); - double etaR = solR(i,0); - double etaL = solL(i,0); - double sq_g_h = sqrt(g/h); - double etaRiemann = (etaL+etaR + (unL-unR)/sq_g_h)/2; - double unRiemann = (unL+unR + (etaL-etaR)*sq_g_h)/2; - double Fun = -g*etaRiemann; - double Feta = -h*unRiemann; - _value(i,0) = Feta; - _value(i,1) = Fun*nx; - _value(i,2) = Fun*ny; + } +}; - _value(i,3) = -_value(i,0); - _value(i,4) = -_value(i,1); - _value(i,5) = -_value(i,2); - } +class dgConservationLawShallowWater2d::source: public dataCacheDouble { + dataCacheDouble &xyz, &solution; + public : + source(dataCacheMap &cacheMap) : + xyz(cacheMap.get("XYZ",this)), + solution(cacheMap.get("Solution",this)) + {} + void _eval () { + int nQP = xyz().size1(); + if(_value.size1() != nQP) + _value = fullMatrix<double>(nQP,3); + for (int i=0; i<nQP; i++) { + double eta = solution(i,0); + double u = solution(i,1); + double v = solution(i,2); + double wind = 0.1*sin(xyz(i,1)/1e6*M_PI)/1000/h; + double f = 1e-4+xyz(i,1)*2e-11; + double gamma = 1e-6; + _value (i,0) = 0; + _value (i,1) = f*v + - gamma*u + wind; + _value (i,2) = -f*u - gamma*v ; } - }; - public: - dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const { - return new advection(cacheMap); - } - dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { - return new riemann(cacheMapLeft, cacheMapRight); - } - dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const { - return 0; - } - dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const { - return new source(cacheMap); } - dgConservationLawShallowWater2d() - { - _nbf = 3; // H U(=Hu) V(=Hv) +}; +class dgConservationLawShallowWater2d::riemann:public dataCacheDouble { + dataCacheDouble &normals, &solL, &solR; + public: + riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight): + normals(cacheMapLeft.get("Normals", this)), + solL(cacheMapLeft.get("Solution", this)), + solR(cacheMapRight.get("Solution", this)) + {}; + void _eval () { + int nQP = solL().size1(); + if(_value.size1() != nQP) + _value = fullMatrix<double>(nQP,6); + for(int i=0; i< nQP; i++) { + double nx = normals(0,i); + double ny = normals(1,i); + double unL = nx*solL(i,1) + ny*solL(i,2); + double unR = nx*solR(i,1) + ny*solR(i,2); + double utL = ny*solL(i,1) - nx*solL(i,2); + double utR = ny*solR(i,1) - nx*solR(i,2); + double etaR = solR(i,0); + double etaL = solL(i,0); + double sq_g_h = sqrt(g/h); + double etaRiemann = (etaL+etaR + (unL-unR)/sq_g_h)/2; + double unRiemann = (unL+unR + (etaL-etaR)*sq_g_h)/2; + double Fun = -g*etaRiemann; + double Feta = -h*unRiemann; + _value(i,0) = Feta; + _value(i,1) = Fun*nx; + _value(i,2) = Fun*ny; + + _value(i,3) = -_value(i,0); + _value(i,4) = -_value(i,1); + _value(i,5) = -_value(i,2); + } } }; -class dgBoundaryConditionShallowWater2dWall : public dgBoundaryCondition { +class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition { class term : public dataCacheDouble { dataCacheDouble &sol,&normals; public: @@ -135,16 +117,36 @@ class dgBoundaryConditionShallowWater2dWall : public dgBoundaryCondition { } }; public: - dgBoundaryConditionShallowWater2dWall(){} + boundaryWall(){} dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const { return new term(cacheMapLeft); } }; -dgBoundaryCondition *dgNewBoundaryConditionShallowWater2dWall() { - return new dgBoundaryConditionShallowWater2dWall(); +dataCacheDouble *dgConservationLawShallowWater2d::newConvectiveFlux( dataCacheMap &cacheMap) const { + return new advection(cacheMap); +} +dataCacheDouble *dgConservationLawShallowWater2d::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { + return new riemann(cacheMapLeft, cacheMapRight); +} +dataCacheDouble *dgConservationLawShallowWater2d::newDiffusiveFlux( dataCacheMap &cacheMap) const { + return 0; +} +dataCacheDouble *dgConservationLawShallowWater2d::newSourceTerm (dataCacheMap &cacheMap) const { + return new source(cacheMap); } -dgConservationLaw *dgNewConservationLawShallowWater2d() { - return new dgConservationLawShallowWater2d(); +dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(){ + return new boundaryWall(); } + +#ifdef HAVE_LUA +#include "Bindings.h" +const char dgConservationLawShallowWater2d::className[]="ShallowWater2d"; +const char dgConservationLawShallowWater2d::parentClassName[]="ConservationLaw"; +methodBinding *dgConservationLawShallowWater2d::methods[]={ + new methodBindingTemplate<dgConservationLawShallowWater2d,dgBoundaryCondition*>("newWallBoundary",&dgConservationLawShallowWater2d::newBoundaryWall), + 0 +}; +constructorBinding *dgConservationLawShallowWater2d::constructorMethod=new constructorBindingTemplate<dgConservationLawShallowWater2d>(); +#endif diff --git a/Solver/dgConservationLawShallowWater2d.h b/Solver/dgConservationLawShallowWater2d.h new file mode 100644 index 0000000000000000000000000000000000000000..c37c9525e2320800933ebaa75bf6ea491a3d8a89 --- /dev/null +++ b/Solver/dgConservationLawShallowWater2d.h @@ -0,0 +1,33 @@ +#ifndef _DG_CONSERVATION_LAW_SHALLOW_WATER_2D_ +#define _DG_CONSERVATION_LAW_SHALLOW_WATER_2D_ +#include "dgConservationLaw.h" + +#ifdef HAVE_LUA +class methodBinding; +class constructorBinding; +#endif +class dataCacheMap; + +class dgConservationLawShallowWater2d : public dgConservationLaw { + class advection; + class source; + class riemann; + class boundaryWall; + public: + dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const; + dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const; + dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const; + dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const; + dgConservationLawShallowWater2d() + { + _nbf = 3; // H U(=Hu) V(=Hv) + } + dgBoundaryCondition *newBoundaryWall(); +#ifdef HAVE_LUA + static const char className[]; + static const char parentClassName[]; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; +#endif +}; +#endif diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp index ab6447ed281d3b2b34a67fc283610e1bfb76dddd..440c8993f7fb89d8770c40d45f914cd5e78db5a5 100644 --- a/Solver/dgConservationLawWaveEquation.cpp +++ b/Solver/dgConservationLawWaveEquation.cpp @@ -1,4 +1,5 @@ #include "dgConservationLaw.h" +#include "dgConservationLawWaveEquation.h" #include "function.h" // dp/dt - rho*c^2 div(u,v) = 0 // du/dt + 1/rho dp/dx = 0 @@ -6,85 +7,81 @@ static double c=1; static double rho=1; -class dgConservationLawWaveEquation : public dgConservationLaw { - int _DIM; - class hyperbolicFlux : public dataCacheDouble { - dataCacheDouble / - const int _DIM,_nbf; - public: - hyperbolicFlux(dataCacheMap &cacheMap,int DIM): - sol(cacheMap.get("Solution",this)),_DIM(DIM),_nbf(DIM+1) - {}; - void _eval () { - int nQP = sol().size1(); - if(_value.size1() != nQP) - _value=fullMatrix<double>(nQP,3*_nbf); - _value.scale(0.); - for(int i=0; i< nQP; i++) { - const double p = sol(i,0); - for (int j=0;j<_DIM;j++){ - _value(i,0 +j*_nbf) = c*c*rho*sol(i,j+1); - _value(i,j+1+j*_nbf) = p/rho; - } +class dgConservationLawWaveEquation::hyperbolicFlux : public dataCacheDouble { + dataCacheDouble / + const int _DIM,_nbf; + public: + hyperbolicFlux(dataCacheMap &cacheMap,int DIM): + sol(cacheMap.get("Solution",this)),_DIM(DIM),_nbf(DIM+1) + {}; + void _eval () { + int nQP = sol().size1(); + if(_value.size1() != nQP) + _value=fullMatrix<double>(nQP,3*_nbf); + _value.scale(0.); + for(int i=0; i< nQP; i++) { + const double p = sol(i,0); + for (int j=0;j<_DIM;j++){ + _value(i,0 +j*_nbf) = c*c*rho*sol(i,j+1); + _value(i,j+1+j*_nbf) = p/rho; } } - }; - class riemann : public dataCacheDouble { - dataCacheDouble &normals, &solL, &solR; - const int _DIM,_nbf; - public: - riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, int DIM): - normals(cacheMapLeft.get("Normals", this)), - solL(cacheMapLeft.get("Solution", this)), - solR(cacheMapRight.get("Solution", this)), - _DIM(DIM),_nbf(DIM+1) - {}; - void _eval () { - int nQP = solL().size1(); - if(_value.size1() != nQP) - _value = fullMatrix<double>(nQP,2*_nbf); - for(int i=0; i< nQP; i++) { - const double n[3] = {normals(0,i),normals(1,i),normals(2,i)}; - double unL=0,unR=0; - for (int j=0;j<_DIM;j++){ - unL += n[j]*solL(i,j+1); - unR += n[j]*solR(i,j+1); - } - double pR = solR(i,0); - double pL = solL(i,0); + } +}; +class dgConservationLawWaveEquation::riemann : public dataCacheDouble { + dataCacheDouble &normals, &solL, &solR; + const int _DIM,_nbf; + public: + riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, int DIM): + normals(cacheMapLeft.get("Normals", this)), + solL(cacheMapLeft.get("Solution", this)), + solR(cacheMapRight.get("Solution", this)), + _DIM(DIM),_nbf(DIM+1) + {}; + void _eval () { + int nQP = solL().size1(); + if(_value.size1() != nQP) + _value = fullMatrix<double>(nQP,2*_nbf); + for(int i=0; i< nQP; i++) { + const double n[3] = {normals(0,i),normals(1,i),normals(2,i)}; + double unL=0,unR=0; + for (int j=0;j<_DIM;j++){ + unL += n[j]*solL(i,j+1); + unR += n[j]*solR(i,j+1); + } + double pR = solR(i,0); + double pL = solL(i,0); - double pRiemann = 0.5 * (pL+pR + (unL-unR)*(rho*c)); - double unRiemann = 0.5 * (unL+unR + (pL-pR)/(rho*c)); + double pRiemann = 0.5 * (pL+pR + (unL-unR)*(rho*c)); + double unRiemann = 0.5 * (unL+unR + (pL-pR)/(rho*c)); - double Fp = -rho*c*c*unRiemann; - double Fun = -pRiemann/rho; + double Fp = -rho*c*c*unRiemann; + double Fun = -pRiemann/rho; - _value(i,0) = Fp; - for (int j=0;j<_DIM;j++) - _value(i,j+1) = Fun*n[j]; - for (int j=0;j<_nbf;j++) - _value(i,j+_nbf) = -_value(i,j); - } + _value(i,0) = Fp; + for (int j=0;j<_DIM;j++) + _value(i,j+1) = Fun*n[j]; + for (int j=0;j<_nbf;j++) + _value(i,j+_nbf) = -_value(i,j); } - }; - public: - dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const { - return new hyperbolicFlux(cacheMap,_DIM); - } - dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { - return new riemann(cacheMapLeft, cacheMapRight,_DIM); - } - dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const { - return 0; - } - dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const { - return 0; - } - dgConservationLawWaveEquation(int DIM) : _DIM(DIM) - { - _nbf = 1 + _DIM; // H U(=Hu) V(=Hv) } }; +dataCacheDouble *dgConservationLawWaveEquation::newConvectiveFlux( dataCacheMap &cacheMap) const { + return new hyperbolicFlux(cacheMap,_DIM); +} +dataCacheDouble *dgConservationLawWaveEquation::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const { + return new riemann(cacheMapLeft, cacheMapRight,_DIM); +} +dataCacheDouble *dgConservationLawWaveEquation::newDiffusiveFlux( dataCacheMap &cacheMap) const { + return 0; +} +dataCacheDouble *dgConservationLawWaveEquation::newSourceTerm (dataCacheMap &cacheMap) const { + return 0; +} +dgConservationLawWaveEquation::dgConservationLawWaveEquation(int DIM) : _DIM(DIM) +{ + _nbf = 1 + _DIM; // H U(=Hu) V(=Hv) +} class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition { int _DIM; @@ -93,19 +90,19 @@ class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition { dataCacheDouble &sol,&normals; public: term(dataCacheMap &cacheMap, int DIM): - sol(cacheMap.get("Solution",this)), - normals(cacheMap.get("Normals",this)), - _DIM(DIM){} + sol(cacheMap.get("Solution",this)), + normals(cacheMap.get("Normals",this)), + _DIM(DIM){} void _eval () { int nQP = sol().size1(); if(_value.size1() != nQP) _value = fullMatrix<double>(nQP,_DIM+1); for(int i=0; i< nQP; i++) { - const double n[3] = {normals(0,i),normals(1,i),normals(2,i)}; + const double n[3] = {normals(0,i),normals(1,i),normals(2,i)}; double p = sol(i,0); _value(i,0) = 0; - for (int j=0;j<_DIM;j++) - _value(i,j+1) = -p/rho*n[j]; + for (int j=0;j<_DIM;j++) + _value(i,j+1) = -p/rho*n[j]; } } }; @@ -115,10 +112,14 @@ class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition { return new term(cacheMapLeft,_DIM); } }; - -dgConservationLaw *dgNewConservationLawWaveEquation(int DIM) { - return new dgConservationLawWaveEquation(DIM); -} -dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall(int DIM) { - return new dgBoundaryConditionWaveEquationWall(DIM); +dgBoundaryCondition *dgConservationLawWaveEquation::newBoundaryWall()const{ + return new dgBoundaryConditionWaveEquationWall(_DIM); } + +#include "Bindings.h" +const char *dgConservationLawWaveEquation::className = "ConservationLawWaveEquation"; +const char *dgConservationLawWaveEquation::parentClassName = "ConservationLaw"; +methodBinding *dgConservationLawWaveEquation::methods[] ={ + new methodBindingTemplate<const dgConservationLawWaveEquation,dgBoundaryCondition*>("newWallBoundary",&dgConservationLawWaveEquation::newBoundaryWall), +0}; +constructorBinding *dgConservationLawWaveEquation::constructorMethod=new constructorBindingTemplate<dgConservationLawWaveEquation,int>(); diff --git a/Solver/dgConservationLawWaveEquation.h b/Solver/dgConservationLawWaveEquation.h new file mode 100644 index 0000000000000000000000000000000000000000..adb3c2d43e94ce8de1d0162c25749bdabb71f880 --- /dev/null +++ b/Solver/dgConservationLawWaveEquation.h @@ -0,0 +1,22 @@ +#ifndef _DG_CONSERVATION_LAW_WAVE_EQUATION_H_ +#define _DG_CONSERVATION_LAW_WAVE_EQUATION_H_ +#include "dgConservationLaw.h" +class methodBinding; +class constructorBinding; +class dgConservationLawWaveEquation : public dgConservationLaw { + int _DIM; + class hyperbolicFlux; + class riemann; + public: + dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const ; + dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const; + dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const ; + dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const ; + dgBoundaryCondition *newBoundaryWall () const; + dgConservationLawWaveEquation(int DIM); + static const char *className; + static const char *parentClassName; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; +}; +#endif diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp index b222417f8a690fae15362657ac25068838c968bf..c1cea4499dab60bd88bca5d3d72c90a7320472b5 100644 --- a/Solver/dgGroupOfElements.cpp +++ b/Solver/dgGroupOfElements.cpp @@ -485,4 +485,10 @@ void dgGroupOfFaces::mapFromInterface ( int nFields, } } } - +/* +const char luaTest::className[]="luaTest"; +methodBinding *luaTest::methods[] = { + new LunaSignature<luaTest,int,int,int,int>("print",&luaTest::print), + 0}; +constructorBinding *luaTest::constructorMethod = new LunaConstructor<luaTest,int>; +*/ diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h index ecf50829ee23f1df497c307bc96b29c7bc6b3dd2..bbea72d73d577ad0b8669351304b13d53cb64bb6 100644 --- a/Solver/dgGroupOfElements.h +++ b/Solver/dgGroupOfElements.h @@ -177,6 +177,4 @@ public: void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight); const polynomialBasis * getPolynomialBasis() const {return _fsFace;} }; - - #endif diff --git a/Solver/dgMainLua.cpp b/Solver/dgMainLua.cpp index cb0ac1845b3172ecdb43fa79ba7a119aee64dbe2..e21c5df73c629d59becec2ea96f6df41b76c3437 100644 --- a/Solver/dgMainLua.cpp +++ b/Solver/dgMainLua.cpp @@ -1,9 +1,19 @@ #include <iostream> #include "Gmsh.h" -#include "LuaBindings_Geo.h" +#include "Bindings.h" #include "dgSystemOfEquations.h" #include "luaFunction.h" #include "function.h" +#include "dgGroupOfElements.h" +#include "dgConservationLawShallowWater2d.h" +#include "dgConservationLawAdvection.h" +#include "dgConservationLawPerfectGas.h" +#include "dgConservationLawWaveEquation.h" +extern "C" { + #include "lua.h" + #include "lualib.h" + #include "lauxlib.h" +} void report_errors(lua_State *L, int status) { if ( status!=0 ) { @@ -18,17 +28,26 @@ int main(int argc, char *argv[]) lua_State *L = lua_open(); luaopen_base(L); luaopen_table(L); + luaopen_os(L); //luaopen_io(L); luaopen_string(L); luaopen_math(L); luaopen_debug(L); - // Register GModel bindings - LuaGModel::Register(L); - dgSystemOfEquations::Register(L); - fullMatrix<double>::Register(L); + // Register Lua bindings + classBinding<GModel>::Register(L); + classBinding<dgSystemOfEquations>::Register(L); + classBinding<dgBoundaryCondition>::Register(L); + classBinding<dgConservationLaw>::Register(L); + classBinding<dgConservationLawShallowWater2d>::Register(L); + classBinding<dgConservationLawAdvection>::Register(L); + classBinding<dgConservationLawWaveEquation>::Register(L); + classBinding<dgPerfectGasLaw2d>::Register(L); + classBinding<fullMatrix<double> >::Register(L); + classBinding<function>::Register(L); + classBinding<functionLua>::Register(L); + classBinding<functionConstant>::Register(L); function::registerDefaultFunctions(); - RegisterFunctions(L); int s = luaL_loadfile(L, argv[1]); diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp index 8ba9b2b7b1f2084bc2b13fea2903389da3793bd1..9035150742661bb60d94358b03570fdbe1dc349a 100644 --- a/Solver/dgSystemOfEquations.cpp +++ b/Solver/dgSystemOfEquations.cpp @@ -1,135 +1,67 @@ #include <stdio.h> #include <stdlib.h> #include "dgSystemOfEquations.h" -#include "LuaBindings_Geo.h" #include "function.h" #include "MElement.h" #include "PView.h" #include "PViewData.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, RK44), - _method(dgSystemOfEquations, exportSolution), - _method(dgSystemOfEquations, setConservationLaw), - _method(dgSystemOfEquations, setOrder), - _method(dgSystemOfEquations, setup), - _method(dgSystemOfEquations, addBoundaryCondition), - _method(dgSystemOfEquations, L2Projection), - {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); -} class dgConservationLawL2Projection : public dgConservationLaw { std::string _functionName; public: - dgConservationLawL2Projection(const std::string & functionName, - dgConservationLaw &_claw) : + dgConservationLawL2Projection(const std::string & functionName, dgConservationLaw &_claw) : _functionName(functionName) { _nbf =_claw.nbFields(); } dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const { - //return new gaussian(cacheMap,0.2,0.3); return &cacheMap.get(_functionName); } }; - -// 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 : _gm->getNumFaces() ? 2 : 1; +dgSystemOfEquations::dgSystemOfEquations(GModel *gm){ + _gm=gm; + _dimension = _gm->getNumRegions() ? 3 : _gm->getNumFaces() ? 2 : 1 ; _solution = 0; } -// set the conservation law as a string (for now) -int dgSystemOfEquations::setConservationLaw(lua_State *L){ - _claw = 0; - //int argc = (int)luaL_checknumber(L,0); - _cLawName = std::string (luaL_checkstring(L, 1)); - if (_cLawName == "WaveEquation") - _claw = dgNewConservationLawWaveEquation(_dimension); - else if (_cLawName == "ShallowWater2d") - _claw = dgNewConservationLawShallowWater2d(); - else if (_cLawName == "PerfectGas2d") - _claw = dgNewPerfectGasLaw2d(); - else if (_cLawName == "AdvectionDiffusion"){ - _claw = dgNewConservationLawAdvection(luaL_checkstring(L,2),luaL_checkstring(L,3)); - } - if (!_claw)throw; - return 0; -} // set the order of interpolation -int dgSystemOfEquations::setOrder(lua_State *L){ - _order = luaL_checkint(L, 1); +void dgSystemOfEquations::setOrder(int o){ + _order = o; } // add a boundary Condition -int dgSystemOfEquations::addBoundaryCondition (lua_State *L){ - std::string physicalName(luaL_checkstring(L, 1)); - std::string bcName(luaL_checkstring(L, 2)); - //generic boundary conditions - if (bcName == "0Flux"){ - _claw->addBoundaryCondition(physicalName,dgBoundaryCondition::new0FluxCondition(*_claw)); - } - else if (bcName == "OutsideValues"){ - _claw->addBoundaryCondition(physicalName,dgBoundaryCondition::newOutsideValueCondition(*_claw,luaL_checkstring(L,3))); - } - //specific boundary conditions - else if (_cLawName == "WaveEquation"){ - if (bcName == "Wall"){ - _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionWaveEquationWall(_dimension)); - } - else throw; - } - else if (_cLawName == "ShallowWater2d"){ - if (bcName == "Wall"){ - _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionShallowWater2dWall()); - } - else throw; - } - else if (_cLawName == "PerfectGas2d"){ - if (bcName == "Wall"){ - _claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionPerfectGasLaw2dWall()); - } - else if (bcName == "FreeStream"){ - std::string freeStreamName(luaL_checkstring(L, 3)); - _claw->addBoundaryCondition(physicalName, - dgNewBoundaryConditionPerfectGasLaw2dFreeStream(freeStreamName)); - } - else throw; - } - else throw; +void dgSystemOfEquations::setConservationLaw (dgConservationLaw *law){ + _claw=law; } +#ifdef HAVE_LUA +#include "Bindings.h" +const char dgSystemOfEquations::className[] = "dgSystemOfEquations"; +const char dgSystemOfEquations::parentClassName[] = ""; +constructorBinding *dgSystemOfEquations::constructorMethod = new constructorBindingTemplate<dgSystemOfEquations,GModel*>(); +methodBinding *dgSystemOfEquations::methods[]={ + new methodBindingTemplate<dgSystemOfEquations,void,int>("setOrder",&dgSystemOfEquations::setOrder), + new methodBindingTemplate<dgSystemOfEquations,void,dgConservationLaw*>("setConservationLaw",&dgSystemOfEquations::setConservationLaw), + new methodBindingTemplate<dgSystemOfEquations,void>("setup",&dgSystemOfEquations::setup), + new methodBindingTemplate<dgSystemOfEquations,void,std::string>("exportSolution",&dgSystemOfEquations::exportSolution), + new methodBindingTemplate<dgSystemOfEquations,void,std::string>("L2Projection",&dgSystemOfEquations::L2Projection), + new methodBindingTemplate<dgSystemOfEquations,double,double>("RK44",&dgSystemOfEquations::RK44), + 0}; + // do a L2 projection -int dgSystemOfEquations::L2Projection (lua_State *L){ - dgConservationLawL2Projection Law(std::string(luaL_checkstring(L, 1)),*_claw); +void dgSystemOfEquations::L2Projection (std::string functionName){ + dgConservationLawL2Projection Law(functionName,*_claw); for (int i=0;i<_elementGroups.size();i++){ _algo->residualVolume(Law,*_elementGroups[i],*_solution->_dataProxys[i],*_rightHandSide->_dataProxys[i]); _algo->multAddInverseMassMatrix(*_elementGroups[i],*_rightHandSide->_dataProxys[i],*_solution->_dataProxys[i]); } - return 0; } // ok, we can setup the groups and create solution vectors -int dgSystemOfEquations::setup(lua_State *L){ +void dgSystemOfEquations::setup(){ if (!_claw) throw; _algo->buildGroups(_gm, _dimension, @@ -143,18 +75,13 @@ int dgSystemOfEquations::setup(lua_State *L){ } -int dgSystemOfEquations::RK44(lua_State *L){ - double dt = luaL_checknumber(L, 1); +double dgSystemOfEquations::RK44(double dt){ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide); - double normSolution = _solution->_data->norm(); - lua_pushnumber (L, normSolution); - return 1; + return _solution->_data->norm(); } -int dgSystemOfEquations::exportSolution(lua_State *L){ - std::string outputFile(luaL_checkstring(L, 1)); +void dgSystemOfEquations::exportSolution(std::string outputFile){ export_solution_as_is(outputFile); - return 0; } #endif // HAVE_LUA diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h index 45c5108a6a9fe5cc363dba3b7ec0094dbabd993f..7671f76d81383bb69f0da1959deef16d42e18084 100644 --- a/Solver/dgSystemOfEquations.h +++ b/Solver/dgSystemOfEquations.h @@ -8,17 +8,6 @@ #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 - struct dgDofContainer { private: dgDofContainer (const dgDofContainer&); @@ -30,6 +19,8 @@ public: ~dgDofContainer (); }; +class methodBinding; +class constructorBinding; class dgSystemOfEquations { // the mesh and the model @@ -54,20 +45,19 @@ class dgSystemOfEquations { std::vector<dgGroupOfFaces*> _boundaryGroups; dgSystemOfEquations(const dgSystemOfEquations &) {} public: - // lua stuff -#if defined(HAVE_LUA) + void setOrder (int order); // set the polynomial order + void setConservationLaw (dgConservationLaw *law); // set the conservationLaw + dgSystemOfEquations(GModel *_gm); + void setup (); // setup the groups and allocate + void exportSolution (std::string filename); // export the solution + double RK44 (double dt); + void L2Projection (std::string functionName); // assign the solution to a given function + static const char className[]; - static Luna<dgSystemOfEquations>::RegType methods[]; - static void Register(lua_State *L); - 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 - int L2Projection (lua_State *L); // assign the solution to a given function - int RK44 (lua_State *L); // assign the solution to a given function - dgSystemOfEquations(lua_State *L); -#endif // HAVE LUA + static const char parentClassName[]; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; + inline const fullMatrix<double> getSolutionProxy (int iGroup, int iElement){ return fullMatrix<double> ( *_solution->_dataProxys [iGroup] , iElement * _claw->nbFields(),_claw->nbFields()); diff --git a/Solver/function.cpp b/Solver/function.cpp index 9fed214400df13e83b0b74162a6bc2b2fd30d908..f092620340b91f1859f14dab80308ae105f54522 100644 --- a/Solver/function.cpp +++ b/Solver/function.cpp @@ -1,6 +1,7 @@ #include "function.h" #include "SPoint3.h" #include "MElement.h" +#include <sstream> // dataCache members @@ -84,11 +85,13 @@ class functionXYZ : public function { private: dataCacheElement &_element; dataCacheDouble &_uvw; + int count; public: data(dataCacheMap *m) : dataCacheDouble(m->getNbEvaluationPoints(),3), _element(m->getElement(this)), _uvw(m->get("UVW", this)) - {} + { + } void _eval() { for(int i = 0; i < _uvw().size1(); i++){ @@ -99,6 +102,8 @@ class functionXYZ : public function { _value(i, 2) = p.z(); } } + ~data(){ + } }; public: dataCacheDouble *newDataCache(dataCacheMap *m) @@ -107,39 +112,52 @@ class functionXYZ : public function { } }; + // constant values copied over each line -class functionConstant : public function { - private : - class data : public dataCacheDouble { - const functionConstant *_function; - public: - data(const functionConstant * function,dataCacheMap *m): - dataCacheDouble(m->getNbEvaluationPoints(),function->_source.size1()){ - _function = function; - } - void _eval() { - for(int i=0;i<_value.size1();i++) - for(int j=0;j<_function->_source.size1();j++) - _value(i,j)=_function->_source(j,0); - } - }; - fullMatrix<double> _source; +class functionConstant::data : public dataCacheDouble { + const functionConstant *_function; public: - dataCacheDouble *newDataCache(dataCacheMap *m) - { - return new data(this,m); - } - functionConstant(const fullMatrix<double> &source){ - _source = source; - } + data(const functionConstant * function,dataCacheMap *m): + dataCacheDouble(m->getNbEvaluationPoints(),function->_source.size1()){ + _function = function; + } + void _eval() { + for(int i=0;i<_value.size1();i++) + for(int j=0;j<_function->_source.size1();j++) + _value(i,j)=_function->_source(j,0); + } }; - -function *function::newFunctionConstant(const fullMatrix<double> &source){ - return new functionConstant(source); +dataCacheDouble *functionConstant::newDataCache(dataCacheMap *m) +{ + return new data(this,m); } +functionConstant::functionConstant(const fullMatrix<double> *source){ + _source = *source; + static int c=0; + std::ostringstream oss; + oss<<"FunctionConstant_"<<c++; + _name = oss.str(); + function::add(_name,this); +} + +#include "Bindings.h" +const char *functionConstant::className="FunctionConstant"; +const char *functionConstant::parentClassName="Function"; +methodBinding *functionConstant::methods[]={0}; +constructorBinding *functionConstant::constructorMethod=new constructorBindingTemplate<functionConstant,fullMatrix<double>*>(); + + + + void function::registerDefaultFunctions() { function::add("XYZ", new functionXYZ); } +const char *function::className="Function"; +const char *function::parentClassName=""; +methodBinding *function::methods[]={new methodBindingTemplate<const function,std::string>("getName",&function::getName), +0}; +constructorBinding *function::constructorMethod=0; + diff --git a/Solver/function.h b/Solver/function.h index 89d775ebc23070aa2820942d150adec34aa40014..666c60686eec06bbe48acc35070b695c2b7fa44d 100644 --- a/Solver/function.h +++ b/Solver/function.h @@ -114,6 +114,7 @@ class dataCacheDouble : public dataCache { class function { private: static std::map<std::string, function*> _allFunctions; + protected: std::string _name; public: static void registerDefaultFunctions(); @@ -122,8 +123,13 @@ class function { virtual dataCacheDouble *newDataCache(dataCacheMap *m) =0; - //we need a parser for this - static function *newFunctionConstant(const fullMatrix<double> &source); + inline std::string getName()const {return _name;} + + static const char *className; + static const char *parentClassName; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; + virtual ~function(){}; }; // A special node in the dependency tree for which all the leafs @@ -149,8 +155,7 @@ class dataCacheMap { std::map<std::string, dataCacheDouble*> _cacheDoubleMap; class providedDataDouble : public dataCacheDouble // for data provided by the algorithm and that does not have an _eval function - // (typically "UVW") this class is not stricly necessary, we could write - // a function for each case + // (typically "UVW") { void _eval() {throw;}; public: @@ -166,4 +171,18 @@ class dataCacheMap { inline int getNbEvaluationPoints(){return _nbEvaluationPoints;} ~dataCacheMap(); }; +class functionConstant : public function { + public: + class data ; + fullMatrix<double> _source; + dataCacheDouble *newDataCache(dataCacheMap *m); + functionConstant(const fullMatrix<double> *source); + static const char *className; + static const char *parentClassName; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; + ~functionConstant(){ + printf("delete fc\n"); + } +}; #endif diff --git a/Solver/luaFunction.cpp b/Solver/luaFunction.cpp index fd553a2a9102eaf7b4a4e2b6955212a286ebdb23..4e7083cfd5f23391bfc0a16f6417045ec20b9260 100644 --- a/Solver/luaFunction.cpp +++ b/Solver/luaFunction.cpp @@ -1,97 +1,54 @@ +#define _FUNCTION_LUA_H #include "luaFunction.h" #if defined(HAVE_LUA) +#include <sstream> #include <string> #include <vector> #include "function.h" +#include "Bindings.h" // function that is defined in Lua -class functionLua : public function { - lua_State *_L; - std::string _luaFunctionName; - std::vector<std::string> _dependenciesName; - int _nbCol; - private: - class data : public dataCacheDouble{ - private: - const functionLua *_function; - std::vector<dataCacheDouble *> _dependencies; - public: - data(const functionLua *f, dataCacheMap *m): - dataCacheDouble(m->getNbEvaluationPoints(),f->_nbCol), - _function(f) - { - _dependencies.resize ( _function->_dependenciesName.size()); - for (int i=0;i<_function->_dependenciesName.size();i++) - _dependencies[i] = &m->get(_function->_dependenciesName[i],this); - } - void _eval() - { - lua_getfield(_function->_L, LUA_GLOBALSINDEX, _function->_luaFunctionName.c_str()); - for (int i=0;i< _dependencies.size();i++){ - const fullMatrix<double> *data = &(*_dependencies[i])(); - lua_pushlightuserdata (_function->_L, (void*) data); - } - lua_pushlightuserdata (_function->_L, &_value); - lua_call(_function->_L,_dependencies.size()+1,0); /* call Lua function */ - } - }; - public: - functionLua (int nbCol, - std::string &luaFunctionName, - std::vector<std::string> &dependenciesName, - lua_State *L) - : _luaFunctionName(luaFunctionName), _dependenciesName(dependenciesName),_L(L),_nbCol(nbCol) +class functionLua::data : public dataCacheDouble{ + private: + const functionLua *_function; + std::vector<dataCacheDouble *> _dependencies; + public: + data(const functionLua *f, dataCacheMap *m): + dataCacheDouble(m->getNbEvaluationPoints(),f->_nbCol), + _function(f) { + _dependencies.resize ( _function->_dependenciesName.size()); + for (int i=0;i<_function->_dependenciesName.size();i++) + _dependencies[i] = &m->get(_function->_dependenciesName[i],this); } - - dataCacheDouble *newDataCache(dataCacheMap *m) + void _eval() { - return new data(this,m); + lua_getfield(_function->_L, LUA_GLOBALSINDEX, _function->_luaFunctionName.c_str()); + for (int i=0;i< _dependencies.size();i++){ + const fullMatrix<double> *data = &(*_dependencies[i])(); + classBinding<const fullMatrix<double> >::push(_function->_L,data,false); + } + classBinding<const fullMatrix<double> >::push(_function->_L,&_value,false); + lua_call(_function->_L,_dependencies.size()+1,0); /* call Lua function */ } }; - -static int newLuaFunction (lua_State *L){ - int n = lua_gettop(L); - std::vector<std::string> dependenciesName; - int nbFields = luaL_checkinteger(L, 1); - std::string luaFunctionName = std::string(luaL_checkstring(L, 2)); - for (int i=3;i<=n;i++) - dependenciesName.push_back(std::string(luaL_checkstring(L, i))); - int ITER = 1; - std::string functionName = luaFunctionName; - while (function::get(functionName, true)){ - char toto[256]; - sprintf(toto,"%d",ITER); - functionName = luaFunctionName + toto; - } - function::add(functionName,new functionLua(nbFields,luaFunctionName,dependenciesName,L)); - lua_pushstring(L, functionName.c_str()); - - return 1; +functionLua::functionLua (int nbCol, std::string &luaFunctionName, std::vector<std::string> &dependencies, lua_State *L) + : _luaFunctionName(luaFunctionName), _dependenciesName(dependencies),_L(L),_nbCol(nbCol) +{ + static int c=0; + std::ostringstream oss; + oss<<"luaFunction_"<<c++; + _name = oss.str(); + function::add(_name,this); } - -static int newConstantFunction (lua_State *L){ - fullMatrix<double> * _ud = Luna<fullMatrix<double> >::check(L, 1); - int ITER = 1; - std::string functionName = "ConstantFunction"; - while (function::get(functionName, true)){ - char toto[256]; - sprintf(toto,"%d",ITER); - functionName = functionName + toto; - } - - function::add(functionName,function::newFunctionConstant(*_ud)); - lua_pushstring(L, functionName.c_str()); - return 1; +dataCacheDouble *functionLua::newDataCache(dataCacheMap *m) +{ + return new data(this,m); } -int RegisterFunctions (lua_State *L) { - luaL_reg fcts[] = {{"lua", newLuaFunction }, - {"constant", newConstantFunction }, - {0,0}}; - luaL_register(L, "createFunction", fcts); -}; - - +const char *functionLua::className="FunctionLua"; +const char *functionLua::parentClassName="Function"; +methodBinding *functionLua::methods[]={0}; +constructorBinding *functionLua::constructorMethod=new constructorBindingTemplate<functionLua,int,std::string,std::vector<std::string> ,lua_State*>(); #endif // HAVE LUA diff --git a/Solver/luaFunction.h b/Solver/luaFunction.h index bc161db373b85a9b23399b9a0b661c93ef3644e8..63e66c937e420ccaf2f229bb854404788d42a68a 100644 --- a/Solver/luaFunction.h +++ b/Solver/luaFunction.h @@ -2,12 +2,24 @@ #ifndef _LUA_FUNCTION_H_ #define _LUA_FUNCTION_H_ #if defined(HAVE_LUA) -// include lua stuff -extern "C" { -#include "lua.h" -#include "lauxlib.h" -#include "lualib.h" -} -int RegisterFunctions (lua_State *L); +#include "function.h" +class lua_State; +#include <string> +#include <vector> +class functionLua : public function { + lua_State *_L; + std::string _luaFunctionName; + std::vector<std::string> _dependenciesName; + int _nbCol; + class data; + public: + functionLua (int nbCol, std::string &luaFunctionName, std::vector<std::string> &dependenciesName, lua_State *L); + + dataCacheDouble *newDataCache(dataCacheMap *m); + static const char *className; + static const char *parentClassName; + static methodBinding *methods[]; + static constructorBinding *constructorMethod; +}; #endif // HAVE LUA #endif // _LUA_FUNCTION_H_