/
test_8.py
49 lines (40 loc) · 1.02 KB
/
test_8.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
from dolfin import *
class MyExpression(Expression):
def __init__(self, a=0, b=0):
self.a = a
self.b = b
def eval(self, value, x):
dx = x[0] - self.a
dy = x[1] - self.b
value[0] = exp(-(dx*dx + dy*dy)/0.02)
value[1] = exp(-(dx*dx + dy*dy)/0.02)
def value_shape(self):
return (2,)
cpp_code = '''
class MyCppExpression : public Expression
{
public:
MyCppExpression() : Expression(2), a(0), b(0) { }
void eval(Array<double>& values, const Array<double>& x) const
{
const double dx = x[0] - a;
const double dy = x[1] - b;
values[0] = exp(-(dx*dx + dy*dy)/0.02);
values[1] = exp(-(dx*dx + dy*dy)/0.02);
}
public:
double a, b;
};'''
mesh = UnitSquareMesh(20, 20)
f = MyExpression(a=0.5, b=0.5)
f.a = 0.5
f.b = 0.25
f_cpp = Expression(cpp_code)
f_cpp.a = 0.5
f_cpp.b = 0.25
# Compare the two
V = VectorFunctionSpace(mesh, 'CG', 1)
F = interpolate(f, V).vector()
F_cpp = interpolate(f_cpp, V).vector()
F.axpy(-1, F_cpp)
print '|F-F_cpp|', F.norm('l2')