← Back to branch summary

~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to test/unit/python/adaptivity/test_error_control.py

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2015-03-17 07:57:11 UTC
  • mfrom: (1.1.18) (19.1.24 experimental)
  • Revision ID: package-import@ubuntu.com-20150317075711-1v207zbty9qmygow
Tags: 1.5.0-1
http://bugs.debian.org/780359
* New upstream release (closes: #780359).
* debian/control:
  - Bump Standards-Version to 3.9.6 (no changes needed).
  - Bump X-Python-Version to >= 2.7.
  - Update package names for new SONAME 1.5 (libdolfin1.4 ->
    libdolfin1.5, libdolfin1.4-dbg -> libdolfin1.5-dbg and
    libdolfin1.4-dev -> libdolfin1.5-dev).
  - Bump minimum required version for python-instant, python-ufl and
    python-ffc to 1.5.0.
  - Add python-sympy and python-six to Depends for binary package
    python-dolfin.
  - Add dh-python to Build-Depends.
  - Remove libcgal-dev from {Build-}Depends.
* Remove CSGCGALMeshGenerator3D-oom.patch since CGAL is no longer used
  by DOLFIN.
* Move debian/libdolfin1.4.install -> debian/libdolfin1.5.install.
* debian/rules: No longer any non DFSG-free stuff to remove, so update
  get-orig-source target (update debian/watch accordingly).
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/unit/python/adaptivity/test_error_control.py
 
1
#!/usr/bin/env py.test
 
2
 
 
3
"""Unit tests for error control"""
 
4
 
 
5
# Copyright (C) 2011 Marie E. Rognes
 
6
#
 
7
# This file is part of DOLFIN.
 
8
#
 
9
# DOLFIN is free software: you can redistribute it and/or modify
 
10
# it under the terms of the GNU Lesser General Public License as published by
 
11
# the Free Software Foundation, either version 3 of the License, or
 
12
# (at your option) any later version.
 
13
#
 
14
# DOLFIN is distributed in the hope that it will be useful,
 
15
# but WITHOUT ANY WARRANTY; without even the implied warranty of
 
16
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 
17
# GNU Lesser General Public License for more details.
 
18
#
 
19
# You should have received a copy of the GNU Lesser General Public License
 
20
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 
21
 
 
22
import pytest
 
23
from ufl.algorithms import replace
 
24
 
 
25
from dolfin import *
 
26
from dolfin.fem.adaptivesolving import *
 
27
from dolfin_utils.test import skip_in_parallel
 
28
 
 
29
# FIXME: Move this to dolfin for user access?
 
30
def reconstruct_refined_form(form, functions, mesh):
 
31
    function_mapping = {}
 
32
    for u in functions:
 
33
        w = Function(u.leaf_node().function_space())
 
34
        w.assign(u.leaf_node())
 
35
        function_mapping[u] = w
 
36
    domain = mesh.leaf_node().ufl_domain()
 
37
    newform = replace_integral_domains(replace(form, function_mapping), domain)
 
38
    return newform, function_mapping
 
39
 
 
40
 
 
41
# This must be scope function, because the tests will modify some of the objects,
 
42
# including the mesh which gets its hierarchial adapted submeshes attached.
 
43
fixt = pytest.fixture(scope="function")
 
44
 
 
45
@fixt
 
46
def mesh():
 
47
    return UnitSquareMesh(8, 8)
 
48
 
 
49
@fixt
 
50
def V(mesh):
 
51
    return FunctionSpace(mesh, "Lagrange", 1)
 
52
 
 
53
@fixt
 
54
def u(V):
 
55
    return Function(V)
 
56
 
 
57
@fixt
 
58
def a(V):
 
59
    v = TestFunction(V)
 
60
    u = TrialFunction(V)
 
61
    return inner(grad(u), grad(v))*dx()
 
62
 
 
63
@fixt
 
64
def L(V):
 
65
    v = TestFunction(V)
 
66
    f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=1)
 
67
    g = Expression("sin(5*x[0])", degree=1)
 
68
    return f*v*dx() + g*v*ds()
 
69
 
 
70
@fixt
 
71
def problem(a, L, u, V):
 
72
    bc = [DirichletBC(V, 0.0, "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS")]
 
73
    return LinearVariationalProblem(a, L, u, bc)
 
74
 
 
75
@fixt
 
76
def goal(u):
 
77
    return u*dx()
 
78
 
 
79
@fixt
 
80
def ec(problem, goal):
 
81
    return generate_error_control(problem, goal)
 
82
 
 
83
 
 
84
@skip_in_parallel
 
85
def test_check_domains(goal, mesh, a, L):
 
86
    # Asserting that domains are ok before trying error control generation
 
87
    msg = "Expecting only the domain from the mesh to get here through u."
 
88
    assert len(goal.domains()) == 1, msg
 
89
    assert goal.domains()[0] == mesh.ufl_domain(), msg
 
90
    assert len(a.domains()) == 1, msg
 
91
    assert a.domains()[0] == mesh.ufl_domain(), msg
 
92
    assert len(L.domains()) == 1, msg
 
93
    assert L.domains()[0] == mesh.ufl_domain(), msg
 
94
 
 
95
 
 
96
@skip_in_parallel
 
97
def test_error_estimation(problem, u, ec):
 
98
 
 
99
    # Solve variational problem once
 
100
    solver = LinearVariationalSolver(problem)
 
101
    solver.solve()
 
102
 
 
103
    # Compute error estimate
 
104
    error_estimate = ec.estimate_error(u, problem.bcs())
 
105
 
 
106
    # Compare estimate with defined reference
 
107
    reference = 0.0011789985750808342
 
108
    assert round(error_estimate - reference, 7) == 0
 
109
 
 
110
 
 
111
@skip_in_parallel
 
112
def test_error_indicators(problem, u, mesh):
 
113
 
 
114
    # Solve variational problem once
 
115
    solver = LinearVariationalSolver(problem)
 
116
    solver.solve()
 
117
 
 
118
    # Compute error indicators
 
119
    indicators = Vector(mesh.mpi_comm(), u.function_space().mesh().num_cells())
 
120
    indicators[0] = 1.0
 
121
    #ec.compute_indicators(indicators, u) #
 
122
 
 
123
    reference = 1.0 # FIXME
 
124
    assert round(indicators.sum() - reference, 7) == 0
 
125
 
 
126
 
 
127
@skip_in_parallel
 
128
def _test_adaptive_solve(problem, goal, u, mesh):
 
129
 
 
130
    # Solve problem adaptively
 
131
    solver = AdaptiveLinearVariationalSolver(problem, goal)
 
132
    tol = 0.00087
 
133
    solver.solve(tol)
 
134
 
 
135
    # Note: This old approach is now broken, as it doesn't change the integration domain:
 
136
    #M = replace(goal, {u: w})
 
137
    # This new approach handles the integration domain properly:
 
138
    M, fm = reconstruct_refined_form(goal, [u], mesh)
 
139
 
 
140
    # Compare computed goal with reference
 
141
    reference = 0.12583303389560166
 
142
    assert round(assemble(M) - reference, 7) == 0