dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #16749
Re: [Fwd: [Branch ~dolfin-core/dolfin/main] Rev 4352: Work on Expression in PyDOLFIN]
Where does this show up?
--
Anders
On Mon, Nov 30, 2009 at 09:44:24AM +0000, Garth N. Wells wrote:
> 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: noreply@xxxxxxxxxxxxx
> Reply-To: noreply@xxxxxxxxxxxxx
> To: Garth Wells <gnw20@xxxxxxxxx>
>
> ------------------------------------------------------------
> revno: 4352
> committer: Johan Hake <johan.hake@xxxxxxxxx>
> 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
>
>
> === 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 (logg@xxxxxxxxx)"
> -__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 : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
Attachment:
signature.asc
Description: Digital signature
Follow ups
References
