dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #16750
Re: [Fwd: [Branch ~dolfin-core/dolfin/main] Rev 4352: Work on Expression in PyDOLFIN]
Anders Logg wrote:
> Where does this show up?
>
Not sure what's up. Demo seems to work now.
Garth
> --
> 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
>
Follow ups
References
