File: basistest.py

package info (click to toggle)
dune-functions 2.10.0-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 1,544 kB
  • sloc: cpp: 14,241; python: 661; makefile: 3
file content (160 lines) | stat: -rw-r--r-- 5,446 bytes parent folder | download
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
# SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later

from array import array

# Check consistency of basis.size(prefix)
def checkBasisSizeConsistency(basis, multiIndexSet):

    # Based on the index tree, build a map that contains all possible
    # prefixes and map each prefix to the size (of the subsequent digit).
    prefixSet = dict()

    for index in multiIndexSet:
        prefix = [];
        for digit in index:
            if not tuple(prefix) in prefixSet:
                prefixSet[tuple(prefix)] = 0
            prefixSet[tuple(prefix)] = max(prefixSet[tuple(prefix)], digit+1);
            prefix.append(digit);

        prefixSet[tuple(prefix)] = 0;


    # Now check for all prefixes, if the size computed from the
    # index tree is consistent with basis.size(prefix).
    for prefix, size in prefixSet.items():

        prefixSize = basis.size(prefix);

        if prefixSize != size:
            raise IndexError("basis.size(" + str(prefix) + ") = " + str(prefixSize) + ", but should be " + str(size))


# Check indices of basis:
# - First store the whole index tree in a set
# - Check if this corresponds to a consistent index tree
# - Check if index tree is consistent with basis.size(prefix) and basis.dimension()
def checkBasisIndices(basis):

    multiIndexSet = set()

    localView = basis.localView()
    gridView = basis.gridView

    for element in gridView.elements:
        localView.bind(element)

        for i in range(len(localView)):
            multiIndex = localView.index(i)
            for digit in multiIndex:
                if digit < 0:
                    raise IndexError(str(element) + ": Global multi-index contains negative entry for shape function " + str(i))

            multiIndexSet.add(tuple(multiIndex))

    # TODO: Implement this!
    # checkBasisIndexTreeConsistency(multiIndexSet);
    checkBasisSizeConsistency(basis, multiIndexSet);

    if (basis.dimension != len(multiIndexSet)):
        raise ValueError("basis.dimension() does not equal the total number of basis functions.")


# Check the methods of individual nodes of a function bases tree
# This method calls itself recursively to traverse the entire tree.
def checkTreeNode(node, localIndices):

    if node.isLeaf:

        if node.size() != node.finiteElement().size():
            raise ValueError("Size of leaf node and finite element are different.")

        for i in range(node.size()):

            if node.localIndex(i) >= len(localIndices):
                raise ValueError("Local index exceeds localView size.")

            # Log that we have seen this local index
            localIndices[node.localIndex(i)] += 1

    elif node.isPower or node.isComposite:

        # TODO: Test features of non-leaf nodes

        # Recursively test the children
        for i in range(node.degree()):
            checkTreeNode(node.child(i), localIndices)

    else:
        raise NotImplementedError("Found an unsupported node type")


def checkLocalView(basis, localView, checkLocalIndexSetCompleteness):

    if (localView.size() > localView.maxSize()):
        raise ValueError("localView.size() is " + str(localView.size()) + " but localView.maxSize() is " + str(localView.maxSize()))

    tree1 = localView.tree()
    tree2 = localView.tree()
    if (tree1 != tree2):
        raise ValueError("Multiple calls to localView.tree() do not hand out the same tree.")

    if (tree1.isLeaf and tree1.isPower):
        raise ValueError("Tree claims to be both 'leaf' and 'power'.")

    if (tree1.isLeaf and tree1.isComposite):
        raise ValueError("Tree claims to be both 'leaf' and 'composite'.")

    # Recursively test all nodes in the tree
    localIndices = array('I', [0] * localView.size())

    checkTreeNode(tree1, localIndices)

    # Check that each local index appeared exactly once.
    if checkLocalIndexSetCompleteness:
        for localIndex in localIndices:

            if localIndex < 1:
                raise ValueError("Local index " + str(localIndex) + " did not appear.")
            elif localIndex > 1:
                raise ValueError("Local index " + str(localIndex) + " appeared multiple times.")


# Perform tests that don't modify the basis
def checkConstBasis(basis, *, checkLocalIndexSetCompleteness):

    # Perform all local tests
    localView = basis.localView()
    gridView = basis.gridView

    for element in gridView.elements:
        localView.bind(element)

        # Check the 'element' method
        boundElement = localView.element()
        if boundElement != element:
            raise ValueError("LocalView object does not yield the element it is bound to.")

        checkLocalView(basis,localView, checkLocalIndexSetCompleteness)

    # Perform global index tests
    checkBasisIndices(basis)

    # Perform continuity check
    # TODO


def checkBasis(basis, *, checkLocalIndexSetCompleteness=True):

    # Check the basis
    checkConstBasis(basis, checkLocalIndexSetCompleteness=checkLocalIndexSetCompleteness)

    # Check whether the basis can be copied
    copiedBasis = basis
    checkConstBasis(copiedBasis, checkLocalIndexSetCompleteness=checkLocalIndexSetCompleteness)

    # Can the 'update' method be called?
    gridView = basis.gridView
    # TODO: 'update' method not exposed by the interface yet!
    #basis.update(gridView)