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/__init__.py b/fuse/__init__.py index dfd083c..7480548 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, ParameterisationKernel diff --git a/fuse/cells.py b/fuse/cells.py index 5840c2e..9cf7544 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 @@ -656,9 +657,9 @@ def to_tikz(self, show=True, scale=3): return "\n".join(tikz_commands) return tikz_commands - def generate_facet_parameterisation(self, facet_num): - #facet = self.d_entities(self.dimension - 1)[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: @@ -667,11 +668,9 @@ def generate_facet_parameterisation(self, facet_num): 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]) - midpoint = np.average(v_coords, axis=0) - stacked = np.c_[np.ones((self.dimension,)), v_coords[:, 0].reshape(self.dimension,1)] + 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 @@ -679,7 +678,6 @@ def generate_facet_parameterisation(self, facet_num): 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: @@ -790,10 +788,12 @@ def attachment(self, source, dst): assert all(np.isclose(val, vals[0]).all() for val in vals) 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): @@ -803,7 +803,7 @@ def cell_attachment(self, dst): def orient(self, o): """ Orientation node is always labelled with -1 """ - #if self.oriented: + # if self.oriented: # o = self.oriented * o oriented_point = copy.deepcopy(self) top_level_node = oriented_point.d_entities_ids( @@ -1103,7 +1103,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 5f9152c..dbb847f 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 @@ -36,6 +36,7 @@ def __call__(self, kernel, v, cell): return v(*kernel.pt) def convert_to_fiat(self, ref_el, dof, interpolant_deg): + raise NotImplementedError("this should be outdated") pt = dof.eval(FuseFunction(lambda *x: x)) # pt1 = dof.tabulate([[1]]) return PointEvaluation(ref_el, pt) @@ -73,27 +74,11 @@ 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): + raise NotImplementedError("this should be outdated") pass def add_entity(self, entity): @@ -110,12 +95,13 @@ def permute(self, g): res = L2Pairing() if self.entity: res.entity = self.entity.orient(g) - #if self.orientation is not None: + # if self.orientation is not None: # g = self.orientation * g res.orientation = g return res def convert_to_fiat(self, ref_el, dof, interpolant_degree): + raise NotImplementedError("this should be outdated") 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()) @@ -128,8 +114,8 @@ def convert_to_fiat(self, ref_el, dof, interpolant_degree): 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 + # (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. @@ -194,10 +180,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): + 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} @@ -224,13 +208,12 @@ 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, symbols=self.syms) def __call__(self, *args): res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) @@ -238,8 +221,8 @@ def __call__(self, *args): return [res] return res - def tabulate(self, Qpts, attachment=None): - return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts): + return Qpts, np.array([wt*self(*pt) 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} @@ -251,6 +234,7 @@ def dict_id(self): def _from_dict(obj_dict): return PolynomialKernel(obj_dict["fn"]) + class ComponentKernel(BaseKernel): def __init__(self, comp): @@ -258,7 +242,7 @@ def __init__(self, comp): super(ComponentKernel, self).__init__() def __repr__(self): - return f"[{self.comp}]" + return f"[{self.comp}]" def degree(self, interpolant_degree): return interpolant_degree @@ -269,9 +253,9 @@ def permute(self, g): def __call__(self, *args): return args[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): + 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} @@ -283,9 +267,11 @@ def dict_id(self): def _from_dict(obj_dict): return ComponentKernel(obj_dict["comp"]) + class ParameterisationKernel(BaseKernel): def __init__(self, g=None): + raise NotImplementedError("Parameterisation kernel not needed, use polynomial") self.g = g self.fn = None self.syms = None @@ -293,7 +279,7 @@ def __init__(self, g=None): def add_context(self, fn=None, syms=None): new_cls = ParameterisationKernel(self.g) - new_cls.fn = fn + new_cls.fn = fn new_cls.syms = syms return new_cls @@ -329,11 +315,10 @@ def __call__(self, *args): res = sympy_to_numpy(fn, self.syms, args[:len(self.syms)]) if not hasattr(res, '__iter__'): res = [res] - #res = np.ones_like(res) return res - def tabulate(self, Qpts, attachment=None): - return np.array([self(*pt) for pt in Qpts]).astype(np.float64) + def evaluate(self, Qpts, Qwts): + return Qpts, np.array([wt*self(*pt) 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} @@ -345,13 +330,14 @@ def dict_id(self): def _from_dict(obj_dict): return ParameterisationKernel(obj_dict["fn"]) + class DOF(): def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=None, g=None, immersed=False, generation=None, sub_id=None, cell=None): 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 @@ -368,7 +354,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) @@ -380,8 +366,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 @@ -391,12 +378,25 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non self.sub_id = generator_id if self.immersed and isinstance(self.kernel, ParameterisationKernel): - fn, syms = self.cell.generate_facet_parameterisation(self.trace_entity.id) + 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): - return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree) + 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,)) + pts, wts, comps = self.kernel.evaluate(Qpts, Qwts) + 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, self.g) + 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)} + return pt_dict def __repr__(self, fn="v"): return str(self.pairing).format(fn=fn, kernel=self.kernel) @@ -428,29 +428,29 @@ 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, self.g) return self.pairing(self.kernel, attached_fn, self.cell) def tabulate(self, Qpts): # modify this to take reference space q pts immersion = self.target_space.tabulate(Qpts, self.pairing.entity, self.g) - res = self.kernel.tabulate(Qpts, self.attachment) + 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_trace_entity.id) + 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): @@ -505,5 +505,3 @@ def _to_dict(self): def dict_id(self): return "Function" - - diff --git a/fuse/groups.py b/fuse/groups.py index 14a5a41..3ddc77c 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -165,12 +165,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): @@ -285,7 +284,6 @@ def get_member(self, perm): 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 56bbce1..422cf05 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -147,8 +147,8 @@ def __init__(self, weights, spaces): self.weights = weights self.spaces = 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] + + 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)]) diff --git a/fuse/traces.py b/fuse/traces.py index 2f77606..e283dfb 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -15,7 +15,7 @@ def __call__(self, trace_entity, g): def plot(self, ax, coord, trace_entity, g, **kwargs): raise NotImplementedError("Trace uninstanitated") - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity, g): raise NotImplementedError("Tabulation uninstantiated") def _to_dict(self): @@ -54,8 +54,8 @@ def plot(self, ax, coord, trace_entity, g, **kwargs): def to_tikz(self, coord, trace_entity, g, 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, g): + return Qwts def __repr__(self): return "H1" @@ -80,7 +80,7 @@ def plot(self, ax, coord, trace_entity, g, **kwargs): vec = self.tabulate([], trace_entity, g).squeeze() ax.quiver(*coord, *vec, **kwargs) - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity, g): entityBasis = np.array(trace_entity.basis_vectors()) cellEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) basis = np.matmul(entityBasis, cellEntityBasis) @@ -117,7 +117,7 @@ def apply(*x): return tuple(result) return apply - def tabulate(self, Qpts, trace_entity, g): + def tabulate(self, Qwts, trace_entity, g): tangent = np.array(trace_entity.basis_vectors()) subEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity)) result = np.matmul(tangent, subEntityBasis) diff --git a/fuse/triples.py b/fuse/triples.py index bab3630..fb4ab21 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,6 +102,7 @@ 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 = [] @@ -118,16 +119,18 @@ 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)) 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 + 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() + else: + self.matrices = entity_perms form_degree = 1 if self.spaces[0].set_shape else 0 # TODO: Change this when Dense case in Firedrake @@ -157,12 +160,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, dof.g, 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, dof.g, scale, color)] if show: tikz_commands += ['\\end{tikzpicture}'] return "\n".join(tikz_commands) @@ -190,7 +193,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, dof.g, color=color) else: ax.scatter(*coord, color=color) ax.text(*coord, dof.id) @@ -212,12 +215,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, dof.g, 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, dof.g, color=color, length=0.2) ax.text(*coord, dof.id) plt.axis('off') ax.get_xaxis().set_visible(False) @@ -251,7 +254,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): @@ -261,7 +264,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) @@ -274,11 +277,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() + # 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): @@ -294,7 +297,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 @@ -310,8 +313,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]) @@ -356,9 +359,9 @@ 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) @@ -378,8 +381,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)) @@ -393,11 +397,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(): @@ -433,6 +436,7 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): 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 @@ -440,8 +444,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()#and dim < self.cell.dim() - if entity_orientations[num_ents + e_id] != 0 : + 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] @@ -535,7 +538,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 3b3e069..fcfa3d9 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -34,7 +34,7 @@ def sympy_to_numpy(array, symbols, values): if len(nparray.shape) > 1: return nparray.squeeze() - + if len(nparray.shape) == 0: return nparray.item() else: diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index de6be11..8b64f6e 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) diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index 128a53e..f405009 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) @@ -594,10 +594,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 d007fd7..7afbbe6 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 @@ -139,45 +141,65 @@ def test_permute_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_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 43b9c36..77a7cb3 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -3,21 +3,23 @@ from firedrake import * from fuse import * import sympy as sp -from test_convert_to_fiat import create_cg1, helmholtz_solve, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg0, create_dg1 +from test_convert_to_fiat import create_cg1, construct_nd, construct_rt, create_cg2_tri, construct_cg3, create_dg1 import os -#with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): + +# with mock.patch.object(ElementTriple, 'make_dof_perms', new=dummy_dof_perms): 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()) + # oriented_mats_by_entity, flat_by_entity = cls._initialise_entity_dicts(cls.generate()) oriented_mats_by_entity = cls.make_entity_dense_matrices(*args) - #for key1, val1 in oriented_mats_by_entity.items(): + # 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 == 1 and key3 == 1: - # oriented_mats_by_entity[key1][key2][key3] = np.array([[0, 1],[1,0]]) + # oriented_mats_by_entity[key1][key2][key3] = np.array([[0, 1],[1,0]]) return oriented_mats_by_entity + def construct_nd2(tri=None): if tri is None: tri = polygon(3) @@ -26,28 +28,26 @@ def construct_nd2(tri=None): x = sp.Symbol("x") y = sp.Symbol("y") - xs = [DOF(L2Pairing(), ParameterisationKernel())] - + xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] 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)] + 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 - + 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) @@ -59,7 +59,6 @@ def construct_nd2_for_fiat(tri=None): 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) @@ -76,7 +75,7 @@ def construct_nd2_for_fiat(tri=None): 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_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) @@ -84,23 +83,21 @@ def construct_nd2_for_fiat(tri=None): 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 - + 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_for_fiat, "N1curl", 2), (construct_nd2, "N1curl", 2)]) def test_surface_const_nd(elem_gen, elem_code, deg): - with mock.patch.object(ElementTriple, 'dummy_function', new=dummy_dof_perms): cell = polygon(3) elem = elem_gen(cell) - ones = as_vector((0,1)) + ones = as_vector((0, 1)) - for n in range(2, 3): + 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()) @@ -109,8 +106,8 @@ def test_surface_const_nd(elem_gen, 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) @@ -118,7 +115,7 @@ def test_surface_const_nd(elem_gen, elem_code, deg): 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) @@ -129,12 +126,13 @@ 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) +@pytest.mark.xfail(reason="orientations") def test_surface_vec(): cell = polygon(3) rt_elem = construct_rt(cell) @@ -145,7 +143,7 @@ def test_surface_vec(): 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) @@ -155,13 +153,13 @@ 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) + # V2 = FunctionSpace(mesh, create_dg0(cell).to_ufl()) + # u = TestFunction(V2) res2 = assemble(dot(vec1, normal) * ds) else: V = FunctionSpace(mesh, "N1curl", 1) @@ -170,13 +168,13 @@ def test_surface_vec(): res2 = assemble(dot(vec1, normal) * ds) print(f"curl {n}: {res2}") - #assert np.allclose(0, res1) + assert np.allclose(0, res1) assert np.allclose(0, res2) -def test_interpolate_vs_project(V): +def interpolate_vs_project(V): mesh = V.mesh() - dim = mesh.geometric_dimension() + dim = mesh.geometric_dimension if dim == 2: x, y = SpatialCoordinate(mesh) elif dim == 3: @@ -207,14 +205,15 @@ def test_interpolate_vs_project(V): def test_degree2_interpolation(): cell = polygon(3) - elem = construct_nd2_for_fiat(cell) - mesh = UnitSquareMesh(1,1) + elem = construct_nd2(cell) + mesh = UnitSquareMesh(1, 1) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V = FunctionSpace(mesh, elem.to_ufl()) else: V = FunctionSpace(mesh, "N1curl", 2) - test_interpolate_vs_project(V) + interpolate_vs_project(V) + @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg2_tri, "CG", 2), (create_cg1, "CG", 1), @@ -224,20 +223,20 @@ def test_degree2_interpolation(): def test_interpolation(elem_gen, elem_code, deg): cell = polygon(3) elem = elem_gen(cell) - mesh = UnitSquareMesh(1,1) + 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) - test_interpolate_vs_project(V) + interpolate_vs_project(V) + def test_create_fiat_nd(): cell = polygon(3) nd = construct_nd2(cell) nd_fv = construct_nd2_for_fiat(cell) ref_el = cell.to_fiat() - sd = ref_el.get_spatial_dimension() deg = 2 from FIAT.nedelec import Nedelec @@ -246,7 +245,7 @@ def test_create_fiat_nd(): for f in fiat_elem.dual_basis(): print(f.pt_dict) - print("fiat: fuse's version") + print("fiat: fuse's version") for d in nd_fv.generate(): print(d.convert_to_fiat(ref_el, deg).pt_dict) @@ -255,18 +254,16 @@ def test_create_fiat_nd(): e = cell.edges()[0] g = S3.add_cell(cell).members()[2] print(g) - nodes = [d.convert_to_fiat(ref_el, deg) for d in dofs] - new_nodes = [d(g).convert_to_fiat(ref_el, deg) if d.trace_entity == e else d.convert_to_fiat(ref_el, deg) for d in dofs] + 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"{dofs[i]}: {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}") + 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}") print() - my_elem = nd.to_fiat() - #my_elem = nd_fv.to_fiat() - #breakpoint() - - + nd.to_fiat() + nd_fv.to_fiat()