This change appears to have broken compiled expressions. Garth
-------- Original Message -------- Subject: [Branch ~dolfin-core/dolfin/main] Rev 4352: Work on Expression in PyDOLFIN Date: Mon, 30 Nov 2009 08:11:09 -0000 From: [email protected] Reply-To: [email protected] To: Garth Wells <[email protected]> ------------------------------------------------------------ revno: 4352 committer: Johan Hake <[email protected]> branch nick: dolfin timestamp: Sun 2009-11-29 23:59:29 -0800 message: Work on Expression in PyDOLFIN - Simplify expression.py * I have split the logic in create_expression_class into on for python Expressions and one for compiled ones * Removed the possibilty of creating compiled classes that is not instantiated. (never used?) - Fixes in some typemaps * director typemaps for std::vector<double> &values now works * I have turned off director for eval(values,x). Could not get it to work It now works through eval(values,data) in python * Got argout typemap for std::vector<double> & values to work. Bottom line we are up to par in features wrt Expression in PyDOLFIN - Plotting works - Subclassing in Python works for both eval and eval_data - CompiledExpressions work: * MultiLine expressions need to set its value shape if it is different from 0 TODO: - Cleanup in the SWIG interface files and docstrings in expression.py - get_row, and set_row in GenericMatrix does not work as intended. Needed to sacrify this one... modified: demo/pde/poisson/python/demo.py dolfin/function/Expression.cpp dolfin/function/Expression.h dolfin/swig/function_post.i dolfin/swig/function_pre.i dolfin/swig/std_vector_typemaps.i site-packages/dolfin/expression.py -- lp:dolfin https://code.launchpad.net/~dolfin-core/dolfin/main Your team DOLFIN Core Team is subscribed to branch lp:dolfin. To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/main/+edit-subscription.
=== modified file 'demo/pde/poisson/python/demo.py' --- demo/pde/poisson/python/demo.py 2009-11-29 21:17:18 +0000 +++ demo/pde/poisson/python/demo.py 2009-11-30 07:59:29 +0000 @@ -13,7 +13,7 @@ """ __author__ = "Anders Logg ([email protected])" -__date__ = "2007-08-16 -- 2009-11-24" +__date__ = "2007-08-16 -- 2009-11-29" __copyright__ = "Copyright (C) 2007-2009 Anders Logg" __license__ = "GNU LGPL Version 2.1" @@ -25,7 +25,7 @@ class Source(Expression): def eval(self, values, x): - values[0] = 4.0*DOLFIN_PI*DOLFIN_PI*DOLFIN_PI*DOLFIN_PI*sin(DOLFIN_PI*x[0])*sin(DOLFIN_PI*x[1]) + values[0] = 10*exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02) # Define Dirichlet boundary (x = 0 or x = 1) def boundary(x): @@ -38,8 +38,8 @@ # Define variational problem v = TestFunction(V) u = TrialFunction(V) -#f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)") -f = Source(V) +f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)") +#f = Source() g = Expression("sin(5*x[0])") a = inner(grad(v), grad(u))*dx L = v*f*dx - v*g*ds @@ -53,4 +53,4 @@ file << u # Plot solution -#plot(u, interactive=True) +plot(u, interactive=True) === modified file 'dolfin/function/Expression.cpp' --- dolfin/function/Expression.cpp 2009-11-29 12:18:27 +0000 +++ dolfin/function/Expression.cpp 2009-11-30 07:59:29 +0000 @@ -2,7 +2,7 @@ // Licensed under the GNU LGPL Version 2.1. // // First added: 2009-09-28 -// Last changed: 2009-10-11 +// Last changed: 2009-11-29 // // Modified by Johan Hake, 2009. @@ -114,8 +114,8 @@ //----------------------------------------------------------------------------- void Expression::eval(std::vector<double>& values, const std::vector<double>& x) const { - eval(&values[0], x); - //error("Missing eval() for Expression (must be overloaded)."); + //eval(&values[0], x); + error("Missing eval() for Expression (must be overloaded)."); } //----------------------------------------------------------------------------- === modified file 'dolfin/function/Expression.h' --- dolfin/function/Expression.h 2009-11-29 21:17:18 +0000 +++ dolfin/function/Expression.h 2009-11-30 07:59:29 +0000 @@ -2,7 +2,7 @@ // Licensed under the GNU LGPL Version 2.1. // // First added: 2009-09-28 -// Last changed: 2009-10-11 +// Last changed: 2009-11-29 #ifndef __EXPRESSION_H #define __EXPRESSION_H @@ -75,13 +75,6 @@ /// Evaluate expression, must be overloaded by user (simple version) virtual void eval(std::vector<double>& values, const std::vector<double>& x) const; - // Tempory fix while figuring out SWIG - virtual void eval(double* values, const std::vector<double>& x) const - { - cout << "In eval " << endl; - error("Missing eval() for Expression (must be overloaded)."); - } - protected: // Value shape === modified file 'dolfin/swig/function_post.i' --- dolfin/swig/function_post.i 2009-11-16 15:13:08 +0000 +++ dolfin/swig/function_post.i 2009-11-30 07:59:29 +0000 @@ -25,16 +25,16 @@ //----------------------------------------------------------------------------- // Extend the Data class with an accessor function for the x coordinates //----------------------------------------------------------------------------- -//%extend dolfin::Data { -// PyObject* x_() { -// npy_intp adims[1]; -// adims[0] = self->cell().dim(); -// PyArrayObject* array = reinterpret_cast<PyArrayObject*>(PyArray_SimpleNewFromData(1, adims, NPY_DOUBLE, (char *)(self->x))); -// if ( array == NULL ) return NULL; -// PyArray_INCREF(array); -// return reinterpret_cast<PyObject*>(array); -// } -//} +%extend dolfin::Data { + PyObject* x_() { + npy_intp adims[1]; + adims[0] = self->cell().dim(); + PyArrayObject* array = reinterpret_cast<PyArrayObject*>(PyArray_SimpleNewFromData(1, adims, NPY_DOUBLE, reinterpret_cast<char *>( &(const_cast<std::vector<double>& >(self->x))[0] ))); + if ( array == NULL ) return NULL; + PyArray_INCREF(array); + return reinterpret_cast<PyObject*>(array); + } +} //----------------------------------------------------------------------------- // Clear director typemaps @@ -42,3 +42,4 @@ %clear const double* x; %clear double* y; %clear const std::vector<double>& x; +%clear std::vector<double>& values; === modified file 'dolfin/swig/function_pre.i' --- dolfin/swig/function_pre.i 2009-11-29 21:17:18 +0000 +++ dolfin/swig/function_pre.i 2009-11-30 07:59:29 +0000 @@ -9,7 +9,7 @@ // Modified by Kent-Andre Mardal, 2009 // // First added: 2007-08-16 -// Last changed: 2009-10-07 +// Last changed: 2009-11-29 // =========================================================================== // SWIG directives for the DOLFIN function kernel module (pre) @@ -60,8 +60,7 @@ //----------------------------------------------------------------------------- %ignore dolfin::Data::x; %rename (x) dolfin::Data::x_(); -%ignore dolfin::Expression::eval(std::vector<double>& values, const std::vector<double>& x) const; -//%ignore dolfin::Expression::eval(double* values, const std::vector<double>& x) const; +//%ignore dolfin::Expression::eval(std::vector<double>& values, const std::vector<double>& x) const; //----------------------------------------------------------------------------- // Modifying the interface of Constant @@ -118,19 +117,18 @@ %feature("novaluewrapper") std::vector<double>; //----------------------------------------------------------------------------- -// Instantiate a dummy std::vector<dolfin::uint> so value wrapper is not used +// Instantiate a dummy std::vector<dolfin::double> so value wrapper is not used //----------------------------------------------------------------------------- %template () std::vector<double>; //----------------------------------------------------------------------------- -// Typemap for std::vector<dolfin::uint> values +// Typemap for std::vector<dolfin::double> values (used in Constant constructor) //----------------------------------------------------------------------------- -//%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY) std::vector<double> values -//{ -// $1 = PyList_Check($input) ? 1 : 0; -//} +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY) std::vector<double> values +{ + $1 = PyList_Check($input) ? 1 : 0; +} -/* %typemap (in) std::vector<double> values { if (PyList_Check($input)) @@ -153,7 +151,7 @@ SWIG_exception(SWIG_TypeError, "expected list of floats"); } } -*/ + //----------------------------------------------------------------------------- // Add director classes //----------------------------------------------------------------------------- @@ -163,28 +161,14 @@ %feature("nodirector") dolfin::Expression::gather; %feature("nodirector") dolfin::Expression::value_dimension; %feature("nodirector") dolfin::Expression::value_rank; +%feature("nodirector") dolfin::Expression::eval(std::vector<double>& values, const std::vector<double>& x) const; //----------------------------------------------------------------------------- // Director typemap for values in Expression //----------------------------------------------------------------------------- -%typemap(directorin) double* values -{ - { - // Compute size of value (number of entries in tensor value) - dolfin::uint size = 1; - for (dolfin::uint i = 0; i < this->value_rank(); i++) - size *= this->value_dimension(i); - - npy_intp dims[1] = {size}; - $input = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, - reinterpret_cast<char*>($1_name)); - } -} - %typemap(directorin) std::vector<double>& values { { - std::cout << "In typemap " << std::endl; // Compute size of x npy_intp dims[1] = {$1_name.size()}; $input = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, @@ -195,18 +179,11 @@ //----------------------------------------------------------------------------- // Director typemap for coordinates in Expression //----------------------------------------------------------------------------- -//%typemap(directorin) const double* x { -// { -// // Compute size of x -// npy_intp dims[1] = {this->geometric_dimension()}; -// $input = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, reinterpret_cast<char *>(const_cast<double*>($1_name))); -// } -//} - // FIXME: Is there a better way to map a std::vector to a numpy array? %typemap(directorin) const std::vector<double>& x { { + std::cout << "In typemap & x" << std::endl; // Compute size of x npy_intp dims[1] = {$1_name.size()}; $input = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, @@ -214,3 +191,43 @@ } } +////----------------------------------------------------------------------------- +//// In typemaps for std::vector<double> _array +////----------------------------------------------------------------------------- +//%typemap(in) std::vector<double> & _array (std::vector<double> vec_tmp) +//{ +// // Check arguments +// if (!PyArray_Check($input)){ +// PyErr_SetString(PyExc_TypeError, "numpy array of 'double' expected. Make sure that the numpy array use dtype='d'."); +// return NULL; +// } +// +// PyArrayObject *xa = reinterpret_cast<PyArrayObject*>(input); +// if (!PyArray_TYPE(xa) == NPY_DOUBLE ){ +// PyErr_SetString(PyExc_TypeError, "numpy array of 'double' expected. Make sure that the numpy array use dtype='d'."); +// return NULL; +// } +// +// // Get size and reserve the tmp vector +// npy_int size = PyArray_Size(xa); +// vec_tmp.reserve(size); +// +// // Get the data +// double * data = static_cast<double*>(PyArray_DATA(xa)); +// for (int i=0, i<size; i++) +// vec_tmp[i] = data[i]; +// +// // Provide the out argument +// $1 = &vec_tmp; +//} +// +//%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY) std::vector<double> & values +//{ +// $1 = PyArray_Check($input) ? 1 : 0; +//} +// +//%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY) const std::vector<double> & x +//{ +// $1 = PyArray_Check($input) ? 1 : 0; +//} +// === modified file 'dolfin/swig/std_vector_typemaps.i' --- dolfin/swig/std_vector_typemaps.i 2009-11-29 21:17:18 +0000 +++ dolfin/swig/std_vector_typemaps.i 2009-11-30 07:59:29 +0000 @@ -3,7 +3,7 @@ // Licensed under the GNU LGPL Version 2.1. // // First added: 2009-08-31 -// Last changed: 2009-09-29 +// Last changed: 2009-11-29 //============================================================================= // In this file we declare what types that should be able to be passed using a @@ -132,6 +132,7 @@ //----------------------------------------------------------------------------- %typemap (in,numinputs=0) std::vector<TYPE>& ARG_NAME (std::vector<TYPE> vec_temp) { + // ARGOUT ARG_NAME $1 = &vec_temp; } @@ -184,8 +185,11 @@ IN_TYPEMAPS_STD_VECTOR_OF_POINTERS(FunctionSpace) ARGOUT_TYPEMAP_STD_VECTOR_OF_PRIMITIVES(dolfin::uint, INT32, cells, NPY_INT) -ARGOUT_TYPEMAP_STD_VECTOR_OF_PRIMITIVES(dolfin::uint, INT32, columns, NPY_INT) -ARGOUT_TYPEMAP_STD_VECTOR_OF_PRIMITIVES(double, DOUBLE, values, NPY_DOUBLE) +// FIXME: We cannot have argout typemaps for columns and values +// FIXME: They work for get_row, but they interfere with eval interface. +// FIXME: They should also _not_ kick in for const std::vector<TYPE>, but they do +//ARGOUT_TYPEMAP_STD_VECTOR_OF_PRIMITIVES(dolfin::uint, INT32, columns, NPY_INT) +//ARGOUT_TYPEMAP_STD_VECTOR_OF_PRIMITIVES(double, DOUBLE, values, NPY_DOUBLE) //----------------------------------------------------------------------------- // NumPy to const std::vector<double>& typemap. @@ -194,28 +198,6 @@ $1 = PyArray_Check($input) ? 1 : 0; } -%typemap(in) std::vector<double>& x (std::vector<double> temp) -{ - { - if (PyArray_Check($input)) - { - PyArrayObject *xa = reinterpret_cast<PyArrayObject*>($input); - if ( PyArray_TYPE(xa) == NPY_DOUBLE ) - { - const unsigned int size = PyArray_DIM(xa, 0); - temp.resize(size); - double* array = static_cast<double*>(PyArray_DATA(xa)); - std::copy(array, array + size, temp.begin()); - $1 = &temp; - } - else - SWIG_exception(SWIG_TypeError, "NumPy array expected"); - } - else - SWIG_exception(SWIG_TypeError, "NumPy array expected"); - } -} - %typemap(in) const std::vector<double>& x (std::vector<double> temp) { { @@ -237,19 +219,64 @@ SWIG_exception(SWIG_TypeError, "NumPy array expected"); } } -//----------------------------------------------------------------------------- -// const std::vector<double>& to NumPy typemap. -//----------------------------------------------------------------------------- -%typemap(argout) const std::vector<double>& x -{ - // Do nothing -} -//----------------------------------------------------------------------------- -// std::vector<double>& to NumPy typemap. -//----------------------------------------------------------------------------- -%typemap(argout) std::vector<double>& x -{ - SWIG_exception(SWIG_TypeError, "std::vector<double> to NumPy (argout) not implemented"); -} -//----------------------------------------------------------------------------- - + +//----------------------------------------------------------------------------- +// std::vector<double>& to NumPy typemap. +//----------------------------------------------------------------------------- +%typemap(in) std::vector<double>& values (std::vector<double> temp, PyArrayObject *out_array = 0) +{ + if (!PyArray_Check($input)) + SWIG_exception(SWIG_TypeError, "numpy array of 'double' expected. Make sure that the numpy array use dtype='d'."); + out_array = reinterpret_cast<PyArrayObject*>($input); + if ( !PyArray_TYPE(out_array) == NPY_DOUBLE ) + SWIG_exception(SWIG_TypeError, "numpy array of 'double' expected. Make sure that the numpy array use dtype='d'."); + + // Use the size of the incomming array to pass a correct sized vector to the function + const unsigned int size = PyArray_DIM(out_array, 0); + temp.resize(size); + $1 = &temp; +} + +%typemap(argout) std::vector<double>& values +{ + // Get the pointer to the data in the passed NumPy array + double* array = static_cast<double*>(PyArray_DATA(out_array$argnum)); + + // Copy the content from the temp array to the NumPy array + std::copy(temp$argnum.begin(), temp$argnum.end(), array); + + // Tell the function to return 'None', which means not return value + $result = Py_None; + +} + +//----------------------------------------------------------------------------- +// std::vector<double>& to NumPy typemap. +//----------------------------------------------------------------------------- +%typemap(in) std::vector<double>& vertex_values (std::vector<double> temp, PyArrayObject *out_array = 0, dolfin::uint init_size = 0) +{ + if (!PyArray_Check($input)) + SWIG_exception(SWIG_TypeError, "numpy array of 'double' expected. Make sure that the numpy array use dtype='d'."); + out_array = reinterpret_cast<PyArrayObject*>($input); + if ( !PyArray_TYPE(out_array) == NPY_DOUBLE ) + SWIG_exception(SWIG_TypeError, "numpy array of 'double' expected. Make sure that the numpy array use dtype='d'."); + + // Use the size of the incomming array to pass a correct sized vector to the function + const unsigned int size = PyArray_DIM(out_array, 0); + init_size = size; + temp.resize(size); + $1 = &temp; +} + +%typemap(argout) std::vector<double>& vertex_values +{ + // Get the pointer to the data in the passed NumPy array + double* array = static_cast<double*>(PyArray_DATA(out_array$argnum)); + + // Copy the content from the temp array to the NumPy array + std::copy(temp$argnum.begin(), temp$argnum.end(), array); + + // Tell the function to return 'None', which means not return value + $result = Py_None; + +} === modified file 'site-packages/dolfin/expression.py' --- site-packages/dolfin/expression.py 2009-11-28 21:18:15 +0000 +++ site-packages/dolfin/expression.py 2009-11-30 07:59:29 +0000 @@ -42,7 +42,6 @@ __license__ = "GNU LGPL Version 2.1" # Modified by Anders Logg, 2008-2009. -# Modified by Garth N. Wells, 2009. __all__ = ["Expression", "Expressions"] @@ -64,96 +63,93 @@ from compile_expressions import compile_expressions from functionspace import * -def create_expression_class(name, - cpp_base, - user_bases = None, - user_dict = None, - dim_needs_to_be_passed = False): +def create_compiled_expression_class(cpp_base): + # Check the cpp_base + assert(isinstance(cpp_base, (types.ClassType, type))) + + def __init__(self, cpparg, defaults=None, element=None, degree=None): + """ Initialize the Expression """ + # Initialize the cpp base class first and extract value_shape + cpp_base.__init__(self) + value_shape = tuple(self.value_dimension(i) \ + for i in range(self.value_rank())) + + # Store the dim + if len(value_shape) == 0: + self._dim = 0 + elif len(value_shape) == 1: + self._dim = value_shape[0] + else: + self._dim = value_shape + + # Select an appropriate element if not specified + if element is None: + element = _auto_select_element_from_shape(value_shape, degree) + else: + # Check that we have an element + if not isinstance(element, ufl.FiniteElementBase): + raise TypeError, "The 'element' argument must be a UFL finite element." + + # Check same value shape of compiled expression and passed element + if not element.value_shape() == value_shape: + raise ValueError, "The value shape of the passed 'element', is not equal to the value shape of the compiled expression." + + # Initialize UFL base class + self._ufl_element = element + ufl.Function.__init__(self, self._ufl_element) + + # Create and return the class + return type("CompiledExpression", (Expression, ufl.Function, cpp_base), {"__init__":__init__}) + +def create_python_derived_expression_class(name, user_bases, user_dict): """Return Expression class This function is used to create all the dynamically created Expression classes. It takes a name, and a compiled cpp.Expression and returns - a dolfin.Expression class. In addition to cpp.Expression and dolfin.Expression - it also inherits from ufl.Function. + a class that inherits the compiled class together with dolfin.Expression + and ufl.Function. @param name: The name of the class - @param cpp_base: - The cpp.Expression base class which the created - Expression class shall inherit. @param user_bases: - Optional user defined bases + User defined bases @param user_dict: - Optional dict with user specified function or attributes - @param dim_needs_to_be_passed: - Optional if a simple expression using cpparg is created with no - information about geometric dimensions + Dict with user specified function or attributes """ - # Check the name + # Check args assert(isinstance(name, str)) - assert(name != "Expression"), "Cannot create a sub class of Expression with the same name as Expression" - - assert(isinstance(cpp_base, (types.ClassType, type))) + assert(isinstance(user_bases, list)) + assert(isinstance(user_dict, dict)) # Define the bases - user_bases = user_bases or [] assert(all([isinstance(t, (types.ClassType, type)) for t in user_bases])) - bases = tuple([Expression, ufl.Function, cpp_base] + user_bases) - - # Define the dictionary of the class - dict_ = user_dict or {} - assert isinstance(dict_, dict) - + bases = tuple([Expression, ufl.Function, cpp.Expression] + user_bases) + # If a user init is not provided create a dummy one - if "__init__" not in dict_: + if "__init__" not in user_dict: def user_init(self, *arg, **kwargs): pass else: - user_init = dict_.pop("__init__") - - # If a user init is not provided create a dummy one - if "dim" not in dict_: - def user_dim(self): return 1 - else: - user_dim = dict_.pop("dim") + user_init = user_dict.pop("__init__") def __init__(self, *args, **kwargs): - # This is called if no user defined init function is provided # Get element and degree - element = kwargs.pop("element", None) - degree = kwargs.pop("degree", None) + element = kwargs.get("element", None) + degree = kwargs.get("degree", None) # Select an appropriate element if not specified - if element is None and len(args) > 0: - element = _auto_select_element_from_cpparg(args[0], degree) - elif element is None: - element = _auto_select_element_from_dim(user_dim(self), degree) - - # Check that we have have an element - if not isinstance(element, ufl.FiniteElementBase): + if element is None: + element = _auto_select_element_from_dim(self.dim(), degree) + elif not isinstance(element, ufl.FiniteElementBase): raise TypeError, "The 'element' argument must be a UFL finite element." - + # Initialize UFL base class self._ufl_element = element ufl.Function.__init__(self, self._ufl_element) # Initialize cpp_base class - - # First check if we are instantiating a user-defined Python class - if "eval" in dict_ or "eval_data" in dict_: - assert cpp_base == cpp.Expression - cpp_base.__init__(self, list(self._ufl_element.value_shape())) - else: - cpp_base.__init__(self) - - # Check that the value_shape of the ufl.FiniteElement corresponds with the - # created cpp.Expression - shape = self._ufl_element.value_shape() - if not (self.value_rank() == len(shape) and - all(dim == self.value_dimension(i) for i, dim in enumerate(shape))): - exp_shape = tuple(self.value_dimension(i) for i in xrange(self.value_rank())) - raise ValueError, "value_shape of passed element does not match value_shape of the Expression: %s != %s"%(str(shape), str(exp_shape)) + cpp.Expression.__init__(self, list(self._ufl_element.value_shape())) # Calling the user defined_init user_init(self, *args, **kwargs) @@ -165,6 +161,7 @@ __init__.__doc__ = """ Initialize the Expression""" # NOTE: Do not prevent the user to overload attributes "reserved" by PyDOLFIN + # NOTE: Why not? ## Collect reserved attributes from both cpp.Function and ufl.Function #reserved_attr = dir(ufl.Function) @@ -180,22 +177,22 @@ # if attr in dict_: # raise TypeError, "The Function attribute '%s' is reserved by PyDOLFIN."%attr - # Fill the dict_ with constructed function - dict_["__init__"] = __init__ + # Add __init__ to the user_dict + user_dict["__init__"] = __init__ # Create the class and return it - return type(name, bases, dict_) + return type(name, bases, user_dict) class ExpressionMetaClass(type): - + """ Meta Class for Expression""" def __new__(cls, name, bases, dict_): - """ Return a new Expression class """ + """ Returns a new Expression class """ assert(isinstance(name, str)), "Expecting a 'str'" assert(isinstance(bases, tuple)), "Expecting a 'tuple'" assert(isinstance(dict_, dict)), "Expecting a 'dict'" - # First check if we are creating the Function class + # First check if we are creating the Expression class if name == "Expression": # Assert that the class is _not_ a subclass of Expression, # i.e., a user have tried to: @@ -209,7 +206,9 @@ # this module return type.__new__(cls, name, bases, dict_) - # If subclassing Expression directly (used in specialfunctions.py) + # If creating a fullfledged derived expression class, i.e, inheriting + # dolfin.Expression, ufl.Function and cpp.Expression (or a subclass) + # then just return the new class. if len(bases) >= 3 and bases[0] == Expression and \ bases[1] == ufl.Function and issubclass(bases[2], cpp.Expression): # Return the instantiated class @@ -221,59 +220,28 @@ # remove Expression, to be added later user_bases.remove(Expression) - # Check the cppcode and eval attributes - if 'cpparg' in dict_ and ('eval' in dict_ or 'eval_data' in dict_) : - raise TypeError, "Cannot create class with both 'eval'/'eval_data' and 'cpparg' attributes defined." - - # If the Expression class is a user defined python class, case 4. from docstring - if 'eval' in dict_ or 'eval_data' in dict_: - # Get name of eval function - eval_name = 'eval' if 'eval' in dict_ else 'eval_data' - - user_eval = dict_[eval_name] - - # Check type and number of arguments of user_eval function - if not isinstance(user_eval, types.FunctionType): - raise TypeError, "'%s' attribute must be a 'function'"%eval_name - if not user_eval.func_code.co_argcount == 3: - raise TypeError, "The overloaded '%s' function must use three arguments"%eval_name - - return create_expression_class(name, cpp.Expression, user_bases, dict_) - - # If cpparg is provided, case 5. from docstring - if 'cpparg' in dict_: - - # Check the handed attributes and return an args tuple - cpparg = dict_.pop('cpparg') - defaults = dict_.pop("defaults",None) - - # Check that the user has not provide any other attributes - # than the allowed ones. - if len(dict_) > 1: - dict_.pop('__module__') - raise TypeError, "Not allowed to provide user defined attributes to a sub class of Expression when the compiled function interface is used. Found: %s"%\ - (", ".join(["'%s'"%key for key in dict_.iterkeys()])) - - # Check arguments - _check_cpparg(cpparg) - _check_defaults(defaults) - - # Compile the cppargs - cpp_base = compile_expressions([cpparg], [defaults])[0] - - # Add back the cpparg as an attribute - cpp_base.cpparg = cpparg - - # If defaults where handed add it back too - if defaults is not None: - cpp_base.defaults = defaults - - # Create the Expression class and return it - return create_expression_class(name, cpp_base, user_bases, - dim_needs_to_be_passed = not _is_complex_expression(cpparg)) - - # If we have reached this stage raise error - raise TypeError, "Error in subclassing Expression. For correct usage see 4. and 5. in Expression docstring." + # Check the user has provided either an eval or eval_data method + if not ('eval' in dict_ or 'eval_data' in dict_): + raise TypeError, "expected an overload 'eval' or 'eval_data' method" + + # Get name of eval function + eval_name = 'eval' if 'eval' in dict_ else 'eval_data' + + user_eval = dict_[eval_name] + + # Check type and number of arguments of user_eval function + if not isinstance(user_eval, types.FunctionType): + raise TypeError, "'%s' attribute must be a 'function'"%eval_name + if not user_eval.func_code.co_argcount == 3: + raise TypeError, "The overloaded '%s' function must use three arguments"%eval_name + # A nice hack to get around some SWIG director problems + # In short: eval_data works but not eval... + if eval_name == 'eval': + def eval_data(self, values, data): + user_eval(self, values, data.x()) + dict_['eval_data'] = eval_data + + return create_python_derived_expression_class(name, user_bases, dict_) #--- The user interface --- @@ -324,7 +292,6 @@ return ufl.Function.__call__(self,*args) # Some help variables - #dim = self.geometric_dimension() value_size = ufl.common.product(self.ufl_element().value_shape()) # If values (return argument) is passed, check the type and length @@ -349,19 +316,14 @@ raise TypeError, "expected a scalar or an iterable as coordinate argument" # Check for scalar x if isinstance(x[0], (int, float)): - #if not dim == 1: - # raise TypeError, "expected a coordinate argument of length %d"%dim x = numpy.fromiter(x, 'd') else: x = x[0] - #if len(x) != dim: - # raise TypeError, "expected an iterable of length %d as coordinate argument"%dim if isinstance(x, (list, tuple)): x = numpy.fromiter(x, 'd') # If several x arguments have been provided else: - #if len(x) != dim or not all(isinstance(v,(int,float)) for v in x): if not all(isinstance(v,(int,float)) for v in x): raise TypeError, "expected different number of scalar arguments for the coordinates" x = numpy.fromiter(x,'d') @@ -550,8 +512,8 @@ Optional quadrature degree element. """ - # If the __new__ function is called because we are instantiating a sub - # class of Expression, then use the object's __new__ function instead + # If the __new__ function is called because we are instantiating a python sub + # class of Expression, then just return a new instant of the passed class if cls.__name__ != "Expression": return object.__new__(cls) @@ -562,85 +524,111 @@ # Compile module and get the cpp.Expression class cpp_base = compile_expressions([cpparg], [defaults])[0] - # Store arguments for later use + # Store compile arguments for later use cpp_base.cpparg = cpparg cpp_base.defaults = defaults - cpp_base.element = element - cpp_base.degree = degree - return object.__new__(create_expression_class("CompiledExpression", cpp_base)) - + + return object.__new__(create_compiled_expression_class(cpp_base)) + + # This method is only included so a user can check what arguments + # one should use in IPython using tab completion + def __init__(self, cpparg=None, defaults=None, element=None, degree=None):pass + + # Reuse the docstring from __new__ + __init__.__doc__ = __new__.__doc__ + def ufl_element(self): " Return the ufl FiniteElement." return self._ufl_element def __str__(self): + "x.__str__() <==> print(x)" # FIXME: We might change this using rank and dimension instead return "<Expression on a %s>" % str(self._ufl_element) def __repr__(self): + "x.__repr__() <==> repr(x)" return ufl.Function.__repr__(self) + # Default value of dim + _dim = 0 + + def dim(self): + """ Returns the dimension of the value""" + return self._dim + + __call__ = expression__call__ + def Expressions(*args, **kwargs): """ Batch-processed user-defined JIT-compiled expressions ------------------------------------------------------- By specifying several cppargs one may compile more than one expression - at a time. These may either be instantiated using a single FunctionSpace - common for all expressions, using the optional kwarg 'V', or with - a separate FunctionSpace for each cpparg: + at a time: >>> f0, f1 = Expressions('sin(x[0]) + cos(x[1])', 'exp(x[1])', degree=3) - >>> f0, f1 = Expressions('sin(x[0]) + cos(x[1])', 'exp(x[1])', element=element) - - Here cppcode is a code snippet, which should be a string of C++ + + >>> f0, f1, f2 = Expressions((('A*sin(x[0])', 'B*cos(x[1])') + ('0','1')), {'A':2.0,'B':3.0}, + code, + (('cos(x[0])','sin(x[1])'), + ('sin(x[0])','cos(x[1])')), element=element) + + Here code is a C++ code snippet, which should be a string of C++ code that implements a class that inherits from dolfin::Expression, see user case 3. in Expression docstring Batch-processing of JIT-compiled expressions may significantly speed up JIT-compilation at run-time. - -""" + """ # Get the element/degree from kwarg + if len(kwargs) > 1: + raise TypeError, "Can only define one kwarg and that can only be 'degree' or 'element'." degree = kwargs.pop("degree", None) element = kwargs.pop("element", None) - if len(kwargs) > 1: - raise TypeError, "Can only define one kwarg and that can only be 'degree' or 'element'." # Iterate over the *args and collect input to compile_expressions cppargs = []; defaults = []; i = 0; - while i < len(args): + while i < nargs: + + # Check type of cppargs if not isinstance(args[i],(tuple, list, str)): raise TypeError, "Expected either a 'list', 'tuple' or 'str' for argument %d"%i + cppargs.append(args[i]) - defaults.append(None) i += 1 + # If we have more args and the next is a dict + if i < nargs and isinstance(args[i], dict): + # Append the dict to defaults + _check_defaults(args[i]) + defaults.append(args[i]) + i += 1 + else: + # If not append None + defaults.append(None) + # Compile the cpp.Expressions cpp_bases = compile_expressions(cppargs, defaults) # Instantiate the return arguments return_expressions = [] - for i, cpp_base in enumerate(cpp_bases): - return_expressions.append(create_expression_class("CompiledExpression", cpp_base) - (cppargs[i], element=element, degree=degree)) + + for cpp_base in cpp_bases: + # If we only want the cpp.Expression + return_expressions.append(create_compiled_expression_class(cpp_base)(\ + degree=degree, + element=element)) # Return the instantiated Expressions return tuple(return_expressions) -# Assign doc string -expression__call__.__doc__ - -# Assign __call__ method -Expression.__call__ = types.MethodType(expression__call__, None, Expression) - #--- Utility functions --- def _check_cpparg(cpparg): "Check that cpparg makes sense" - if cpparg is None: return - # Check that we get a string expression or nested expression if not isinstance(cpparg, (str, tuple, list)): raise TypeError, "Please provide a 'str', 'tuple' or 'list' for the 'cpparg' argument." @@ -665,14 +653,6 @@ "Check if cpparg is a complex expression" return isinstance(cpparg, str) and "class" in cpparg and "Expression" in cpparg -def _auto_select_element_from_cpparg(cpparg, degree=None): - "Automatically select an appropriate element from cpparg." - - # Use numpy to get the shape - shape = numpy.shape(cpparg) - - return _auto_select_element_from_shape(shape, degree) - def _auto_select_element_from_dim(dim, degree=None): "Automatically select an appropriate element from dim." @@ -708,3 +688,12 @@ cpp.info("Automatic selection of expression element: " + str(element)) return element + +def _check_name_and_base(name, cpp_base): + # Check the name + assert(isinstance(name, str)) + assert(name != "Expression"), "Cannot create a sub class of Expression with the same name as Expression" + + assert(isinstance(cpp_base, (types.ClassType, type))) + +
_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
