fuse/test/test_dofs.py at d19632043bf8ce22123dfd1257abcbda96fd9b34 · firedrakeproject/fuse · GitHub

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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
from fuse import *
from test_convert_to_fiat import create_cg1, create_dg1, construct_cg3, construct_rt, construct_nd
from test_orientations import construct_nd2

import sympy as sp
import numpy as np


def test_permute_dg1():
    cell = Point(1, [Point(0), Point(0)], vertex_num=2)

    dg1 = create_dg1(cell)

    for dof in dg1.generate():
        print(dof)

    for g in dg1.cell.group.members():
        print("g", g)
        for dof in dg1.generate():
            print(dof, "->", dof(g))


def test_permute_cg1():
    cell = Point(1, [Point(0), Point(0)], vertex_num=2)

    cg1 = create_cg1(cell)

    for dof in cg1.generate():
        print(dof)

    for g in cg1.cell.group.members():
        print("g", g)
        for dof in cg1.generate():
            print(dof, "->", dof(g))


def test_permute_cg3():
    cell = polygon(3)

    cg3 = construct_cg3(cell)

    for dof in cg3.generate():
        print(dof)

    for g in cg3.cell.group.members():
        print(g)
        for dof in cg3.generate():
            print(dof, "->", dof(g))


def test_permute_rt():
    cell = polygon(3)

    rt = construct_rt(cell)
    x = sp.Symbol("x")
    y = sp.Symbol("y")
    func = FuseFunction(sp.Matrix([x, -1/3 + 2*y]), symbols=(x, y))

    for dof in rt.generate():
        print(dof)

    for g in rt.cell.group.members():
        print(g.numeric_rep())
        for dof in rt.generate():
            print(dof, "->", dof(g), "eval, ", dof(g).eval(func))


def test_permute_nd():
    cell = polygon(3)

    nd = construct_nd(cell)
    x = sp.Symbol("x")
    y = sp.Symbol("y")
    func = FuseFunction(sp.Matrix([x, -1/3 + 2*y]), symbols=(x, y))

    # for dof in nd.generate():
    #     print(dof)

    for g in nd.cell.group.members():
        print("g:", g, g.numeric_rep())
        for dof in nd.generate():
            print(dof(g).convert_to_fiat(cell.to_fiat(), 0).pt_dict)
            print(dof, "->", dof(g), "eval, ", dof(g).eval(func))


def test_permute_nd2():
    cell = polygon(3)

    nd = construct_nd2(cell)
    x = sp.Symbol("x")
    y = sp.Symbol("y")
    func = FuseFunction(sp.Matrix([x, -1/3 + 2*y]), symbols=(x, y))

    # for dof in nd.generate():
    #     print(dof)

    for g in nd.cell.group.members():
        print("g:", g.numeric_rep())
        if g.numeric_rep() == 2:
            for i, dof in enumerate(nd.generate()):
                if 0 < i < 2:
                    print(dof(g).convert_to_fiat(cell.to_fiat(), 2, (2,)).pt_dict)
                    print(dof, "->", dof(g), "eval, ", dof(g).eval(func))


def test_permute_nd_old():
    cell = polygon(3)

    nd = construct_nd(cell)
    x = sp.Symbol("x")
    y = sp.Symbol("y")
    # func = FuseFunction(sp.Matrix([x, -1/3 + 2*y]), symbols=(x, y))

    # phi_0 = FuseFunction(sp.Matrix([-0.333333333333333*y - 0.192450089729875, 0.333333333333333*x + 0.333333333333333]), symbols=(x, y))
    # phi_1 = FuseFunction(sp.Matrix([0.333333333333333*y + 0.192450089729875, 0.333333333333333 - 0.333333333333333*x]), symbols=(x, y))

    # # original dofs
    phi_2 = FuseFunction(sp.Matrix([1/3 - (np.sqrt(3)/6)*y, (np.sqrt(3)/6)*x]), symbols=(x, y))
    phi_0 = FuseFunction(sp.Matrix([-1/6 - (np.sqrt(3)/6)*y, (-np.sqrt(3)/6) + (np.sqrt(3)/6)*x]), symbols=(x, y))
    phi_1 = FuseFunction(sp.Matrix([-1/6 - (np.sqrt(3)/6)*y,
                                    (np.sqrt(3)/6) + (np.sqrt(3)/6)*x]), symbols=(x, y))

    for g in nd.cell.group.members():
        if g.numeric_rep() == 0 or g.numeric_rep() == 1:
            print(g)
            for dof in nd.generate():
                print(dof, "->", dof(g), dof(g).convert_to_fiat(cell.to_fiat(), 1).pt_dict)
                print(dof, "->", dof(g), "eval p2 ", dof(g).eval(phi_2), "eval p0 ", dof(g).eval(phi_0), "eval p1 ", dof(g).eval(phi_1))

    # reflected dofs
    phi_2 = FuseFunction(sp.Matrix([0.288675134594813*y - 0.333333333333333, -0.288675134594813*x]), symbols=(x, y))
    phi_0 = FuseFunction(sp.Matrix([0.288675134594813*y + 0.166666666666667, -0.288675134594813*x - 0.288675134594813]), symbols=(x, y))
    phi_1 = FuseFunction(sp.Matrix([0.288675134594813*y + 0.166666666666667, 0.288675134594813 - 0.288675134594813*x]), symbols=(x, y))
    reflect = nd.cell.group.get_member([0, 1, 2])
    print(nd.cell.permute_entities(reflect, 1))
    reflect = nd.cell.group.get_member([2, 0, 1])
    print(nd.cell.permute_entities(reflect, 1))
    # print(reflect)
    print(nd.cell.get_topology())
    # nd.cell.plot(filename="test_perms.png")
    for g in nd.cell.group.members():
        if g.numeric_rep() == 0 or g.numeric_rep() == 1:
            print(g)
            for dof in nd.generate():
                print(dof, "->", dof(g), dof(g).convert_to_fiat(cell.to_fiat(), 1).pt_dict)
                print(dof, "->", dof(g), "eval p2 ", dof(g).eval(phi_2), "eval p0 ", dof(g).eval(phi_0), "eval p1 ", dof(g).eval(phi_1))
    #     # print(dof.convert_to_fiat(cell.to_fiat(), 1)(lambda x: np.array([1/3 - (np.sqrt(3)/6)*x[1], (np.sqrt(3)/6)*x[0]])))

    # for g in nd.cell.group.members():
    #     print(g)
    #     print(nd.cell.permute_entities(g, 0))
    #     print(nd.cell.permute_entities(g, 1))


def test_permute_nodes():
    cell = polygon(3)
    cg1 = create_cg1(cell)
    degree = cg1.spaces[0].degree()
    nodes = [d.convert_to_fiat(cell.to_fiat(), degree) for d in cg1.generate()]
    print([n.pt_dict for n in nodes])
    for g in cg1.cell.group.members():
        print(g, [d(g).convert_to_fiat(cell.to_fiat(), degree).pt_dict for d in cg1.generate()])


# def test_edge_parametrisation():
#     tri = polygon(3)
#     for i in tri.d_entities(1):
#         print(tri.generate_facet_parameterisation(i.id))
#     from fuse.dof import ParameterisationKernel
#
#     deg = 2
#     edge = tri.edges()[0]
#     x = sp.Symbol("x")
#     y = sp.Symbol("y")
#
#     xs = [DOF(L2Pairing(), ParameterisationKernel())]
#
#     dofs = DOFGenerator(xs, S2, S2)
#     int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs)
#
#     xs = [DOF(L2Pairing(), ComponentKernel((0,))),
#           DOF(L2Pairing(), ComponentKernel((1,)))]
#     center_dofs = DOFGenerator(xs, S1, S3)
#     xs = [immerse(tri, int_ned1, TrHCurl)]
#     tri_dofs = DOFGenerator(xs, C3, S1)
#
#     vec_Pk = PolynomialSpace(deg - 1, set_shape=True)
#     Pk = PolynomialSpace(deg - 1)
#     M = sp.Matrix([[y, -x]])
#     nd = vec_Pk + (Pk.restrict(deg-2, deg-1))*M
#
#     ned = ElementTriple(tri, (nd, CellHCurl, C0), [tri_dofs, center_dofs])
#     for n in ned.generate():
#         print(n)
#
#     from test_orientations import construct_nd2_for_fiat
#
#     ned_fiat = construct_nd2_for_fiat(tri)
#
#     print("fiat")
#     for n in ned_fiat.generate():
#         print(n)


def test_generate_quadrature():
    cell = polygon(3)
    degree = 2
    # elem = create_dg1(cell)
    # elem = create_cg1(cell)
    # elem = construct_nd(cell)
    elem = construct_nd2(cell)
    # elem = construct_nd2_for_fiat(cell)
    from FIAT.nedelec import Nedelec
    fiat_elem = Nedelec(cell.to_fiat(), degree)
    # from FIAT.lagrange import Lagrange
    # fiat_elem = Lagrange(cell.to_fiat(), degree)
    degree = elem.spaces[0].degree()
    print(degree)
    for d in fiat_elem.dual_basis():
        print("fiat", d.pt_dict)
    print()
    for d in elem.generate():
        print("fuse", d.to_quadrature(degree, value_shape=(2,)))

    elem.to_fiat()