Skip to content

Commit

Permalink
new interior dofs
Browse files Browse the repository at this point in the history
  • Loading branch information
@pbrubeck
pbrubeck committed Oct 27, 2024
1 parent 053895f commit 950671d
Show file tree
Hide file tree
Showing 7 changed files with 64 additions and 50 deletions.
14 changes: 7 additions & 7 deletions FIAT/arnold_winther.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,11 @@
from FIAT.reference_element import TRIANGLE
from FIAT.quadrature_schemes import create_quadrature
from FIAT.functional import (ComponentPointEvaluation,
IntegralMoment,
TensorBidirectionalIntegralMoment,
IntegralMomentOfTensorDivergence,
IntegralLegendreNormalNormalMoment,
IntegralLegendreNormalTangentialMoment)


import numpy


Expand All @@ -28,7 +27,6 @@ def __init__(self, ref_el, degree=2):
raise ValueError("Nonconforming Arnold-Winther elements are only defined for degree 2.")
top = ref_el.get_topology()
sd = ref_el.get_spatial_dimension()
shp = (sd, sd)
entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
nodes = []

Expand All @@ -45,9 +43,10 @@ def __init__(self, ref_el, degree=2):

# internal dofs: constant moments of three unique components
cur = len(nodes)
n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
Q = create_quadrature(ref_el, degree)
phi = numpy.ones(Q.get_weights().shape)
nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp)
phi = numpy.full(Q.get_weights().shape, 1/ref_el.volume())
nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
for i in range(sd) for j in range(i, sd))
entity_ids[2][0].extend(range(cur, len(nodes)))

Expand Down Expand Up @@ -102,10 +101,11 @@ def __init__(self, ref_el, degree=3):
entity_ids[1][entity].extend(range(cur, len(nodes)))

# internal dofs: moments of unique components against P_{k-3}
n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
Q = create_quadrature(ref_el, 2*(degree-1))
P = polynomial_set.ONPolynomialSet(ref_el, degree-3)
P = polynomial_set.ONPolynomialSet(ref_el, degree-3, scale="L2 piola")
phis = P.tabulate(Q.get_points())[(0,)*sd]
nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp)
nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
for phi in phis for i in range(sd) for j in range(i, sd))

# constraint dofs: moments of divergence against P_{k-1} \ P_{k-2}
Expand Down
13 changes: 7 additions & 6 deletions FIAT/expansions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -279,16 +279,19 @@ def __init__(self, ref_el, scale=None, variant=None):
self._dmats_cache = {}
self._cell_node_map_cache = {}

def get_scale(self, cell=0):
def get_scale(self, n, cell=0):
scale = self.scale
sd = self.ref_el.get_spatial_dimension()
if isinstance(scale, str):
sd = self.ref_el.get_spatial_dimension()
vol = self.ref_el.volume_of_subcomplex(sd, cell)
scale = scale.lower()
if scale == "orthonormal":
scale = math.sqrt(1.0 / vol)
elif scale == "l2 piola":
scale = 1.0 / vol
elif n == 0 and sd > 1 and len(self.affine_mappings) == 1:
# return 1 for n=0 to make regression tests pass
scale = 1
return scale

def get_num_members(self, n):
Expand All @@ -310,9 +313,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
ref_pts = numpy.add(numpy.dot(pts, A.T), b).T
Jinv = A if direction is None else numpy.dot(A, direction)[:, None]
sd = self.ref_el.get_spatial_dimension()

# Always return 1 for n=0 to make regression tests pass
scale = 1.0 if n == 0 and len(self.affine_mappings) == 1 else self.get_scale(cell=cell)
scale = self.get_scale(n, cell=cell)
phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
scale, variant=self.variant)
if self.continuity == "C0":
Expand Down Expand Up @@ -549,7 +550,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
Jinv = A[0, 0] if direction is None else numpy.dot(A, direction)
xs = numpy.add(numpy.dot(pts, A.T), b)
results = {}
scale = self.get_scale(cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
scale = self.get_scale(n, cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
for k in range(order+1):
v = numpy.zeros((n + 1, len(xs)), xs.dtype)
if n >= k:
Expand Down
2 changes: 1 addition & 1 deletion FIAT/functional.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -625,7 +625,7 @@ def __init__(self, ref_el, v, w, pt):
super().__init__(ref_el, shp, pt_dict, {}, "PointwiseInnerProductEval")


class TensorBidirectionalMomentInnerProductEvaluation(FrobeniusIntegralMoment):
class TensorBidirectionalIntegralMoment(FrobeniusIntegralMoment):
r"""
This is a functional on symmetric 2-tensor fields. Let u be such a
field, f a function tabulated at points, and v,w be vectors. This implements the evaluation
Expand Down
11 changes: 7 additions & 4 deletions FIAT/hu_zhang.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,12 @@
from FIAT.quadrature_schemes import create_quadrature
from FIAT.functional import (ComponentPointEvaluation,
PointwiseInnerProductEvaluation,
IntegralMoment,
TensorBidirectionalIntegralMoment,
IntegralLegendreNormalNormalMoment,
IntegralLegendreNormalTangentialMoment)

import numpy


class HuZhangDual(dual_set.DualSet):
def __init__(self, ref_el, degree, variant, qdegree):
Expand Down Expand Up @@ -63,10 +65,11 @@ def __init__(self, ref_el, degree, variant, qdegree):

elif variant == "integral":
# Moments of unique components against a basis for P_{k-2}
Q = create_quadrature(ref_el, qdegree + degree-2)
P = polynomial_set.ONPolynomialSet(ref_el, degree-2)
n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
Q = create_quadrature(ref_el, 2*degree-2)
P = polynomial_set.ONPolynomialSet(ref_el, degree-2, scale="L2 piola")
phis = P.tabulate(Q.get_points())[(0,)*sd]
nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp)
nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
for phi in phis for i in range(sd) for j in range(i, sd))

entity_ids[2][0].extend(range(cur, len(nodes)))
Expand Down
24 changes: 11 additions & 13 deletions FIAT/johnson_mercier.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
from FIAT import finite_element, dual_set, macro, polynomial_set
from FIAT.functional import IntegralMoment, FrobeniusIntegralMoment
from FIAT.functional import TensorBidirectionalIntegralMoment
from FIAT.quadrature import FacetQuadratureRule
from FIAT.quadrature_schemes import create_quadrature
import numpy
Expand All @@ -18,31 +18,29 @@ def __init__(self, ref_complex, degree, variant=None):
nodes = []

# Face dofs: bidirectional (nn and nt) Legendre moments
R = numpy.array([[0, 1], [-1, 0]])
dim = sd - 1
R = numpy.array([[0, 1], [-1, 0]])
ref_facet = ref_el.construct_subelement(dim)
Qref = create_quadrature(ref_facet, 2*degree)
P = polynomial_set.ONPolynomialSet(ref_facet, degree)
phis = P.tabulate(Qref.get_points())[(0,) * dim]
for facet in sorted(top[dim]):
for f in sorted(top[dim]):
cur = len(nodes)
Q = FacetQuadratureRule(ref_el, dim, facet, Qref)
thats = ref_el.compute_tangents(dim, facet)
Q = FacetQuadratureRule(ref_el, dim, f, Qref)
thats = ref_el.compute_tangents(dim, f)
nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats)
normal = nhat / Q.jacobian_determinant()

uvecs = (nhat, *thats)
comps = [numpy.outer(normal, uvec) for uvec in uvecs]
nodes.extend(FrobeniusIntegralMoment(ref_el, Q, comp[:, :, None] * phi[None, None, :])
for phi in phis for comp in comps)
entity_ids[dim][facet].extend(range(cur, len(nodes)))
nodes.extend(TensorBidirectionalIntegralMoment(ref_el, normal, comp, Q, phi)
for phi in phis for comp in (nhat, *thats))
entity_ids[dim][f].extend(range(cur, len(nodes)))

cur = len(nodes)
# Interior dofs: moments for each independent component
n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
Q = create_quadrature(ref_complex, 2*degree-1)
P = polynomial_set.ONPolynomialSet(ref_el, degree-1)
P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola")
phis = P.tabulate(Q.get_points())[(0,) * sd]
nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(i, j))
nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
for phi in phis for i in range(sd) for j in range(i, sd))

entity_ids[sd][0].extend(range(cur, len(nodes)))
Expand Down
28 changes: 17 additions & 11 deletions test/unit/test_awc.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
import numpy as np
from FIAT import ufc_simplex, ArnoldWinther, make_quadrature, expansions
from FIAT import ufc_simplex, ArnoldWinther, create_quadrature, expansions


def test_dofs():
Expand All @@ -20,7 +20,7 @@ def test_dofs():
assert np.allclose(vert_vals[3*i+j, :, :, (i+k) % 3], np.zeros((2, 2)))

# check edge moments
Qline = make_quadrature(line, 6)
Qline = create_quadrature(line, 6)

linebfs = expansions.LineExpansionSet(line)
linevals = linebfs.tabulate(1, Qline.pts)
Expand Down Expand Up @@ -73,15 +73,21 @@ def test_dofs():
assert np.allclose(ntmoments[bf, :], np.zeros(2))

# check internal dofs
Q = make_quadrature(T, 6)
ns = list(map(T.compute_scaled_normal, range(3)))
Q = create_quadrature(T, 3)
qpvals = AW.tabulate(0, Q.pts)[(0, 0)]
const_moms = qpvals @ Q.wts
assert np.allclose(const_moms[:21], np.zeros((21, 2, 2)))
assert np.allclose(const_moms[24:], np.zeros((6, 2, 2)))
assert np.allclose(const_moms[21:24, 0, 0], np.asarray([1, 0, 0]))
assert np.allclose(const_moms[21:24, 0, 1], np.asarray([0, 1, 0]))
assert np.allclose(const_moms[21:24, 1, 0], np.asarray([0, 1, 0]))
assert np.allclose(const_moms[21:24, 1, 1], np.asarray([0, 0, 1]))
const_moms = qpvals @ Q.wts / T.volume()
nn_moms = const_moms.copy()
for j in range(2):
for i in range(2):
comp = np.outer(ns[i+1], ns[j+1])
nn_moms[:, i, j] = np.tensordot(const_moms, comp, ((1, 2), (0, 1)))
assert np.allclose(nn_moms[:21], np.zeros((21, 2, 2)))
assert np.allclose(nn_moms[24:], np.zeros((6, 2, 2)))
assert np.allclose(nn_moms[21:24, 0, 0], np.asarray([1, 0, 0]))
assert np.allclose(nn_moms[21:24, 0, 1], np.asarray([0, 1, 0]))
assert np.allclose(nn_moms[21:24, 1, 0], np.asarray([0, 1, 0]))
assert np.allclose(nn_moms[21:24, 1, 1], np.asarray([0, 0, 1]))


def frob(a, b):
Expand All @@ -94,7 +100,7 @@ def test_projection():

AW = ArnoldWinther(T, 3)

Q = make_quadrature(T, 4)
Q = create_quadrature(T, 6)
qpts = np.asarray(Q.pts)
qwts = np.asarray(Q.wts)
nqp = len(Q.wts)
Expand Down
22 changes: 14 additions & 8 deletions test/unit/test_awnc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
def test_dofs():
line = ufc_simplex(1)
T = ufc_simplex(2)
T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)])
T.vertices = ((0.0, 0.0), (1.0, 0.25), (-0.75, 1.1))
AW = ArnoldWintherNC(T, 2)

Qline = make_quadrature(line, 6)
Expand Down Expand Up @@ -61,12 +61,18 @@ def test_dofs():
assert np.allclose(ntmoments[bf, :], np.zeros(2), atol=1.e-7)

# check internal dofs
ns = list(map(T.compute_scaled_normal, range(3)))
Q = make_quadrature(T, 6)
qpvals = AW.tabulate(0, Q.pts)[(0, 0)]
const_moms = qpvals @ Q.wts
assert np.allclose(const_moms[:12], np.zeros((12, 2, 2)))
assert np.allclose(const_moms[15:], np.zeros((3, 2, 2)))
assert np.allclose(const_moms[12:15, 0, 0], np.asarray([1, 0, 0]))
assert np.allclose(const_moms[12:15, 0, 1], np.asarray([0, 1, 0]))
assert np.allclose(const_moms[12:15, 1, 0], np.asarray([0, 1, 0]))
assert np.allclose(const_moms[12:15, 1, 1], np.asarray([0, 0, 1]))
const_moms = qpvals @ Q.wts / T.volume()
nn_moms = const_moms.copy()
for j in range(2):
for i in range(2):
comp = np.outer(ns[i+1], ns[j+1])
nn_moms[:, i, j] = np.tensordot(const_moms, comp, ((1, 2), (0, 1)))
assert np.allclose(nn_moms[:12], np.zeros((12, 2, 2)))
assert np.allclose(nn_moms[15:], np.zeros((3, 2, 2)))
assert np.allclose(nn_moms[12:15, 0, 0], np.asarray([1, 0, 0]))
assert np.allclose(nn_moms[12:15, 0, 1], np.asarray([0, 1, 0]))
assert np.allclose(nn_moms[12:15, 1, 0], np.asarray([0, 1, 0]))
assert np.allclose(nn_moms[12:15, 1, 1], np.asarray([0, 0, 1]))

0 comments on commit 950671d

Please sign in to comment.