← Back to branch summary

~registry/fenics-book/main

« back to all changes in this revision

Viewing changes to chapters/langtangen/src/stationary/poisson/d6_p2D.py

  • Committer: Anders Logg
  • Date: 2011-10-26 20:12:17 UTC
  • Revision ID: logg@simula.no-20111026201217-vt684n4qmx4nfbv2
Update code for Chapter 1

Show diffs side-by-side

added added

removed removed

Lines of Context:
chapters/langtangen/src/stationary/poisson/d6_p2D.py
21
21
 
22
22
    # Exact solution
23
23
    omega = 1.0
24
 
    u_exact = Expression('sin(omega*pi*x[0])*sin(omega*pi*x[1])',
25
 
                         omega=omega)
 
24
    u_e = Expression('sin(omega*pi*x[0])*sin(omega*pi*x[1])',
 
25
                     omega=omega)
26
26
 
27
27
    # Define variational problem
28
28
    u = TrialFunction(V)
30
30
    #f = Function(V,
31
31
    #    '2*pow(pi,2)*pow(omega,2)*sin(omega*pi*x[0])*sin(omega*pi*x[1])',
32
32
    #    {'omega': omega})
33
 
    f = 2*pi**2*omega**2*u_exact
 
33
    f = 2*pi**2*omega**2*u_e
34
34
    a = inner(nabla_grad(u), nabla_grad(v))*dx
35
35
    L = f*v*dx
36
36
 
47
47
    # Compute error norm
48
48
 
49
49
    # Function - Expression
50
 
    error = (u - u_exact)**2*dx
 
50
    error = (u - u_e)**2*dx
51
51
    E1 = sqrt(assemble(error))
52
52
 
53
 
    # Explicit interpolation of u_exact onto the same space as u:
54
 
    u_e = interpolate(u_exact, V)
55
 
    error = (u - u_e)**2*dx
 
53
    # Explicit interpolation of u_e onto the same space as u:
 
54
    u_e_V = interpolate(u_e, V)
 
55
    error = (u - u_e_V)**2*dx
56
56
    E2 = sqrt(assemble(error))
57
57
 
58
 
    # Explicit interpolation of u_exact to higher-order elements,
 
58
    # Explicit interpolation of u_e to higher-order elements,
59
59
    # u will also be interpolated to the space Ve before integration
60
60
    Ve = FunctionSpace(mesh, 'Lagrange', degree=5)
61
 
    u_e = interpolate(u_exact, Ve)
62
 
    error = (u - u_e)**2*dx
 
61
    u_e_Ve = interpolate(u_e, Ve)
 
62
    error = (u - u_e_Ve)**2*dx
63
63
    E3 = sqrt(assemble(error))
64
64
 
65
 
    # errornorm interpolates u and u_exact to a space with
 
65
    # errornorm interpolates u and u_e to a space with
66
66
    # given degree, and creates the error field by subtracting
67
67
    # the degrees of freedom, then the error field is integrated
68
 
    # TEMPORARY BUG - doesn't accept Expression for u_exact
69
 
    #E4 = errornorm(u_exact, u, normtype='l2', degree=3)
 
68
    # TEMPORARY BUG - doesn't accept Expression for u_e
 
69
    #E4 = errornorm(u_e, u, normtype='l2', degree=3)
70
70
    # Manual implementation
71
 
    def errornorm(u_exact, u, Ve):
 
71
    def errornorm(u_e, u, Ve):
72
72
        u_Ve = interpolate(u, Ve)
73
 
        u_e_Ve = interpolate(u_exact, Ve)
 
73
        u_e_Ve = interpolate(u_e, Ve)
74
74
        e_Ve = Function(Ve)
75
75
        # Subtract degrees of freedom for the error field
76
 
        e_Ve.vector()[:] = u_e_Ve.vector().array() - \
77
 
                           u_Ve.vector().array()
 
76
        e_Ve.vector()[:] = u_e_Ve.vector().array() - u_Ve.vector().array()
78
77
        # More efficient computation (avoids the rhs array result above)
79
78
        #e_Ve.assign(u_e_Ve)                      # e_Ve = u_e_Ve
80
79
        #e_Ve.vector().axpy(-1.0, u_Ve.vector())  # e_Ve += -1.0*u_Ve
81
80
        error = e_Ve**2*dx
82
81
        return sqrt(assemble(error, mesh=Ve.mesh())), e_Ve
83
 
    E4, e_Ve = errornorm(u_exact, u, Ve)
 
82
    E4, e_Ve = errornorm(u_e, u, Ve)
84
83
 
85
84
    # Infinity norm based on nodal values
86
 
    u_e = interpolate(u_exact, V)
87
 
    E5 = abs(u_e.vector().array() - u.vector().array()).max()
 
85
    u_e_V = interpolate(u_e, V)
 
86
    E5 = abs(u_e_V.vector().array() - u.vector().array()).max()
88
87
    print 'E2:', E2
89
88
    print 'E3:', E3
90
89
    print 'E4:', E4
95
94
    E6 = sqrt(assemble(error))
96
95
 
97
96
    # Collect error measures in a dictionary with self-explanatory keys
98
 
    errors = {'u - u_exact': E1,
99
 
              'u - interpolate(u_exact,V)': E2,
100
 
              'interpolate(u,Ve) - interpolate(u_exact,Ve)': E3,
 
97
    errors = {'u - u_e': E1,
 
98
              'u - interpolate(u_e,V)': E2,
 
99
              'interpolate(u,Ve) - interpolate(u_e,Ve)': E3,
101
100
              'error field': E4,
102
101
              'infinity norm (of dofs)': E5,
103
102
              'grad(error field)': E6}
104
 
    
 
103
 
105
104
    return errors
106
105
 
107
106
# Perform experiments
123
122
        r = ln(Ei/Eim1)/ln(h[i]/h[i-1])
124
123
        print 'h=%8.2E E=%8.2E r=%.2f' % (h[i], Ei, r)
125
124
 
126
 
    
 
125
 
127
126