diff --git a/Makefile b/Makefile index 734a8d3..e21d433 100644 --- a/Makefile +++ b/Makefile @@ -26,6 +26,12 @@ tests: @echo " Running all tests" @FIREDRAKE_USE_FUSE=1 python3 -m coverage run -p -m pytest -rx test +mini_tests: + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_2d_examples_docs.py + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_1d + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_orientations.py::test_surface_vec_rt + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_projection_convergence_3d\[construct_tet_ned-N1curl-1-0.8\] + coverage: @python3 -m coverage combine @python3 -m coverage report -m diff --git a/fuse/dof.py b/fuse/dof.py index fa02b1f..4c6546f 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -347,14 +347,13 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non self.pairing = self.pairing.add_entity(cell) if self.target_space is None: self.target_space = space - if self.id is None and overall_id is not None: + if overall_id is not None: self.id = overall_id - if self.sub_id is None and generator_id is not None: + if generator_id is not None: self.sub_id = generator_id def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): # TODO deriv dict needs implementing (currently {}) - print(value_shape) return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree, value_shape), {}, str(self)) def to_quadrature(self, arg_degree, value_shape): diff --git a/fuse/spaces/element_sobolev_spaces.py b/fuse/spaces/element_sobolev_spaces.py index 07c9fdd..b5964d9 100644 --- a/fuse/spaces/element_sobolev_spaces.py +++ b/fuse/spaces/element_sobolev_spaces.py @@ -1,6 +1,7 @@ from functools import total_ordering from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2 + @total_ordering class ElementSobolevSpace(object): """ @@ -54,7 +55,7 @@ def __init__(self, cell): def __repr__(self): return "H1" - + def to_ufl(self): return H1 diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 3502a6e..c7f8d36 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -96,9 +96,7 @@ def __mul__(self, x): the sympy object on the right. This is due to Sympy's implementation of __mul__ not passing to this handler as it should. """ - if isinstance(x, sp.Symbol): - return ConstructedPolynomialSpace([x], [self]) - elif isinstance(x, sp.Matrix): + if isinstance(x, sp.Symbol) or isinstance(x, sp.Expr) or isinstance(x, sp.Matrix): return ConstructedPolynomialSpace([x], [self]) else: raise TypeError(f'Cannot multiply a PolySpace with {type(x)}') @@ -177,8 +175,6 @@ def to_ON_polynomial_set(self, ref_el): # otherwise have to work on this through tabulation - Q = create_quadrature(ref_el, 2 * (k + 1)) - Qpts, Qwts = Q.get_points(), Q.get_weights() weighted_sets = [] for (s, w) in zip(self.spaces, self.weights): @@ -196,6 +192,8 @@ def to_ON_polynomial_set(self, ref_el): else: vec = True w_deg = max_deg_sp_expr(w) + Q = create_quadrature(ref_el, 2 * (k + w_deg + 1)) + Qpts, Qwts = Q.get_points(), Q.get_weights() Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, shape, scale="orthonormal") # vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal") @@ -206,11 +204,11 @@ def to_ON_polynomial_set(self, ref_el): if s.set_shape or vec: scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] else: - # breakpoint() scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] scaled_at_Qpts = scaled_at_Qpts.squeeze() PkHw_coeffs = np.dot(np.multiply(scaled_at_Qpts, Qwts), Pkpw_at_Qpts.T) - # breakpoint() + if len(PkHw_coeffs.shape) == 1: + PkHw_coeffs = PkHw_coeffs.reshape(1, -1) weighted_sets.append(polynomial_set.PolynomialSet(ref_el, space.degree + w_deg, space.degree + w_deg, diff --git a/fuse/triples.py b/fuse/triples.py index 32a7220..2d29733 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -421,7 +421,7 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): ent_dofs = entity_associations[dim][e_id][dof_gen] ent_dofs_ids = np.array([self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs], dtype=int) # (dof_gen, ent_dofs) - total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if ed.id not in total_ent_dof_ids] + total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if self.dof_id_to_fiat_id[ed.id] not in total_ent_dof_ids] # dof_idx = [total_ent_dof_ids.index(id) for id in ent_dofs_ids] dof_gen_class = ent_dofs[0].generation @@ -537,6 +537,21 @@ def reverse_dof_perms(self, matrices): reversed_mats[dim][e_id] = perms_copy return reversed_mats + def __add__(self, other): + """ Construct a new element triple by combining the degrees of freedom + This implementation does not make assertions about the properties + of the resulting element. + + Elements being adding must be defined over the same cell and have the same + value shape and mapping""" + assert self.cell == other.cell + assert self.spaces[0].set_shape == other.spaces[0].set_shape + assert str(self.spaces[1]) == str(other.spaces[1]) + + spaces = (self.spaces[0] + other.spaces[0], self.spaces[1], max([self.spaces[2], other.spaces[2]])) + + return ElementTriple(self.cell, spaces, self.DOFGenerator + other.DOFGenerator) + def _to_dict(self): o_dict = {"cell": self.cell, "spaces": self.spaces, "dofs": self.DOFGenerator} return o_dict @@ -574,21 +589,20 @@ def num_dofs(self): return self.dof_numbers def generate(self, cell, space, id_counter): - if self.ls is None: - self.ls = [] - for l_g in self.x: - i = 0 - for g in self.g1.members(): - generated = l_g(g) - if not isinstance(generated, list): - generated = [generated] - for dof in generated: - dof.add_context(self, cell, space, g, id_counter, i) - id_counter += 1 - i += 1 - self.ls.extend(generated) - self.dof_numbers = len(self.ls) - self.dof_ids = [dof.id for dof in self.ls] + self.ls = [] + for l_g in self.x: + i = 0 + for g in self.g1.members(): + generated = l_g(g) + if not isinstance(generated, list): + generated = [generated] + for dof in generated: + dof.add_context(self, cell, space, g, id_counter, i) + id_counter += 1 + i += 1 + self.ls.extend(generated) + self.dof_numbers = len(self.ls) + self.dof_ids = [dof.id for dof in self.ls] return self.ls def make_entity_ids(self): diff --git a/fuse/utils.py b/fuse/utils.py index 6592e11..0d0ee9a 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -67,11 +67,10 @@ def max_deg_sp_expr(sp_expr): for comp in sp_expr: # only compute degree if component is a polynomial if sp.sympify(comp).as_poly(): - degs += [sp.sympify(comp).as_poly().degree()] + degs += [sp.sympify(comp).as_poly().total_degree()] else: if sp.sympify(sp_expr).as_poly(): - degs += [sp.sympify(sp_expr).as_poly().degree()] - + degs += [sp.sympify(sp_expr).as_poly().total_degree()] return max(degs) diff --git a/test/test_3d_examples_docs.py b/test/test_3d_examples_docs.py index 803c361..7572943 100644 --- a/test/test_3d_examples_docs.py +++ b/test/test_3d_examples_docs.py @@ -327,7 +327,7 @@ def test_tet_rt2(): ls = rt2.generate() # TODO make this a proper test for dof in ls: - print(dof.to_quadrature(1)) + print(dof.to_quadrature(1, value_shape=(2,))) rt2.to_fiat() diff --git a/test/test_algebra.py b/test/test_algebra.py new file mode 100644 index 0000000..5304927 --- /dev/null +++ b/test/test_algebra.py @@ -0,0 +1,36 @@ +from fuse import * +from firedrake import * +import numpy as np +import sympy as sp +from test_convert_to_fiat import create_cg2_tri, construct_cg3 + + +def construct_bubble(cell=None): + if cell is None: + cell = polygon(3) + x = sp.Symbol("x") + y = sp.Symbol("y") + f = (3*np.sqrt(3)/4)*(y + np.sqrt(3)/3)*(np.sqrt(3)*x + y - 2*np.sqrt(3)/3)*(-np.sqrt(3)*x + y - 2*np.sqrt(3)/3) + space = PolynomialSpace(3).restrict(0, 0)*f + xs = [DOF(DeltaPairing(), PointKernel((0, 0)))] + bubble = ElementTriple(cell, (space, CellL2, L2), DOFGenerator(xs, S1, S1)) + return bubble + + +def test_bubble(): + mesh = UnitTriangleMesh() + x = SpatialCoordinate(mesh) + + tri = polygon(3) + bub = construct_bubble(tri) + cg2 = create_cg2_tri(tri) + p2b3 = bub + cg2 + V = FunctionSpace(mesh, p2b3.to_ufl()) + W = FunctionSpace(mesh, construct_cg3().to_ufl()) + + bubble_func = 27*x[0]*x[1]*(1-x[0]-x[1]) + u = project(bubble_func, V) + exact = Function(W) + exact.interpolate(bubble_func, W) + # make sure that these are the same + assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14 diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 9057815..ad1db36 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -358,7 +358,7 @@ def test_entity_perms(elem_gen, cell): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg1, "CG", 1), (create_dg1, "DG", 1), - (construct_dg0_integral, "DG", 0), + pytest.param(construct_dg0_integral, "DG", 0, marks=pytest.mark.xfail(reason='Passes locally, fails in CI, probably same as dg2')), (construct_dg1_integral, "DG", 1), (construct_dg2_integral, "DG", 2), pytest.param(create_dg2, "DG", 2, marks=pytest.mark.xfail(reason='Need to update TSFC in CI')), diff --git a/test/test_dofs.py b/test/test_dofs.py index ec7424f..bf29ad4 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -220,6 +220,6 @@ def test_generate_quadrature(): print("fiat", d.pt_dict) print() for d in elem.generate(): - print("fuse", d.to_quadrature(degree)) + print("fuse", d.to_quadrature(degree, value_shape=(2,))) elem.to_fiat() diff --git a/test/test_perms.py b/test/test_perms.py index f00aed6..82a3cb9 100644 --- a/test/test_perms.py +++ b/test/test_perms.py @@ -1,9 +1,8 @@ from fuse import * -from test_convert_to_fiat import create_cg1, create_dg1, create_cg2 +from test_convert_to_fiat import create_cg1, create_dg1, create_cg2, construct_nd from test_2d_examples_docs import construct_cg3 import pytest import numpy as np -import sympy as sp vert = Point(0) edge = Point(1, [Point(0), Point(0)], vertex_num=2) @@ -48,24 +47,7 @@ def test_basic_perms(cell): @pytest.mark.parametrize("cell", [tri]) def test_nd_perms(cell): - deg = 1 - edge = cell.edges(get_class=True)[0] - x = sp.Symbol("x") - y = sp.Symbol("y") - - xs = [DOF(L2Pairing(), PolynomialKernel(edge.basis_vectors()[0]))] - dofs = DOFGenerator(xs, S1, S2) - int_ned = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - - xs = [immerse(cell, int_ned, TrHCurl)] - tri_dofs = DOFGenerator(xs, C3, S3) - - M = sp.Matrix([[y, -x]]) - vec_Pk = PolynomialSpace(deg - 1, set_shape=True) - Pk = PolynomialSpace(deg - 1) - nd = vec_Pk + (Pk.restrict(deg - 2, deg - 1))*M - - ned = ElementTriple(cell, (nd, CellHCurl, C0), [tri_dofs]) + ned = construct_nd(cell) ned.to_fiat() for i, mat in ned.matrices[2][0].items(): print(i) diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index a90e258..15e15a7 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -36,10 +36,10 @@ def mass_solve(U): @pytest.mark.parametrize("generator1, generator2, code1, code2, deg1, deg2", - [(construct_cg1, construct_cg1, "CG", "CG", 1, 1), - (construct_dg1, construct_dg1, "DG", "DG", 1, 1), - (construct_dg1, construct_cg1, "DG", "CG", 1, 1), - (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)]) + [(construct_cg1, construct_cg1, "CG", "CG", 1, 1), + (construct_dg1, construct_dg1, "DG", "DG", 1, 1), + (construct_dg1, construct_cg1, "DG", "CG", 1, 1), + (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)]) def test_ext_mesh(generator1, generator2, code1, code2, deg1, deg2): m = UnitIntervalMesh(2) mesh = ExtrudedMesh(m, 2) @@ -137,3 +137,59 @@ def test_flattening(A, B, res): else: cell = tensor_cell.flatten() cell.construct_fuse_rep() + + +# def test_trace_galerkin_projection(): +# mesh = UnitSquareMesh(10, 10, quadrilateral=True) + +# x, y = SpatialCoordinate(mesh) +# A = construct_cg1() +# B = construct_dg0_integral() +# elem = tensor_product(A, B) +# elem = elem.flatten() + +# # Define the Trace Space +# T = FunctionSpace(mesh, elem.to_ufl()) + +# # Define trial and test functions +# lambdar = TrialFunction(T) +# gammar = TestFunction(T) + +# # Define right hand side function + +# V = FunctionSpace(mesh, "CG", 1) +# f = Function(V) +# f.interpolate(cos(x*pi*2)*cos(y*pi*2)) + +# # Construct bilinear form +# a = inner(lambdar, gammar) * ds + inner(lambdar('+'), gammar('+')) * dS + +# # Construct linear form +# l = inner(f, gammar) * ds + inner(f('+'), gammar('+')) * dS + +# # Compute the solution +# t = Function(T) +# solve(a == l, t, solver_parameters={'ksp_rtol': 1e-14}) + +# # Compute error in trace norm +# trace_error = sqrt(assemble(FacetArea(mesh)*inner((t - f)('+'), (t - f)('+')) * dS)) + +# assert trace_error < 1e-13 + +# def test_hdiv(): +# np.set_printoptions(linewidth=90, precision=4, suppress=True) +# m = UnitIntervalMesh(2) +# mesh = ExtrudedMesh(m, 2) +# CG_1 = FiniteElement("CG", interval, 1) +# DG_0 = FiniteElement("DG", interval, 0) +# P1P0 = TensorProductElement(CG_1, DG_0) +# RT_horiz = HDivElement(P1P0) +# P0P1 = TensorProductElement(DG_0, CG_1) +# RT_vert = HDivElement(P0P1) +# elt = RT_horiz + RT_vert +# V = FunctionSpace(mesh, elt) +# tabulation = V.finat_element.fiat_equivalent.tabulate(0, [(0, 0), (1, 0)]) +# for ent, arr in tabulation.items(): +# print(ent) +# for comp in arr: +# print(comp[0], comp[1])