5839.1.7
by Johan Hake
Rename python subdirectories avoiding name mangling |
1 |
"""This module handles the Expression class in Python.
|
2 |
||
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.
|
|
7 |
||
8 |
The resulting Expression class may thus act both as a variable in a
|
|
9 |
UFL form expression and as a DOLFIN C++ Expression.
|
|
10 |
||
11 |
This module make heavy use of creation of Expression classes and
|
|
12 |
instantiation of these dynamically at runtime.
|
|
13 |
||
14 |
The whole logic behind this somewhat magic behaviour is handle by:
|
|
15 |
||
16 |
1) function __new__ in the Expression class
|
|
17 |
2) meta class ExpressionMetaClass
|
|
18 |
3) function compile_expressions from the compiledmodule/expression
|
|
19 |
module
|
|
20 |
4) functions create_compiled_expression_class and
|
|
21 |
create_python_derived_expression_class
|
|
22 |
||
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.
|
|
26 |
||
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.
|
|
30 |
||
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
|
|
35 |
"""
|
|
36 |
||
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" |
|
41 |
||
42 |
# Modified by Anders Logg, 2008-2009.
|
|
43 |
||
44 |
__all__ = ["Expression", "Expressions"] |
|
45 |
||
46 |
# FIXME: Make all error messages uniform according to the following template:
|
|
47 |
#
|
|
48 |
# if not isinstance(foo, Foo):
|
|
49 |
# raise TypeError, "Illegal argument for creation of Bar, not a Foo: " + str(foo)
|
|
50 |
||
51 |
# Python imports
|
|
52 |
import types |
|
53 |
||
54 |
# Import UFL and SWIG-generated extension module (DOLFIN C++)
|
|
55 |
import ufl |
|
56 |
import dolfin.cpp as cpp |
|
57 |
import numpy |
|
58 |
||
59 |
# Local imports
|
|
60 |
from dolfin.compilemodules.expressions import compile_expressions |
|
61 |
||
62 |
def create_compiled_expression_class(cpp_base): |
|
63 |
# Check the cpp_base
|
|
64 |
assert(isinstance(cpp_base, (types.ClassType, type))) |
|
65 |
||
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())) |
|
72 |
||
73 |
# Store the value_shape
|
|
74 |
self._value_shape = value_shape |
|
75 |
||
76 |
# Select an appropriate element if not specified
|
|
77 |
if element is None: |
|
78 |
element = _auto_select_element_from_shape(value_shape, degree) |
|
79 |
else: |
|
80 |
# Check that we have an element
|
|
81 |
if not isinstance(element, ufl.FiniteElementBase): |
|
82 |
raise TypeError, "The 'element' argument must be a UFL"\ |
|
83 |
" finite element."
|
|
84 |
||
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 "\
|
|
89 |
"expression."
|
|
90 |
||
91 |
# Initialize UFL base class
|
|
92 |
self._ufl_element = element |
|
93 |
ufl.Coefficient.__init__(self, self._ufl_element) |
|
94 |
||
95 |
# Create and return the class
|
|
96 |
return type("CompiledExpression", (Expression, ufl.Coefficient, cpp_base),\ |
|
97 |
{"__init__":__init__}) |
|
98 |
||
99 |
def create_python_derived_expression_class(name, user_bases, user_dict): |
|
100 |
"""Return Expression class
|
|
101 |
||
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
|
|
105 |
and ufl.Coefficient.
|
|
106 |
||
107 |
*Arguments*
|
|
108 |
name
|
|
109 |
The name of the class
|
|
110 |
user_bases
|
|
111 |
User defined bases
|
|
112 |
user_dict
|
|
113 |
Dict with user specified function or attributes
|
|
114 |
"""
|
|
115 |
||
116 |
# Check args
|
|
117 |
assert(isinstance(name, str)) |
|
118 |
assert(isinstance(user_bases, list)) |
|
119 |
assert(isinstance(user_dict, dict)) |
|
120 |
||
121 |
# Define the bases
|
|
122 |
assert(all([isinstance(t, (types.ClassType, type)) for t in user_bases])) |
|
123 |
bases = tuple([Expression, ufl.Coefficient, cpp.Expression] + user_bases) |
|
124 |
||
125 |
user_init = user_dict.pop("__init__", None) |
|
126 |
||
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"
|
|
131 |
||
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) |
|
137 |
||
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'"
|
|
148 |
if len(args) != 0: |
|
149 |
raise TypeError, "expected no 'args'" |
|
150 |
||
151 |
# Select an appropriate element if not specified
|
|
152 |
if element is None: |
|
153 |
element = _auto_select_element_from_shape(self.value_shape(), |
|
154 |
degree) |
|
155 |
elif isinstance(element, ufl.FiniteElementBase): |
|
156 |
pass
|
|
157 |
else: |
|
158 |
raise TypeError, "The 'element' argument must be a UFL finite"\ |
|
159 |
" element."
|
|
160 |
||
161 |
# Initialize UFL base class
|
|
162 |
self._ufl_element = element |
|
163 |
ufl.Coefficient.__init__(self, element) |
|
164 |
||
165 |
# Initialize cpp_base class
|
|
166 |
cpp.Expression.__init__(self, list(element.value_shape())) |
|
167 |
||
168 |
# Calling the user defined_init
|
|
169 |
if user_init is not None: |
|
170 |
user_init(self, *args, **kwargs) |
|
171 |
||
172 |
# NOTE: Do not prevent the user to overload attributes "reserved" by PyDOLFIN
|
|
173 |
# NOTE: Why not?
|
|
174 |
||
175 |
## Collect reserved attributes from both cpp.Function and ufl.Coefficient
|
|
176 |
#reserved_attr = dir(ufl.Coefficient)
|
|
177 |
#reserved_attr.extend(dir(cpp.Function))
|
|
178 |
#
|
|
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)
|
|
183 |
#
|
|
184 |
## Check the dict_ for reserved attributes
|
|
185 |
#for attr in reserved_attr:
|
|
186 |
# if attr in dict_:
|
|
187 |
# raise TypeError, "The Function attribute '%s' is reserved by PyDOLFIN."%attr
|
|
188 |
||
189 |
# Add __init__ to the user_dict
|
|
190 |
user_dict["__init__"] = __init__ |
|
191 |
||
192 |
# Create the class and return it
|
|
193 |
return type(name, bases, user_dict) |
|
194 |
||
195 |
class ExpressionMetaClass(type): |
|
196 |
"""Meta Class for Expression"""
|
|
197 |
def __new__(mcs, name, bases, dict_): |
|
198 |
"""Returns a new Expression class"""
|
|
199 |
||
200 |
assert(isinstance(name, str)), "Expecting a 'str'" |
|
201 |
assert(isinstance(bases, tuple)), "Expecting a 'tuple'" |
|
202 |
assert(isinstance(dict_, dict)), "Expecting a 'dict'" |
|
203 |
||
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:
|
|
208 |
#
|
|
209 |
# class Expression(Expression):
|
|
210 |
# ...
|
|
211 |
if len(bases) > 1 and bases[0] != object: |
|
212 |
raise TypeError, "Cannot name a subclass of Expression:"\ |
|
213 |
" 'Expression'"
|
|
214 |
||
215 |
# Return the new class, which just is the original Expression defined in
|
|
216 |
# this module
|
|
217 |
return type.__new__(mcs, name, bases, dict_) |
|
218 |
||
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], \ |
|
224 |
cpp.Expression): |
|
225 |
# Return the instantiated class
|
|
226 |
return type.__new__(mcs, name, bases, dict_) |
|
227 |
||
228 |
# Handle any user provided base classes
|
|
229 |
user_bases = list(bases) |
|
230 |
||
231 |
# remove Expression, to be added later
|
|
232 |
user_bases.remove(Expression) |
|
233 |
||
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" |
|
237 |
||
238 |
# Get name of eval function
|
|
239 |
eval_name = 'eval' if 'eval' in dict_ else 'eval_cell' |
|
240 |
||
241 |
user_eval = dict_[eval_name] |
|
242 |
||
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 |
|
253 |
||
254 |
return create_python_derived_expression_class(name, user_bases, dict_) |
|
255 |
||
256 |
#--- The user interface ---
|
|
257 |
||
258 |
# Places here so it can be reused in Function
|
|
259 |
def expression__call__(self, *args, **kwargs): |
|
260 |
"""
|
|
261 |
Evaluates the Expression
|
|
262 |
||
263 |
*Example*
|
|
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])
|
|
269 |
>>> v0 = fs(x0)
|
|
270 |
>>> v1 = fs(x1)
|
|
271 |
>>> v2 = fs(x2)
|
|
272 |
||
273 |
2) Using multiple scalar args for x, interpreted as a point coordinate
|
|
274 |
>>> v0 = f(1.,0.5,0.5)
|
|
275 |
||
276 |
3) Passing return array
|
|
277 |
>>> fv = Expression(("sin(x[0])*cos(x[1])*sin(x[3])",
|
|
278 |
"2.0","0.0"))
|
|
279 |
>>> x0 = numpy.array([1.,0.5,0.5])
|
|
280 |
>>> v0 = numpy.zeros(3)
|
|
281 |
>>> fv(x0, values = v0)
|
|
282 |
||
283 |
.. note::
|
|
284 |
||
285 |
A longer values array may be passed. In this way one can fast
|
|
286 |
fill up an array with different evaluations.
|
|
287 |
||
288 |
>>> values = numpy.zeros(9)
|
|
289 |
>>> for i in xrange(0,10,3):
|
|
290 |
... fv(x[i:i+3], values = values[i:i+3])
|
|
291 |
||
292 |
"""
|
|
293 |
if len(args)==0: |
|
294 |
raise TypeError, "expected at least 1 argument" |
|
295 |
||
296 |
# Test for ufl restriction
|
|
297 |
if len(args) == 1 and args[0] in ('+','-'): |
|
298 |
return ufl.Coefficient.__call__(self, *args) |
|
299 |
||
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) |
|
303 |
||
304 |
# Some help variables
|
|
305 |
value_size = ufl.common.product(self.ufl_element().value_shape()) |
|
306 |
||
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 |
|
317 |
else: |
|
318 |
values_provided = False |
|
319 |
values = numpy.zeros(value_size, dtype='d') |
|
320 |
||
321 |
# Assume all args are x argument
|
|
322 |
x = args |
|
323 |
||
324 |
# If only one x argument has been provided, unpack it if it's an iterable
|
|
325 |
if len(x) == 1: |
|
326 |
if hasattr(x[0], '__iter__'): |
|
327 |
x = x[0] |
|
328 |
||
329 |
# Convert it to an 1D numpy array
|
|
330 |
try: |
|
331 |
x = numpy.fromiter(x, 'd') |
|
332 |
except (TypeError, ValueError, AssertionError): |
|
333 |
raise TypeError, "expected scalar arguments for the coordinates" |
|
334 |
||
335 |
if len(x) == 0: |
|
336 |
raise TypeError, "coordinate argument too short" |
|
337 |
||
338 |
# The actual evaluation
|
|
339 |
self.eval(values, x) |
|
340 |
||
341 |
# If scalar return statement, return scalar value.
|
|
342 |
if value_size == 1 and not values_provided: |
|
343 |
return values[0] |
|
344 |
||
345 |
return values |
|
346 |
||
347 |
class Expression(object): |
|
348 |
"""
|
|
349 |
This class represents a user-defined expression.
|
|
350 |
||
351 |
Expressions can be used as coefficients in variational forms or
|
|
352 |
interpolated into finite element spaces.
|
|
353 |
||
354 |
*Arguments*
|
|
355 |
cppcode
|
|
356 |
C++ argument, see below
|
|
357 |
defaults
|
|
358 |
Optional C++ argument, see below
|
|
359 |
element
|
|
360 |
Optional element argument
|
|
361 |
||
362 |
*1. Simple user-defined JIT-compiled expressions*
|
|
363 |
One may alternatively specify a C++ code for evaluation of the
|
|
364 |
Expression as follows:
|
|
365 |
||
366 |
>>> f0 = Expression('sin(x[0]) + cos(x[1])')
|
|
367 |
>>> f1 = Expression(('cos(x[0])', 'sin(x[1])'), element = V.ufl_element())
|
|
368 |
||
369 |
Here, f0 is is scalar and f1 is vector-valued.
|
|
370 |
||
371 |
Tensor expressions of rank 2 (matrices) may also be created:
|
|
372 |
||
373 |
>>> f2 = Expression((('exp(x[0])','sin(x[1])'),
|
|
374 |
('sin(x[0])','tan(x[1])')))
|
|
375 |
||
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).
|
|
379 |
||
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
|
|
383 |
the user.
|
|
384 |
||
385 |
Expression parameters can be included as follows:
|
|
386 |
||
387 |
>>> f = Expression('A*sin(x[0]) + B*cos(x[1])')
|
|
388 |
>>> f.A = 2.0
|
|
389 |
>>> f.B = 4.0
|
|
390 |
||
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':
|
|
393 |
||
394 |
>>> f = Expression('A*sin(x[0]) + B*cos(x[1])',
|
|
395 |
defaults = {'A': 2.0,'B': 4.0})
|
|
396 |
||
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.
|
|
401 |
||
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.
|
|
405 |
||
406 |
Note the use of the 'data' parameter.
|
|
407 |
||
408 |
>>> code = '''
|
|
409 |
... class MyFunc : public Expression
|
|
410 |
... {
|
|
411 |
... public:
|
|
412 |
...
|
|
413 |
... MeshFunction<uint> *cell_data;
|
|
414 |
...
|
|
415 |
... MyFunc() : Expression(2), cell_data(0)
|
|
416 |
... {
|
|
417 |
... }
|
|
418 |
...
|
|
419 |
... void eval(Array<double>& values, const Data& data) const
|
|
420 |
... {
|
|
421 |
... assert(cell_data);
|
|
422 |
... switch ((*cell_data)[data.cell()])
|
|
423 |
... {
|
|
424 |
... case 0:
|
|
425 |
... values[0] = exp(-data.x[0]);
|
|
426 |
... break;
|
|
427 |
... case 1:
|
|
428 |
... values[0] = exp(-data.x[2]);
|
|
429 |
... break;
|
|
430 |
... default:
|
|
431 |
... values[0] = 0.0;
|
|
432 |
... }
|
|
433 |
... }
|
|
434 |
... };'''
|
|
435 |
||
436 |
>>> cell_data = MeshFunction('uint', V.mesh(), 2)
|
|
437 |
>>> f = Expression(code)
|
|
438 |
>>> f.cell_data = cell_data
|
|
439 |
||
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.
|
|
444 |
||
445 |
>>> class MyExpression0(Expression):
|
|
446 |
... def eval(self, value, x):
|
|
447 |
... dx = x[0] - 0.5
|
|
448 |
... dy = x[1] - 0.5
|
|
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):
|
|
452 |
... return (2,)
|
|
453 |
>>> f0 = MyExpression0()
|
|
454 |
||
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:
|
|
458 |
||
459 |
>>> V = FunctionSpace(mesh, "BDM", 1)
|
|
460 |
>>> f1 = MyExpression0(element = V.ufl_element())
|
|
461 |
||
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.
|
|
465 |
||
466 |
>>> class MyExpression1(Expression):
|
|
467 |
... def eval_cell(self, value, ufc_cell):
|
|
468 |
... if ufc_cell.index > 10:
|
|
469 |
... value[0] = 1.0
|
|
470 |
... else:
|
|
471 |
... value[0] = -1.0
|
|
472 |
||
473 |
>>> f2 = MyExpression1()
|
|
474 |
||
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__:
|
|
478 |
||
479 |
>>> class MyExpression1(Expression):
|
|
480 |
... def __init__(self, mesh, domain):
|
|
481 |
... self._mesh = mesh
|
|
482 |
... self._domain = domain
|
|
483 |
... def eval(self, values, x):
|
|
484 |
... ...
|
|
485 |
||
486 |
>>> f3 = MyExpression1(mesh=mesh, domain=domain)
|
|
487 |
||
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.
|
|
491 |
"""
|
|
492 |
||
493 |
# Set the meta class
|
|
494 |
__metaclass__ = ExpressionMetaClass |
|
495 |
||
496 |
def __new__(cls, cppcode=None, defaults=None, element=None, cell=None, \ |
|
497 |
degree=None, **kwargs): |
|
498 |
"""
|
|
499 |
Instantiate a new Expression
|
|
500 |
||
501 |
*Arguments*
|
|
502 |
cppcode
|
|
503 |
C++ argument.
|
|
504 |
defaults
|
|
505 |
Optional C++ argument.
|
|
506 |
element
|
|
507 |
Optional ufl.FiniteElement argument
|
|
508 |
degree
|
|
509 |
Optional element degree when element is not given
|
|
510 |
"""
|
|
511 |
||
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) |
|
516 |
||
517 |
# Check arguments
|
|
518 |
_check_cppcode(cppcode) |
|
519 |
_check_defaults(defaults) |
|
520 |
||
521 |
# Compile module and get the cpp.Expression class
|
|
522 |
cpp_base = compile_expressions([cppcode], [defaults])[0] |
|
523 |
||
524 |
# Store compile arguments for later use
|
|
525 |
cpp_base.cppcode = cppcode |
|
526 |
cpp_base.defaults = defaults |
|
527 |
||
528 |
return object.__new__(create_compiled_expression_class(cpp_base)) |
|
529 |
||
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 |
|
534 |
||
535 |
# Reuse the docstring from __new__
|
|
536 |
__init__.__doc__ = __new__.__doc__ |
|
537 |
||
538 |
def ufl_element(self): |
|
539 |
"Return the ufl FiniteElement."
|
|
540 |
return self._ufl_element |
|
541 |
||
542 |
def __str__(self): |
|
543 |
"x.__str__() <==> print(x)"
|
|
544 |
return "<Expression %s>" % str(self._value_shape) |
|
545 |
||
546 |
def __repr__(self): |
|
547 |
"x.__repr__() <==> repr(x)"
|
|
548 |
return ufl.Coefficient.__repr__(self) |
|
549 |
||
550 |
# Default value for the value shape
|
|
551 |
_value_shape = () |
|
552 |
||
553 |
def value_shape(self): |
|
554 |
"""Returns the value shape of the expression"""
|
|
555 |
return self._value_shape |
|
556 |
||
557 |
__call__ = expression__call__ |
|
558 |
||
559 |
def Expressions(*args, **kwargs): |
|
560 |
"""
|
|
561 |
Batch-processed user-defined JIT-compiled expressions
|
|
562 |
||
563 |
By specifying several cppcodes one may compile more than one expression
|
|
564 |
at a time:
|
|
565 |
||
566 |
>>> f0, f1 = Expressions('sin(x[0]) + cos(x[1])', 'exp(x[1])', degree=3)
|
|
567 |
||
568 |
>>> f0, f1, f2 = Expressions((('A*sin(x[0])', 'B*cos(x[1])')
|
|
569 |
('0','1')), {'A':2.0,'B':3.0},
|
|
570 |
code,
|
|
571 |
(('cos(x[0])','sin(x[1])'),
|
|
572 |
('sin(x[0])','cos(x[1])')), element=element)
|
|
573 |
||
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
|
|
577 |
||
578 |
Batch-processing of JIT-compiled expressions may significantly speed up
|
|
579 |
JIT-compilation at run-time.
|
|
580 |
"""
|
|
581 |
||
582 |
# Get the element, cell and degree degree from kwarg
|
|
583 |
if len(kwargs) > 1: |
|
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) |
|
588 |
||
589 |
# Iterate over the *args and collect input to compile_expressions
|
|
590 |
cppcodes = []; defaults = []; i = 0; |
|
591 |
while i < len(args): |
|
592 |
||
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 |
|
596 |
||
597 |
cppcodes.append(args[i]) |
|
598 |
i += 1 |
|
599 |
||
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]) |
|
605 |
i += 1 |
|
606 |
else: |
|
607 |
# If not append None
|
|
608 |
defaults.append(None) |
|
609 |
||
610 |
# Compile the cpp.Expressions
|
|
611 |
cpp_bases = compile_expressions(cppcodes, defaults) |
|
612 |
||
613 |
# Instantiate the return arguments
|
|
614 |
return_expressions = [] |
|
615 |
||
616 |
for cppcode, cpp_base in zip(cppcodes, cpp_bases): |
|
617 |
return_expressions.append(create_compiled_expression_class(cpp_base)(\ |
|
618 |
cppcode, |
|
619 |
degree=degree, |
|
620 |
element=element)) |
|
621 |
||
622 |
# Return the instantiated Expressions
|
|
623 |
return tuple(return_expressions) |
|
624 |
||
625 |
#--- Utility functions ---
|
|
626 |
||
627 |
def _check_cppcode(cppcode): |
|
628 |
"Check that cppcode makes sense"
|
|
629 |
||
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." |
|
633 |
||
634 |
def _check_defaults(defaults): |
|
635 |
"Check that defaults makes sense"
|
|
636 |
||
637 |
if defaults is None: return |
|
638 |
||
639 |
# Check that we get a dictionary
|
|
640 |
if not isinstance(defaults, dict): |
|
641 |
raise TypeError, "Please provide a 'dict' for the 'defaults' argument." |
|
642 |
||
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." |
|
649 |
||
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 |
|
653 |
||
654 |
def _auto_select_element_from_shape(shape, degree=None, cell=None): |
|
655 |
"Automatically select an appropriate element from cppcode."
|
|
656 |
||
657 |
# Default element, change to quadrature when working
|
|
658 |
Family = "Lagrange" |
|
659 |
||
660 |
# Check if scalar, vector or tensor valued
|
|
661 |
if len(shape) == 0: |
|
662 |
element = ufl.FiniteElement(Family, cell, degree) |
|
663 |
elif len(shape) == 1: |
|
664 |
element = ufl.VectorElement(Family, cell, degree, dim=shape[0]) |
|
665 |
else: |
|
666 |
element = ufl.TensorElement(Family, cell, degree, shape=shape) |
|
667 |
||
668 |
cpp.debug("Automatic selection of expression element: " + str(element)) |
|
669 |
||
670 |
return element |
|
671 |
||
672 |
def _check_name_and_base(name, cpp_base): |
|
673 |
# Check the name
|
|
674 |
assert(isinstance(name, str)) |
|
675 |
assert(name != "Expression"), "Cannot create a sub class of Expression with the same name as Expression" |
|
676 |
||
677 |
assert(isinstance(cpp_base, (types.ClassType, type))) |