dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #16751
Re: [Fwd: [Branch ~dolfin-core/dolfin/main] Rev 4352: Work on Expression in PyDOLFIN]
On Monday 30 November 2009 10:03:40 Garth N. Wells wrote:
> Anders Logg wrote:
> > Where does this show up?
>
> Not sure what's up. Demo seems to work now.
I thought you where aiming at Expressions with an 's'. That was broken until
it was fixed simultaneously by me and Garth.
Johan
> 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(DOLF
> >>IN_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
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>
References
