Skip to content
Snippets Groups Projects
Commit 2b2e585f authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dg : move function (and numpy dependency) from gmsh/Solver to projects/dg

parent fcceb044
Branches
Tags
No related merge requests found
......@@ -41,7 +41,6 @@ option(ENABLE_MPEG_ENCODE "Enable built-in MPEG encoder" ON)
option(ENABLE_MPI "Enable MPI parallelization" OFF)
option(ENABLE_MSVC_STATIC_RUNTIME "Use static Visual C++ runtime" OFF)
option(ENABLE_NATIVE_FILE_CHOOSER "Enable native file chooser in GUI" ON)
option(ENABLE_NUMPY "Enable Numpy function in solvers (requires SWIG)" ON)
option(ENABLE_NETGEN "Enable Netgen mesh generator" ON)
option(ENABLE_OCC "Enable Open CASCADE geometrical models" ON)
option(ENABLE_OSMESA "Use OSMesa for offscreen rendering" OFF)
......@@ -810,12 +809,6 @@ if(ENABLE_SWIG)
message("WARNING: Python bindings require SWIG >= 2: disabling Python")
else(SWIG_MAJOR_VERSION EQUAL 1)
set_config_option(HAVE_SWIG "Swig")
if(ENABLE_NUMPY)
find_path(NUMPY_INC "numpyconfig.h" PATH_SUFFIXES include/numpy)
if(NUMPY_INC)
set_config_option(HAVE_NUMPY "Numpy")
endif(NUMPY_INC)
endif(ENABLE_NUMPY)
endif(SWIG_MAJOR_VERSION EQUAL 1)
endif(SWIG_FOUND AND PYTHONLIBS_FOUND)
endif(ENABLE_SWIG)
......
......@@ -38,7 +38,6 @@
#cmakedefine HAVE_NETGEN
#cmakedefine HAVE_NO_SOCKLEN_T
#cmakedefine HAVE_NO_VSNPRINTF
#cmakedefine HAVE_NUMPY
#cmakedefine HAVE_OCC
#cmakedefine HAVE_OPENGL
#cmakedefine HAVE_OSMESA
......
......@@ -14,8 +14,6 @@ set(SRC
SElement.cpp
eigenSolver.cpp
multiscaleLaplace.cpp
function.cpp
functionDerivator.cpp
functionSpace.cpp
filters.cpp
sparsityPattern.cpp
......
......@@ -22,7 +22,6 @@
#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
#include "function.h"
#endif
static void printLinearSystem(linearSystemCSRTaucs<double> *lsys)
......
......@@ -12,7 +12,6 @@
#include "dofManager.h"
#include "simpleFunction.h"
#include "functionSpace.h"
#include "function.h"
class GModel;
class PView;
......
This diff is collapsed.
#ifndef _FUNCTION_H_
#define _FUNCTION_H_
#include "fullMatrix.h"
#include <map>
#include <set>
#include <list>
#include <string>
#include <vector>
class dataCacheDouble;
class dataCacheMap;
class dgDataCacheMap;
class function;
class functionConstant;
class functionReplace;
class functionReplaceCache;
class MElement;
// An abstract interface to functions
class function {
public:
// additionnal types
class dependency {
public:
unsigned iMap;
const function *f;
dependency(int iMap_, const function *f_) {
iMap = iMap_;
f = f_;
}
bool operator < (const dependency &d) const {
return (d.iMap <iMap || d.f < f);
}
};
class argument {
public:
int iMap; // id of the dataCacheMap, e.g. on interfaces
const function *f;
fullMatrix<double> *val;
argument(fullMatrix<double> &v, int iMap_, const function *f_) {
val = &v;
iMap = iMap_;
f = f_;
}
};
// data
int _nbCol;
bool _invalidatedOnElement;
std::vector<functionReplace*> _functionReplaces;
std::vector<argument> arguments;
std::set<dependency> dependencies;
private:
static functionConstant *_timeFunction;
static functionConstant *_dtFunction;
static functionConstant *_dtSubFunction;
public:
function(int nbCol, bool invalidatedOnElement = true)
: _nbCol(nbCol), _invalidatedOnElement(invalidatedOnElement) {}
virtual ~function(){}
void addFunctionReplace(functionReplace &fr);
void setArgument(fullMatrix<double> &v, const function *f, int iMap = 0);
void setArgumentWrapped(fullMatrix<double> &v, const function *f, int expectedNbCol, std::string funcName="a function", int iMap=0);
virtual void call(dataCacheMap *m, fullMatrix<double> &res) = 0;
virtual void registerInDataCacheMap(dataCacheMap *m, dataCacheDouble *d) {}
// get or print information
inline bool isInvalitedOnElement() { return _invalidatedOnElement; }
inline int getNbCol() const {if(_nbCol ==0) Msg::Error("Cannot ask nbCol of functionSolution"); return _nbCol; }
static functionConstant *getTime();
static functionConstant *getDT();
static functionConstant *getSubDT();
static function *getSolution();
static function *getCoordinates();
static function *getSolutionGradient();
static function *getNormals();
void printDep()
{
for (std::set<dependency>::iterator it = dependencies.begin();
it != dependencies.end(); it++)
printf("%i %p\n", it->iMap, it->f);
}
};
// Additionnal class to get solution
class functionSolution : public function {
static functionSolution *_instance;
// constructor is private only 1 instance can exists, call get to
// access the instance
functionSolution() : function(0) {};
public:
void call(dataCacheMap *m, fullMatrix<double> &sol)
{
Msg::Error("A function requires the solution but this algorithm does "
"not provide the solution");
throw;
}
void setNbFields(int nbFields){ _nbCol = nbFields; }
static functionSolution *get()
{
if(!_instance)
_instance = new functionSolution();
return _instance;
}
};
class functionSolutionGradient : public function {
static functionSolutionGradient *_instance;
// constructor is private only 1 instance can exists, call get to
// access the instance
functionSolutionGradient();
public:
void call(dataCacheMap *m, fullMatrix<double> &sol);
static function *get();
};
class functionReplaceCache {
public:
dataCacheMap *map;
std::vector <dataCacheDouble*> toReplace;
std::vector <dataCacheDouble*> toCompute;
};
class functionReplace {
friend class dataCacheMap;
friend class dataCacheDouble;
public:
int _nChildren;
function *_master;
functionReplaceCache *currentCache;
std::set <function::dependency> _replaced;
std::set <function::dependency> _fromParent;
std::vector <function::argument> _toReplace;
std::vector <function::argument> _toCompute;
functionReplace();
void get(fullMatrix<double> &v, const function *, int iMap = 0);
void replace(fullMatrix<double> &v, const function *, int iMap = 0);
void compute();
void addChild();
};
// A dataCache when the value is a matrix of double. The user should
// provide the number of rows by evaluating points and the number of
// columns; then the size of the matrix is automatically adjusted
class dataCacheDouble {
friend class dataCacheMap;
friend class dgDataCacheMap;
friend class functionReplace;
std::vector<dataCacheDouble*> _directDependencies;
function *_function;
dataCacheMap &_cacheMap;
std::vector<functionReplaceCache> functionReplaceCaches;
protected:
// pointers to all of the dataCache depending on me
std::set<dataCacheDouble*> _dependOnMe;
std::set<dataCacheDouble*> _iDependOn;
fullMatrix<double> _value;
bool _valid;
dataCacheDouble(dataCacheMap *,function *f);
// do the actual computation and put the result into _value
void _eval();
void resize(int nrow);
public:
inline function * getFunction() {return _function;}
//set the value (without computing it by _eval) and invalidate the
// dependencies this function is needed to be able to pass the
// _value to functions like gemm or mult but you cannot keep the
// reference to the _value, you should always use the set function
// to modify the _value take care if you use this to set a proxy you
// must ensure that the value pointed to are not modified without
// further call to set because the dependencies won't be invalidate
inline fullMatrix<double> &set()
{
if(_valid) {
for(std::set<dataCacheDouble*>::iterator it = _dependOnMe.begin();
it!=_dependOnMe.end(); it++)
(*it)->_valid=false;
}
_valid = true;
return _value;
}
// access _value and compute it if necessary
inline const fullMatrix<double> &get()
{
if(!_valid)
_eval();
return _value;
}
// dataCacheMap is the only one supposed to call this
inline bool somethingDependOnMe() { return !_dependOnMe.empty(); }
inline bool doIDependOn(dataCacheDouble &other)
{
return (_iDependOn.find(&other)!=_iDependOn.end());
}
inline std::vector<dataCacheDouble*> & getDirectDependencies() { return _directDependencies; }
};
class dataCacheMap {
const function *_functionSolution, *_functionSolutionGradient, *_functionCoordinates, *_containerSolution, *_containerSolutionGradient;
//handle function solution and funciton solution gradient
const function * _translate (const function *) const;
public:
dataCacheMap *_parent;
std::list<dataCacheMap*> _children;
std::vector<dataCacheMap*> _secondaryCaches;
int _nbEvaluationPoints;
std::map<const function*, dataCacheDouble*> _cacheDoubleMap;
std::set<dataCacheDouble*> _allDataCaches;
std::vector<dataCacheDouble*> _toInvalidateOnElement;
MElement *_element;
dataCacheMap() {
_functionSolution = _functionSolutionGradient = _functionCoordinates = _containerSolution = _containerSolutionGradient = NULL;
_nbEvaluationPoints = 0;
_parent = NULL;
}
~dataCacheMap();
void addDataCacheDouble(dataCacheDouble *data, bool invalidatedOnElement)
{
_allDataCaches.insert(data);
if(invalidatedOnElement)
_toInvalidateOnElement.push_back(data);
}
virtual dgDataCacheMap *asDgDataCacheMap()
{
Msg::Error("I'm not a dgDataCacheMap");
return NULL;
}
dataCacheMap *getSecondaryCache(int i)
{
if (i==0)
return this;
return _secondaryCaches[i-1];
}
void addSecondaryCache(dataCacheMap *s)
{
_secondaryCaches.push_back(s);
}
dataCacheDouble *get(const function *f, dataCacheDouble *caller=0, bool createIfNotPresent = true);
virtual void setElement(MElement *element)
{
_element=element;
for(std::vector<dataCacheDouble*>::iterator it=_toInvalidateOnElement.begin();
it!= _toInvalidateOnElement.end(); it++) {
(*it)->_valid=false;
}
for(std::list<dataCacheMap*>::iterator it=_children.begin(); it!=_children.end(); it++) {
(*it)->setElement(element);
}
}
inline MElement *getElement() { return _element; }
virtual dataCacheMap *newChild()
{
dataCacheMap *m = new dataCacheMap();
m->_parent = this;
_children.push_back(m);
m->_nbEvaluationPoints = 0;
return m;
}
inline void setFunctionCoordinates(const function *functionCoordinates) {
_functionCoordinates = functionCoordinates;
}
inline void setSolutionFunction(const function *functionSolution, const function *functionSolutionGradient) {
_functionSolution = functionSolution;
_functionSolutionGradient = functionSolutionGradient;
for(std::vector<dataCacheMap*>::iterator it = _secondaryCaches.begin(); it != _secondaryCaches.end(); it++) {
(*it)->setSolutionFunction(functionSolution, functionSolutionGradient);
}
}
inline void setReferenceSolutionFunction(const function *functionSolution, const function *functionSolutionGradient) {
_containerSolution = functionSolution;
_containerSolutionGradient = functionSolutionGradient;
for(std::list<dataCacheMap*>::iterator it = _children.begin(); it != _children.end(); it++) {
(*it)->setReferenceSolutionFunction(functionSolution, functionSolutionGradient);
}
for(std::vector<dataCacheMap*>::iterator it = _secondaryCaches.begin(); it != _secondaryCaches.end(); it++) {
(*it)->setReferenceSolutionFunction(functionSolution, functionSolutionGradient);
}
}
void setNbEvaluationPoints(int nbEvaluationPoints);
inline int getNbEvaluationPoints() { return _nbEvaluationPoints; }
};
// functionConstant (constant values copied over each line)
class functionConstant : public function {
public:
fullMatrix<double> _source;
void call(dataCacheMap *m, fullMatrix<double> &val);
functionConstant(std::vector<double> source);
functionConstant(double source);
void set(double val);
};
class functionC : public function {
std::vector<fullMatrix<double> > args;
void (*callback)(void);
public:
static void buildLibrary(std::string code, std::string filename) ;
static void buildLibraryFromFile(const std::string cfilename, const std::string libfilename);
void call (dataCacheMap *m, fullMatrix<double> &val) ;
functionC (std::string file, std::string symbol, int nbCol,
std::vector<const function *> dependencies);
};
function *functionLevelsetNew (const function *f0, const double valMin, const double valPlus);
function *functionLevelsetSmoothNew (const function *f0, const double valMin, const double valPlus, const double E);
function *functionSumNew (const function *f0, const function *f1);
function *functionMinusNew (const function *f0, const function *f1);
function *functionProdNew (const function *f0, const function *f1);
function *functionQuotientNew (const function *f0, const function *f1);
function *functionScaleNew (const function *f0, const double s);
function *functionAbsNew (const function *f0);
function *functionExtractCompNew (const function *f0, const int iComp);
function *functionCatCompNew(std::vector<const function *> fArray);
function *functionMeanP1New(const function *f, const function *df);
#endif
#include "functionDerivator.h"
#include "function.h"
const fullMatrix<double> &functionDerivator::compute() {
const int nbColX = _x.get().size2();
const double o_eps = 1. / _epsilon;
_fRef = _f.get();
_dfdx.resize(_fRef.size1(),_fRef.size2() * nbColX, false);
for (int j = 0; j < nbColX; ++j) {
fullMatrix<double> &x =_x.set();
for (int i=0;i<_fRef.size1();i++)
x(i,j) += _epsilon;
const fullMatrix<double> &f = _f.get();
for (int k=0; k<_fRef.size2(); k++)
for (int i=0; i<_fRef.size1(); i++)
_dfdx(i, k * nbColX + j) = (f(i,k)-_fRef(i,k)) * o_eps;
for (int i=0;i<_fRef.size1();i++)
x(i,j) -= _epsilon;
}
_x.set();
return _dfdx;
};
#ifndef _FUNCTION_DERIVATOR_H_
#define _FUNCTION_DERIVATOR_H_
#include "fullMatrix.h"
#include "function.h"
class functionDerivator {
dataCacheDouble &_f,&_x;
fullMatrix<double> _fRef, _xRef, _dfdx;
double _epsilon;
public:
functionDerivator (dataCacheDouble &f, dataCacheDouble &x, double epsilon):
_f(f),
_x(x),
_epsilon(epsilon)
{}
const fullMatrix<double> &compute();
inline double get(int iQP, int iF, int iX)
{
return _dfdx(iQP, iF*_xRef.size2()+iX);
}
};
#endif
#ifndef _FUNCTION_NUMPY_H_
#define _FUNCTION_NUMPY_H_
#include "GmshConfig.h"
#ifdef HAVE_NUMPY
#include "Python.h"
#include "numpy/arrayobject.h"
class functionNumpy : public function {
PyObject *_pycallback;
std::vector<fullMatrix<double> > args;
public:
static PyObject *pyArrayFromFullMatrix(fullMatrix<double> &m) {
long int n[2] = {m.size1(), m.size2()};
return PyArray_New(&PyArray_Type, 2, n, NPY_DOUBLE, NULL, &m(0, 0), 0, NPY_FARRAY, NULL);
}
void call (dataCacheMap *m, fullMatrix<double> &res)
{
PyObject *swigR;
PyObject *pyargs;
std::vector<PyObject*> swigA(args.size());
swigR = pyArrayFromFullMatrix(res);
for (int i = 0; i < args.size(); i++) {
swigA[i] = pyArrayFromFullMatrix(args[i]);
}
switch(args.size()) {
case 0 : pyargs = Py_BuildValue("(N)", swigR); break;
case 1 : pyargs = Py_BuildValue("(NN)", swigR, swigA[0]); break;
case 2 : pyargs = Py_BuildValue("(NNN)", swigR, swigA[0], swigA[1]); break;
case 3 : pyargs = Py_BuildValue("(NNNN)", swigR, swigA[0], swigA[1], swigA[2]); break;
case 4 : pyargs = Py_BuildValue("(NNNNN)", swigR, swigA[0], swigA[1], swigA[2], swigA[3]); break;
case 5 : pyargs = Py_BuildValue("(NNNNNN)", swigR, swigA[0], swigA[1], swigA[2], swigA[3], swigA[4]); break;
case 6 : pyargs = Py_BuildValue("(NNNNNNN)", swigR, swigA[0], swigA[1], swigA[2], swigA[3], swigA[4], swigA[5]); break;
case 7 : pyargs = Py_BuildValue("(NNNNNNNN)", swigR, swigA[0], swigA[1], swigA[2], swigA[3], swigA[4], swigA[5], swigA[6]); break;
case 8 : pyargs = Py_BuildValue("(NNNNNNNNN)", swigR, swigA[0], swigA[1], swigA[2], swigA[3], swigA[4], swigA[5], swigA[6], swigA[7]); break;
default:Msg::Error("python function not implemented for more than 8 arguments");
}
PyObject *result = PyEval_CallObject(_pycallback, pyargs);
if (result) {
Py_DECREF(result);
}
else {
PyErr_Print();
Msg::Fatal("An error occurs in the python function.");
}
/* for (int i = 0; i < args.size(); i++) {
Py_DECREF(swigA[i]);
}
Py_DECREF(swigR);*/
}
functionNumpy (int nbCol, PyObject *callback, std::vector<const function*> dependencies)
: function(nbCol), _pycallback(callback)
{
args.resize(dependencies.size());
for (int i = 0; i < dependencies.size(); i++) {
setArgument(args[i], dependencies[i]);
}
static bool _arrayImported = false;
if (! _arrayImported){
_import_array();
_arrayImported = true;
}
}
functionNumpy (int nbCol, PyObject *callback, std::vector<std::pair<const function*, int> > dependencies)
: function(nbCol), _pycallback(callback)
{
args.resize(dependencies.size());
for (int i = 0; i < dependencies.size(); i++) {
setArgument(args[i], dependencies[i].first, dependencies[i].second);
}
printf("import array !!!\n");
static bool _arrayImported = false;
if (! _arrayImported) {
printf("import array !!!\n");
_import_array();
_arrayImported = true;
}
}
};
#endif
#endif
#ifndef _FUNCTION_PYTHON_H_
#define _FUNCTION_PYTHON_H_
#include "Python.h"
class functionPythonReturnMatrix : public fullMatrix<double> {
dataCacheMap *_cacheMap;
public:
void setDataCacheMap(dataCacheMap *cacheMap) {_cacheMap = cacheMap;}
dataCacheMap *getDataCacheMap() const {return _cacheMap;}
};
class functionPython : public function {
PyObject *_pycallback, *_pyargs;
PyObject *_swigR, *_swigCacheMap;
std::vector<PyObject*> _swigA;
std::string _luaFunctionName;
std::vector<fullMatrix<double> > args;
functionPythonReturnMatrix R;
void _init() {
_swigA.resize(args.size());
_swigR = SWIG_NewPointerObj((void*)&R,SWIGTYPE_p_functionPythonReturnMatrix, 0);
for (int i = 0; i < args.size(); i++) {
_swigA[i] = SWIG_NewPointerObj((void*)&args[i],SWIGTYPE_p_fullMatrixT_double_t, 0);
}
switch(args.size()) {
case 0 : _pyargs = Py_BuildValue("(O)", _swigR); break;
case 1 : _pyargs = Py_BuildValue("(OO)", _swigR, _swigA[0]); break;
case 2 : _pyargs = Py_BuildValue("(OOO)", _swigR, _swigA[0], _swigA[1]); break;
case 3 : _pyargs = Py_BuildValue("(OOOO)", _swigR, _swigA[0], _swigA[1], _swigA[2]); break;
case 4 : _pyargs = Py_BuildValue("(OOOOO)", _swigR, _swigA[0], _swigA[1], _swigA[2], _swigA[3]); break;
case 5 : _pyargs = Py_BuildValue("(OOOOOO)", _swigR, _swigA[0], _swigA[1], _swigA[2], _swigA[3], _swigA[4]); break;
case 6 : _pyargs = Py_BuildValue("(OOOOOOO)", _swigR, _swigA[0], _swigA[1], _swigA[2], _swigA[3], _swigA[4], _swigA[5]); break;
case 7 : _pyargs = Py_BuildValue("(OOOOOOOO)", _swigR, _swigA[0], _swigA[1], _swigA[2], _swigA[3], _swigA[4], _swigA[5], _swigA[6]); break;
case 8 : _pyargs = Py_BuildValue("(OOOOOOOOO)", _swigR, _swigA[0], _swigA[1], _swigA[2], _swigA[3], _swigA[4], _swigA[5], _swigA[6], _swigA[7]); break;
default:Msg::Error("python function not implemented for more than 8 arguments");
}
}
public:
void call (dataCacheMap *m, fullMatrix<double> &res)
{
R.setAsProxy(res);
R.setDataCacheMap(m);
PyObject *result = PyEval_CallObject(_pycallback, _pyargs);
if (result) {
Py_DECREF(result);
}
else {
PyErr_Print();
Msg::Fatal("An error occurs in the python function.");
}
}
functionPython (int nbCol, PyObject *callback, std::vector<const function*> dependencies)
: function(nbCol), _pycallback(callback)
{
args.resize(dependencies.size());
for (int i = 0; i < dependencies.size(); i++) {
setArgument(args[i], dependencies[i]);
}
_init();
}
functionPython (int nbCol, PyObject *callback, std::vector<std::pair<const function*, int> > dependencies)
: function(nbCol), _pycallback(callback)
{
args.resize(dependencies.size());
for (int i = 0; i < dependencies.size(); i++) {
setArgument(args[i], dependencies[i].first, dependencies[i].second);
}
_init();
}
};
#endif
......@@ -58,16 +58,12 @@ ENDMACRO(SWIG_GET_WRAPPER_DEPENDENCIES)
option(ENABLE_PYTHON_LIB_API "Export all C header files needed to build the python library" OFF)
if(ENABLE_PYTHON_LIB_API)
set(GMSH_API ${GMSH_API} Geo/Curvature.h Mesh/Generator.h Mesh/highOrderTools.h Mesh/meshGFaceLloyd.h Numeric/DivideAndConquer.h Post/PViewFactory.h Solver/function.h Solver/functionDerivator.h Solver/functionPython.h Solver/functionNumpy.h Solver/linearSystemPETSc.h Fltk/FlGui.h Solver/functionSpace.h Solver/STensor43.h Solver/sparsityPattern.h Solver/SElement.h Solver/groupOfElements.h PARENT_SCOPE)
set(GMSH_API ${GMSH_API} Geo/Curvature.h Mesh/Generator.h Mesh/highOrderTools.h Mesh/meshGFaceLloyd.h Numeric/DivideAndConquer.h Post/PViewFactory.h Solver/linearSystemPETSc.h Fltk/FlGui.h Solver/functionSpace.h Solver/STensor43.h Solver/sparsityPattern.h Solver/SElement.h Solver/groupOfElements.h PARENT_SCOPE)
endif(ENABLE_PYTHON_LIB_API)
if(HAVE_SWIG)
include(${SWIG_USE_FILE})
include_directories(${PYTHON_INCLUDE_DIR})
if(HAVE_NUMPY)
include_directories(${NUMPY_INC})
add_definitions(-DHAVE_NUMPY)
endif(HAVE_NUMPY)
foreach(module ${SWIG_MODULES})
set_source_files_properties(${module}.i PROPERTIES CPLUSPLUS ON)
......
......@@ -7,32 +7,19 @@
%{
#include "GmshConfig.h"
#include "dofManager.h"
#include "elasticitySolver.h"
#include "function.h"
#include "functionDerivator.h"
#include "functionPython.h"
#include "functionNumpy.h"
#include "linearSystem.h"
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
#include "linearSystemPETSc.h"
%}
namespace std {
%template(VectorFunctionConst) vector<const function*, std::allocator<const function*> >;
}
%include "GmshConfig.h"
%include "dofManager.h"
%template(dofManagerDouble) dofManager<double>;
%include "elasticitySolver.h"
%include "function.h"
%include "functionDerivator.h"
%include "functionPython.h"
%include "functionNumpy.h"
%include "linearSystem.h"
%template(linearSystemDouble) linearSystem<double>;
%template(linearSystemFullMatrixDouble) linearSystem<fullMatrix<double> >;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment