dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #16735
[Fwd: [Branch ~dolfin-core/dolfin/main] Rev 4352: Work on Expression in PyDOLFIN]
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
--
lp:dolfin
https://code.launchpad.net/~dolfin-core/dolfin/main
Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
To unsubscribe from this branch go to
https://code.launchpad.net/~dolfin-core/dolfin/main/+edit-subscription.
=== modified file 'demo/pde/poisson/python/demo.py'
--- demo/pde/poisson/python/demo.py 2009-11-29 21:17:18 +0000
+++ demo/pde/poisson/python/demo.py 2009-11-30 07:59:29 +0000
@@ -13,7 +13,7 @@
"""
__author__ = "Anders Logg (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)))
+
+
Follow ups
