File: test_multimesh_cell_types.py

package info (click to toggle)
dolfin 2018.1.0.post1-16
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 28,764 kB
  • sloc: xml: 104,040; cpp: 98,856; python: 22,511; makefile: 204; sh: 182
file content (98 lines) | stat: -rw-r--r-- 3,182 bytes parent folder | download | duplicates (3)
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
"""Unit tests for MultiMesh cell types"""

# Copyright (C) 2016 Magne Nordaas
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# First added:  2016-11-28
# Last changed: 2016-11-28

from dolfin import *
import pytest
from dolfin_utils.test import skip_in_parallel

# test case with interface-edge overlap
def case_1_impl(M, N):
    multimesh = MultiMesh()
    mesh0 = UnitSquareMesh(M, M)
    mesh1 = RectangleMesh(Point(0.25, 0.25), Point(0.75, 0.75), N, N)
    multimesh.add(mesh0)
    multimesh.add(mesh1)
    multimesh.build()
    return multimesh

# test case with squares on the diagonal
def case_2_impl(width, offset, Nx):

    # Mesh width (must be less than 1)
    assert width < 1

    # Mesh placement (must be less than the width)
    assert offset < width

    # Background mesh
    mesh_0 = UnitSquareMesh(Nx, Nx)

    # Create multimesh
    multimesh = MultiMesh()
    multimesh.add(mesh_0)

    # Now we have num_parts = 1
    num_parts = multimesh.num_parts()

    while num_parts*offset + width < 1:
        a = num_parts*offset
        b = a + width
        mesh_top = RectangleMesh(Point(a,a), Point(b,b), Nx, Nx)
        multimesh.add(mesh_top)
        num_parts = multimesh.num_parts()

    multimesh.build()
    return multimesh

@pytest.fixture
def case_1(M, N):
    return case_1_impl(M,N)

@pytest.fixture
def case_2(width, offset, Nx):
    return case_2_impl(width, offset, Nx)

test_cases = [case_1_impl(4,3),
              case_2_impl(DOLFIN_PI/5, 0.1111, 3)]

@skip_in_parallel
@pytest.mark.parametrize("multimesh", test_cases)
def test_cut_cell_has_quadrature(multimesh):
    # Test that every cut cell has a nontrivial interface quadrature rule
    for part in range(multimesh.num_parts()):
        for cell in multimesh.cut_cells(part):
            assert multimesh.quadrature_rules_interface(part, cell)

@skip_in_parallel
@pytest.mark.parametrize("multimesh", test_cases)
def test_multimesh_cell_types(multimesh):
    # Test that every cell in the multimesh is either cut, uncut, or covered
    for part in range(multimesh.num_parts()):
        cells = set(range(multimesh.part(part).num_cells()))
        cut_cells = set(multimesh.cut_cells(part))
        uncut_cells = set(multimesh.uncut_cells(part))
        covered_cells = set(multimesh.covered_cells(part))

        assert cut_cells.union(uncut_cells).union(covered_cells) == cells
        assert cut_cells.intersection(uncut_cells) == set()
        assert cut_cells.intersection(covered_cells) == set()
        assert uncut_cells.intersection(covered_cells) == set()