diff --git a/.github/workflows/setup_repos.sh b/.github/workflows/setup_repos.sh index f986e1f..ea5d1fe 100644 --- a/.github/workflows/setup_repos.sh +++ b/.github/workflows/setup_repos.sh @@ -16,12 +16,12 @@ git checkout indiamai/integrate_fuse git status python3 -m pip install --break-system-packages -e . -/usr/bin/git config --global --add safe.directory ~ -cd ~ -git clone https://github.com/firedrakeproject/ufl.git -/usr/bin/git config --global --add safe.directory ~/ufl -cd ufl -git fetch -git checkout indiamai/integrate-fuse -git status -python3 -m pip install --break-system-packages -e . \ No newline at end of file +#/usr/bin/git config --global --add safe.directory ~ +#cd ~ +#git clone https://github.com/firedrakeproject/ufl.git +#/usr/bin/git config --global --add safe.directory ~/ufl +#cd ufl +#git fetch +#git checkout indiamai/integrate-fuse +#git status +#python3 -m pip install --break-system-packages -e . diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index a265da9..88b3d0f 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -14,7 +14,8 @@ jobs: # Run on the Github hosted runner runs-on: ubuntu-latest container: - image: firedrakeproject/firedrake-vanilla-default:latest + #image: firedrakeproject/firedrake-vanilla-default:latest + image: firedrakeproject/firedrake-vanilla-default:dev-main # Steps represent a sequence of tasks that will be executed as # part of the jobs steps: @@ -44,17 +45,6 @@ jobs: git checkout indiamai/integrate_fuse git status python3 -m pip install --break-system-packages -e . - - name: Checkout correct UFL branch - run: | - /usr/bin/git config --global --add safe.directory ~ - cd ~ - git clone https://github.com/firedrakeproject/ufl.git - /usr/bin/git config --global --add safe.directory ~/ufl - cd ufl - git fetch - git checkout indiamai/integrate-fuse - git status - python3 -m pip install --break-system-packages -e . - name: Run tests run: | pip list @@ -104,4 +94,4 @@ jobs: message: ${{ env.total }}% minColorRange: 50 maxColorRange: 90 - valColorRange: ${{ env.total }} \ No newline at end of file + valColorRange: ${{ env.total }} diff --git a/Makefile b/Makefile index 969fd02..734a8d3 100644 --- a/Makefile +++ b/Makefile @@ -19,12 +19,12 @@ lint: test_examples: @echo " Running examples" - @python3 -m pytest test/test_2d_examples_docs.py - @python3 -m pytest test/test_3d_examples_docs.py + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_2d_examples_docs.py + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_3d_examples_docs.py tests: @echo " Running all tests" - @python3 -m coverage run -p -m pytest -rx test + @FIREDRAKE_USE_FUSE=1 python3 -m coverage run -p -m pytest -rx test coverage: @python3 -m coverage combine @@ -34,13 +34,13 @@ coverage: test_cells: @echo " Running all cell comparison tests" @firedrake-clean - @python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect0] + @FIREDRAKE_USE_FUSE=1 python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect0] @firedrake-clean - @python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect1] + @FIREDRAKE_USE_FUSE=1 python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect1] clean: @(cd docs/ && make clean) prepush: lint tests @rm .coverage.* - make clean docs \ No newline at end of file + make clean docs diff --git a/conftest.py b/conftest.py index c8c4e00..acf59ca 100644 --- a/conftest.py +++ b/conftest.py @@ -1,9 +1,7 @@ -import pytest - def pytest_addoption(parser): parser.addoption( "--run-cleared", action="store_true", default=False, help="Run tests that require a cleared cache", - ) \ No newline at end of file + ) diff --git a/docs/source/conf.py b/docs/source/conf.py index 154171e..a8781a8 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -87,8 +87,11 @@ # so a file named "default.css" will overwrite the builtin "default.css". html_static_path = ['_static'] # html_css_files = ["additional.css"] + + def setup(app): - app.add_css_file('additional.css') + app.add_css_file('additional.css') + html_theme_options = { 'navbar_links': [ @@ -109,11 +112,11 @@ def setup(app): 'globaltoc_depth': 2, } -html_sidebars = {'api/*': ['localtoc.html'], +html_sidebars = {'api/*': ['localtoc.html'], 'api': ['localtoc.html'], 'about': [], 'install': [], 'index': ['localtoc.html'], 'manual': ['localtoc.html'], 'manual/*': ['localtoc.html'], - '_generated/*': ['custom-sidebar.html'], } \ No newline at end of file + '_generated/*': ['custom-sidebar.html'], } diff --git a/fuse/.triples.py.swp b/fuse/.triples.py.swp new file mode 100644 index 0000000..c6ec9c8 Binary files /dev/null and b/fuse/.triples.py.swp differ diff --git a/fuse/__init__.py b/fuse/__init__.py index c0ffec1..80c5970 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,4 +1,3 @@ - from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, tri_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, PolynomialKernel, ComponentKernel diff --git a/fuse/cells.py b/fuse/cells.py index fb9694a..171b273 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -12,6 +12,7 @@ from sympy.combinatorics.named_groups import SymmetricGroup from fuse.utils import sympy_to_numpy, fold_reduce, numpy_to_str_tuple, orientation_value from FIAT.reference_element import Simplex, TensorProductCell as FiatTensorProductCell, Hypercube +from FIAT.quadrature_schemes import create_quadrature from ufl.cell import Cell, TensorProductCell from functools import cache @@ -94,7 +95,7 @@ def compute_scaled_verts(d, n): :param: n: number of vertices """ if d == 2: - source = np.array([0, 1]) + source = np.array([-np.sqrt(3)/2, -1/2]) rot_coords = [source for i in range(0, n)] rot_mat = np.array([[np.cos((2*np.pi)/n), -np.sin((2*np.pi)/n)], [np.sin((2*np.pi)/n), np.cos((2*np.pi)/n)]]) @@ -544,7 +545,7 @@ def get_node(self, node, return_coords=False): if return_coords: top_level_node = self.d_entities_ids(self.graph_dim())[0] if self.dimension == 0: - return [()] + return () return self.attachment(top_level_node, node)() return self.G.nodes.data("point_class")[node] @@ -623,7 +624,11 @@ def basis_vectors(self, return_coords=True, entity=None): if self.dimension == 0: # return [[] raise ValueError("Dimension 0 entities cannot have Basis Vectors") - top_level_node = self_levels[0][0] + if self.oriented: + # ordered_vertices() handles the orientation so we want to drop the orientation node + top_level_node = self_levels[1][0] + else: + top_level_node = self_levels[0][0] v_0 = vertices[0] if return_coords: v_0_coords = self.attachment(top_level_node, v_0)() @@ -656,6 +661,27 @@ def to_tikz(self, show=True, scale=3): return "\n".join(tikz_commands) return tikz_commands + def generate_facet_parameterisation(self, facet_num): + raise NotImplementedError("Facet Parameterisation can be expressed using polynomials") + # facet = self.d_entities(self.dimension - 1)[facet_num] + facet = self.get_node(facet_num) + facet_dim = facet.dimension + if facet_dim != self.dimension - 1: + raise ValueError(f"Supplied node {facet_num} is not a facet") + if facet_dim > 1: + raise NotImplementedError("Facet parameterisation is not implemented for dimensions greater than 1") + verts = facet.vertices() + v_coords = np.array([self.get_node(v.id, return_coords=True) for v in verts]) + stacked = np.c_[np.ones((self.dimension,)), v_coords[:, 0].reshape(self.dimension, 1)] + b = np.array([0, 1]) + coeffs = np.linalg.solve(stacked, b) + symbol_names = ["x", "y", "z"] + symbols = [1] + [sp.Symbol(symbol_names[d]) for d in range(facet_dim)] + res = 0 + for d in range(facet_dim + 1): + res += coeffs[d] * symbols[d] + return res, symbols[1:] + def plot(self, show=True, plain=False, ax=None, filename=None): """ for now into 2 dimensional space """ if self.dimension == 3: @@ -767,6 +793,12 @@ def attachment(self, source, dst): return lambda *x: fold_reduce(attachments[0], *x) + def quadrature(self, degree): + fiat_el = self.to_fiat() + Q = create_quadrature(fiat_el, degree) + pts, wts = Q.get_points(), Q.get_weights() + return pts, wts + def cell_attachment(self, dst): if not isinstance(dst, int): raise ValueError @@ -775,6 +807,8 @@ def cell_attachment(self, dst): def orient(self, o): """ Orientation node is always labelled with -1 """ + if self.oriented: + o = self.oriented * o oriented_point = copy.deepcopy(self) top_level_node = oriented_point.d_entities_ids( oriented_point.dimension)[0] @@ -838,10 +872,7 @@ def __call__(self, *x): if hasattr(self.attachment, '__iter__'): res = [] for attach_comp in self.attachment: - if len(attach_comp.atoms(sp.Symbol)) == len(x): - res.append(sympy_to_numpy(attach_comp, syms, x)) - else: - res.append(attach_comp.subs({syms[i]: x[i] for i in range(len(x))})) + res.append(sympy_to_numpy(attach_comp, syms, x)) return tuple(res) return sympy_to_numpy(self.attachment, syms, x) return x @@ -885,7 +916,6 @@ def get_spatial_dimension(self): def get_sub_entities(self): self.A.get_sub_entities() self.B.get_sub_entities() - breakpoint() def dimension(self): return tuple(self.A.dimension, self.B.dimension) @@ -1076,7 +1106,7 @@ def __init__(self, cell, name=None): super(CellComplexToUFL, self).__init__(name) def to_fiat(self): - return self.cell_complex.to_fiat(name=self.cellname()) + return self.cell_complex.to_fiat(name=self.cellname) def __repr__(self): return super(CellComplexToUFL, self).__repr__() diff --git a/fuse/dof.py b/fuse/dof.py index 52afcb5..fa3dfa6 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -1,6 +1,6 @@ from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule -from FIAT.functional import PointEvaluation, IntegralMoment, FrobeniusIntegralMoment +from FIAT.functional import PointEvaluation, IntegralMoment, FrobeniusIntegralMoment, Functional from fuse.utils import sympy_to_numpy import numpy as np import sympy as sp @@ -35,11 +35,8 @@ def __call__(self, kernel, v, cell): assert isinstance(kernel, PointKernel) return v(*kernel.pt) - def convert_to_fiat(self, ref_el, dof, interpolant_deg): - pt = dof.eval(FuseFunction(lambda *x: x)) - # pt1 = dof.tabulate([[1]]) - return PointEvaluation(ref_el, pt) - # return PointEvaluation(ref_el, tuple(pt1[0])) + def tabulate(self): + return 1 def add_entity(self, entity): res = DeltaPairing() @@ -51,7 +48,10 @@ def add_entity(self, entity): def permute(self, g): res = DeltaPairing() if self.entity: - res.entity = self.entity.orient(g) + res.entity = self.entity + # res.entity = self.entity.orient(g) + if self.orientation is not None: + g = g * self.orientation res.orientation = g return res @@ -71,30 +71,18 @@ class L2Pairing(Pairing): def __init__(self): super(L2Pairing, self).__init__() - def __call__(self, kernel, v, cell): - # TODO get degree of v - # if cell == self.entity: - # ref_el = self.entity.to_fiat() - # # print("evaluating", kernel, v, "on", self.entity) - # Q = create_quadrature(self.entity.to_fiat(), 5) - # # need quadrature here too - therefore need the information from the triple. - # else: - # ref_el = cell.to_fiat() - # ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] - # entity_ref = ref_el.construct_subelement(self.entity.dim()) - # entity = ref_el.construct_subelement(self.entity.dim(), ent_id, self.orientation) - # Q_ref = create_quadrature(entity, 5) - # Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref, self.orientation) - Q = create_quadrature(self.entity.to_fiat(), 5) - - def kernel_dot(x): - return np.dot(kernel(*x), v(*x)) - - return Q.integrate(kernel_dot) + Qpts, Qwts = self.entity.quadrature(5) + return sum([wt*np.dot(kernel(*pt), v(*pt)) for pt, wt in zip(Qpts, Qwts)]) def tabulate(self): - pass + + bvs = np.array(self.entity.basis_vectors()) + if self.orientation: + new_bvs = np.array(self.entity.orient(self.orientation).basis_vectors()) + basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) + return basis_change + return np.eye(bvs.shape[0]) def add_entity(self, entity): res = L2Pairing() @@ -109,35 +97,15 @@ def permute(self, g): res = L2Pairing() if self.entity: - res.entity = self.entity.orient(g) + res.entity = self.entity + if self.orientation is not None: + g = g * self.orientation res.orientation = g return res - def convert_to_fiat(self, ref_el, dof, interpolant_degree): - total_deg = dof.kernel.degree(interpolant_degree) - ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()] - # entity_ref = ref_el.construct_subelement(self.entity.dim()) - entity = ref_el.construct_subelement(self.entity.dim(), ent_id, self.orientation) - Q_ref = create_quadrature(entity, total_deg) - # pts_ref, wts_ref = Q_ref.get_points(), Q_ref.get_weights() - - # pts, wts, J = map_quadrature(pts_ref, wts_ref, Q_ref.ref_el, entity_ref, jacobian=True) - Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref, self.orientation) - Jdet = Q.jacobian_determinant() - # Jdet = pseudo_determinant(J) - qpts, _ = Q.get_points(), Q.get_weights() - # (np.sqrt(2)) * - f_at_qpts = dof.tabulate(qpts).T / Jdet - comp = None - if hasattr(dof.kernel, "comp"): - # TODO Value shape should be a variable. - comp = dof.kernel.comp - functional = IntegralMoment(ref_el, Q, f_at_qpts, comp=comp, shp=(2,)) - else: - functional = FrobeniusIntegralMoment(ref_el, Q, f_at_qpts) - return functional - def __repr__(self): + if self.orientation is not None: + return "integral_{}({})({{kernel}} * {{fn}}) dx) ".format(str(self.orientation), str(self.entity)) return "integral_{}({{kernel}} * {{fn}}) dx) ".format(str(self.entity)) def dict_id(self): @@ -152,8 +120,13 @@ def _from_dict(obj_dict): class BaseKernel(): def __init__(self): + self.cell = None + self.entity = None self.attachment = False + def add_context(self, cell, entity): + return self + def permute(self, g): raise NotImplementedError("This method should be implemented by the subclass") @@ -185,10 +158,8 @@ def permute(self, g): def __call__(self, *args): return self.pt - def tabulate(self, Qpts, attachment=None): - #if attachment: - # return np.array([attachment(*self.pt) for _ in Qpts]).astype(np.float64) - return np.array([self.pt for _ in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts, basis_change): + return np.array([self.pt for _ in Qpts]).astype(np.float64), np.ones_like(Qwts), [[tuple()] for pt in Qpts] def _to_dict(self): o_dict = {"pt": self.pt} @@ -203,10 +174,11 @@ def _from_dict(obj_dict): class PolynomialKernel(BaseKernel): - def __init__(self, fn, symbols=[]): + def __init__(self, fn, g=None, symbols=[]): if len(symbols) != 0 and not sp.sympify(fn).as_poly(): raise ValueError("Function argument must be able to be interpreted as a sympy polynomial") self.fn = sp.sympify(fn) + self.g = g self.syms = symbols super(PolynomialKernel, self).__init__() @@ -215,25 +187,28 @@ def __repr__(self): def degree(self, interpolant_degree): if len(self.fn.free_symbols) == 0: - return interpolant_degree + return interpolant_degree return self.fn.as_poly().total_degree() * interpolant_degree def permute(self, g): - return self - # new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) - # return PolynomialKernel(new_fn, symbols=self.syms) + new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) + return PolynomialKernel(new_fn, g=g, symbols=self.syms) def __call__(self, *args): res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) - if not hasattr(res, '__iter__'): - return [res] + #if not hasattr(res, '__iter__'): + # return [res] + #if self.g: + # print(self.g, self.g(res)) + # return self.g(res) return res - def tabulate(self, Qpts, attachment=None): - # TODO do we need to attach qpts - # if attachment: - # return np.array([self(*attachment(*pt)) for pt in Qpts]).astype(np.float64) - return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts, basis_change): + # convert basis change to right format + #breakpoint() + # TODO i don't think this makes sense + return Qpts, np.array([wt*self(*pt)*basis_change for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(len(pt) + 1)] for pt in Qpts] + #return Qpts, np.array([wt*self(*(np.matmul(pt, basis_change))) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), [[(i,) for i in range(len(pt) + 1)] for pt in Qpts] def _to_dict(self): o_dict = {"fn": self.fn} @@ -245,14 +220,16 @@ def dict_id(self): def _from_dict(obj_dict): return PolynomialKernel(obj_dict["fn"]) + class ComponentKernel(BaseKernel): def __init__(self, comp): + self.comp = comp super(ComponentKernel, self).__init__() def __repr__(self): - return f"[{self.comp}]" + return f"[{self.comp}]" def degree(self, interpolant_degree): return interpolant_degree @@ -261,11 +238,12 @@ def permute(self, g): return self def __call__(self, *args): - return args[self.comp] + return tuple(args[i] if i in self.comp else 0 for i in range(len(args))) +# return tuple(args[c] for c in self.comp) - def tabulate(self, Qpts, attachment=None): - return np.array([1 for _ in Qpts]).astype(np.float64) - #return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts, basis_change): + return Qpts, Qwts, [[self.comp] for pt in Qpts] + # return Qpts, np.array([self(*pt) for pt in Qpts]).astype(np.float64) def _to_dict(self): o_dict = {"comp": self.comp} @@ -285,7 +263,7 @@ def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=N self.pairing = pairing self.kernel = kernel self.immersed = immersed - self.trace_entity = entity + self.cell_defined_on = entity self.attachment = attachment self.target_space = target_space self.g = g @@ -302,7 +280,7 @@ def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=N def __call__(self, g): new_generation = self.generation.copy() - return DOF(self.pairing.permute(g), self.kernel.permute(g), self.trace_entity, self.attachment, self.target_space, g, self.immersed, new_generation, self.sub_id, self.cell) + return DOF(self.pairing.permute(g), self.kernel.permute(g), self.cell_defined_on, self.attachment, self.target_space, g, self.immersed, new_generation, self.sub_id, self.cell) def eval(self, fn, pullback=True): return self.pairing(self.kernel, fn, self.cell) @@ -314,8 +292,9 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non # For some of these, we only want to store the first instance of each self.generation[cell.dim()] = dof_gen self.cell = cell - if self.trace_entity is None: + if self.cell_defined_on is None: self.trace_entity = cell + self.cell_defined_on = cell self.pairing = self.pairing.add_entity(cell) if self.target_space is None: self.target_space = space @@ -324,8 +303,36 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non if self.sub_id is None and generator_id is not None: self.sub_id = generator_id - def convert_to_fiat(self, ref_el, interpolant_degree): - return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree) + #if self.immersed and isinstance(self.kernel, ParameterisationKernel): + # fn, syms = self.cell.generate_facet_parameterisation(self.cell_defined_on.id) + # self.kernel = self.kernel.add_context(fn, syms) + + def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): + # TODO deriv dict needs implementing (currently {}) + return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree), {}, str(self)) + + def to_quadrature(self, arg_degree): + Qpts, Qwts = self.cell_defined_on.quadrature(arg_degree) + Qwts = Qwts.reshape(Qwts.shape + (1,)) + + if self.cell_defined_on.get_spatial_dimension() > 0 and self.pairing.orientation: + bvs = np.array(self.cell_defined_on.basis_vectors()) + new_bvs = np.array(self.cell_defined_on.orient(self.pairing.orientation).basis_vectors()) + basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) + else: + basis_change = np.eye(self.cell_defined_on.get_spatial_dimension()) + pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change) + if self.immersed: + # need to compute jacobian from attachment. + pts = [self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts] + immersion = self.target_space.tabulate(wts, self.pairing.entity) + wts = np.outer(wts, immersion) + + # pt dict is { pt: (weight, component)} + pt_dict = {tuple(pt): [(w, c) for w, c in zip(wt, cp)] for pt, wt, cp in zip(pts, wts, comps)} + #if isinstance(self.kernel, ComponentKernel): + # breakpoint() + return pt_dict def __repr__(self, fn="v"): return str(self.pairing).format(fn=fn, kernel=self.kernel) @@ -357,25 +364,30 @@ def eval(self, fn, pullback=True): attached_fn = fn.attach(self.attachment) if pullback: - attached_fn = self.target_space(attached_fn, self.trace_entity, self.g) + attached_fn = self.target_space(attached_fn, self.cell_defined_on) return self.pairing(self.kernel, attached_fn, self.cell) def tabulate(self, Qpts): - immersion = self.target_space.tabulate(Qpts, self.trace_entity, self.g) - res = self.kernel.tabulate(Qpts, self.attachment) + # modify this to take reference space q pts + immersion = self.target_space.tabulate(Qpts, self.pairing.entity) + res, _ = self.kernel.tabulate(Qpts, self.attachment) return immersion*res def __call__(self, g): - index_trace = self.cell.d_entities_ids(self.trace_entity.dim()).index(self.trace_entity.id) - permuted_e, permuted_g = self.cell.permute_entities(g, self.trace_entity.dim())[index_trace] - new_trace_entity = self.cell.get_node(permuted_e).orient(permuted_g) + index_trace = self.cell.d_entities_ids(self.cell_defined_on.dim()).index(self.cell_defined_on.id) + permuted_e, permuted_g = self.cell.permute_entities(g, self.cell_defined_on.dim())[index_trace] + new_cell_defined_on = self.cell.get_node(permuted_e) + #.orient(permuted_g) + #if self.immersed and isinstance(self.kernel, ParameterisationKernel): + # fn, syms = self.cell.generate_facet_parameterisation(new_cell_defined_on.id) + # self.kernel = self.kernel.add_context(fn, syms) new_attach = lambda *x: g(self.attachment(*x)) - return ImmersedDOF(self.pairing.permute(permuted_g), self.kernel.permute(permuted_g), new_trace_entity, + return ImmersedDOF(self.pairing.permute(permuted_g), self.kernel.permute(permuted_g), new_cell_defined_on, new_attach, self.target_space, g, self.triple, self.generation, self.sub_id, self.cell) def __repr__(self): - fn = "tr_{1}_{0}(v)".format(str(self.trace_entity), str(self.target_space)) + fn = "tr_{1}_{0}(v)".format(str(self.cell_defined_on), str(self.target_space)) return super(ImmersedDOF, self).__repr__(fn) def immerse(self, entity, attachment, trace, g): diff --git a/fuse/groups.py b/fuse/groups.py index a08d96e..d6e7457 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -9,12 +9,12 @@ def perm_matrix_to_perm_array(p_mat): summed = np.sum(p_mat, axis=0) - if np.all(summed == np.zeros_like(summed)): + if np.allclose(summed, np.zeros_like(summed)): return list(np.zeros_like(summed)) - assert np.all(summed == np.ones_like(summed)) + assert np.allclose(summed, np.ones_like(summed)) res = [] for row in p_mat: - indices = list(row).index(1) + indices = list(np.isclose(row, 1)).index(True) res += [indices] return res @@ -78,7 +78,11 @@ def __hash__(self): def __mul__(self, x): assert isinstance(x, GroupMemberRep) - return self.group.get_member(self.perm * x.perm) + max_size = max(self.perm.size, x.perm.size) + larger_group = self.group if len(self.group.members()) >= len(x.group.members()) else x.group + x_perm = Permutation(x.perm, size=max_size) + self_perm = Permutation(self.perm, size=max_size) + return larger_group.get_member(self_perm * x_perm) def __invert__(self): return self.group.get_member(~self.perm) @@ -165,12 +169,11 @@ def get_member(self, perm): return m raise ValueError("Permutation not a member of group") - def get_member_by_val(self, val): for m in self.members(): if m.numeric_rep() == val: return m - + raise ValueError("Value does not represent a group member") def compute_num_reps(self, base_val=0): @@ -218,6 +221,8 @@ def __init__(self, base_group, cell=None): if cell is not None: self.cell = cell vertices = cell.vertices(return_coords=True) + verts = cell.ordered_vertices() + vertices = [cell.get_node(v, return_coords=True) for v in verts] self._members = [] counter = 0 @@ -280,12 +285,13 @@ def transform_between_perms(self, perm1, perm2): assert perm2 in member_perms return ~self.get_member(Permutation(perm1)) * self.get_member(Permutation(perm2)) - def get_member(self, perm): - for m in self.members(): - if m.perm == perm: - return m - raise ValueError("Permutation not a member of group") - + #def get_member(self, perm): + # if not isinstance(perm, Permutation): + # perm = Permutation.from_sequence(perm) + # for m in self.members(): + # if m.perm == perm: + # return m + # raise ValueError("Permutation not a member of group") # def compute_reps(self, g, path, remaining_members): # # breadth first search to find generator representations of all members diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index c05844a..422cf05 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -148,8 +148,10 @@ def __init__(self, weights, spaces): self.weights = weights self.spaces = spaces - maxdegree = max([space.maxdegree for space in spaces]) - mindegree = min([space.mindegree for space in spaces]) + weight_degrees = [0 if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)) else max_deg_sp_mat(w) for w in self.weights] + + maxdegree = max([space.maxdegree + w_deg for space, w_deg in zip(spaces, weight_degrees)]) + mindegree = min([space.mindegree + w_deg for space, w_deg in zip(spaces, weight_degrees)]) vec = any([s.set_shape for s in spaces]) super(ConstructedPolynomialSpace, self).__init__(maxdegree, -1, mindegree, set_shape=vec) diff --git a/fuse/traces.py b/fuse/traces.py index 2f77606..b400bfd 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -9,13 +9,13 @@ class Trace(): def __init__(self, cell): self.domain = cell - def __call__(self, trace_entity, g): + def __call__(self, trace_entity): raise NotImplementedError("Trace uninstanitated") - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): raise NotImplementedError("Trace uninstanitated") - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity): raise NotImplementedError("Tabulation uninstantiated") def _to_dict(self): @@ -45,17 +45,17 @@ class TrH1(Trace): def __init__(self, cell): super(TrH1, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): return v - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): ax.scatter(*coord, **kwargs) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): + def to_tikz(self, coord, trace_entity, scale, color="black"): return f"\\filldraw[{color}] {numpy_to_str_tuple(coord, scale)} circle (2pt) node[anchor = south] {{}};" - def tabulate(self, Qpts, trace_entity, g): - return np.ones_like(Qpts) + def tabulate(self, Qwts, trace_entity): + return Qwts def __repr__(self): return "H1" @@ -66,21 +66,21 @@ class TrHDiv(Trace): def __init__(self, cell): super(TrHDiv, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): def apply(*x): - result = np.dot(self.tabulate(None, trace_entity, g), np.array(v(*x)).squeeze()) + result = np.dot(self.tabulate(None, trace_entity), np.array(v(*x)).squeeze()) if isinstance(result, np.float64): # todo: might always be a float return (result,) return tuple(result) return apply - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): # plot dofs of the type associated with this space - vec = self.tabulate([], trace_entity, g).squeeze() + vec = self.tabulate([], trace_entity).squeeze() ax.quiver(*coord, *vec, **kwargs) - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity): entityBasis = np.array(trace_entity.basis_vectors()) cellEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) basis = np.matmul(entityBasis, cellEntityBasis) @@ -94,8 +94,8 @@ def tabulate(self, Qpts, trace_entity, g): return result - def to_tikz(self, coord, trace_entity, g, scale, color="black"): - vec = self.tabulate([], trace_entity, g).squeeze() + def to_tikz(self, coord, trace_entity, scale, color="black"): + vec = self.tabulate([], trace_entity).squeeze() end_point = [coord[i] + 0.25*vec[i] for i in range(len(coord))] arw = "-{Stealth[length=3mm, width=2mm]}" return f"\\draw[thick, {color}, {arw}] {numpy_to_str_tuple(coord, scale)} -- {numpy_to_str_tuple(end_point, scale)};" @@ -109,26 +109,26 @@ class TrHCurl(Trace): def __init__(self, cell): super(TrHCurl, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): def apply(*x): - result = np.dot(self.tabulate(None, trace_entity, g), np.array(v(*x)).squeeze()) + result = np.dot(self.tabulate(None, trace_entity), np.array(v(*x)).squeeze()) if isinstance(result, np.float64): return (result,) return tuple(result) return apply - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity): tangent = np.array(trace_entity.basis_vectors()) subEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) result = np.matmul(tangent, subEntityBasis) return result - def plot(self, ax, coord, trace_entity, g, **kwargs): - vec = self.tabulate([], trace_entity, g).squeeze() + def plot(self, ax, coord, trace_entity, **kwargs): + vec = self.tabulate([], trace_entity).squeeze() ax.quiver(*coord, *vec, **kwargs) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): - vec = self.tabulate([], trace_entity, g).squeeze() + def to_tikz(self, coord, trace_entity, scale, color="black"): + vec = self.tabulate([], trace_entity).squeeze() end_point = [coord[i] + 0.25*vec[i] for i in range(len(coord))] arw = "-{Stealth[length=3mm, width=2mm]}" return f"\\draw[thick, {color}, {arw}] {numpy_to_str_tuple(coord, scale)} -- {numpy_to_str_tuple(end_point, scale)};" @@ -142,8 +142,9 @@ class TrGrad(Trace): def __init__(self, cell): super(TrGrad, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): # Compute grad v and then dot with tangent rotated according to the group member + raise NotImplementedError("Gradient immersions are under development") tangent = np.array(g(np.array(self.domain.basis_vectors())[0])) def apply(*x): @@ -159,11 +160,11 @@ def apply(*x): return tuple(result) return apply - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): circle1 = plt.Circle(coord, 0.075, fill=False, **kwargs) ax.add_patch(circle1) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): + def to_tikz(self, coord, trace_entity, scale, color="black"): return f"\\draw[{color}] {numpy_to_str_tuple(coord, scale)} circle (4pt) node[anchor = south] {{}};" def __repr__(self): @@ -175,7 +176,8 @@ class TrHess(Trace): def __init__(self, cell): super(TrHess, self).__init__(cell) - def __call__(self, v, trace_entity, g): + def __call__(self, v, trace_entity): + raise NotImplementedError("Hessian trace needs reviewing") b0, b1 = self.domain.basis_vectors() tangent0 = np.array(g(b0)) tangent1 = np.array(g(b1)) @@ -192,11 +194,11 @@ def apply(*x): return tuple(result) return apply - def plot(self, ax, coord, trace_entity, g, **kwargs): + def plot(self, ax, coord, trace_entity, **kwargs): circle1 = plt.Circle(coord, 0.15, fill=False, **kwargs) ax.add_patch(circle1) - def to_tikz(self, coord, trace_entity, g, scale, color="black"): + def to_tikz(self, coord, trace_entity, scale, color="black"): return f"\\draw[{color}] {numpy_to_str_tuple(coord, scale)} circle (6pt) node[anchor = south] {{}};" def __repr__(self): diff --git a/fuse/triples.py b/fuse/triples.py index 3d4e3f7..14cb01a 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -77,16 +77,16 @@ def degree(self): def get_dof_info(self, dof, tikz=True): colours = {False: {0: "b", 1: "r", 2: "g", 3: "b"}, True: {0: "blue", 1: "red", 2: "green", 3: "black"}} - if dof.trace_entity.dimension == 0: - center = self.cell.cell_attachment(dof.trace_entity.id)() - elif dof.trace_entity.dimension == 1: - center = self.cell.cell_attachment(dof.trace_entity.id)(0) - elif dof.trace_entity.dimension == 2: - center = self.cell.cell_attachment(dof.trace_entity.id)(0, 0) + if dof.cell_defined_on.dimension == 0: + center = self.cell.cell_attachment(dof.cell_defined_on.id)() + elif dof.cell_defined_on.dimension == 1: + center = self.cell.cell_attachment(dof.cell_defined_on.id)(0) + elif dof.cell_defined_on.dimension == 2: + center = self.cell.cell_attachment(dof.cell_defined_on.id)(0, 0) else: center = list(sum(np.array(self.cell.vertices(return_coords=True)))) - return center, colours[tikz][dof.trace_entity.dimension] + return center, colours[tikz][dof.cell_defined_on.dimension] def get_value_shape(self): # TODO Shape should be specificed somewhere else probably @@ -102,14 +102,13 @@ def to_fiat(self): ref_el = self.cell.to_fiat() dofs = self.generate() degree = self.spaces[0].degree() + value_shape = self.get_value_shape() entity_ids = {} entity_perms = {} nodes = [] top = ref_el.get_topology() min_ids = self.cell.get_starter_ids() poly_set = self.spaces[0].to_ON_polynomial_set(ref_el) - #from FIAT.nedelec import Nedelec - #poly_set = Nedelec(ref_el, 2).poly_set for dim in sorted(top): entity_ids[dim] = {i: [] for i in top[dim]} @@ -120,16 +119,23 @@ def to_fiat(self): for entity in entities: dim = entity[0] for i in range(len(dofs)): - if entity[1] == dofs[i].trace_entity.id - min_ids[dim]: - entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].append(counter) - nodes.append(dofs[i].convert_to_fiat(ref_el, degree)) + if entity[1] == dofs[i].cell_defined_on.id - min_ids[dim]: + entity_ids[dim][dofs[i].cell_defined_on.id - min_ids[dim]].append(counter) + nodes.append(dofs[i].convert_to_fiat(ref_el, degree, value_shape)) + print(nodes[-1].pt_dict) counter += 1 self.matrices_by_entity = self.make_entity_dense_matrices(ref_el, entity_ids, nodes, poly_set) + self.matrices_by_entity = self.dummy_function(ref_el, entity_ids, nodes, poly_set) # TODO remove mat_perms, entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set) - if not pure_perm: - self.matrices = mat_perms - self.reverse_dof_perms() - self.orient_mat_perms() + #if entity_perms is None: + self.matrices = mat_perms + self.reverse_dof_perms() + entity_perms = None + #if not pure_perm: + # self.matrices = mat_perms + #else: + # self.matrices = entity_perms + form_degree = 1 if self.spaces[0].set_shape else 0 # TODO: Change this when Dense case in Firedrake @@ -139,6 +145,9 @@ def to_fiat(self): dual = DualSet(nodes, ref_el, entity_ids) return CiarletElement(poly_set, dual, degree, form_degree) + def dummy_function(self, ref_el, entity_ids, nodes, poly_set): + return self.matrices_by_entity + def to_tikz(self, show=True, scale=3): """Generates tikz code for the element diagram @@ -156,12 +165,12 @@ def to_tikz(self, show=True, scale=3): if isinstance(dof.pairing, DeltaPairing): coord = dof.eval(identity, pullback=False) if isinstance(dof.target_space, Trace): - tikz_commands += [dof.target_space.to_tikz(coord, dof.trace_entity, dof.g, scale, color)] + tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, scale, color)] else: tikz_commands += [f"\\filldraw[{color}] {numpy_to_str_tuple(coord, scale)} circle (2pt) node[anchor = south] {{}};"] elif isinstance(dof.pairing, L2Pairing): coord = center - tikz_commands += [dof.target_space.to_tikz(coord, dof.trace_entity, dof.g, scale, color)] + tikz_commands += [dof.target_space.to_tikz(coord, dof.cell_defined_on, scale, color)] if show: tikz_commands += ['\\end{tikzpicture}'] return "\n".join(tikz_commands) @@ -189,7 +198,7 @@ def plot(self, filename="temp.png"): if len(coord) == 1: coord = (coord[0], 0) if isinstance(dof.target_space, Trace): - dof.target_space.plot(ax, coord, dof.trace_entity, dof.g, color=color) + dof.target_space.plot(ax, coord, dof.cell_defined_on, color=color) else: ax.scatter(*coord, color=color) ax.text(*coord, dof.id) @@ -211,12 +220,12 @@ def plot(self, filename="temp.png"): if isinstance(dof.pairing, DeltaPairing): coord = dof.eval(identity, pullback=False) if isinstance(dof.target_space, Trace): - dof.target_space.plot(ax, coord, dof.trace_entity, dof.g, color=color) + dof.target_space.plot(ax, coord, dof.cell_defined_on, color=color) else: ax.scatter(*coord, color=color) elif isinstance(dof.pairing, L2Pairing): coord = center - dof.target_space.plot(ax, center, dof.trace_entity, dof.g, color=color, length=0.2) + dof.target_space.plot(ax, center, dof.cell_defined_on, color=color, length=0.2) ax.text(*coord, dof.id) plt.axis('off') ax.get_xaxis().set_visible(False) @@ -250,7 +259,7 @@ def compute_dense_matrix(self, ref_el, entity_ids, nodes, poly_set): def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): degree = self.spaces[0].degree() min_ids = self.cell.get_starter_ids() - nodes = [d.convert_to_fiat(ref_el, degree) for d in self.generate()] + nodes = [d.convert_to_fiat(ref_el, degree, self.get_value_shape()) for d in self.generate()] sub_ents = [] res_dict = {} for d in range(0, self.cell.dim() + 1): @@ -260,7 +269,7 @@ def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): dim = e.dim() e_id = e.id - min_ids[dim] res_dict[dim][e_id] = {} - dof_ids = [d.id for d in self.generate() if d.trace_entity == e] + dof_ids = [d.id for d in self.generate() if d.cell_defined_on == e] res_dict[dim][e_id][0] = np.eye(len(dof_ids)) original_V, original_basis = self.compute_dense_matrix(ref_el, entity_ids, nodes, poly_set) @@ -273,9 +282,11 @@ def make_entity_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): elif g.perm.is_Identity: res_dict[dim][e_id][val] = np.eye(len(dof_ids)) else: - new_nodes = [d(g).convert_to_fiat(ref_el, degree) if d.trace_entity == e else d.convert_to_fiat(ref_el, degree) for d in self.generate()] + new_nodes = [d(g).convert_to_fiat(ref_el, degree, self.get_value_shape()) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, degree, self.get_value_shape()) for d in self.generate()] transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_ids, new_nodes, poly_set) res_dict[dim][e_id][val] = np.matmul(transformed_basis, original_V.T)[np.ix_(dof_ids, dof_ids)] + #if dim == 1 and permuted_g.perm.array_form == [1,0]: + # breakpoint() return res_dict def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): @@ -291,7 +302,7 @@ def make_overall_dense_matrices(self, ref_el, entity_ids, nodes, poly_set): if g.perm.is_Identity: res_dict[dim][e_id][val] = np.eye(len(nodes)) else: - new_nodes = [d(g).convert_to_fiat(ref_el, degree) for d in self.generate()] + new_nodes = [d(g).convert_to_fiat(ref_el, degree, self.get_value_shape()) for d in self.generate()] transformed_V, transformed_basis = self.compute_dense_matrix(ref_el, entity_ids, new_nodes, poly_set) res_dict[dim][e_id][val] = np.matmul(transformed_basis, original_V.T) return res_dict @@ -307,8 +318,8 @@ def _entity_associations(self, dofs): # construct mapping of entities to the dof generators and the dofs they generate for d in dofs: - sub_dim = d.trace_entity.dim() - sub_dict = entity_associations[sub_dim][d.trace_entity.id - min_ids[sub_dim]] + sub_dim = d.cell_defined_on.dim() + sub_dict = entity_associations[sub_dim][d.cell_defined_on.id - min_ids[sub_dim]] for dim in set([sub_dim, cell_dim]): dof_gen = str(d.generation[dim]) @@ -353,9 +364,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): dofs = self.generate() min_ids = self.cell.get_starter_ids() entity_associations, pure_perm, sub_pure_perm = self._entity_associations(dofs) - if pure_perm is False: - # TODO think about where this call goes - return self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set), None, pure_perm + # if pure_perm is False: + # TODO think about where this call goes + # return self.make_overall_dense_matrices(ref_el, entity_ids, nodes, poly_set), None, pure_perm + # return self.matrices_by_entity, None, pure_perm oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(dofs) # for each entity, look up generation on that entity and permute the @@ -374,8 +386,9 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): # (dof_gen, ent_dofs) total_ent_dof_ids += [ed.id for ed in ent_dofs if ed.id not in total_ent_dof_ids] dof_gen_class = ent_dofs[0].generation - if not len(dof_gen_class[dim].g2.members()) == 1: + if not len(dof_gen_class[dim].g2.members()) == 1 and dim == min(dof_gen_class.keys()): # if DOFs on entity are not perms, get the matrix + # only get this if they are defined on the current dimension oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = self.matrices_by_entity[dim][e_id][val] elif g.perm.is_Identity or (pure_perm and len(ent_dofs_ids) == 1): oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = np.eye(len(ent_dofs_ids)) @@ -389,11 +402,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() else: # TODO what if an orientation is not in G1 - #warnings.warn("FUSE: should probably be using equivalent members of g in g1 for this") - #sub_mat = g.matrix_form() - #oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() - - #raise NotImplementedError(f"Orientation {g} is not in group {dof_gen_class[dim].g1.members()}") + warnings.warn("FUSE: orientation case not covered") + # sub_mat = g.matrix_form() + # oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() + # raise NotImplementedError(f"Orientation {g} is not in group {dof_gen_class[dim].g1.members()}") pass if len(dof_gen_class.keys()) == 2 and dim == self.cell.dim(): @@ -421,14 +433,15 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): # remove immersed DOFs from flat form oriented_mats_overall = oriented_mats_by_entity[dim][0] - for val, mat in oriented_mats_overall.items(): - cell_dofs = entity_ids[dim][0] - flat_by_entity[dim][e_id][val] = perm_matrix_to_perm_array(mat[np.ix_(cell_dofs, cell_dofs)]) if pure_perm and sub_pure_perm: + for val, mat in oriented_mats_overall.items(): + cell_dofs = entity_ids[dim][0] + flat_by_entity[dim][e_id][val] = perm_matrix_to_perm_array(mat[np.ix_(cell_dofs, cell_dofs)]) return oriented_mats_by_entity, flat_by_entity, True return oriented_mats_by_entity, None, False def orient_mat_perms(self): + raise NotImplementedError("This should not be necessary") min_ids = self.cell.get_starter_ids() entity_orientations = compare_topologies(ufc_cell(self.cell.to_ufl().cellname()).get_topology(), self.cell.get_topology()) num_ents = 0 @@ -436,8 +449,7 @@ def orient_mat_perms(self): ents = self.cell.d_entities(dim) for e in ents: e_id = e.id - min_ids[dim] - members = e.group.members() - if entity_orientations[num_ents + e_id] != 0 and dim < self.cell.dim(): + if entity_orientations[num_ents + e_id] != 0: modifier = self.matrices[dim][e_id][entity_orientations[num_ents+e_id]] reverse_modifier_val = (~e.group.get_member_by_val(entity_orientations[num_ents+e_id])).numeric_rep() reverse_modifier = self.matrices[dim][e_id][reverse_modifier_val] @@ -531,7 +543,7 @@ def make_entity_ids(self): entity_ids[dim] = {i: [] for i in top[dim]} for i in range(len(dofs)): - entity = dofs[i].trace_entity + entity = dofs[i].cell_defined_on dim = entity.dim() entity_ids[dim][entity.id - min_ids[dim]].append(i) return entity_ids diff --git a/fuse/utils.py b/fuse/utils.py index 63d1e48..fcfa3d9 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -17,7 +17,8 @@ def fold_reduce(func_list, *prev): def sympy_to_numpy(array, symbols, values): """ - Convert a sympy array to a numpy array. + TODO: rename this function + Evaluate symbols at values, then convert to numpy if all have been replaced :param: array: sympy array :param: symbols: array of symbols contained in the sympy exprs @@ -27,13 +28,17 @@ def sympy_to_numpy(array, symbols, values): is greater than 1 dimension to remove extra dimensions """ substituted = array.subs({symbols[i]: values[i] for i in range(len(values))}) - nparray = np.array(substituted).astype(np.float64) - if len(nparray.shape) > 1: - return nparray.squeeze() + if len(array.atoms(sp.Symbol)) == len(values) and all(not isinstance(v, sp.Expr) for v in values): + nparray = np.array(substituted).astype(np.float64) - if len(nparray.shape) == 0: - return nparray.item() + if len(nparray.shape) > 1: + return nparray.squeeze() + + if len(nparray.shape) == 0: + return nparray.item() + else: + nparray = substituted return nparray diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index de6be11..2ed8f5e 100644 --- a/test/test_2d_examples_docs.py +++ b/test/test_2d_examples_docs.py @@ -165,7 +165,7 @@ def construct_nd(tri=None): y = sp.Symbol("y") # xs = [DOF(L2Pairing(), PointKernel(edge.basis_vectors()[0]))] - #xs = [DOF(L2Pairing(), PointKernel((np.sqrt(2),)))] + # xs = [DOF(L2Pairing(), PointKernel((np.sqrt(2),)))] xs = [DOF(L2Pairing(), PolynomialKernel((1,)))] dofs = DOFGenerator(xs, S1, S2) @@ -245,6 +245,7 @@ def test_rt_example(): for dof in rt.generate(): assert [np.allclose(1, dof.eval(basis_func).flatten()) for basis_func in basis_funcs].count(True) == 1 assert [np.allclose(0, dof.eval(basis_func).flatten()) for basis_func in basis_funcs].count(True) == 2 + rt.to_fiat() def construct_hermite(): @@ -270,19 +271,19 @@ def construct_hermite(): [v_dofs, v_derv_dofs, v_derv2_dofs, i_dofs]) return her - -def test_hermite_example(): - her = construct_hermite() - - # TODO improve this test - x = sp.Symbol("x") - y = sp.Symbol("y") - phi_0 = FuseFunction(x**2 + 3*y**3 + 4*x*y, symbols=(x, y)) - ls = her.generate() - print("num dofs ", her.num_dofs()) - for dof in ls: - print(dof) - print("dof eval", dof.eval(phi_0)) +# draft of hermite test, immersions need work +#def test_hermite_example(): +# her = construct_hermite() +# +# # TODO improve this test +# x = sp.Symbol("x") +# y = sp.Symbol("y") +# phi_0 = FuseFunction(x**2 + 3*y**3 + 4*x*y, symbols=(x, y)) +# ls = her.generate() +# print("num dofs ", her.num_dofs()) +# for dof in ls: +# print(dof) +# print("dof eval", dof.eval(phi_0)) def test_square_cg(): diff --git a/test/test_3d_examples_docs.py b/test/test_3d_examples_docs.py index 9123a5a..1d1d5dd 100644 --- a/test/test_3d_examples_docs.py +++ b/test/test_3d_examples_docs.py @@ -1,4 +1,6 @@ from fuse import * +import pytest +from sympy.combinatorics import Permutation import sympy as sp import numpy as np @@ -93,7 +95,7 @@ def construct_tet_rt(cell=None): Pd = PolynomialSpace(deg - 1) rt_space = vec_Pd + (Pd.restrict(deg - 2, deg - 1))*M - xs = [DOF(L2Pairing(), PolynomialKernel((1, 1, 1)))] + xs = [DOF(L2Pairing(), PolynomialKernel(1))] dofs = DOFGenerator(xs, S1, S3) face_vec = ElementTriple(face, (rt_space, CellHDiv, "C0"), dofs) @@ -106,11 +108,42 @@ def construct_tet_rt(cell=None): return rt1 +def construct_tet_ned(cell=None): + deg = 1 + tri = make_tetrahedron() + edge = tri.edges()[0] + + x = sp.Symbol("x") + y = sp.Symbol("y") + z = sp.Symbol("z") + + M1 = sp.Matrix([[0, z, -y]]) + M2 = sp.Matrix([[z, 0, -x]]) + M3 = sp.Matrix([[y, -x, 0]]) + + vec_Pd = PolynomialSpace(deg - 1, set_shape=True) + Pd = PolynomialSpace(deg - 1) + nd_space = vec_Pd + (Pd.restrict(deg - 2, deg - 1))*M1 + (Pd.restrict(deg - 2, deg - 1))*M2 + (Pd.restrict(deg - 2, deg - 1))*M3 + + # [test_tet_ned 0] + xs = [DOF(L2Pairing(), PolynomialKernel(1))] + dofs = DOFGenerator(xs, S1, S2) + + edges = ElementTriple(edge, (vec_Pd, CellHCurl, L2), dofs) + xs = [immerse(tri, edges, TrHCurl)] + tet_edges = PermutationSetRepresentation([Permutation([0, 1, 2, 3]), Permutation([1, 2, 3, 0]), + Permutation([2, 3, 0, 1]), Permutation([1, 3, 0, 2]), + Permutation([2, 0, 1, 3]), Permutation([3, 0, 1, 2])]) + edge_dofs = DOFGenerator(xs, tet_edges, S1) + # [test_tet_ned 1] + + return ElementTriple(tri, (nd_space, CellHCurl, L2), [edge_dofs]) + + def plot_tet_rt(): rt = construct_tet_rt() rt.plot() - def test_tet_rt(): tetra = make_tetrahedron() rt1 = construct_tet_rt(tetra) @@ -118,24 +151,18 @@ def test_tet_rt(): # TODO make this a proper test for dof in ls: print(dof) + rt1.to_fiat() - -def construct_tet_ned(): +@pytest.mark.xfail(reason='DOFs not forming full rank matrix') +def test_tet_nd(): tetra = make_tetrahedron() - edge = tetra.edges()[0] - # [test_tet_ned 0] - xs = [DOF(L2Pairing(), PolynomialKernel(1))] - dofs = DOFGenerator(xs, S1, S2) - int_ned = ElementTriple(edge, (P1, CellHCurl, "C0"), dofs) - - im_xs = [immerse(tetra, int_ned, TrHCurl)] - edge = DOFGenerator(im_xs, tet_edges, Z4) - - ned = ElementTriple(tetra, (P1, CellHCurl, "C0"), - [edge]) - # [test_tet_ned 1] - - return ned + nd1 = construct_tet_ned(tetra) + ls = nd1.generate() + # TODO make this a proper test + for dof in ls: + print(dof) + nd1.to_fiat() + def plot_tet_ned(): @@ -143,9 +170,9 @@ def plot_tet_ned(): ned.plot() -def test_tet_ned(): - ned = construct_tet_ned() - ls = ned.generate() - # TODO make this a proper test - for dof in ls: - print(dof) +#def test_tet_ned(): +# ned = construct_tet_ned() +# ls = ned.generate() +# # TODO make this a proper test +# for dof in ls: +# print(dof) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 128a53e..ee284c4 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -470,7 +470,7 @@ def test_helmholtz_3d(elem_gen, elem_code, deg, conv_rate): def helmholtz_solve(V, mesh): # Define variational problem - dim = mesh.ufl_cell().topological_dimension() + dim = mesh.ufl_cell().topological_dimension u = TrialFunction(V) v = TestFunction(V) f = Function(V) @@ -493,10 +493,11 @@ def helmholtz_solve(V, mesh): # l_a = assemble(L) # elem = V.finat_element.fiat_equivalent - # W = VectorFunctionSpace(mesh, V.ufl_element()) - # X = assemble(interpolate(mesh.coordinates, W)) - # print(X.dat.data) - # print(assemble(a).M.values) + #W = VectorFunctionSpace(mesh, V.ufl_element()) + #X = assemble(interpolate(mesh.coordinates, W)) + #print(X.dat.data) + #print(assemble(a).M.values) + #breakpoint() # Compute solution sol = Function(V) @@ -594,10 +595,10 @@ def test_project_vec(elem_gen, elem_code, deg): mesh = UnitTriangleMesh() U = FunctionSpace(mesh, elem_code, deg) - assert np.allclose(project(U, mesh, as_vector((1,1))), 0, rtol=1e-5) + assert np.allclose(project(U, mesh, as_vector((1, 1))), 0, rtol=1e-5) U = FunctionSpace(mesh, elem.to_ufl()) - assert np.allclose(project(U, mesh, as_vector((1,1))), 0, rtol=1e-5) + assert np.allclose(project(U, mesh, as_vector((1, 1))), 0, rtol=1e-5) @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_dg1_tet, "DG", 1)]) diff --git a/test/test_dofs.py b/test/test_dofs.py index 5b56f43..72c0838 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -1,5 +1,7 @@ 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 @@ -75,11 +77,32 @@ def test_permute_nd(): # print(dof) for g in nd.cell.group.members(): - print("g:", g.numeric_rep()) + 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)) + breakpoint() + + def test_permute_nd_old(): cell = polygon(3) @@ -138,3 +161,66 @@ def test_permute_nodes(): 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)) + + elem.to_fiat() diff --git a/test/test_orientations.py b/test/test_orientations.py index de81546..23d7c2e 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -1,60 +1,100 @@ import unittest.mock as mock +import pytest from firedrake import * from fuse import * import sympy as sp -from test_convert_to_fiat import create_cg1, helmholtz_solve, construct_nd, construct_rt +from test_convert_to_fiat import create_cg1, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg1 import os -def dummy_dof_perms(cls, *args, **kwargs): - # return -1s of right shape here - oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) - for key1, val1 in oriented_mats_by_entity.items(): - for key2, val2 in oriented_mats_by_entity[key1].items(): - for key3, val3 in oriented_mats_by_entity[key1][key2].items(): - if key1 == 2: - oriented_mats_by_entity[key1][key2][key3] = 1 * np.identity(val3.shape[0]) - #oriented_mats_by_entity[key1][key2][key3][0] = key1*100 + key2*10 + key3 - if key3 == 5: - - oriented_mats_by_entity[key1][key2][key3] = 100 * np.identity(val3.shape[0]) - return oriented_mats_by_entity, False, None - - -def test_interpolation_values(): - cell = polygon(3) - elem = construct_nd(cell) - print() - #with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): - mesh = UnitSquareMesh(3, 3) - ones = as_vector((0,1)) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): - V = FunctionSpace(mesh, elem.to_ufl()) - else: - V = FunctionSpace(mesh, "N1curl", 1) - print(V.cell_node_list) - u = TestFunction(V) - res1= assemble(interpolate(ones, V)) - for i in range(len(res1.dat.data)): - print(f"{i}: {res1.dat.data[i]}") +def construct_nd2(tri=None): + if tri is None: + tri = polygon(3) + deg = 2 + edge = tri.edges()[0] + x = sp.Symbol("x") + y = sp.Symbol("y") + + xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] + + dofs = DOFGenerator(xs, S2, S2) + int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) + v_2 = np.array(tri.get_node(tri.ordered_vertices()[2], return_coords = True)) + v_1 = np.array(tri.get_node(tri.ordered_vertices()[1], return_coords = True)) + xs = [DOF(L2Pairing(), PolynomialKernel((v_2 - v_1)/2))] + + center_dofs = DOFGenerator(xs, S2, 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]) + return ned + + +def construct_nd2_for_fiat(tri=None): + if tri is None: + tri = polygon(3) + deg = 2 + edge = tri.edges()[0] + x = sp.Symbol("x") + y = sp.Symbol("y") + + xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), + DOF(L2Pairing(), PolynomialKernel((3*x*np.sqrt(2)/np.sqrt(3)), symbols=[x]))] + dofs = DOFGenerator(xs, S1, S2) + int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) -def test_surface_const_nd(): + xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), + DOF(L2Pairing(), PolynomialKernel((-np.sqrt(2)/np.sqrt(3))*(6*x + 3), symbols=[x]))] + dofs = DOFGenerator(xs, S1, S2) + int_ned2 = ElementTriple(tri.edges()[0], (P1, CellHCurl, C0), dofs) + + xs = [DOF(L2Pairing(), PolynomialKernel(np.sqrt(2))), + DOF(L2Pairing(), PolynomialKernel((-np.sqrt(2)/np.sqrt(3))*(6*x - 3), symbols=[x]))] + dofs = DOFGenerator(xs, S1, S2) + int_ned3 = ElementTriple(tri.edges()[0], (P1, CellHCurl, C0), dofs) + + xs = [DOF(L2Pairing(), ComponentKernel((0,))), + DOF(L2Pairing(), ComponentKernel((1,)))] + center_dofs = DOFGenerator(xs, S1, S3) + xs = [immerse(tri, int_ned2, TrHCurl, node=0)] + xs += [immerse(tri, int_ned1, TrHCurl, node=1)] + xs += [immerse(tri, int_ned3, TrHCurl, node=2)] + tri_dofs = DOFGenerator(xs, S1, 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]) + return ned + + +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), + (construct_nd2, "N1curl", 2)]) +def test_surface_const_nd(elem_gen, elem_code, deg): cell = polygon(3) - elem = construct_nd(cell) - ones = as_vector((0,1)) + elem = elem_gen(cell) + ones = as_vector((0, 1)) - for n in range(1, 6): + for n in range(2, 6): mesh = UnitSquareMesh(n, n) - if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: - V = FunctionSpace(mesh, "N1curl", 1) + V = FunctionSpace(mesh, elem_code, deg) + print(V.cell_node_list) normal = FacetNormal(mesh) ones1 = interpolate(ones, V) - res1= assemble(dot(ones1, normal) * ds) - + res1 = assemble(dot(ones1, normal) * ds) + print(f"{n}: {res1}") assert np.allclose(res1, 0) @@ -62,7 +102,7 @@ def test_surface_const_nd(): def test_surface_const_rt(): cell = polygon(3) elem = construct_rt(cell) - ones = as_vector((1,0)) + ones = as_vector((1, 0)) for n in range(1, 6): mesh = UnitSquareMesh(n, n) @@ -73,23 +113,25 @@ def test_surface_const_rt(): V = FunctionSpace(mesh, "RT", 1) normal = FacetNormal(mesh) ones1 = interpolate(ones, V) - res1= assemble(dot(ones1, normal) * ds) - + res1 = assemble(dot(ones1, normal) * ds) + print(f"{n}: {res1}") assert np.allclose(res1, 0) -def test_surface_vec(): +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_nd, "N1curl", 1), + (construct_nd2, "N1curl", 2)]) +def test_surface_vec(elem_gen,elem_code,deg): cell = polygon(3) rt_elem = construct_rt(cell) - nd_elem = construct_nd(cell) + nd_elem = elem_gen(cell) for n in range(1, 6): mesh = UnitSquareMesh(n, n) x, y = SpatialCoordinate(mesh) normal = FacetNormal(mesh) - test_vec = as_vector((-y,x)) + test_vec = as_vector((-y, x)) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, rt_elem.to_ufl()) vec1 = interpolate(test_vec, V) @@ -99,13 +141,17 @@ def test_surface_vec(): vec2 = interpolate(test_vec, V) res1 = assemble(dot(vec2, normal) * ds) print(f"div {n}: {res1}") - + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, nd_elem.to_ufl()) + print(V.cell_node_list) vec1 = interpolate(test_vec, V) + # V2 = FunctionSpace(mesh, create_dg0(cell).to_ufl()) + # u = TestFunction(V2) res2 = assemble(dot(vec1, normal) * ds) else: - V = FunctionSpace(mesh, "N1curl", 1) + V = FunctionSpace(mesh, elem_code, deg) + print(V.cell_node_list) vec1 = interpolate(test_vec, V) res2 = assemble(dot(vec1, normal) * ds) print(f"curl {n}: {res2}") @@ -114,9 +160,9 @@ def test_surface_vec(): assert np.allclose(0, res2) -def test_interpolate_vs_project(V): +def get_expression(V): mesh = V.mesh() - dim = mesh.geometric_dimension() + dim = mesh.geometric_dimension if dim == 2: x, y = SpatialCoordinate(mesh) elif dim == 3: @@ -125,103 +171,226 @@ def test_interpolate_vs_project(V): shape = V.value_shape if dim == 2: if len(shape) == 0: - expression = x + y + exact = Function(FunctionSpace(mesh, 'CG', 5)) + expression = x elif len(shape) == 1: - expression = as_vector([x, y]) - elif len(shape) == 2: - expression = as_tensor(([x, y], [x, y])) + exact = Function(VectorFunctionSpace(mesh, 'CG', 5)) + expr = x + y + expression = as_vector([expr, expr]) elif dim == 3: if len(shape) == 0: + exact = Function(FunctionSpace(mesh, 'CG', 5)) expression = x + y + z elif len(shape) == 1: + exact = Function(FunctionSpace(mesh, 'CG', 5)) expression = as_vector([x, y, z]) - elif len(shape) == 2: - expression = as_tensor(([x, y, z], [x, y, z], [x, y, z])) + return expression, exact +def interpolate_vs_project(V, expression, exact): f = assemble(interpolate(expression, V)) - expect = project(expression, V) + #expect = project(expression, V) + exact.interpolate(expression) + print(exact.dat.data) + #return sqrt(assemble(inner((expect - exact), (expect - exact)) * dx)), print(f.dat.data) - print(expect.dat.data) - assert np.allclose(f.dat.data, expect.dat.data, atol=1e-06) + return 0, sqrt(assemble(inner((f - exact), (f - exact)) * dx)), -def construct_nd2(tri=None): - if tri is None: - tri = polygon(3) - deg = 2 - edge = tri.edges()[0] - x = sp.Symbol("x") - y = sp.Symbol("y") - #xs = [DOF(L2Pairing(), PointKernel((-np.sqrt(2),)))] - xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x+1), symbols=[x])), - DOF(L2Pairing(), PolynomialKernel((1/2)*(1-x), symbols=[x]))] - - - dofs = DOFGenerator(xs, S1, S2) - int_ned = 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_ned, 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]) - return ned -def test_degree2_interpolation(): +@pytest.mark.parametrize("elem_gen,elem_code,deg,conv_rate", [(construct_nd2, "N1curl", 2, 1.8)]) +def test_convergence(elem_gen,elem_code,deg,conv_rate): cell = polygon(3) - #elem = construct_nd2(cell) - mesh = UnitSquareMesh(1,1) - - #print("fuse") - #V = FunctionSpace(mesh, elem.to_ufl()) - #test_interpolate_vs_project(V) + elem = elem_gen(cell) + scale_range = range(2,4) + diff_proj = [0 for i in scale_range] + diff_inte = [0 for i in scale_range] + for n in scale_range: + mesh = UnitSquareMesh(2**n, 2**n) + + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + x,y = SpatialCoordinate(mesh) + expr = cos(x*pi*2)*sin(y*pi*2) + expr = as_vector([expr,expr]) + _, exact = get_expression(V) + diff_proj[n-min(scale_range)], diff_inte[n-min(scale_range)] = interpolate_vs_project(V, expr, exact) + + #print("projection l2 error norms:", diff_proj) + #diff_proj = np.array(diff_proj) + #conv1 = np.log2(diff_proj[:-1] / diff_proj[1:]) + #print("convergence order:", conv1) + #assert all([c > conv_rate for c in conv1]) + print("interpolation l2 error norms:", diff_inte) + diff_inte = np.array(diff_inte) + conv1 = np.log2(diff_inte[:-1] / diff_inte[1:]) + print("convergence order:", conv1) + assert all([c > conv_rate for c in conv1]) + + +@pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), + (create_cg1, "CG", 1), + (create_dg1, "DG", 1), + (construct_cg3, "CG", 3), + (construct_nd2, "N1curl", 2)]) +def test_interpolation(elem_gen, elem_code, deg): + cell = polygon(3) + elem = elem_gen(cell) + #mesh = UnitTriangleMesh(0) + mesh = UnitSquareMesh(2,2) - print("firedrake") - V = FunctionSpace(mesh, "N1curl", 2) - test_interpolate_vs_project(V) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + expression, _ = get_expression(V) + expect = project(expression, V) + f = assemble(interpolate(expression, V)) + assert np.allclose(f.dat.data, expect.dat.data) + +@pytest.mark.xfail(reason="issues with tets") +def test_projection_convergence(elem_gen, elem_code, deg, conv_rate): + cell = make_tetrahedron() + elem = elem_gen(cell) + function = lambda x: cos((3/4)*pi*x[0]) + + scale_range = range(1, 6) + diff = [0 for i in scale_range] + for i in scale_range: + mesh = UnitSquareMesh(2 ** i, 2 ** i) + x, y = SpatialCoordinate(mesh) + expr = cos(x*pi*2)*sin(y*pi*2) + expr = as_vector([expr, expr]) + expression = as_vector([x, y]) + exact = Function(VectorFunctionSpace(mesh, 'CG', 5)) + exact.interpolate(expr) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V2 = FunctionSpace(mesh, elem.to_ufl()) + res2 = project(function(x), V2) + diff[i - 1] = res2 + else: + V = FunctionSpace(mesh, elem_code, deg) + #res = project(function(x), V) + res = project(expr, V, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'}) + diff[i - 1] = sqrt(assemble(inner((res - exact), (res - exact)) * dx)) -def test_create_fiat_nd(): - cell = polygon(3) - nd = construct_nd2(cell) - ref_el = cell.to_fiat() - sd = ref_el.get_spatial_dimension() - deg = 2 - from FIAT.nedelec import Nedelec - fiat_elem = Nedelec(ref_el, deg) - print("fiat") - for f in fiat_elem.dual_basis(): - print(f.pt_dict) - - print("fuse") - for d in nd.generate(): - print(d.convert_to_fiat(ref_el, deg).pt_dict) - print(d) - breakpoint() - my_elem = nd.to_fiat() + #assert np.allclose(res1, res2) + print("l2 error norms:", diff) + diff = np.array(diff) + conv1 = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv1) - Q = create_quadrature(ref_el, 2*(deg+1)) - Qpts, _ = Q.get_points(), Q.get_weights() + #print("fuse l2 error norms:", diff) + #diff = np.array(diff) + #conv2 = np.log2(diff[:-1] / diff[1:]) + #print("fuse convergence order:", conv2) - fiat_vals = fiat_elem.tabulate(0, Qpts) - my_vals = my_elem.tabulate(0, Qpts) + assert (np.array(conv1) > conv_rate).all() + #assert (np.array(conv2) > conv_rate).all() - fiat_vals = flatten(fiat_vals[(0,) * sd]) - my_vals = flatten(my_vals[(0,) * sd]) - (x, res, _, _) = np.linalg.lstsq(fiat_vals.T, my_vals.T) - x1 = np.linalg.inv(x) - assert np.allclose(np.linalg.norm(my_vals.T - fiat_vals.T @ x), 0) - assert np.allclose(np.linalg.norm(fiat_vals.T - my_vals.T @ x1), 0) - assert np.allclose(res, 0) +@pytest.mark.parametrize("elem_gen,elem_gen2,elem_code,deg,deg2", [ + # (create_cg2_tri, "CG", 2), + (create_cg1, create_cg2_tri, "CG", 1, 2), + # (create_dg1, "DG", 1), + # (construct_cg3, "CG", 3), + (construct_nd, construct_nd2, "N1curl", 1, 2), + ]) +def test_two_form(elem_gen, elem_gen2, elem_code, deg, deg2): + cell = polygon(3) + elem = elem_gen(cell) + mesh = UnitSquareMesh(1, 1) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V = FunctionSpace(mesh, elem.to_ufl()) + else: + V = FunctionSpace(mesh, elem_code, deg) + if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): + V2 = FunctionSpace(mesh, elem_gen2(cell).to_ufl()) + else: + V2 = FunctionSpace(mesh, elem_code, deg2) + v = TestFunction(V) + u = TrialFunction(V2) + res =assemble(inner(u, v)*dx).M.values + for row in res: + print(row) + +#def test_create_fiat_nd(): +# cell = polygon(3) +# nd = construct_nd2(cell) +# nd_fv = construct_nd2_for_fiat(cell) +# ref_el = cell.to_fiat() +# deg = 2 +# +# from FIAT.nedelec import Nedelec +# fiat_elem = Nedelec(ref_el, deg) +# print("fiat") +# for f in fiat_elem.dual_basis(): +# print(f.pt_dict) +# +# #print("fiat: fuse's version") +# #for d in nd_fv.generate(): +# # print(d.convert_to_fiat(ref_el, deg).pt_dict) +# +# print("fuse") +# dofs = nd.generate() +# e = cell.edges()[0] +# s3 = S3.add_cell(cell) +# g = s3.get_member([1,2,0]) +# print(dofs[-1].to_quadrature(2)) +# print(dofs[-1](g).to_quadrature(2)) +# for d in nd.generate(): +# print(d.convert_to_fiat(ref_el, deg, (2,)).pt_dict) +# print(g) +# for d in nd.generate(): +# print(d(g).convert_to_fiat(ref_el, deg, (2,)).pt_dict) +# #nodes = [d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] +# #new_nodes = [d(g).convert_to_fiat(ref_el, deg, (2,)) if d.cell_defined_on == e else d.convert_to_fiat(ref_el, deg, (2,)) for d in dofs] +# #for i in range(len(new_nodes)): +# # print(f"{nodes[i].pt_dict}") +# # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") +# # print(f"{dofs[i]}: {new_nodes[i].pt_dict}") +# # for g in S2.add_cell(cell).members(): +# # print(d(g)) +# # print(f"{g} {d(g).convert_to_fiat(ref_el, deg).pt_dict}") +# +# nd1= construct_nd(cell) +# print("nd1") +# for d in nd1.generate(): +# print(d.convert_to_fiat(ref_el, deg).pt_dict) +# print(g) +# for d in nd1.generate(): +# print(d(g).convert_to_fiat(ref_el, deg).pt_dict) +# nd1.to_fiat() +# +# nd.to_fiat() +# #nd_fv.to_fiat() +# +# +#def test_cg3(): +# tri = polygon(3) +# elem = construct_cg3(tri) +# for d in elem.generate(): +# print(d.to_quadrature(1)) +# breakpoint() +# +# +#def test_int_nd(): +# tri = polygon(3) +# deg = 2 +# edge = tri.edges()[0] +# x = sp.Symbol("x") +# y = sp.Symbol("y") +# +# xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] +# +# dofs = DOFGenerator(xs, S2, S2) +# int_ned1 = ElementTriple(edge, (P1, CellHCurl, C0), dofs) +# dofs = int_ned1.generate() +# for d in dofs: +# print(d.to_quadrature(1)) +# breakpoint() diff --git a/test/test_polynomial_space.py b/test/test_polynomial_space.py index cb1c769..0da6486 100644 --- a/test/test_polynomial_space.py +++ b/test/test_polynomial_space.py @@ -152,6 +152,41 @@ def flatten(coeffs): nc = np.reshape(coeffs, newshape) return nc +@pytest.mark.parametrize("deg", [1, 2, 3, 4]) +def test_3d_nd_construction(deg): + cell = make_tetrahedron() + ref_el = cell.to_fiat() + sd = ref_el.get_spatial_dimension() + x = sp.Symbol("x") + y = sp.Symbol("y") + z = sp.Symbol("z") + M1 = sp.Matrix([[0, z, -y]]) + M2 = sp.Matrix([[z, 0, -x]]) + M3 = sp.Matrix([[y, -x, 0]]) + + vec_Pd = PolynomialSpace(deg - 1, set_shape=True) + Pd = PolynomialSpace(deg - 1) + composite = vec_Pd + (Pd.restrict(deg - 2, deg - 1))*M1 + (Pd.restrict(deg - 2, deg - 1))*M2 + (Pd.restrict(deg - 2, deg - 1))*M3 + + assert composite.set_shape + assert isinstance(composite, ConstructedPolynomialSpace) + on_set = composite.to_ON_polynomial_set(ref_el) + + from FIAT.nedelec import NedelecSpace3D + fiat_nd_space = NedelecSpace3D(ref_el, deg) + + Q = create_quadrature(ref_el, 2*(deg + 1)) + Qpts, _ = Q.get_points(), Q.get_weights() + fiat_vals = fiat_nd_space.tabulate(Qpts)[(0,) * sd] + my_vals = on_set.tabulate(Qpts)[(0,) * sd] + fiat_vals = flatten(fiat_vals) + my_vals = flatten(my_vals) + + (x, res, _, _) = np.linalg.lstsq(fiat_vals.T, my_vals.T) + x1 = np.linalg.inv(x) + assert np.allclose(np.linalg.norm(my_vals.T - fiat_vals.T @ x), 0) + assert np.allclose(np.linalg.norm(fiat_vals.T - my_vals.T @ x1), 0) + assert np.allclose(res, 0) @pytest.mark.parametrize("deg", [1, 2, 3, 4]) def test_3d_rt_construction(deg):