Skip to content

Commit

Permalink
Merge branch 'master' into pbrubeck/h1curl
Browse files Browse the repository at this point in the history
  • Loading branch information
@pbrubeck
pbrubeck authored Nov 19, 2024
2 parents aca9c84 + c6369c7 commit ba09d58
Show file tree
Hide file tree
Showing 4 changed files with 157 additions and 62 deletions.
20 changes: 16 additions & 4 deletions FIAT/kong_mulder_veldhuizen.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# Copyright (C) 2020 Robert C. Kirby (Baylor University)
#
# contributions by Keith Roberts (University of São Paulo)
# and Alexandre Olender (University of São Paulo)
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
Expand Down Expand Up @@ -77,6 +78,14 @@ def _get_entity_ids(ref_el, degree):
1: dict((i, etop[i]) for i in range(3)),
2: {0: [i for i in range(15, 30)]},
}
elif degree == 6:
if sd == 2:
etop = [[6, 12, 3, 13, 7], [9, 15, 4, 14, 8], [10, 16, 5, 17, 11]]
entity_ids = {
0: dict((i, [i]) for i in range(3)),
1: dict((i, etop[i]) for i in range(3)),
2: {0: [i for i in range(18, 39)]},
}
return entity_ids


Expand All @@ -89,7 +98,7 @@ def bump(T, deg):
if sd == 2:
if deg < 5:
return (1, 1)
elif deg == 5:
elif deg == 5 or deg == 6:
return (2, 2)
else:
raise ValueError("Degree not supported")
Expand Down Expand Up @@ -160,13 +169,16 @@ class KongMulderVeldhuizen(finite_element.CiarletElement):
NEW HIGHER-ORDER MASS-LUMPED TETRAHEDRAL ELEMENTS
S. GEEVERS, W.A. MULDER, AND J.J.W. VAN DER VEGT
"""
More Continuous Mass-Lumped Triangular Finite Elements
W. A. MULDER
"""

def __init__(self, ref_el, degree):
if ref_el != TRIANGLE and ref_el != TETRAHEDRON:
raise ValueError("KMV is only valid for triangles and tetrahedrals")
if degree > 5 and ref_el == TRIANGLE:
raise NotImplementedError("Only P < 6 for triangles are implemented.")
if degree > 6 and ref_el == TRIANGLE:
raise NotImplementedError("Only P < 7 for triangles are implemented.")
if degree > 3 and ref_el == TETRAHEDRON:
raise NotImplementedError("Only P < 4 for tetrahedrals are implemented.")
S = KongMulderVeldhuizenSpace(ref_el, degree)
Expand Down
186 changes: 131 additions & 55 deletions FIAT/quadrature_schemes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,8 +253,10 @@ def _kmv_lump_scheme(ref_el, degree):
elif degree == 4:
if sd == 2:
alpha = 0.2113248654051871 # 0.2113248654051871
beta1 = 0.4247639617258106 # 0.4247639617258106
beta2 = 0.130791593829745 # 0.130791593829745
betas = [
0.4247639617258106,
0.130791593829745
]
x = list(ref_el.vertices)
for e in range(3):
x.extend(ref_el.make_points(1, e, 2)) # edge midpoints
Expand All @@ -268,68 +270,56 @@ def _kmv_lump_scheme(ref_el, degree):
(1 - alpha, 0.0),
] # edge points
)
x.extend(
[(beta1, beta1), (1 - 2 * beta1, beta1), (beta1, 1 - 2 * beta1)]
) # points in center of cell
x.extend(
[(beta2, beta2), (1 - 2 * beta2, beta2), (beta2, 1 - 2 * beta2)]
) # points in center of cell
for beta in betas:
x.extend(
[(beta, beta), (1 - 2 * beta, beta), (beta, 1 - 2 * beta)]
) # points in center of cell
w = numpy.arange(18, dtype=numpy.float64)
w[0:3] = 0.003174603174603175 # chk
w[3:6] = 0.0126984126984127 # chk 0.0126984126984127
w[6:12] = 0.01071428571428571 # chk 0.01071428571428571
w[12:15] = 0.07878121446939182 # chk 0.07878121446939182
w[15:18] = 0.05058386489568756 # chk 0.05058386489568756
w[0:3] = 0.003174603174603175
w[3:6] = 0.0126984126984127
w[6:12] = 0.01071428571428571
w[12:15] = 0.07878121446939182
w[15:18] = 0.05058386489568756
else:
raise ValueError("Dimension not supported")

elif degree == 5:
if sd == 2:
alpha1 = 0.3632980741536860e-00
alpha2 = 0.1322645816327140e-00
beta1 = 0.4578368380791611e-00
beta2 = 0.2568591072619591e-00
beta3 = 0.5752768441141011e-01
gamma1 = 0.7819258362551702e-01
delta1 = 0.2210012187598900e-00
alphas = [
0.3632980741536860e-00,
0.1322645816327140e-00
]
betas = [
0.4578368380791611e-00,
0.2568591072619591e-00,
0.5752768441141011e-01
]
gamma = 0.7819258362551702e-01
delta = 0.2210012187598900e-00
x = list(ref_el.vertices)
for alpha in alphas:
x.extend(
[
(1 - alpha, alpha),
(alpha, 1 - alpha),
(0.0, 1 - alpha),
(0.0, alpha),
(alpha, 0.0),
(1 - alpha, 0.0),
] # edge points
)
for beta in betas:
x.extend(
[(beta, beta), (1 - 2 * beta, beta), (beta, 1 - 2 * beta)]
) # points in center of cell
x.extend(
[
(1 - alpha1, alpha1),
(alpha1, 1 - alpha1),
(0.0, 1 - alpha1),
(0.0, alpha1),
(alpha1, 0.0),
(1 - alpha1, 0.0),
] # edge points
)
x.extend(
[
(1 - alpha2, alpha2),
(alpha2, 1 - alpha2),
(0.0, 1 - alpha2),
(0.0, alpha2),
(alpha2, 0.0),
(1 - alpha2, 0.0),
] # edge points
)
x.extend(
[(beta1, beta1), (1 - 2 * beta1, beta1), (beta1, 1 - 2 * beta1)]
) # points in center of cell
x.extend(
[(beta2, beta2), (1 - 2 * beta2, beta2), (beta2, 1 - 2 * beta2)]
) # points in center of cell
x.extend(
[(beta3, beta3), (1 - 2 * beta3, beta3), (beta3, 1 - 2 * beta3)]
) # points in center of cell
x.extend(
[
(gamma1, delta1),
(1 - gamma1 - delta1, delta1),
(gamma1, 1 - gamma1 - delta1),
(delta1, gamma1),
(1 - gamma1 - delta1, gamma1),
(delta1, 1 - gamma1 - delta1),
(gamma, delta),
(1 - gamma - delta, delta),
(gamma, 1 - gamma - delta),
(delta, gamma),
(1 - gamma - delta, gamma),
(delta, 1 - gamma - delta),
] # edge points
)
w = numpy.arange(30, dtype=numpy.float64)
Expand All @@ -343,6 +333,92 @@ def _kmv_lump_scheme(ref_el, degree):
else:
raise ValueError("Dimension not supported")

elif degree == 6:
if sd == 2:
x = list(ref_el.vertices)
episilon = 5.00000000000000e-1
alphas = [
8.29411811106452e-2,
2.68649695592714e-1,
]
betas = [
4.68059729056814e-1,
7.93088545089875e-2,
3.92931636618867e-1,
]
gammas = [
2.48172758709406e-1,
1.56582066033687e-1,
]
deltas = [
6.99812197147049e-1,
2.43089592364562e-1,
]
weights = [
5.35113520281665e-4,
4.29435346026293e-3,
3.02990950926060e-3,
3.16396316646563e-3,
2.43035184285235e-2,
1.66312091329395e-2,
3.42178857644876e-2,
1.73480160090330e-2,
1.98004044953264e-2,
]

x.extend(
[
(episilon, episilon),
(0.0, episilon),
(episilon, 0.0),
]
) # edge midpoints (class 2 points)
for alpha in alphas:
x.extend(
[
(1 - alpha, alpha),
(alpha, 1 - alpha),
(0.0, 1 - alpha),
(0.0, alpha),
(alpha, 0.0),
(1 - alpha, 0.0),
] # edge points (sets of class 3 points)
)

for beta in betas:
x.extend(
[
(beta, beta),
(1 - 2 * beta, beta),
(beta, 1 - 2 * beta)
] # interior points on bisector (sets of class 5 points)
)

for gamma, delta in zip(gammas, deltas):
x.extend(
[
(gamma, delta),
(1 - gamma - delta, delta),
(gamma, 1 - gamma - delta),
(delta, gamma),
(1 - gamma - delta, gamma),
(delta, 1 - gamma - delta),
] # interior points (sets of class 6 points)
)

w = numpy.arange(39, dtype=numpy.float64)
w[0:3] = weights[0] # class 1 points (vertices)
w[3:6] = weights[1] # class 2 points
w[6:12] = weights[2] # class 3 points
w[12:18] = weights[3]
w[18:21] = weights[4] # class 5 points
w[21:24] = weights[5]
w[24:27] = weights[6]
w[27:33] = weights[7] # class 6 points
w[33:39] = weights[8]

else:
raise ValueError("Dimension not supported")
# Return scheme
return QuadratureRule(T, x, w)

Expand Down
7 changes: 7 additions & 0 deletions test/unit/test_fiat.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@
from FIAT.bubble import Bubble
from FIAT.enriched import EnrichedElement # noqa: F401
from FIAT.nodal_enriched import NodalEnrichedElement
from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen # noqa: F401

P = Point()
I = UFCInterval() # noqa: E741
Expand Down Expand Up @@ -339,6 +340,12 @@ def __init__(self, a, b):
"HuZhang(T, 4)",
"HuZhang(T, 3, 'point')",
"HuZhang(T, 4, 'point')",
"KongMulderVeldhuizen(T,1)",
"KongMulderVeldhuizen(T,2)",
"KongMulderVeldhuizen(T,3)",
"KongMulderVeldhuizen(T,4)",
"KongMulderVeldhuizen(T,5)",
"KongMulderVeldhuizen(T,6)",

# Macroelements
"Lagrange(T, 1, 'iso')",
Expand Down
6 changes: 3 additions & 3 deletions test/unit/test_kong_mulder_veldhuizen.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ def test_kmv_quad_tet_schemes(p_d): # noqa: W503
)


@pytest.mark.parametrize("p_d", [(1, 1), (2, 3), (3, 5), (4, 7), (5, 9)])
@pytest.mark.parametrize("p_d", [(1, 1), (2, 3), (3, 5), (4, 7), (5, 9), (6, 11)])
def test_kmv_quad_tri_schemes(p_d):
fct = math.factorial
p, d = p_d
Expand All @@ -44,7 +44,7 @@ def test_kmv_quad_tri_schemes(p_d):

@pytest.mark.parametrize(
"element_degree",
[(KMV(T, 1), 1), (KMV(T, 2), 2), (KMV(T, 3), 3), (KMV(T, 4), 4), (KMV(T, 5), 5)],
[(KMV(T, 1), 1), (KMV(T, 2), 2), (KMV(T, 3), 3), (KMV(T, 4), 4), (KMV(T, 5), 5), (KMV(T, 6), 6)],
)
def test_Kronecker_property_tris(element_degree):
"""
Expand Down Expand Up @@ -101,7 +101,7 @@ def test_edge_degree(degree):

@pytest.mark.parametrize(
"element_degree",
[(KMV(T, 1), 1), (KMV(T, 2), 2), (KMV(T, 3), 3), (KMV(T, 4), 4), (KMV(T, 5), 5)],
[(KMV(T, 1), 1), (KMV(T, 2), 2), (KMV(T, 3), 3), (KMV(T, 4), 4), (KMV(T, 5), 5), (KMV(T, 6), 6)],
)
def test_interpolate_monomials_tris(element_degree):
element, degree = element_degree
Expand Down

0 comments on commit ba09d58

Please sign in to comment.