1
"""This module handles the Expression class in Python.
3
The Expression class needs special handling and is not mapped directly
4
by SWIG from the C++ interface. Instead, a new Expression class is
5
created which inherits both from the DOLFIN C++ Expression class and
6
the ufl Coefficient class.
8
The resulting Expression class may thus act both as a variable in a
9
UFL form expression and as a DOLFIN C++ Expression.
11
This module make heavy use of creation of Expression classes and
12
instantiation of these dynamically at runtime.
14
The whole logic behind this somewhat magic behaviour is handle by:
16
1) function __new__ in the Expression class
17
2) meta class ExpressionMetaClass
18
3) function compile_expressions from the compiledmodule/expression
20
4) functions create_compiled_expression_class and
21
create_python_derived_expression_class
23
The __new__ method in the Expression class take cares of the logic
24
when the class Expression is used to create an instance of Expression,
25
see use cases 1-4 in the docstring of Expression.
27
The meta class ExpressionMetaClass take care of the logic when a user
28
subclasses Expression to create a user-defined Expression, see use
29
cases 3 in the docstring of Expression.
31
The function compile_expression is a JIT compiler. It compiles and
32
returns different kinds of cpp.Expression classes, depending on the
33
arguments. These classes are sent to the
34
create_compiled_expression_class
37
__author__ = "Johan Hake (hake@simula.no)"
38
__date__ = "2008-11-03 -- 2009-12-16"
39
__copyright__ = "Copyright (C) 2008-2009 Johan Hake"
40
__license__ = "GNU LGPL Version 2.1"
42
# Modified by Anders Logg, 2008-2009.
44
__all__ = ["Expression", "Expressions"]
46
# FIXME: Make all error messages uniform according to the following template:
48
# if not isinstance(foo, Foo):
49
# raise TypeError, "Illegal argument for creation of Bar, not a Foo: " + str(foo)
54
# Import UFL and SWIG-generated extension module (DOLFIN C++)
56
import dolfin.cpp as cpp
60
from dolfin.compilemodules.expressions import compile_expressions
62
def create_compiled_expression_class(cpp_base):
64
assert(isinstance(cpp_base, (types.ClassType, type)))
66
def __init__(self, cppcode, defaults=None, element=None, degree=None):
67
"""Initialize the Expression """
68
# Initialize the cpp base class first and extract value_shape
69
cpp_base.__init__(self)
70
value_shape = tuple(self.value_dimension(i) \
71
for i in range(self.value_rank()))
73
# Store the value_shape
74
self._value_shape = value_shape
76
# Select an appropriate element if not specified
78
element = _auto_select_element_from_shape(value_shape, degree)
80
# Check that we have an element
81
if not isinstance(element, ufl.FiniteElementBase):
82
raise TypeError, "The 'element' argument must be a UFL"\
85
# Check same value shape of compiled expression and passed element
86
if not element.value_shape() == value_shape:
87
raise ValueError, "The value shape of the passed 'element',"\
88
" is not equal to the value shape of the compiled "\
91
# Initialize UFL base class
92
self._ufl_element = element
93
ufl.Coefficient.__init__(self, self._ufl_element)
95
# Create and return the class
96
return type("CompiledExpression", (Expression, ufl.Coefficient, cpp_base),\
97
{"__init__":__init__})
99
def create_python_derived_expression_class(name, user_bases, user_dict):
100
"""Return Expression class
102
This function is used to create all the dynamically created Expression
103
classes. It takes a name, and a compiled cpp.Expression and returns
104
a class that inherits the compiled class together with dolfin.Expression
109
The name of the class
113
Dict with user specified function or attributes
117
assert(isinstance(name, str))
118
assert(isinstance(user_bases, list))
119
assert(isinstance(user_dict, dict))
122
assert(all([isinstance(t, (types.ClassType, type)) for t in user_bases]))
123
bases = tuple([Expression, ufl.Coefficient, cpp.Expression] + user_bases)
125
user_init = user_dict.pop("__init__", None)
127
# Check for deprecated dim and rank methods
128
if "dim" in user_dict or "rank" in user_dict:
129
raise DeprecationWarning, "'rank' and 'dim' are depcrecated, overload"\
130
" 'value_shape' instead"
132
def __init__(self, *args, **kwargs):
133
"""Initialize the Expression"""
134
# Get element, cell and degree
135
element = kwargs.get("element", None)
136
degree = kwargs.get("degree", None)
138
# Check if user has passed too many arguments if no
139
# user_init is provided
140
if user_init is None:
141
from operator import add
142
# First count how many valid kwargs is passed
143
num_valid_kwargs = reduce(add, [kwarg is not None \
144
for kwarg in [element, degree]])
145
if len(kwargs) != num_valid_kwargs:
146
raise TypeError, "expected only 'kwargs' from one of "\
147
"'element', or 'degree'"
149
raise TypeError, "expected no 'args'"
151
# Select an appropriate element if not specified
153
element = _auto_select_element_from_shape(self.value_shape(),
155
elif isinstance(element, ufl.FiniteElementBase):
158
raise TypeError, "The 'element' argument must be a UFL finite"\
161
# Initialize UFL base class
162
self._ufl_element = element
163
ufl.Coefficient.__init__(self, element)
165
# Initialize cpp_base class
166
cpp.Expression.__init__(self, list(element.value_shape()))
168
# Calling the user defined_init
169
if user_init is not None:
170
user_init(self, *args, **kwargs)
172
# NOTE: Do not prevent the user to overload attributes "reserved" by PyDOLFIN
175
## Collect reserved attributes from both cpp.Function and ufl.Coefficient
176
#reserved_attr = dir(ufl.Coefficient)
177
#reserved_attr.extend(dir(cpp.Function))
179
## Remove attributes that will be set by python
180
#for attr in ["__module__"]:
181
# while attr in reserved_attr:
182
# reserved_attr.remove(attr)
184
## Check the dict_ for reserved attributes
185
#for attr in reserved_attr:
187
# raise TypeError, "The Function attribute '%s' is reserved by PyDOLFIN."%attr
189
# Add __init__ to the user_dict
190
user_dict["__init__"] = __init__
192
# Create the class and return it
193
return type(name, bases, user_dict)
195
class ExpressionMetaClass(type):
196
"""Meta Class for Expression"""
197
def __new__(mcs, name, bases, dict_):
198
"""Returns a new Expression class"""
200
assert(isinstance(name, str)), "Expecting a 'str'"
201
assert(isinstance(bases, tuple)), "Expecting a 'tuple'"
202
assert(isinstance(dict_, dict)), "Expecting a 'dict'"
204
# First check if we are creating the Expression class
205
if name == "Expression":
206
# Assert that the class is _not_ a subclass of Expression,
207
# i.e., a user have tried to:
209
# class Expression(Expression):
211
if len(bases) > 1 and bases[0] != object:
212
raise TypeError, "Cannot name a subclass of Expression:"\
215
# Return the new class, which just is the original Expression defined in
217
return type.__new__(mcs, name, bases, dict_)
219
# If creating a fullfledged derived expression class, i.e, inheriting
220
# dolfin.Expression, ufl.Coefficient and cpp.Expression (or a subclass)
221
# then just return the new class.
222
if len(bases) >= 3 and bases[0] == Expression and \
223
bases[1] == ufl.Coefficient and issubclass(bases[2], \
225
# Return the instantiated class
226
return type.__new__(mcs, name, bases, dict_)
228
# Handle any user provided base classes
229
user_bases = list(bases)
231
# remove Expression, to be added later
232
user_bases.remove(Expression)
234
# Check the user has provided either an eval or eval_cell method
235
if not ('eval' in dict_ or 'eval_cell' in dict_):
236
raise TypeError, "expected an overload 'eval' or 'eval_cell' method"
238
# Get name of eval function
239
eval_name = 'eval' if 'eval' in dict_ else 'eval_cell'
241
user_eval = dict_[eval_name]
243
# Check type and number of arguments of user_eval function
244
if not isinstance(user_eval, types.FunctionType):
245
raise TypeError, "'%s' attribute must be a 'function'"%eval_name
246
if eval_name == "eval" and not user_eval.func_code.co_argcount == 3:
247
raise TypeError, "The overloaded '%s' function must use "\
248
"three arguments"%eval_name
249
if eval_name == "eval_cell" and \
250
not user_eval.func_code.co_argcount == 4:
251
raise TypeError, "The overloaded '%s' function must "\
252
"use three arguments"%eval_name
254
return create_python_derived_expression_class(name, user_bases, dict_)
256
#--- The user interface ---
258
# Places here so it can be reused in Function
259
def expression__call__(self, *args, **kwargs):
261
Evaluates the Expression
264
1) Using an iterable as x:
265
>>> fs = Expression("sin(x[0])*cos(x[1])*sin(x[3])")
266
>>> x0 = (1.,0.5,0.5)
267
>>> x1 = [1.,0.5,0.5]
268
>>> x2 = numpy.array([1.,0.5,0.5])
273
2) Using multiple scalar args for x, interpreted as a point coordinate
274
>>> v0 = f(1.,0.5,0.5)
276
3) Passing return array
277
>>> fv = Expression(("sin(x[0])*cos(x[1])*sin(x[3])",
279
>>> x0 = numpy.array([1.,0.5,0.5])
280
>>> v0 = numpy.zeros(3)
281
>>> fv(x0, values = v0)
285
A longer values array may be passed. In this way one can fast
286
fill up an array with different evaluations.
288
>>> values = numpy.zeros(9)
289
>>> for i in xrange(0,10,3):
290
... fv(x[i:i+3], values = values[i:i+3])
294
raise TypeError, "expected at least 1 argument"
296
# Test for ufl restriction
297
if len(args) == 1 and args[0] in ('+','-'):
298
return ufl.Coefficient.__call__(self, *args)
300
# Test for ufl mapping
301
if len(args) == 2 and isinstance(args[1], dict) and self in args[1]:
302
return ufl.Coefficient.__call__(self, *args)
304
# Some help variables
305
value_size = ufl.common.product(self.ufl_element().value_shape())
307
# If values (return argument) is passed, check the type and length
308
values = kwargs.get("values", None)
309
if values is not None:
310
if not isinstance(values, numpy.ndarray):
311
raise TypeError, "expected a NumPy array for 'values'"
312
if len(values) != value_size or \
313
not numpy.issubdtype(values.dtype, 'd'):
314
raise TypeError, "expected a double NumPy array of length"\
315
" %d for return values."%value_size
316
values_provided = True
318
values_provided = False
319
values = numpy.zeros(value_size, dtype='d')
321
# Assume all args are x argument
324
# If only one x argument has been provided, unpack it if it's an iterable
326
if hasattr(x[0], '__iter__'):
329
# Convert it to an 1D numpy array
331
x = numpy.fromiter(x, 'd')
332
except (TypeError, ValueError, AssertionError):
333
raise TypeError, "expected scalar arguments for the coordinates"
336
raise TypeError, "coordinate argument too short"
338
# The actual evaluation
341
# If scalar return statement, return scalar value.
342
if value_size == 1 and not values_provided:
347
class Expression(object):
349
This class represents a user-defined expression.
351
Expressions can be used as coefficients in variational forms or
352
interpolated into finite element spaces.
356
C++ argument, see below
358
Optional C++ argument, see below
360
Optional element argument
362
*1. Simple user-defined JIT-compiled expressions*
363
One may alternatively specify a C++ code for evaluation of the
364
Expression as follows:
366
>>> f0 = Expression('sin(x[0]) + cos(x[1])')
367
>>> f1 = Expression(('cos(x[0])', 'sin(x[1])'), element = V.ufl_element())
369
Here, f0 is is scalar and f1 is vector-valued.
371
Tensor expressions of rank 2 (matrices) may also be created:
373
>>> f2 = Expression((('exp(x[0])','sin(x[1])'),
374
('sin(x[0])','tan(x[1])')))
376
In general, a single string expression will be interpreted as a
377
scalar, a tuple of strings as a tensor of rank 1 (a vector) and a
378
tuple of tuples of strings as a tensor of rank 2 (a matrix).
380
The expressions may depend on x[0], x[1], and x[2] which carry
381
information about the coordinates where the expression is
382
evaluated. All math functions defined in <cmath> are available to
385
Expression parameters can be included as follows:
387
>>> f = Expression('A*sin(x[0]) + B*cos(x[1])')
391
The parameters can only be scalars, and are all initialized to 0.0. The
392
parameters can also be given default values, using the argument 'defaults':
394
>>> f = Expression('A*sin(x[0]) + B*cos(x[1])',
395
defaults = {'A': 2.0,'B': 4.0})
397
*2. Complex user-defined JIT-compiled Expressions*
398
One may also define a Expression using more complicated logic with
399
the 'cppcode'. This argument should be a string of C++
400
code that implements a class that inherits from dolfin::Expression.
402
The following code illustrates how to define a Expression that depends
403
on material properties of the cells in a Mesh. A MeshFunction is
404
used to mark cells with different properties.
406
Note the use of the 'data' parameter.
409
... class MyFunc : public Expression
413
... MeshFunction<uint> *cell_data;
415
... MyFunc() : Expression(2), cell_data(0)
419
... void eval(Array<double>& values, const Data& data) const
421
... assert(cell_data);
422
... switch ((*cell_data)[data.cell()])
425
... values[0] = exp(-data.x[0]);
428
... values[0] = exp(-data.x[2]);
436
>>> cell_data = MeshFunction('uint', V.mesh(), 2)
437
>>> f = Expression(code)
438
>>> f.cell_data = cell_data
440
*3. User-defined expressions by subclassing*
441
The user can subclass Expression and overload the 'eval' function. The
442
value_shape of such an Expression will default to 0. If a user want a vector
443
or tensor Expression, he needs to over load the value_shape method.
445
>>> class MyExpression0(Expression):
446
... def eval(self, value, x):
449
... value[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02)
450
... value[1] = 250.0*exp(-(dx*dx + dy*dy)/0.01)
451
... def value_shape(self):
453
>>> f0 = MyExpression0()
455
If a user want to use the Expression in a UFL form, and have more controll in
456
which Finite element room the Expression shall be interpolated in, he can pass
457
this information using the element kwarg:
459
>>> V = FunctionSpace(mesh, "BDM", 1)
460
>>> f1 = MyExpression0(element = V.ufl_element())
462
The user can also subclass Expression overloading the eval_cell function. By
463
this the user get access to the more powerfull Data structure, with e.g., cell,
464
facet and normal information, during assemble.
466
>>> class MyExpression1(Expression):
467
... def eval_cell(self, value, ufc_cell):
468
... if ufc_cell.index > 10:
473
>>> f2 = MyExpression1()
475
The user can customize initialization of derived Expressions, however because of
476
magic behind the sceens a user need pass optional arguments to __init__ using
477
``**kwargs``, and _not_ calling the base class __init__:
479
>>> class MyExpression1(Expression):
480
... def __init__(self, mesh, domain):
481
... self._mesh = mesh
482
... self._domain = domain
483
... def eval(self, values, x):
486
>>> f3 = MyExpression1(mesh=mesh, domain=domain)
488
Note that subclassing may be significantly slower than using JIT-compiled
489
expressions. This is because a callback from C++ to Python will be involved
490
each time a Expression needs to be evaluated during assemble.
494
__metaclass__ = ExpressionMetaClass
496
def __new__(cls, cppcode=None, defaults=None, element=None, cell=None, \
497
degree=None, **kwargs):
499
Instantiate a new Expression
505
Optional C++ argument.
507
Optional ufl.FiniteElement argument
509
Optional element degree when element is not given
512
# If the __new__ function is called because we are instantiating a python sub
513
# class of Expression, then just return a new instant of the passed class
514
if cls.__name__ != "Expression":
515
return object.__new__(cls)
518
_check_cppcode(cppcode)
519
_check_defaults(defaults)
521
# Compile module and get the cpp.Expression class
522
cpp_base = compile_expressions([cppcode], [defaults])[0]
524
# Store compile arguments for later use
525
cpp_base.cppcode = cppcode
526
cpp_base.defaults = defaults
528
return object.__new__(create_compiled_expression_class(cpp_base))
530
# This method is only included so a user can check what arguments
531
# one should use in IPython using tab completion
532
def __init__(self, cppcode=None, defaults=None, element=None, cell=None,
533
degree=None, **kwargs): pass
535
# Reuse the docstring from __new__
536
__init__.__doc__ = __new__.__doc__
538
def ufl_element(self):
539
"Return the ufl FiniteElement."
540
return self._ufl_element
543
"x.__str__() <==> print(x)"
544
return "<Expression %s>" % str(self._value_shape)
547
"x.__repr__() <==> repr(x)"
548
return ufl.Coefficient.__repr__(self)
550
# Default value for the value shape
553
def value_shape(self):
554
"""Returns the value shape of the expression"""
555
return self._value_shape
557
__call__ = expression__call__
559
def Expressions(*args, **kwargs):
561
Batch-processed user-defined JIT-compiled expressions
563
By specifying several cppcodes one may compile more than one expression
566
>>> f0, f1 = Expressions('sin(x[0]) + cos(x[1])', 'exp(x[1])', degree=3)
568
>>> f0, f1, f2 = Expressions((('A*sin(x[0])', 'B*cos(x[1])')
569
('0','1')), {'A':2.0,'B':3.0},
571
(('cos(x[0])','sin(x[1])'),
572
('sin(x[0])','cos(x[1])')), element=element)
574
Here code is a C++ code snippet, which should be a string of C++
575
code that implements a class that inherits from dolfin::Expression,
576
see user case 3. in Expression docstring
578
Batch-processing of JIT-compiled expressions may significantly speed up
579
JIT-compilation at run-time.
582
# Get the element, cell and degree degree from kwarg
584
raise TypeError, "Can only define one kwarg and that can only be 'degree' or 'element'."
585
element = kwargs.pop("element", None)
586
cell = kwargs.pop("cell", None)
587
degree = kwargs.pop("degree", None)
589
# Iterate over the *args and collect input to compile_expressions
590
cppcodes = []; defaults = []; i = 0;
593
# Check type of cppcodes
594
if not isinstance(args[i],(tuple, list, str)):
595
raise TypeError, "Expected either a 'list', 'tuple' or 'str' for argument %d"%i
597
cppcodes.append(args[i])
600
# If we have more args and the next is a dict
601
if i < len(args) and isinstance(args[i], dict):
602
# Append the dict to defaults
603
_check_defaults(args[i])
604
defaults.append(args[i])
608
defaults.append(None)
610
# Compile the cpp.Expressions
611
cpp_bases = compile_expressions(cppcodes, defaults)
613
# Instantiate the return arguments
614
return_expressions = []
616
for cppcode, cpp_base in zip(cppcodes, cpp_bases):
617
return_expressions.append(create_compiled_expression_class(cpp_base)(\
622
# Return the instantiated Expressions
623
return tuple(return_expressions)
625
#--- Utility functions ---
627
def _check_cppcode(cppcode):
628
"Check that cppcode makes sense"
630
# Check that we get a string expression or nested expression
631
if not isinstance(cppcode, (str, tuple, list)):
632
raise TypeError, "Please provide a 'str', 'tuple' or 'list' for the 'cppcode' argument."
634
def _check_defaults(defaults):
635
"Check that defaults makes sense"
637
if defaults is None: return
639
# Check that we get a dictionary
640
if not isinstance(defaults, dict):
641
raise TypeError, "Please provide a 'dict' for the 'defaults' argument."
643
# Check types of the values in the dict
644
for key, val in defaults.iteritems():
645
if not isinstance(key, str):
646
raise TypeError, "All keys in 'defaults' must be a 'str'."
647
if not numpy.isscalar(val):
648
raise TypeError, "All values in 'defaults' must be scalars."
650
def _is_complex_expression(cppcode):
651
"Check if cppcode is a complex expression"
652
return isinstance(cppcode, str) and "class" in cppcode and "Expression" in cppcode
654
def _auto_select_element_from_shape(shape, degree=None, cell=None):
655
"Automatically select an appropriate element from cppcode."
657
# Default element, change to quadrature when working
660
# Check if scalar, vector or tensor valued
662
element = ufl.FiniteElement(Family, cell, degree)
663
elif len(shape) == 1:
664
element = ufl.VectorElement(Family, cell, degree, dim=shape[0])
666
element = ufl.TensorElement(Family, cell, degree, shape=shape)
668
cpp.debug("Automatic selection of expression element: " + str(element))
672
def _check_name_and_base(name, cpp_base):
674
assert(isinstance(name, str))
675
assert(name != "Expression"), "Cannot create a sub class of Expression with the same name as Expression"
677
assert(isinstance(cpp_base, (types.ClassType, type)))
added
removed
site-packages/dolfin/functions/expression.py
