diff --git a/Makefile b/Makefile index 734a8d3..e21d433 100644 --- a/Makefile +++ b/Makefile @@ -26,6 +26,12 @@ tests: @echo " Running all tests" @FIREDRAKE_USE_FUSE=1 python3 -m coverage run -p -m pytest -rx test +mini_tests: + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_2d_examples_docs.py + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_1d + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_orientations.py::test_surface_vec_rt + @FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_projection_convergence_3d\[construct_tet_ned-N1curl-1-0.8\] + coverage: @python3 -m coverage combine @python3 -m coverage report -m diff --git a/fuse/__init__.py b/fuse/__init__.py index e7ed061..f7e89d3 100644 --- a/fuse/__init__.py +++ b/fuse/__init__.py @@ -1,6 +1,7 @@ -from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex + +from fuse.cells import Point, Edge, polygon, line, make_tetrahedron, constructCellComplex, TensorProductPoint from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, diff_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group -from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, VectorKernel, PolynomialKernel, ComponentKernel +from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, VectorKernel, BarycentricPolynomialKernel, PolynomialKernel, ComponentKernel from fuse.triples import ElementTriple, DOFGenerator, immerse from fuse.traces import TrH1, TrGrad, TrHess, TrHCurl, TrHDiv from fuse.tensor_products import tensor_product diff --git a/fuse/cells.py b/fuse/cells.py index ab06eb9..8b1942c 100644 --- a/fuse/cells.py +++ b/fuse/cells.py @@ -151,6 +151,13 @@ def compute_scaled_verts(d, n): raise ValueError("Dimension {} not supported".format(d)) +def line(): + """ + Constructs the default 1D interval + """ + return Point(1, [Point(0), Point(0)], vertex_num=2) + + def polygon(n): """ Constructs the 2D default cell with n sides/vertices @@ -315,8 +322,8 @@ def __init__(self, d, edges=[], vertex_num=None, oriented=False, group=None, edg self.group = group if not group: self.group = self.compute_cell_group() - self.group = self.group.add_cell(self) + self.basis_group = self.compute_basis_group() def compute_attachments(self, n, points, orientations=None): """ @@ -373,6 +380,7 @@ def compute_cell_group(self): """ verts = self.ordered_vertices() v_coords = [self.get_node(v, return_coords=True) for v in verts] + n = len(verts) max_group = SymmetricGroup(n) edges = [edge.ordered_vertices() for edge in self.edges()] @@ -388,6 +396,19 @@ def compute_cell_group(self): break return fuse_groups.PermutationSetRepresentation(list(accepted_perms)) + def compute_basis_group(self): + if self.dimension == 0: + return fuse_groups.S1.add_cell(self) + bvs = self.basis_vectors() + bv_0 = bvs[0] + basis_group_members = [] + for bv in bvs: + for g in self.group.members(): + if np.allclose(np.array(g(bv_0)), np.array(bv)): + basis_group_members += [g.perm] + break + return fuse_groups.PermutationSetRepresentation(list(basis_group_members)).add_cell(self) + def get_spatial_dimension(self): return self.dimension @@ -475,12 +496,8 @@ def _subentity_traversal(self, sub_ents, min_ids): sub_ents = self.get_node(p)._subentity_traversal(sub_ents, min_ids) if dim > 1: connections = [(c.point.id, c.point.group.identity) for c in self.connections] - # if self.oriented: - # connections = self.permute_entities(self.oriented, dim - 1) # if self.dimension == 2: # connections = [connections[-1]] + connections[:-1] - # print([self.get_node(c[0]).id - min_ids[1] for c in connections]) - # print([c.point.id - min_ids[1] for c in self.connections]) for e, o in connections: p = self.get_node(e).orient(o) p_dim = p.get_spatial_dimension() @@ -502,6 +519,15 @@ def get_starter_ids(self): min_ids = [min(dimension) for dimension in structure] return min_ids + def local_id(self, node): + structure = [sorted(generation) for generation in nx.topological_generations(self.G)] + structure.reverse() + min_id = self.get_starter_ids() + for d in range(len(structure)): + if node.id in structure[d]: + return node.id - min_id[d] + raise ValueError("Node not found in cell") + def graph_dim(self): if self.oriented: dim = self.dimension + 1 @@ -544,6 +570,9 @@ def ordered_vertices(self, get_class=False): return self.oriented.permute(verts) return verts + def ordered_vertex_coords(self): + return [self.get_node(o, return_coords=True) for o in self.ordered_vertices()] + def d_entities_ids(self, d): return self.d_entities(d, get_class=False) @@ -639,13 +668,12 @@ def permute_entities(self, g, d): return reordered_entities - def basis_vectors(self, return_coords=True, entity=None, order=False): + def basis_vectors(self, return_coords=True, entity=None, order=False, norm=True): if not entity: entity = self self_levels = [generation for generation in nx.topological_generations(self.G)] vertices = entity.ordered_vertices() if self.dimension == 0: - # return [[] raise ValueError("Dimension 0 entities cannot have Basis Vectors") if self.oriented: # ordered_vertices() handles the orientation so we want to drop the orientation node @@ -659,7 +687,10 @@ def basis_vectors(self, return_coords=True, entity=None, order=False): for v in vertices[1:]: if return_coords: v_coords = self.attachment(top_level_node, v)() - sub = normalise(np.subtract(v_coords, v_0_coords)) + if norm: + sub = normalise(np.subtract(v_coords, v_0_coords)) + else: + sub = np.subtract(v_coords, v_0_coords) if not hasattr(sub, "__iter__"): basis_vecs.append((sub,)) else: @@ -816,12 +847,35 @@ def attachment(self, source, dst): return lambda *x: fold_reduce(attachments[0], *x) + def attachment_J_det(self, source, dst): + attachment = self.attachment(source, dst) + symbol_names = ["x", "y", "z"] + symbols = [] + if self.dim_of_node(dst) == 0: + return 1 + for i in range(self.dim_of_node(dst)): + symbols += [sp.Symbol(symbol_names[i])] + J = sp.Matrix(attachment(*symbols)).jacobian(sp.Matrix(symbols)) + return np.sqrt(abs(float(sp.det(J.T * J)))) + 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 cartesian_to_barycentric(self, pts): + verts = np.array(self.ordered_vertex_coords()) + v_0 = self.ordered_vertex_coords()[0] + bvs = np.array(self.basis_vectors(norm=False)) + bary_coords = [] + for pt in pts: + res = np.matmul(np.linalg.inv(bvs.T), np.array(pt - v_0)) + assert np.allclose(sum(bvs[i]*res[i] for i in range(len(bvs))), np.array(pt - v_0)) + bary_coords += [(1 - sum(res),) + tuple(res[i] for i in range(len(res)))] + assert np.allclose(np.array(sum(bary_coords[-1][i]*verts[i] for i in range(len(verts)))), pt) + return bary_coords + def cell_attachment(self, dst): if not isinstance(dst, int): raise ValueError @@ -874,6 +928,13 @@ def dict_id(self): def _from_dict(o_dict): return Point(o_dict["dim"], o_dict["edges"], oriented=o_dict["oriented"], cell_id=o_dict["id"]) + def equivalent(self, other): + if self.dimension != other.dimension: + return False + if set(self.ordered_vertex_coords()) != set(other.ordered_vertex_coords()): + return False + return self.get_topology() == other.get_topology() + class Edge(): """ @@ -897,7 +958,12 @@ def __call__(self, *x): if hasattr(self.attachment, '__iter__'): res = [] for attach_comp in self.attachment: - res.append(sympy_to_numpy(attach_comp, syms, x)) + if len(attach_comp.atoms(sp.Symbol)) <= len(x): + res.append(sympy_to_numpy(attach_comp, syms, x)) + else: + res_val = attach_comp.subs({syms[i]: x[i] for i in range(len(x))}) + res.append(res_val) + return tuple(res) return sympy_to_numpy(self.attachment, syms, x) return x @@ -929,11 +995,11 @@ def _from_dict(o_dict): class TensorProductPoint(): - def __init__(self, A, B, flat=False): + def __init__(self, A, B): self.A = A self.B = B self.dimension = self.A.dimension + self.B.dimension - self.flat = flat + self.flat = False def ordered_vertices(self): return self.A.ordered_vertices() + self.B.ordered_vertices() @@ -945,8 +1011,8 @@ def get_sub_entities(self): self.A.get_sub_entities() self.B.get_sub_entities() - def dimension(self): - return tuple(self.A.dimension, self.B.dimension) + def dim(self): + return (self.A.dimension, self.B.dimension) def d_entities(self, d, get_class=True): return self.A.d_entities(d, get_class) + self.B.d_entities(d, get_class) @@ -961,17 +1027,91 @@ def vertices(self, get_class=True, return_coords=False): return verts def to_ufl(self, name=None): - if self.flat: - return CellComplexToUFL(self, "quadrilateral") return TensorProductCell(self.A.to_ufl(), self.B.to_ufl()) def to_fiat(self, name=None): - if self.flat: - return CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name)) return CellComplexToFiatTensorProduct(self, name) def flatten(self): - return TensorProductPoint(self.A, self.B, True) + assert self.A.equivalent(self.B) + return FlattenedPoint(self.A, self.B) + + +class FlattenedPoint(Point, TensorProductPoint): + + def __init__(self, A, B): + self.A = A + self.B = B + self.dimension = self.A.dimension + self.B.dimension + self.flat = True + fuse_edges = self.construct_fuse_rep() + super().__init__(self.dimension, fuse_edges) + + def to_ufl(self, name=None): + return CellComplexToUFL(self, "quadrilateral") + + def to_fiat(self, name=None): + # TODO this should check if it actually is a hypercube + fiat = CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name)) + return fiat + + def construct_fuse_rep(self): + sub_cells = [self.A, self.B] + dims = (self.A.dimension, self.B.dimension) + + points = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells} + attachments = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells} + + for d in range(max(dims) + 1): + for cell in sub_cells: + if d <= cell.dimension: + sub_ent = cell.d_entities(d, get_class=True) + points[cell][d].extend(sub_ent) + for s in sub_ent: + attachments[cell][d].extend(s.connections) + + prod_points = list(itertools.product(*[points[cell][0] for cell in sub_cells])) + # temp = prod_points[1] + # prod_points[1] = prod_points[2] + # prod_points[2] = temp + point_cls = [Point(0) for i in range(len(prod_points))] + edges = [] + + # generate edges of tensor product result + for a in prod_points: + for b in prod_points: + # of all combinations of point, take those where at least one changes and at least one is the same + if any(a[i] == b[i] for i in range(len(a))) and any(a[i] != b[i] for i in range(len(sub_cells))): + # ensure if they change, that edge exists in the existing topology + if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i]._topology[1].values()) for i in range(len(sub_cells))]): + edges.append((a, b)) + # hasse level 1 + edge_cls1 = {e: None for e in edges} + for i in range(len(sub_cells)): + for (a, b) in edges: + a_idx = prod_points.index(a) + b_idx = prod_points.index(b) + if a[i] != b[i]: + a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] + b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0] + edge_cls1[(a, b)] = Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o), + Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)]) + edge_cls2 = [] + # hasse level 2 + for i in range(len(sub_cells)): + for (a, b) in edges: + if a[i] == b[i]: + x = sp.Symbol("x") + a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0] + if i == 0: + attach = (x,) + a_edge.attachment + else: + attach = a_edge.attachment + (x,) + edge_cls2.append(Edge(edge_cls1[(a, b)], attach, a_edge.o)) + return edge_cls2 + + def flatten(self): + return self class CellComplexToFiatSimplex(Simplex): @@ -1161,9 +1301,8 @@ def constructCellComplex(name): return polygon(3).to_ufl(name) # return ufc_triangle().to_ufl(name) elif name == "quadrilateral": - interval = Point(1, [Point(0), Point(0)], vertex_num=2) - return TensorProductPoint(interval, interval).flatten().to_ufl(name) - # return ufc_quad().to_ufl(name) + return TensorProductPoint(line(), line()).flatten().to_ufl(name) + # return firedrake_quad().to_ufl(name) # return polygon(4).to_ufl(name) elif name == "tetrahedron": # return ufc_tetrahedron().to_ufl(name) diff --git a/fuse/dof.py b/fuse/dof.py index 7c07f1c..3f0b6fb 100644 --- a/fuse/dof.py +++ b/fuse/dof.py @@ -152,7 +152,7 @@ def permute(self, g): def __call__(self, *args): return self.pt - def evaluate(self, Qpts, Qwts, basis_change, immersed, dim): + def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): return np.array([self.pt for _ in Qpts]).astype(np.float64), np.ones_like(Qwts), [[tuple()] for pt in Qpts] def _to_dict(self): @@ -189,10 +189,18 @@ def permute(self, g): def __call__(self, *args): return self.pt - def evaluate(self, Qpts, Qwts, basis_change, immersed, dim): + def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): + if len(value_shape) == 0: + comps = [[tuple()] for pt in Qpts] + else: + comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] + if isinstance(self.pt, int): - return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] + return Qpts, np.array([wt*self.pt for wt in Qwts]).astype(np.float64), comps + if not np.allclose(np.matmul(basis_change, self.pt), self.g(self.pt)): + breakpoint() if not immersed: + return Qpts, np.array([wt*np.matmul(self.pt, basis_change)for wt in Qwts]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] return Qpts, np.array([wt*immersed(np.matmul(self.pt, basis_change))for wt in Qwts]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] @@ -207,15 +215,75 @@ def _from_dict(obj_dict): return VectorKernel(tuple(obj_dict["pt"])) +class BarycentricPolynomialKernel(BaseKernel): + + def __init__(self, fn, g=None, symbols=[]): + if hasattr(fn, "__iter__"): + # if len(symbols) != 0 and any(not sp.sympify(fn[i]).as_poly() for i in range(len(fn))): + # raise ValueError("Function components must be able to be interpreted as a sympy polynomial") + self.fn = [sp.sympify(fn[i]).as_poly() for i in range(len(fn))] + self.shape = len(fn) + else: + if len(symbols) != 0 and not sp.sympify(fn).as_poly(): + raise ValueError("Function must be able to be interpreted as a sympy polynomial") + self.fn = sp.sympify(fn) + self.shape = 0 + self.g = g + self.syms = symbols + super(BarycentricPolynomialKernel, self).__init__() + + def __repr__(self): + return str(self.fn) + + def degree(self, interpolant_degree): + if self.shape != 0: + return max([self.fn[i].as_poly().total_degree() for i in range(self.shape)]) + interpolant_degree + if len(self.fn.free_symbols) == 0: # this should probably be removed + return interpolant_degree + return self.fn.as_poly().total_degree() + interpolant_degree + + def permute(self, g): + new_fn = self.fn + return BarycentricPolynomialKernel(new_fn, g=g, symbols=self.syms) + + def __call__(self, *args): + if self.shape == 0: + res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) + else: + res = [sympy_to_numpy(self.fn[i], self.syms, args[:len(self.syms)]) for i in range(self.shape)] + return res + + def evaluate(self, Qpts, bary_pts, Qwts, basis_change, immersed, dim, value_shape): + if len(value_shape) == 0: + comps = [[tuple()] for pt in Qpts] + else: + comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] + if not immersed or self.shape == 0: + return Qpts, np.array([wt*self(*self.g.permute(pt)) for pt, wt in zip(bary_pts, Qwts)]).astype(np.float64), [[(i,) for i in range(dim)] for pt in Qpts] + return Qpts, np.array([wt*immersed(self(*self.g.permute(pt))) for pt, wt in zip(bary_pts, Qwts)]).astype(np.float64), comps + + def _to_dict(self): + o_dict = {"fn": self.fn} + return o_dict + + def dict_id(self): + return "BarycentricPolynomialKernel" + + def _from_dict(obj_dict): + return BarycentricPolynomialKernel(obj_dict["fn"]) + + class PolynomialKernel(BaseKernel): - def __init__(self, fn, g=None, symbols=[], shape=0): - if len(symbols) != 0 and (shape != 0 and any(not sp.sympify(fn[i]).as_poly() for i in range(shape))) and not sp.sympify(fn).as_poly(): - raise ValueError("Function argument or its components must be able to be interpreted as a sympy polynomial") - if shape != 0: - self.fn = [sp.sympify(fn[i]).as_poly() for i in range(shape)] - self.shape = shape + def __init__(self, fn, g=None, symbols=[]): + if hasattr(fn, "__iter__"): + if len(symbols) != 0 and any(not sp.sympify(fn[i]).as_poly() for i in range(len(fn))): + raise ValueError("Function components must be able to be interpreted as a sympy polynomial") + self.fn = [sp.sympify(fn[i]).as_poly() for i in range(len(fn))] + self.shape = len(fn) else: + if len(symbols) != 0 and not sp.sympify(fn).as_poly(): + raise ValueError("Function must be able to be interpreted as a sympy polynomial") self.fn = sp.sympify(fn) self.shape = 0 self.g = g @@ -235,19 +303,23 @@ def degree(self, interpolant_degree): def permute(self, g): # new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))}) new_fn = self.fn - return PolynomialKernel(new_fn, g=g, symbols=self.syms, shape=self.shape) + return PolynomialKernel(new_fn, g=g, symbols=self.syms) def __call__(self, *args): if self.shape == 0: res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)]) else: - res = [] - for i in range(self.shape): - res += [sympy_to_numpy(self.fn[i], self.syms, args[:len(self.syms)])] + res = [sympy_to_numpy(self.fn[i], self.syms, args[:len(self.syms)]) for i in range(self.shape)] return res - def evaluate(self, Qpts, Qwts, basis_change, immersed, dim): - 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(dim)] for pt in Qpts] + def evaluate(self, Qpts, Qwts, basis_change, immersed, dim, value_shape): + if len(value_shape) == 0: + comps = [[tuple()] for pt in Qpts] + else: + comps = [[(i,) for v in value_shape for i in range(v)] for pt in Qpts] + if not immersed or self.shape == 0: + return Qpts, np.array([wt*self(*(np.matmul(pt, basis_change))) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), comps + return Qpts, np.array([wt*immersed(self(*(np.matmul(pt, basis_change)))) for pt, wt in zip(Qpts, Qwts)]).astype(np.float64), comps def _to_dict(self): o_dict = {"fn": self.fn} @@ -277,11 +349,9 @@ def permute(self, g): def __call__(self, *args): 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 evaluate(self, Qpts, Qwts, basis_change, immersed, dim): + def evaluate(self, Qpts, Qwts, immersed, dim): 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} @@ -336,16 +406,16 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non self.pairing = self.pairing.add_entity(cell) if self.target_space is None: self.target_space = space - if self.id is None and overall_id is not None: + if overall_id is not None: self.id = overall_id - if self.sub_id is None and generator_id is not None: + if generator_id is not None: self.sub_id = generator_id def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()): # TODO deriv dict needs implementing (currently {}) - return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree), {}, str(self)) + return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree, value_shape), {}, str(self)) - def to_quadrature(self, arg_degree): + def to_quadrature(self, arg_degree, value_shape): Qpts, Qwts = self.cell_defined_on.quadrature(self.kernel.degree(arg_degree)) Qwts = Qwts.reshape(Qwts.shape + (1,)) dim = self.cell_defined_on.get_spatial_dimension() @@ -355,7 +425,8 @@ def to_quadrature(self, arg_degree): basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) else: basis_change = np.eye(dim) - if self.immersed and isinstance(self.kernel, VectorKernel): + + if self.immersed and (isinstance(self.kernel, VectorKernel) or isinstance(self.kernel, BarycentricPolynomialKernel) or isinstance(self.kernel, PolynomialKernel)): def immersed(pt): basis = np.array(self.cell_defined_on.basis_vectors()).T basis_coeffs = np.matmul(np.linalg.inv(basis), np.array(pt)) @@ -363,30 +434,35 @@ def immersed(pt): return np.matmul(basis_coeffs, immersed_basis) else: immersed = self.immersed - pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change, immersed, self.cell.dimension) + + if isinstance(self.kernel, BarycentricPolynomialKernel): + bary_pts = self.cell_defined_on.cartesian_to_barycentric(Qpts) + pts, wts, comps = self.kernel.evaluate(Qpts, bary_pts, Qwts, basis_change, immersed, self.cell.dimension, value_shape) + else: + pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change, immersed, self.cell.dimension, value_shape) if self.immersed: - # need to compute jacobian from attachment. pts = np.array([self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in pts]) - # if self.pairing.orientation: - # immersion = self.target_space.tabulate(wts, self.pairing.entity.orient(self.pairing.orientation))[0] - # else: + J_det = self.cell.attachment_J_det(self.cell.id, self.cell_defined_on.id) + if not np.allclose(J_det, 1): + raise ValueError("Jacobian Determinant is not 1 did you do something wrong") immersion = self.target_space.tabulate(pts, self.cell_defined_on) - # Special case - force evaluation on different orientation of entity for construction of matrix transforms - if self.entity_o: - immersion = self.target_space.tabulate(wts, self.pairing.entity.orient(self.entity_o)) if isinstance(self.target_space, TrH1): - new_wts = wts + new_wts = wts * J_det else: - new_wts = np.outer(wts, immersion) + new_wts = np.outer(wts * J_det, immersion) + # shape is wrong for 2d face on tet + # if isinstance(self.kernel, BarycentricPolynomialKernel) and self.kernel.shape > 1: + # new_wts = np.array([self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in new_wts]) else: new_wts = wts # 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, new_wts, comps)} - if self.cell_defined_on.dimension == 2: - np.set_printoptions(linewidth=90, precision=4, suppress=True) - for key, val in pt_dict.items(): - print(np.array(key), ":", np.array([v[0] for v in val])) + # if self.cell_defined_on.dimension >= 2: + # print(self) + # np.set_printoptions(linewidth=90, precision=4, suppress=True) + # for key, val in pt_dict.items(): + # print(np.array(key), ":", np.array([v[0] for v in val])) return pt_dict def __repr__(self, fn="v"): diff --git a/fuse/groups.py b/fuse/groups.py index e63a9d3..fd8217e 100644 --- a/fuse/groups.py +++ b/fuse/groups.py @@ -65,9 +65,15 @@ def compute_perm(self, base_val=None): return val, val_list def numeric_rep(self): + """ Uses a standard formula to number permutations in the group. + For the case where this doesn't automatically number from 0..n (ie the group is not the full symmetry group), + a mapping is constructed on group creation""" identity = self.group.identity.perm.array_form m_array = self.perm.array_form - return orientation_value(identity, m_array) + val = orientation_value(identity, m_array) + if self.group.group_rep_numbering is not None: + return self.group.group_rep_numbering[val] + return val def __eq__(self, x): assert isinstance(x, GroupMemberRep) @@ -123,7 +129,8 @@ def __init__(self, perm_list, 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 @@ -144,9 +151,39 @@ def __init__(self, perm_list, cell=None): counter += 1 # self._members = sorted(self._members, key=lambda g: g.numeric_rep()) + self.group_rep_numbering = None + numeric_reps = [m.numeric_rep() for m in self.members()] + if sorted(numeric_reps) != list(range(len(numeric_reps))): + self.group_rep_numbering = {a: b for a, b in zip(sorted(numeric_reps), list(range(len(numeric_reps))))} + def add_cell(self, cell): return PermutationSetRepresentation(self.perm_list, cell=cell) + def conjugacy_class(self, g): + conj_class = set() + for x in self.members(): + res = ~x * g * x + conj_class.add(res) + return conj_class + + def cosets(self, subset): + # Divides current group by given subset + # can be modified to allow members of given subset not to exist in group self + seen = self.members().copy() + cosets = [] + while len(seen) > 0: + g = seen[0] + coset = [] + for h in subset.members(): + try: + coset += [g*h] + seen.remove(g*h) + except ValueError: + # member of subset not a member of superset + pass + cosets += [coset] + return cosets + def members(self, perm=False): if self.cell is None: raise ValueError("Group does not have a domain - members have not been calculated") @@ -221,7 +258,6 @@ def __init__(self, base_group, cell=None): self.generators = [] 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] @@ -229,10 +265,8 @@ def __init__(self, base_group, cell=None): counter = 0 for g in self.base_group.elements: if len(vertices) > g.size: - temp_perm = Permutation(g, size=len(vertices)) - reordered = temp_perm(vertices) - else: - reordered = g(vertices) + g = Permutation(g, size=len(vertices)) + reordered = g(vertices) A = np.c_[np.array(vertices, dtype=float), np.ones(len(vertices))] b = np.array(reordered, dtype=float) @@ -243,6 +277,11 @@ def __init__(self, base_group, cell=None): self.identity = p_rep counter += 1 + self.group_rep_numbering = None + numeric_reps = [m.numeric_rep() for m in self.members()] + if sorted(numeric_reps) != list(range(len(numeric_reps))): + self.group_rep_numbering = {a: b for a, b in zip(sorted(numeric_reps), list(range(len(numeric_reps))))} + # this order produces simpler generator lists # self.generators.reverse() @@ -257,13 +296,6 @@ def __init__(self, base_group, cell=None): else: self.cell = None - def conjugacy_class(self, g): - conj_class = set() - for x in self.members(): - res = ~x * g * x - conj_class.add(res) - return conj_class - def add_cell(self, cell): return GroupRepresentation(self.base_group, cell=cell) @@ -272,29 +304,6 @@ def size(self): assert len(self._members) == self.base_group.order() return self.base_group.order() - def members(self, perm=False): - if self.cell is None: - raise ValueError("Group does not have a domain - members have not been calculated") - if perm: - return [m.perm for m in self._members] - return self._members - - def transform_between_perms(self, perm1, perm2): - member_perms = self.members(perm=True) - perm1 = Permutation.from_sequence(perm1) - perm2 = Permutation.from_sequence(perm2) - assert perm1 in member_perms - assert perm2 in member_perms - return ~self.get_member(Permutation(perm1)) * self.get_member(Permutation(perm2)) - - # 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 # if len(remaining_members) == 0: diff --git a/fuse/spaces/element_sobolev_spaces.py b/fuse/spaces/element_sobolev_spaces.py index 2a3aa4f..b5964d9 100644 --- a/fuse/spaces/element_sobolev_spaces.py +++ b/fuse/spaces/element_sobolev_spaces.py @@ -1,4 +1,5 @@ from functools import total_ordering +from ufl.sobolevspace import H1, HDiv, HCurl, L2, H2 @total_ordering @@ -55,6 +56,9 @@ def __init__(self, cell): def __repr__(self): return "H1" + def to_ufl(self): + return H1 + class CellHDiv(ElementSobolevSpace): @@ -64,6 +68,9 @@ def __init__(self, cell): def __repr__(self): return "HDiv" + def to_ufl(self): + return HDiv + class CellHCurl(ElementSobolevSpace): @@ -73,6 +80,9 @@ def __init__(self, cell): def __repr__(self): return "HCurl" + def to_ufl(self): + return HCurl + class CellH2(ElementSobolevSpace): @@ -82,6 +92,9 @@ def __init__(self, cell): def __repr__(self): return "H2" + def to_ufl(self): + return H2 + class CellL2(ElementSobolevSpace): @@ -90,3 +103,6 @@ def __init__(self, cell): def __repr__(self): return "L2" + + def to_ufl(self): + return L2 diff --git a/fuse/spaces/polynomial_spaces.py b/fuse/spaces/polynomial_spaces.py index 422cf05..c7f8d36 100644 --- a/fuse/spaces/polynomial_spaces.py +++ b/fuse/spaces/polynomial_spaces.py @@ -1,13 +1,16 @@ from FIAT.polynomial_set import ONPolynomialSet +from FIAT.expansions import morton_index2, morton_index3 from FIAT.quadrature_schemes import create_quadrature from FIAT.reference_element import cell_to_simplex from FIAT import expansions, polynomial_set, reference_element from itertools import chain -from fuse.utils import tabulate_sympy, max_deg_sp_mat +from fuse.utils import tabulate_sympy, max_deg_sp_expr import sympy as sp import numpy as np from functools import total_ordering +morton_index = {2: morton_index2, 3: morton_index3} + @total_ordering class PolynomialSpace(object): @@ -47,7 +50,6 @@ def degree(self): return self.maxdegree def to_ON_polynomial_set(self, ref_el, k=None): - # how does super/sub degrees work here if not isinstance(ref_el, reference_element.Cell): ref_el = ref_el.to_fiat() sd = ref_el.get_spatial_dimension() @@ -56,18 +58,25 @@ def to_ON_polynomial_set(self, ref_el, k=None): shape = (sd,) else: shape = tuple() + base_ON = ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal") + indices = None if self.mindegree > 0: - base_ON = ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal") dimPmin = expansions.polynomial_dimension(ref_el, self.mindegree) dimPmax = expansions.polynomial_dimension(ref_el, self.maxdegree) if self.set_shape: indices = list(chain(*(range(i * dimPmin, i * dimPmax) for i in range(sd)))) else: indices = list(range(dimPmin, dimPmax)) - restricted_ON = base_ON.take(indices) - return restricted_ON - return ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal") + + if self.contains != self.maxdegree and self.contains != -1: + indices = [morton_index[sd](p, q) for p in range(self.contains + 1) for q in range(self.contains + 1)] + + if indices is None: + return base_ON + + restricted_ON = base_ON.take(indices) + return restricted_ON def __repr__(self): res = "" @@ -87,9 +96,7 @@ def __mul__(self, x): the sympy object on the right. This is due to Sympy's implementation of __mul__ not passing to this handler as it should. """ - if isinstance(x, sp.Symbol): - return ConstructedPolynomialSpace([x], [self]) - elif isinstance(x, sp.Matrix): + if isinstance(x, sp.Symbol) or isinstance(x, sp.Expr) or isinstance(x, sp.Matrix): return ConstructedPolynomialSpace([x], [self]) else: raise TypeError(f'Cannot multiply a PolySpace with {type(x)}') @@ -148,7 +155,7 @@ 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_expr(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)]) @@ -163,37 +170,49 @@ def to_ON_polynomial_set(self, ref_el): if not isinstance(ref_el, reference_element.Cell): ref_el = ref_el.to_fiat() k = max([s.maxdegree for s in self.spaces]) - space_poly_sets = [s.to_ON_polynomial_set(ref_el) for s in self.spaces] sd = ref_el.get_spatial_dimension() ref_el = cell_to_simplex(ref_el) - if all([w == 1 for w in self.weights]): - weighted_sets = space_poly_sets - # otherwise have to work on this through tabulation - Q = create_quadrature(ref_el, 2 * (k + 1)) - Qpts, Qwts = Q.get_points(), Q.get_weights() weighted_sets = [] - for (space, w) in zip(space_poly_sets, self.weights): + for (s, w) in zip(self.spaces, self.weights): + space = s.to_ON_polynomial_set(ref_el) + if s.set_shape: + shape = (sd,) + else: + shape = tuple() if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)): weighted_sets.append(space) else: - w_deg = max_deg_sp_mat(w) - Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, scale="orthonormal") - vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal") + if isinstance(w, sp.Expr): + w = sp.Matrix([[w]]) + vec = False + else: + vec = True + w_deg = max_deg_sp_expr(w) + Q = create_quadrature(ref_el, 2 * (k + w_deg + 1)) + Qpts, Qwts = Q.get_points(), Q.get_weights() + Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, shape, scale="orthonormal") + # vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal") space_at_Qpts = space.tabulate(Qpts)[(0,) * sd] Pkpw_at_Qpts = Pkpw.tabulate(Qpts)[(0,) * sd] tabulated_expr = tabulate_sympy(w, Qpts).T - scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] + if s.set_shape or vec: + scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] + else: + scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :] + scaled_at_Qpts = scaled_at_Qpts.squeeze() PkHw_coeffs = np.dot(np.multiply(scaled_at_Qpts, Qwts), Pkpw_at_Qpts.T) + if len(PkHw_coeffs.shape) == 1: + PkHw_coeffs = PkHw_coeffs.reshape(1, -1) weighted_sets.append(polynomial_set.PolynomialSet(ref_el, space.degree + w_deg, space.degree + w_deg, - vec_Pkpw.get_expansion_set(), + Pkpw.get_expansion_set(), PkHw_coeffs)) combined_sets = weighted_sets[0] for i in range(1, len(weighted_sets)): diff --git a/fuse/tensor_products.py b/fuse/tensor_products.py index cf16e24..4e0a486 100644 --- a/fuse/tensor_products.py +++ b/fuse/tensor_products.py @@ -43,8 +43,6 @@ def to_ufl(self): if self.flat: return FuseElement(self, self.cell.flatten().to_ufl()) ufl_sub_elements = [e.to_ufl() for e in self.sub_elements()] - # self.setup_matrices() - # breakpoint() return TensorProductElement(*ufl_sub_elements, cell=self.cell.to_ufl()) def flatten(self): diff --git a/fuse/traces.py b/fuse/traces.py index b5f9722..69e5f73 100644 --- a/fuse/traces.py +++ b/fuse/traces.py @@ -139,6 +139,7 @@ def tabulate(self, Qwts, trace_entity): # return result def manipulate_basis(self, basis): + breakpoint() return basis[0] def plot(self, ax, coord, trace_entity, **kwargs): diff --git a/fuse/triples.py b/fuse/triples.py index 32a7220..3444ff5 100644 --- a/fuse/triples.py +++ b/fuse/triples.py @@ -85,6 +85,13 @@ def setup_ids_and_nodes(self): def setup_matrices(self): # self.matrices_by_entity = self.make_entity_dense_matrices(self.ref_el, self.entity_ids, self.nodes, self.poly_set) matrices, entity_perms, pure_perm = self.make_dof_perms(self.ref_el, self.entity_ids, self.nodes, self.poly_set) + # new_matrices = matrices.copy() + # if not pure_perm: + # for j in range(4): + # for i, k in zip([0, 3, 4], [0, 3, 4]): + # new_matrices[2][j][i] = matrices[2][j][k] + # # new_matrices[2][j][k] = matrices[2][j][i] + # matrices = new_matrices reversed_matrices = self.reverse_dof_perms(matrices) if self.perm: self.pure_perm = pure_perm @@ -421,9 +428,10 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): ent_dofs = entity_associations[dim][e_id][dof_gen] ent_dofs_ids = np.array([self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs], dtype=int) # (dof_gen, ent_dofs) - total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if ed.id not in total_ent_dof_ids] + total_ent_dof_ids += [self.dof_id_to_fiat_id[ed.id] for ed in ent_dofs if self.dof_id_to_fiat_id[ed.id] not in total_ent_dof_ids] # dof_idx = [total_ent_dof_ids.index(id) for id in ent_dofs_ids] dof_gen_class = ent_dofs[0].generation + g_to_ent_id = {str(sub_g.perm.array_form): ent_id for ent_id, sub_g in zip(ent_dofs_ids, dof_gen_class[dim].g1.members())} 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 @@ -431,36 +439,66 @@ def make_dof_perms(self, ref_el, entity_ids, nodes, poly_set): bvs = np.array(e.basis_vectors()) new_bvs = np.array(e.orient(~g).basis_vectors()) basis_change = np.matmul(new_bvs, np.linalg.inv(bvs)) - if len(ent_dofs_ids) == basis_change.shape[0]: - sub_mat = basis_change - elif len(dof_gen_class[dim].g2.members()) == 2 and len(ent_dofs_ids) == 1: - # equivalently g1 trivial - sub_mat = ent_dofs[0].target_space.manipulate_basis(basis_change) + + cosets = dof_gen_class[dim].g1.cosets(e.basis_group) + if len(ent_dofs_ids) == len(basis_change): + value_change = basis_change + elif len(cosets[0]) < len(bvs): + if e.basis_group.size() > 2: + raise NotImplementedError("This logic is only valid for facets of 2 dimensions or less") + value_change = None + for coset in cosets: + if g in coset: + value_change = 1 + elif g*e.basis_group.members()[1] in coset: + value_change = -1 + if value_change is None: + raise ValueError("Basis group and generation group do not form the full symmetry group") + else: + value_change = basis_change + if len(cosets) == 1 or len(ent_dofs_ids) == len(basis_change): + sub_mat = value_change else: - # len(dof_gen_class[dim].g2.members()) == 2: - # case where value change is a restriction of the full transformation of the basis - value_change = ent_dofs[0].target_space.manipulate_basis(basis_change) + # check if this should be the cyclic variant sub_mat = np.kron((~g).matrix_form(), value_change) + new_ent_dofs_ids = [int(g_to_ent_id[str(sub_g.perm.array_form)]) for coset in cosets for sub_g in coset] + if not np.allclose(new_ent_dofs_ids, ent_dofs_ids): + # print("original", sub_mat) + # ent_dofs_ids = new_ent_dofs_ids + oriented_mats_by_entity[dim][e_id][val][np.ix_(new_ent_dofs_ids, new_ent_dofs_ids)] = sub_mat.copy() + else: + oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() + # print("added", oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)]) + # if len(ent_dofs_ids) == basis_change.shape[0]: + # sub_mat = basis_change + # elif len(dof_gen_class[dim].g2.members()) == 2 and len(ent_dofs_ids) == 1: + # # equivalently g1 trivial + # sub_mat = ent_dofs[0].target_space.manipulate_basis(basis_change) + # else: + # # len(dof_gen_class[dim].g2.members()) == 2: + # # case where value change is a restriction of the full transformation of the basis + # value_change = ent_dofs[0].target_space.manipulate_basis(basis_change) + # sub_mat = np.kron((~g).matrix_form(), value_change) # sub_mat = (~g).matrix_form() # elif len(ent_dofs_ids) != 1:# more dofs than dimension of g? # case for transforms where the basis vector is already included in the dof # sub_mat = np.kron((~g).matrix_form(), basis_change) # else: - # raise NotImplementedError("Unconsidered permuation case") - oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() - + # raise NotImplementedError("Unconsidered permutation case") 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)) elif dim < self.cell.dim(): # g in dof_gen_class[dim].g1.members() and # Permutation of DOF on the entity they are defined on sub_mat = (~g).matrix_form() oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() - elif len(dof_gen_class.keys()) == 1 and dim == self.cell.dim(): - # case for dofs defined on the cell and not immersed - sub_mat = (~g).matrix_form() - oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() + elif len(dof_gen_class.keys()) == 1 and dim == self.cell.dim() and len(ent_dofs_ids) == g.perm.size: + # case for dofs defined on the cell and not immersed - since they are interior the orientations don't matter + # sub_mat = (~g).matrix_form() + # oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = sub_mat.copy() + oriented_mats_by_entity[dim][e_id][val][np.ix_(ent_dofs_ids, ent_dofs_ids)] = np.eye(len(ent_dofs_ids)) else: # TODO what if an orientation is not in G1 + # also the case of 3 dofs inside a 3d shape 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() @@ -537,6 +575,21 @@ def reverse_dof_perms(self, matrices): reversed_mats[dim][e_id] = perms_copy return reversed_mats + def __add__(self, other): + """ Construct a new element triple by combining the degrees of freedom + This implementation does not make assertions about the properties + of the resulting element. + + Elements being adding must be defined over the same cell and have the same + value shape and mapping""" + assert self.cell == other.cell + assert self.spaces[0].set_shape == other.spaces[0].set_shape + assert str(self.spaces[1]) == str(other.spaces[1]) + + spaces = (self.spaces[0] + other.spaces[0], self.spaces[1], max([self.spaces[2], other.spaces[2]])) + + return ElementTriple(self.cell, spaces, self.DOFGenerator + other.DOFGenerator) + def _to_dict(self): o_dict = {"cell": self.cell, "spaces": self.spaces, "dofs": self.DOFGenerator} return o_dict @@ -574,21 +627,20 @@ def num_dofs(self): return self.dof_numbers def generate(self, cell, space, id_counter): - if self.ls is None: - self.ls = [] - for l_g in self.x: - i = 0 - for g in self.g1.members(): - generated = l_g(g) - if not isinstance(generated, list): - generated = [generated] - for dof in generated: - dof.add_context(self, cell, space, g, id_counter, i) - id_counter += 1 - i += 1 - self.ls.extend(generated) - self.dof_numbers = len(self.ls) - self.dof_ids = [dof.id for dof in self.ls] + self.ls = [] + for l_g in self.x: + i = 0 + for g in self.g1.members(): + generated = l_g(g) + if not isinstance(generated, list): + generated = [generated] + for dof in generated: + dof.add_context(self, cell, space, g, id_counter, i) + id_counter += 1 + i += 1 + self.ls.extend(generated) + self.dof_numbers = len(self.ls) + self.dof_ids = [dof.id for dof in self.ls] return self.ls def make_entity_ids(self): diff --git a/fuse/utils.py b/fuse/utils.py index fcfa3d9..0d0ee9a 100644 --- a/fuse/utils.py +++ b/fuse/utils.py @@ -29,7 +29,7 @@ def sympy_to_numpy(array, symbols, values): """ substituted = array.subs({symbols[i]: values[i] for i in range(len(values))}) - if len(array.atoms(sp.Symbol)) == len(values) and all(not isinstance(v, sp.Expr) for v in values): + 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) > 1: @@ -47,7 +47,7 @@ def tabulate_sympy(expr, pts): # expr: sp matrix expression in x,y,z for components of R^d # pts: n values in R^d # returns: evaluation of expr at pts - res = np.array(pts) + res = np.zeros((pts.shape[0],) + (expr.shape[-1],)) i = 0 syms = ["x", "y", "z"] for pt in pts: @@ -57,16 +57,20 @@ def tabulate_sympy(expr, pts): subbed = np.array(subbed).astype(np.float64) res[i] = subbed[0] i += 1 - final = res.squeeze() - return final + # final = res.squeeze() + return res -def max_deg_sp_mat(sp_mat): +def max_deg_sp_expr(sp_expr): degs = [] - for comp in sp_mat: - # only compute degree if component is a polynomial - if sp.sympify(comp).as_poly(): - degs += [sp.sympify(comp).as_poly().degree()] + if isinstance(sp_expr, sp.Matrix): + for comp in sp_expr: + # only compute degree if component is a polynomial + if sp.sympify(comp).as_poly(): + degs += [sp.sympify(comp).as_poly().total_degree()] + else: + if sp.sympify(sp_expr).as_poly(): + degs += [sp.sympify(sp_expr).as_poly().total_degree()] return max(degs) diff --git a/test/test_2d_examples_docs.py b/test/test_2d_examples_docs.py index a711741..11aa3f1 100644 --- a/test/test_2d_examples_docs.py +++ b/test/test_2d_examples_docs.py @@ -30,6 +30,32 @@ def construct_dg1(): return dg1 +def construct_dg0_integral(cell=None): + edge = Point(1, [Point(0), Point(0)], vertex_num=2) + xs = [DOF(L2Pairing(), VectorKernel(1))] + dg0 = ElementTriple(edge, (P0, CellL2, C0), DOFGenerator(xs, S1, S1)) + return dg0 + + +def construct_dg1_integral(cell=None): + edge = Point(1, [Point(0), Point(0)], vertex_num=2) + x = sp.Symbol("x") + xs = [DOF(L2Pairing(), PolynomialKernel((1/2)*(x + 1), symbols=(x,)))] + dg1 = ElementTriple(edge, (P1, CellL2, C0), DOFGenerator(xs, S2, S1)) + return dg1 + + +def construct_dg2_integral(cell=None): + edge = Point(1, [Point(0), Point(0)], vertex_num=2) + x = sp.Symbol("x") + xs = [DOF(L2Pairing(), PolynomialKernel((x/2)*(x + 1), symbols=(x,)))] + centre = [DOF(L2Pairing(), PolynomialKernel((1 - x**2), symbols=(x,)))] + + dofs = [DOFGenerator(xs, S2, S1), DOFGenerator(centre, S1, S1)] + dg2 = ElementTriple(edge, (PolynomialSpace(2), CellL2, C0), dofs) + return dg2 + + def plot_dg1(): dg1 = construct_dg1() dg1.plot() @@ -164,8 +190,6 @@ def construct_nd(tri=None): x = sp.Symbol("x") y = sp.Symbol("y") - # xs = [DOF(L2Pairing(), PointKernel(edge.basis_vectors()[0]))] - # xs = [DOF(L2Pairing(), PointKernel((1,)))] xs = [DOF(L2Pairing(), VectorKernel(1))] dofs = DOFGenerator(xs, S1, S2) @@ -203,6 +227,31 @@ def construct_nd_2nd_kind(tri=None): return ned +def construct_nd2_2nd_kind(tri=None): + if tri is None: + tri = polygon(3) + deg = 2 + edge = tri.edges()[0] + + s_0 = sp.Symbol("s_0") + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel(s_0*(2*s_0 - 1), symbols=(s_0,)))] + centre = [DOF(L2Pairing(), BarycentricPolynomialKernel(4*s_0*(1 - s_0), symbols=(s_0,)))] + dofs = [DOFGenerator(xs, S2, S2), DOFGenerator(centre, S1, S2)] + int_ned = ElementTriple(edge, (PolynomialSpace(deg, set_shape=True), CellHCurl, C0), dofs) + + xs = [immerse(tri, int_ned, TrHCurl)] + edge_dofs = DOFGenerator(xs, C3, S1) + + s_1 = sp.Symbol("s_1") + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel([-s_0, 1 - s_1], symbols=(s_0, s_1)))] + face_dofs = DOFGenerator(xs, C3, S2) + + nd = PolynomialSpace(deg, set_shape=True) + + ned = ElementTriple(tri, (nd, CellHCurl, C0), [edge_dofs, face_dofs]) + return ned + + def construct_bdm(tri=None): if tri is None: tri = polygon(3) @@ -223,13 +272,32 @@ def construct_bdm(tri=None): return rt +def construct_bdm_bary(tri=None): + if tri is None: + tri = polygon(3) + deg = 1 + edge = tri.edges()[0] + s_0 = sp.Symbol("s_0") + + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel(s_0, symbols=(s_0,)))] + dofs = DOFGenerator(xs, S2, S2) + int_rt = ElementTriple(edge, (PolynomialSpace(1, set_shape=True), CellHDiv, C0), dofs) + + xs = [immerse(tri, int_rt, TrHDiv)] + tri_dofs = DOFGenerator(xs, C3, S1) + + nd = PolynomialSpace(deg, set_shape=True) + + rt = ElementTriple(tri, (nd, CellHDiv, C0), [tri_dofs]) + return rt + + def construct_bdm2(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((x/2)*(x + 1), symbols=(x,)))] centre = [DOF(L2Pairing(), PolynomialKernel((1 - x**2), symbols=(x,)))] @@ -240,18 +308,17 @@ def construct_bdm2(tri=None): xs = [immerse(tri, int_rt, TrHDiv)] tri_dofs = DOFGenerator(xs, C3, S1) - phi_0 = [-1/6 - (np.sqrt(3)/6)*y, (-np.sqrt(3)/6) + (np.sqrt(3)/6)*x] - phi_1 = [-1/6 - (np.sqrt(3)/6)*y, (np.sqrt(3)/6) + (np.sqrt(3)/6)*x] - phi_2 = [1/3 - (np.sqrt(3)/6)*y, (np.sqrt(3)/6)*x] - xs = [DOF(L2Pairing(), PolynomialKernel(phi_0, symbols=(x, y), shape=2)), - DOF(L2Pairing(), PolynomialKernel(phi_1, symbols=(x, y), shape=2)), - DOF(L2Pairing(), PolynomialKernel(phi_2, symbols=(x, y), shape=2))] - interior = DOFGenerator(xs, S1, S1) + s_1 = sp.Symbol("s_1") + s_0 = sp.Symbol("s_0") + phi_0 = [1 - s_1, s_0] + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel(phi_0, symbols=(s_0, s_1)))] + interior = DOFGenerator(xs, C3, S1) - nd = PolynomialSpace(deg, set_shape=True) + space = PolynomialSpace(deg, set_shape=True) - rt = ElementTriple(tri, (nd, CellHDiv, C0), [tri_dofs, interior]) - return rt + bdm2 = ElementTriple(tri, (space, CellHDiv, C0), [tri_dofs, interior]) + dofs = bdm2.generate() + return bdm2 def test_nd_example(): @@ -271,6 +338,7 @@ def test_nd_example(): for dof in ned.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 + ned.to_fiat() def construct_rt(tri=None): diff --git a/test/test_3d_examples_docs.py b/test/test_3d_examples_docs.py index 803c361..d1739af 100644 --- a/test/test_3d_examples_docs.py +++ b/test/test_3d_examples_docs.py @@ -172,27 +172,58 @@ def construct_tet_rt2(cell=None, perm=None): rt_space = vec_Pd + (Pd.restrict(deg - 2, deg - 1))*M xs = [DOF(L2Pairing(), PolynomialKernel(1/3 - (1/2)*x + y/(2*np.sqrt(3)), symbols=(x, y)))] - dofs = DOFGenerator(xs, C3, S2) + dofs = DOFGenerator(xs, C3, S3) face_vec = ElementTriple(face, (rt_space, CellHDiv, "C0"), dofs) im_xs = [immerse(cell, face_vec, TrHDiv)] faces = DOFGenerator(im_xs, tet_faces, S1) - v_0 = np.array(cell.get_node(cell.ordered_vertices()[0], return_coords=True)) - v_1 = np.array(cell.get_node(cell.ordered_vertices()[1], return_coords=True)) - v_2 = np.array(cell.get_node(cell.ordered_vertices()[2], return_coords=True)) - v_3 = np.array(cell.get_node(cell.ordered_vertices()[3], return_coords=True)) - xs = [DOF(L2Pairing(), VectorKernel((v_2 - v_0)/2)), - DOF(L2Pairing(), VectorKernel((v_2 - v_1)/2)), - DOF(L2Pairing(), VectorKernel((v_2 - v_3)/2))] - interior = DOFGenerator(xs, S1, S4) + xs = [DOF(L2Pairing(), VectorKernel(cell.basis_vectors()[0]))] + interior = DOFGenerator(xs, cell.basis_group, S4) rt2 = ElementTriple(cell, (rt_space, CellHDiv, "C0"), [faces, interior]) return rt2 -def construct_tet_bdm(cell=None, perm=None): +def construct_tet_rt3(cell=None, perm=None): + if cell is None: + cell = make_tetrahedron() + face = cell.d_entities(2, get_class=True)[0] + deg = 3 + x = sp.Symbol("x") + y = sp.Symbol("y") + z = sp.Symbol("z") + M = sp.Matrix([[x, y, z]]) + + vec_Pd = PolynomialSpace(deg - 1, set_shape=True) + Pd = PolynomialSpace(deg - 1) + rt_space = vec_Pd + (Pd.restrict(deg - 2, deg - 1))*M + + s_0 = sp.Symbol("s_0") + s_1 = sp.Symbol("s_1") + vertex_dofs = [DOF(L2Pairing(), BarycentricPolynomialKernel(2*s_0**2 + 4*s_0*s_1 - 3*s_0 + 2*s_1**2 - 3*s_1 + 1, symbols=(s_0, s_1)))] + edge_dofs = [DOF(L2Pairing(), BarycentricPolynomialKernel(4*s_0*(1 - s_0 - s_1), symbols=(s_0, s_1)))] + dofs = [DOFGenerator(vertex_dofs, C3, S3), DOFGenerator(edge_dofs, C3, S3)] + face_vec = ElementTriple(face, (rt_space, CellHDiv, "C0"), dofs) + + im_xs = [immerse(cell, face_vec, TrHDiv)] + faces = DOFGenerator(im_xs, tet_faces, S1) + + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel([1 - s_0 - s_1, sp.Poly(0, s_0)], symbols=(s_0, s_1))), + DOF(L2Pairing(), BarycentricPolynomialKernel([sp.Poly(0, s_0), 1 - s_0 - s_1], symbols=(s_0, s_1)))] + interior = DOFGenerator(xs, C3, S4) + + rt2 = ElementTriple(cell, (rt_space, CellHDiv, "C0"), + [faces, interior]) + return rt2 + +# def test_rt3(): +# rt3 = construct_tet_rt3() +# breakpoint() + + +def construct_tet_bdm_old(cell=None, perm=None): if cell is None: cell = make_tetrahedron() face = cell.d_entities(2, get_class=True)[0] @@ -200,16 +231,132 @@ def construct_tet_bdm(cell=None, perm=None): x = sp.Symbol("x") y = sp.Symbol("y") - rt_space = PolynomialSpace(deg, set_shape=True) + space = PolynomialSpace(deg, set_shape=True) xs = [DOF(L2Pairing(), PolynomialKernel(1/3 - (1/2)*x + y/(2*np.sqrt(3)), symbols=(x, y)))] dofs = DOFGenerator(xs, C3, S2) - face_vec = ElementTriple(face, (rt_space, CellHDiv, "C0"), dofs) + face_vec = ElementTriple(face, (space, CellHDiv, "C0"), dofs) im_xs = [immerse(cell, face_vec, TrHDiv)] faces = DOFGenerator(im_xs, tet_faces, S1) - rt2 = ElementTriple(cell, (rt_space, CellHDiv, "C0"), [faces]) - return rt2 + bdm = ElementTriple(cell, (space, CellHDiv, "C0"), [faces]) + return bdm + + +def construct_tet_bdm(cell=None, perm=None): + if cell is None: + cell = make_tetrahedron() + face = cell.d_entities(2, get_class=True)[0] + deg = 1 + + space = PolynomialSpace(deg, set_shape=True) + + s_0 = sp.Symbol("s_0") + s_1 = sp.Symbol("s_1") + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel(1 - s_0 - s_1, symbols=(s_0, s_1)))] + dofs = DOFGenerator(xs, C3, S2) + face_vec = ElementTriple(face, (space, CellHDiv, "C0"), dofs) + im_xs = [immerse(cell, face_vec, TrHDiv)] + faces = DOFGenerator(im_xs, tet_faces, S1) + + bdm = ElementTriple(cell, (space, CellHDiv, "C0"), [faces]) + return bdm + +# def test_compare_bdms(): +# print("NORMAL") +# bdm1 = construct_tet_bdm() +# bdm1.to_fiat() +# print("BARY") +# bdm2 = construct_tet_bdm_bary() +# bdm2.to_fiat() + + +def construct_tet_bdm2(cell=None, perm=None): + if cell is None: + cell = make_tetrahedron() + face = cell.d_entities(2, get_class=True)[0] + deg = 2 + + space = PolynomialSpace(deg, set_shape=True) + + s_0 = sp.Symbol("s_0") + s_1 = sp.Symbol("s_1") + symbols = (s_0, s_1) + vertex_basis = 2*s_0**2 + 4*s_0*s_1 - 3*s_0 + 2*s_1**2 - 3*s_1 + 1 + edge_basis = 4*s_0*(-s_0 - s_1 + 1) + xs = [DOF(L2Pairing(), PolynomialKernel(vertex_basis, symbols=symbols))] + xs1 = [DOF(L2Pairing(), PolynomialKernel(edge_basis, symbols=symbols))] + dofs = [DOFGenerator(xs, C3, S3), DOFGenerator(xs1, C3, S3)] + face_vec = ElementTriple(face, (space, CellHDiv, "C0"), dofs) + im_xs = [immerse(cell, face_vec, TrHDiv)] + faces = DOFGenerator(im_xs, tet_faces, S1) + + bdm = ElementTriple(cell, (space, CellHDiv, "C0"), [faces]) + return bdm + + +def construct_tet_bdm(cell=None, perm=None): + if cell is None: + cell = make_tetrahedron() + face = cell.d_entities(2, get_class=True)[0] + deg = 1 + + space = PolynomialSpace(deg, set_shape=True) + + s_0 = sp.Symbol("s_0") + s_1 = sp.Symbol("s_1") + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel(1 - s_0 - s_1, symbols=(s_0, s_1)))] + dofs = DOFGenerator(xs, C3, S2) + face_vec = ElementTriple(face, (space, CellHDiv, "C0"), dofs) + im_xs = [immerse(cell, face_vec, TrHDiv)] + faces = DOFGenerator(im_xs, tet_faces, S1) + + bdm = ElementTriple(cell, (space, CellHDiv, "C0"), [faces]) + return bdm + +# def test_compare_bdms(): +# print("NORMAL") +# bdm1 = construct_tet_bdm() +# bdm1.to_fiat() +# print("BARY") +# bdm2 = construct_tet_bdm_bary() +# bdm2.to_fiat() + + +def construct_tet_bdm2(cell=None, perm=None): + if cell is None: + cell = make_tetrahedron() + face = cell.d_entities(2, get_class=True)[0] + deg = 2 + + space = PolynomialSpace(deg, set_shape=True) + + s_0 = sp.Symbol("s_0") + s_1 = sp.Symbol("s_1") + symbols = (s_0, s_1) + vertex_basis = 2*s_0**2 + 4*s_0*s_1 - 3*s_0 + 2*s_1**2 - 3*s_1 + 1 + edge_basis = 4*s_0*(-s_0 - s_1 + 1) + xs = [DOF(L2Pairing(), PolynomialKernel(vertex_basis, symbols=symbols))] + xs1 = [DOF(L2Pairing(), PolynomialKernel(edge_basis, symbols=symbols))] + dofs = [DOFGenerator(xs, C3, S3), DOFGenerator(xs1, C3, S3)] + face_vec = ElementTriple(face, (space, CellHDiv, "C0"), dofs) + im_xs = [immerse(cell, face_vec, TrHDiv)] + faces = DOFGenerator(im_xs, tet_faces, S1) + + s_0 = sp.Symbol("s_0") + s_1 = sp.Symbol("s_1") + s_2 = sp.Symbol("s_2") + zero = sp.Poly(0, (s_0, s_1, s_2)) + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel([-s_1 - s_2 + 1, s_0, s_0], symbols=(s_0, s_1, s_2))), + DOF(L2Pairing(), BarycentricPolynomialKernel([s_1, -s_0 - s_2 + 1, s_1], symbols=(s_0, s_1, s_2))), + DOF(L2Pairing(), BarycentricPolynomialKernel([s_2, s_2, -s_0 - s_1 + 1], symbols=(s_0, s_1, s_2))), + DOF(L2Pairing(), BarycentricPolynomialKernel([-s_1, s_0, zero], symbols=(s_0, s_1, s_2))), + DOF(L2Pairing(), BarycentricPolynomialKernel([-s_2, zero, s_0], symbols=(s_0, s_1, s_2))), + DOF(L2Pairing(), BarycentricPolynomialKernel([zero, - s_2, s_1], symbols=(s_0, s_1, s_2)))] + interior = DOFGenerator(xs, S1, S1) + + bdm2 = ElementTriple(cell, (space, CellHDiv, "C0"), [faces, interior]) + return bdm2 def construct_tet_ned(cell=None): @@ -244,7 +391,7 @@ def construct_tet_ned(cell=None): return ElementTriple(tet, (nd_space, CellHCurl, L2), [edge_dofs]) -def construct_tet_ned_2nd_kind(tet=None, perm=None): +def construct_tet_ned_2nd_kind_old(tet=None, perm=None): if tet is None: tet = make_tetrahedron() deg = 1 @@ -259,9 +406,27 @@ def construct_tet_ned_2nd_kind(tet=None, perm=None): int_ned1 = ElementTriple(edge, (PolynomialSpace(1, set_shape=True), CellHCurl, C0), dofs) xs = [immerse(tet, int_ned1, 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) + + ned = ElementTriple(tet, (nd_space, CellHCurl, C0), [edge_dofs]) + return ned + + +def construct_tet_ned_2nd_kind(tet=None, perm=None): + if tet is None: + tet = make_tetrahedron() + deg = 1 + edge = tet.edges()[0] + + vec_Pd = PolynomialSpace(deg, set_shape=True) + nd_space = vec_Pd + + s_0 = sp.Symbol("s_0") + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel(1 - s_0, symbols=(s_0,)))] + dofs = DOFGenerator(xs, S2, S2) + int_ned1 = ElementTriple(edge, (PolynomialSpace(1, set_shape=True), CellHCurl, C0), dofs) + + xs = [immerse(tet, int_ned1, TrHCurl)] edge_dofs = DOFGenerator(xs, tet_edges, S1) ned = ElementTriple(tet, (nd_space, CellHCurl, C0), [edge_dofs]) @@ -289,13 +454,11 @@ def construct_tet_ned2(tet=None, perm=None): dofs = DOFGenerator(xs, S2, S2) int_ned1 = ElementTriple(edge, (PolynomialSpace(1, set_shape=True), CellHCurl, C0), dofs) - v_0 = np.array(face.get_node(face.ordered_vertices()[0], return_coords=True)) - # v_1 = np.array(face.get_node(face.ordered_vertices()[1], return_coords=True)) + # v_0 = np.array(face.get_node(face.ordered_vertices()[0], return_coords=True)) + v_1 = np.array(face.get_node(face.ordered_vertices()[1], return_coords=True)) v_2 = np.array(face.get_node(face.ordered_vertices()[2], return_coords=True)) - xs = [DOF(L2Pairing(), VectorKernel((v_2 - v_0)/2))] - # breakpoint() - # xs = [DOF(L2Pairing(), VectorKernel((v_2 - v_1)/2))] - # xs = [DOF(L2Pairing(), VectorKernel((v_1 - v_0)/2)), DOF(L2Pairing(), VectorKernel((v_2 - v_0)/2)),] + + xs = [DOF(L2Pairing(), VectorKernel((v_2 - v_1)/2))] center_dofs = DOFGenerator(xs, S2, S3) face_vec = ElementTriple(face, (P1, CellHCurl, C0), center_dofs) im_xs = [immerse(tet, face_vec, TrH1)] @@ -310,6 +473,81 @@ def construct_tet_ned2(tet=None, perm=None): ned = ElementTriple(tet, (nd_space, CellHCurl, C0), [edge_dofs, face_dofs]) return ned +def construct_tet_ned_2nd_kind_2(tet=None, both=False): + if tet is None: + tet = make_tetrahedron() + deg = 2 + edge = tet.edges()[0] + face = tet.d_entities(2, get_class=True)[0] + vec_Pd = PolynomialSpace(deg, set_shape=True) + nd_space = vec_Pd + + s_0 = sp.Symbol("s_0") + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel(s_0*(2*s_0 - 1), symbols=(s_0,)))] + centre = [DOF(L2Pairing(), BarycentricPolynomialKernel(4*s_0*(1 - s_0), symbols=(s_0,)))] + dofs = [DOFGenerator(xs, S2, S2), DOFGenerator(centre, S1, S2)] + int_ned1 = ElementTriple(edge, (PolynomialSpace(1, set_shape=True), CellHCurl, C0), dofs) + edge_dofs = DOFGenerator([immerse(tet, int_ned1, TrHCurl)], tet_edges, S1) + + s_1 = sp.Symbol("s_1") + xs = [DOF(L2Pairing(), BarycentricPolynomialKernel([-s_0, -s_1], symbols=(s_0, s_1)))] + face_vec = ElementTriple(face, (P1, CellHCurl, C0), DOFGenerator(xs, C3, S2)) + + face_dofs = DOFGenerator([immerse(tet, face_vec, TrH1)], tet_faces, S1) + + ned = ElementTriple(tet, (nd_space, CellHCurl, C0), [edge_dofs, face_dofs]) + return ned + + +def construct_tet_ned3(tet=None, both=False): + if tet is None: + tet = make_tetrahedron() + deg = 3 + edge = tet.edges()[0] + face = tet.d_entities(2, get_class=True)[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 + + xs = [DOF(L2Pairing(), PolynomialKernel((x/2)*(x + 1), symbols=(x,)))] + centre = [DOF(L2Pairing(), PolynomialKernel((1 - x**2), symbols=(x,)))] + dofs = [DOFGenerator(xs, S2, S2), DOFGenerator(centre, S1, S2)] + int_ned1 = ElementTriple(edge, (PolynomialSpace(1, set_shape=True), CellHCurl, C0), dofs) + edge_dofs = DOFGenerator([immerse(tet, int_ned1, TrHCurl)], tet_edges, S1) + + phi_0 = [-1/6 - (np.sqrt(3)/6)*y, (-np.sqrt(3)/6) + (np.sqrt(3)/6)*x] + + # phi_1 = [-1/6 - (np.sqrt(3)/6)*y, (np.sqrt(3)/6) + (np.sqrt(3)/6)*x] + # phi_2 = [1/3 - (np.sqrt(3)/6)*y, (np.sqrt(3)/6)*x] + xs = [DOF(L2Pairing(), PolynomialKernel(phi_0, symbols=(x, y)))] + # DOF(L2Pairing(), PolynomialKernel(phi_1, symbols=(x, y), shape=2)), + # DOF(L2Pairing(), PolynomialKernel(phi_2, symbols=(x, y), shape=2))] + face_vec = ElementTriple(face, (P1, CellHCurl, C0), DOFGenerator(xs, C3, S1)) + # phi0 = [1/3 - (1/2)*x + y/(2*np.sqrt(3)), sp.Poly(0, (x, y))] + # xs = [DOF(L2Pairing(), PolynomialKernel(phi0, symbols=(x, y), shape=2))] + # if both: + # phi1 = [sp.Poly(0, (x, y)), 1/3 - (1/2)*x + y/(2*np.sqrt(3))] + # xs += [DOF(L2Pairing(), PolynomialKernel(phi1, symbols=(x, y), shape=2))] + # if both: + # face_vec = ElementTriple(face, (P1, CellHCurl, C0), DOFGenerator(xs, C3, S3)) + # else: + # face_vec = ElementTriple(face, (P1, CellHCurl, C0), DOFGenerator(xs, S3, S3)) + face_dofs = DOFGenerator([immerse(tet, face_vec, TrH1)], tet_faces, S1) + + xs = [DOF(L2Pairing(), VectorKernel([1, 0, 0])), + DOF(L2Pairing(), VectorKernel([0, 1, 0])), + DOF(L2Pairing(), VectorKernel([0, 0, 1]))] + int_dofs = DOFGenerator(xs, S1, S1) + + ned = ElementTriple(tet, (nd_space, CellHCurl, C0), [edge_dofs, face_dofs, int_dofs]) + return ned def test_plot_tet_ned2(): nd = construct_tet_ned2() @@ -327,7 +565,7 @@ def test_tet_rt2(): ls = rt2.generate() # TODO make this a proper test for dof in ls: - print(dof.to_quadrature(1)) + print(dof.to_quadrature(1, value_shape=(2,))) rt2.to_fiat() @@ -352,6 +590,38 @@ def test_tet_nd(): nd1.to_fiat() +def construct_V(elem): + nodes = elem.nodes + ref_el = elem.ref_el + entity_ids = elem.entity_ids + poly_set = elem.poly_set + from FIAT.dual_set import DualSet + + dual = DualSet(nodes, ref_el, entity_ids) + + old_coeffs = poly_set.get_coeffs() + dualmat = dual.to_riesz(poly_set) + + shp = dualmat.shape + A = dualmat.reshape((shp[0], -1)) + B = old_coeffs.reshape((shp[0], -1)) + V = np.dot(A, np.transpose(B)) + return V + + +def test_tet_nd3(): + nd3 = construct_tet_ned3() + nd3.to_fiat() + # nd3 = construct_tet_ned_2nd_kind_2() + # bdm2 = construct_tet_bdm2() + # nd3.to_ufl() + # nd3.to_fiat() + # V_nd = construct_V(nd3) + # bdm2.to_ufl() + # bdm2.to_fiat() + # V_bdm = construct_V(bdm2) + + def plot_tet_ned(): ned = construct_tet_ned() ned.plot(filename="tet_nd2_fiat.png") diff --git a/test/test_algebra.py b/test/test_algebra.py new file mode 100644 index 0000000..5304927 --- /dev/null +++ b/test/test_algebra.py @@ -0,0 +1,36 @@ +from fuse import * +from firedrake import * +import numpy as np +import sympy as sp +from test_convert_to_fiat import create_cg2_tri, construct_cg3 + + +def construct_bubble(cell=None): + if cell is None: + cell = polygon(3) + x = sp.Symbol("x") + y = sp.Symbol("y") + f = (3*np.sqrt(3)/4)*(y + np.sqrt(3)/3)*(np.sqrt(3)*x + y - 2*np.sqrt(3)/3)*(-np.sqrt(3)*x + y - 2*np.sqrt(3)/3) + space = PolynomialSpace(3).restrict(0, 0)*f + xs = [DOF(DeltaPairing(), PointKernel((0, 0)))] + bubble = ElementTriple(cell, (space, CellL2, L2), DOFGenerator(xs, S1, S1)) + return bubble + + +def test_bubble(): + mesh = UnitTriangleMesh() + x = SpatialCoordinate(mesh) + + tri = polygon(3) + bub = construct_bubble(tri) + cg2 = create_cg2_tri(tri) + p2b3 = bub + cg2 + V = FunctionSpace(mesh, p2b3.to_ufl()) + W = FunctionSpace(mesh, construct_cg3().to_ufl()) + + bubble_func = 27*x[0]*x[1]*(1-x[0]-x[1]) + u = project(bubble_func, V) + exact = Function(W) + exact.interpolate(bubble_func, W) + # make sure that these are the same + assert sqrt(assemble((u-exact)*(u-exact)*dx)) < 1e-14 diff --git a/test/test_cells.py b/test/test_cells.py index ced8e36..10ce831 100644 --- a/test/test_cells.py +++ b/test/test_cells.py @@ -1,9 +1,9 @@ from fuse import * from firedrake import * -from fuse.cells import ufc_triangle, ufc_tetrahedron +from fuse.cells import ufc_triangle import pytest import numpy as np -from FIAT.reference_element import default_simplex, ufc_simplex +from FIAT.reference_element import default_simplex from test_convert_to_fiat import helmholtz_solve @@ -35,6 +35,29 @@ def test_basis_vectors(C): assert len(bv_ids) == len(bv_coords) +def test_basis_group(C): + if C.dimension == 0: + assert C.basis_group.size() == 1 + else: + bv_coords = C.basis_vectors(return_coords=True) + bv_0 = bv_coords[0] + for i, g in enumerate(C.basis_group.members()): + assert np.allclose(np.array(bv_coords[i]), np.array(g(bv_0))) + if C.dimension == 2: + for i, g in enumerate(C.basis_group.members()): + bvs = np.array(C.basis_vectors()) + new_bvs = np.array(C.orient(g).basis_vectors()) + basis_change = np.matmul(np.linalg.inv(new_bvs), bvs) + assert np.allclose(np.array(bv_coords[i]), np.array(np.matmul(basis_change, bv_0))) + + +def test_cosets(): + cell = polygon(3) + # g = cell.basis_group.members()[1] + conj = cell.group.cosets(cell.basis_group) + print(conj) + + def test_sub_basis_vectors(): cell = polygon(3) @@ -179,117 +202,11 @@ def test_comparison(): # print(tensor_product1 >= tensor_product1) -@pytest.mark.parametrize(["cell"], [(ufc_triangle(),), (polygon(3),)]) -def test_connectivity(cell): - cell = cell.to_fiat() - for dim0 in range(cell.get_spatial_dimension()+1): - connectivity = cell.get_connectivity()[(dim0, 0)] - topology = cell.get_topology()[dim0] - assert len(connectivity) == len(topology) - - assert all(connectivity[i] == t for i, t in topology.items()) - - -def test_tensor_connectivity(): - from test_2d_examples_docs import construct_cg1 - A = construct_cg1() - B = construct_cg1() - cell = tensor_product(A, B).cell - cell = cell.to_fiat() - for dim0 in [(0, 0), (1, 0), (0, 1), (1, 1)]: - connectivity = cell.get_connectivity()[(dim0, (0, 0))] - topology = cell.get_topology()[dim0] - assert len(connectivity) == len(topology) - - assert all(connectivity[i] == t for i, t in topology.items()) - - -@pytest.mark.parametrize(["cell"], [(ufc_triangle(),), (polygon(3),), (make_tetrahedron(), ), (make_tetrahedron(), )]) -def test_new_connectivity(cell): - cell = cell.to_fiat() - for dim0 in range(cell.get_dimension() + 1): - connectivity = cell.get_connectivity()[(dim0, 0)] - topology = cell.get_topology()[dim0] - assert len(connectivity) == len(topology) - for i, t in topology.items(): - print(connectivity[i]) - print(t) - assert all(connectivity[i] == t for i, t in topology.items()) - - -def test_compare_tris(): - fuse_tet = polygon(3) - ufc_tet = ufc_triangle() - fiat_tet = ufc_simplex(2) - - print(fiat_tet.get_topology()) - print(fuse_tet.get_topology()) - print(ufc_tet.get_topology()) - fiat_connectivity = fiat_tet.get_connectivity() - fuse_connectivity = fuse_tet.to_fiat().get_connectivity() - ufc_connectivity = ufc_tet.to_fiat().get_connectivity() - _dim = fiat_tet.get_dimension() - print("fiat") - print(make_entity_cone_lists(fiat_tet)) - for dim0 in range(_dim): - connectivity = fiat_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse") - print(make_entity_cone_lists(fuse_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = fuse_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse ufc") - print(make_entity_cone_lists(ufc_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = ufc_connectivity[(dim0+1, dim0)] - print(connectivity) - - -def test_compare_tets(): - tet = make_tetrahedron() - # perm = tet.group.get_member([1, 2, 0, 3]) - fuse_tet = tet - ufc_tet = ufc_tetrahedron() - fiat_tet = ufc_simplex(3) - # breakpoint() - print(fiat_tet.get_topology()) - print(fuse_tet.get_topology()) - print(ufc_tet.get_topology()) - fiat_connectivity = fiat_tet.get_connectivity() - fuse_connectivity = fuse_tet.to_fiat().get_connectivity() - ufc_connectivity = ufc_tet.to_fiat().get_connectivity() - _dim = fiat_tet.get_dimension() - print("fiat") - print(make_entity_cone_lists(fiat_tet)) - for dim0 in range(_dim): - connectivity = fiat_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse") - print(make_entity_cone_lists(fuse_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = fuse_connectivity[(dim0+1, dim0)] - print(connectivity) - print("fuse ufc") - print(make_entity_cone_lists(ufc_tet.to_fiat())) - for dim0 in range(_dim): - connectivity = ufc_connectivity[(dim0+1, dim0)] - print(connectivity) - - -def make_entity_cone_lists(fiat_cell): - _dim = fiat_cell.get_dimension() - _connectivity = fiat_cell.connectivity - _list = [] - _offset_list = [0 for _ in _connectivity[(0, 0)]] # vertices have no cones - _offset = 0 - _n = 0 # num. of entities up to dimension = _d - for _d in range(_dim): - _n1 = len(_offset_list) - for _conn in _connectivity[(_d + 1, _d)]: - _list += [_c + _n for _c in _conn] # These are indices into cell_closure[some_cell] - _offset_list.append(_offset) - _offset += len(_conn) - _n = _n1 - _offset_list.append(_offset) - return _list, _offset_list +def test_self_equality(C): + assert C == C + + +@pytest.mark.parametrize(["A", "B", "res"], [(ufc_triangle(), polygon(3), False), + (line(), line(), True),]) +def test_equivalence(A, B, res): + assert A.equivalent(B) == res diff --git a/test/test_convert_to_fiat.py b/test/test_convert_to_fiat.py index bee8650..07fde7b 100644 --- a/test/test_convert_to_fiat.py +++ b/test/test_convert_to_fiat.py @@ -5,8 +5,8 @@ from firedrake import * from sympy.combinatorics import Permutation from FIAT.quadrature_schemes import create_quadrature -from test_2d_examples_docs import construct_cg1, construct_nd, construct_rt, construct_cg3 -from test_3d_examples_docs import construct_tet_rt, construct_tet_rt2, construct_tet_ned, construct_tet_ned_2nd_kind, construct_tet_bdm, construct_tet_ned2, construct_tet_cg4 +from test_2d_examples_docs import construct_cg1, construct_nd, construct_rt, construct_cg3, construct_dg0_integral, construct_dg1_integral, construct_dg2_integral +from test_3d_examples_docs import construct_tet_rt, construct_tet_rt2, construct_tet_ned, construct_tet_ned_2nd_kind, construct_tet_ned_2nd_kind_2, construct_tet_bdm, construct_tet_bdm2, construct_tet_ned2, construct_tet_cg4 from test_polynomial_space import flatten from element_examples import CR_n import os @@ -110,13 +110,11 @@ def create_cg1(cell): def create_cg1_quad(): deg = 1 - cell = polygon(4) - # cell = constructCellComplex("quadrilateral").cell_complex - - vert_dg = create_dg0(cell.vertices()[0]) + cell = TensorProductPoint(line(), line()).flatten() + print(cell, type(cell)) + vert_dg = create_dg1(cell.vertices()[0]) xs = [immerse(cell, vert_dg, TrH1)] - - Pk = PolynomialSpace(deg, deg + 1) + Pk = PolynomialSpace(deg + 1, deg) cg = ElementTriple(cell, (Pk, CellL2, C0), DOFGenerator(xs, get_cyc_group(len(cell.vertices())), S1)) return cg @@ -360,6 +358,9 @@ def test_entity_perms(elem_gen, cell): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(create_cg1, "CG", 1), (create_dg1, "DG", 1), + pytest.param(construct_dg0_integral, "DG", 0, marks=pytest.mark.xfail(reason='Passes locally, fails in CI, probably same as dg2')), + (construct_dg1_integral, "DG", 1), + (construct_dg2_integral, "DG", 2), pytest.param(create_dg2, "DG", 2, marks=pytest.mark.xfail(reason='Need to update TSFC in CI')), (create_cg2, "CG", 2) ]) @@ -503,6 +504,7 @@ def helmholtz_solve(V, mesh): def poisson_solve(r, elem, parameters={}, quadrilateral=False): # Create mesh and define function space m = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) + x = SpatialCoordinate(m) V = FunctionSpace(m, elem) @@ -523,6 +525,29 @@ def poisson_solve(r, elem, parameters={}, quadrilateral=False): return sqrt(assemble(inner(u - f, u - f) * dx)) +def run_test_original(r, elem_code, deg, parameters={}, quadrilateral=False): + # Create mesh and define function space + m = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) + + x = SpatialCoordinate(m) + V = FunctionSpace(m, elem_code, deg) + # Define variational problem + u = Function(V) + v = TestFunction(V) + a = inner(grad(u), grad(v)) * dx + + bcs = [DirichletBC(V, Constant(0), 3), + DirichletBC(V, Constant(42), 4)] + + # Compute solution + solve(a == 0, u, solver_parameters=parameters, bcs=bcs) + + f = Function(V) + f.interpolate(42*x[1]) + + return sqrt(assemble(inner(u - f, u - f) * dx)) + + @pytest.mark.parametrize(['params', 'elem_gen'], [(p, d) for p in [{}, {'snes_type': 'ksponly', 'ksp_type': 'preonly', 'pc_type': 'lu'}] @@ -534,7 +559,9 @@ def test_poisson_analytic(params, elem_gen): @pytest.mark.parametrize(['elem_gen'], - [(create_cg1_quad_tensor,), pytest.param(create_cg1_quad, marks=pytest.mark.xfail(reason='Need to allow generation on tensor product quads'))]) + [(create_cg1_quad_tensor,), + pytest.param(create_cg1_quad, marks=pytest.mark.xfail(reason='Issue with cell/mesh')) + ]) def test_quad(elem_gen): elem = elem_gen() r = 0 @@ -542,10 +569,6 @@ def test_quad(elem_gen): assert (poisson_solve(r, ufl_elem, parameters={}, quadrilateral=True) < 1.e-9) -def test_non_tensor_quad(): - create_cg1_quad() - - def project(U, mesh, func): f = assemble(interpolate(func, U)) @@ -604,7 +627,9 @@ def test_project_3d(elem_gen, elem_code, deg): (construct_tet_ned, "N1curl", 1, 0.8), (construct_tet_ned2, "N1curl", 2, 1.8), (construct_tet_ned_2nd_kind, "N2curl", 1, 1.8), - (construct_tet_bdm, "BDM", 1, 1.8)]) + (construct_tet_ned_2nd_kind_2, "N2curl", 2, 2.8), + (construct_tet_bdm, "BDM", 1, 1.8), + (construct_tet_bdm2, "BDM", 2, 2.8)]) def test_projection_convergence_3d(elem_gen, elem_code, deg, conv_rate): cell = make_tetrahedron() elem = elem_gen(cell) @@ -675,7 +700,9 @@ def test_const_vec(elem_gen, elem_code, deg, conv_rate): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_tet_ned2, "N1curl", 2), - (construct_tet_rt2, "RT", 2)]) + (construct_tet_rt2, "RT", 2), + (construct_tet_bdm2, "BDM", 2), + (construct_tet_ned_2nd_kind_2, "N2curl", 2)]) def test_linear_vec(elem_gen, elem_code, deg): cell = make_tetrahedron() elem = elem_gen(cell) @@ -711,7 +738,9 @@ def test_linear_vec(elem_gen, elem_code, deg): @pytest.mark.parametrize("elem_gen,elem_code,deg", [(construct_tet_rt, "RT", 1), (construct_tet_rt2, "RT", 2), (construct_tet_ned, "N1curl", 1), - (construct_tet_ned2, "N1curl", 2)]) + (construct_tet_ned2, "N1curl", 2), + (construct_tet_ned_2nd_kind_2, "N2curl", 2) + ]) def test_vec_two_tet(elem_gen, elem_code, deg): cell = make_tetrahedron() elem = elem_gen(cell) @@ -732,17 +761,19 @@ def vec(mesh, array=False): error_gs = [] error_row_lists = [] for g in group: + # mesh = UnitTetrahedronMesh() mesh = TwoTetMesh(perm=g) print(mesh.entity_orientations) if bool(os.environ.get("FIREDRAKE_USE_FUSE", 0)): V2 = FunctionSpace(mesh, elem.to_ufl()) + # print(elem.matrices[2][0][mesh.entity_orientations[1][10]][np.ix_([12, 13], [12, 13])]) res2 = assemble(interpolate(vec(mesh), V2)) CG3 = VectorFunctionSpace(mesh, create_cg3_tet(cell).to_ufl()) res3 = assemble(interpolate(res2, CG3)) error_rows = [] for i in range(res3.dat.data.shape[0]): if not np.allclose(res3.dat.data[i], vec(mesh, array=True)): - print(res3.dat.data[i]) + # print(res3.dat.data[i]) error_gs += [g] error_rows += [i] error_row_lists += [error_rows] @@ -760,7 +791,7 @@ def vec(mesh, array=False): (construct_tet_cg4, "CG", 4, 1e-13), (construct_tet_rt2, "RT", 2, 1e-13), (construct_tet_ned2, "N1curl", 2, 1e-13)]) -def test_const_two_tet(elem_gen, elem_code, deg, max_err): +def test_project_two_tet(elem_gen, elem_code, deg, max_err): cell = make_tetrahedron() # elem_perms = elem_gen(cell, perm=True) elem_mats = elem_gen(cell, perm=False) diff --git a/test/test_dofs.py b/test/test_dofs.py index ec7424f..bf29ad4 100644 --- a/test/test_dofs.py +++ b/test/test_dofs.py @@ -220,6 +220,6 @@ def test_generate_quadrature(): print("fiat", d.pt_dict) print() for d in elem.generate(): - print("fuse", d.to_quadrature(degree)) + print("fuse", d.to_quadrature(degree, value_shape=(2,))) elem.to_fiat() diff --git a/test/test_groups.py b/test/test_groups.py index 1b26787..0c20703 100644 --- a/test/test_groups.py +++ b/test/test_groups.py @@ -117,3 +117,10 @@ def test_perm_mat_conversion(): mat_form = g.matrix_form() array_form = perm_matrix_to_perm_array(mat_form) assert np.allclose(g.perm.array_form, array_form) + + +def test_numeric_reps(): + cell = polygon(4) + rot4 = get_cyc_group(4).add_cell(cell) + + assert sorted([m.numeric_rep() for m in rot4.members()]) == list(range(len(rot4.members()))) diff --git a/test/test_orientations.py b/test/test_orientations.py index b82edbe..319b7c3 100644 --- a/test/test_orientations.py +++ b/test/test_orientations.py @@ -3,7 +3,7 @@ from fuse import * import numpy as np import sympy as sp -from test_2d_examples_docs import construct_cg3, construct_nd, construct_rt, construct_nd_2nd_kind, construct_bdm, construct_bdm2 +from test_2d_examples_docs import construct_cg3, construct_nd, construct_rt, construct_nd_2nd_kind, construct_nd2_2nd_kind, construct_bdm, construct_bdm2, construct_bdm_bary from test_convert_to_fiat import create_cg1, create_cg2_tri, create_dg1 import os @@ -280,8 +280,10 @@ def test_convergence(elem_gen, elem_code, deg, conv_rate): (construct_nd_2nd_kind, "N2curl", 1, 1.8), (construct_rt2, "RT", 2, 1.8), (construct_bdm, "BDM", 1, 1.8), + (construct_bdm_bary, "BDM", 1, 1.8), (construct_bdm2, "BDM", 2, 2.8), - (construct_nd2, "N1curl", 2, 1.8)]) + (construct_nd2, "N1curl", 2, 1.8), + (construct_nd2_2nd_kind, "N2curl", 2, 2.8),]) def test_convergence_vector(elem_gen, elem_code, deg, conv_rate): cell = polygon(3) elem = elem_gen(cell) diff --git a/test/test_perms.py b/test/test_perms.py index f00aed6..82a3cb9 100644 --- a/test/test_perms.py +++ b/test/test_perms.py @@ -1,9 +1,8 @@ from fuse import * -from test_convert_to_fiat import create_cg1, create_dg1, create_cg2 +from test_convert_to_fiat import create_cg1, create_dg1, create_cg2, construct_nd from test_2d_examples_docs import construct_cg3 import pytest import numpy as np -import sympy as sp vert = Point(0) edge = Point(1, [Point(0), Point(0)], vertex_num=2) @@ -48,24 +47,7 @@ def test_basic_perms(cell): @pytest.mark.parametrize("cell", [tri]) def test_nd_perms(cell): - deg = 1 - edge = cell.edges(get_class=True)[0] - x = sp.Symbol("x") - y = sp.Symbol("y") - - xs = [DOF(L2Pairing(), PolynomialKernel(edge.basis_vectors()[0]))] - dofs = DOFGenerator(xs, S1, S2) - int_ned = ElementTriple(edge, (P1, CellHCurl, C0), dofs) - - xs = [immerse(cell, int_ned, TrHCurl)] - tri_dofs = DOFGenerator(xs, C3, S3) - - M = sp.Matrix([[y, -x]]) - vec_Pk = PolynomialSpace(deg - 1, set_shape=True) - Pk = PolynomialSpace(deg - 1) - nd = vec_Pk + (Pk.restrict(deg - 2, deg - 1))*M - - ned = ElementTriple(cell, (nd, CellHCurl, C0), [tri_dofs]) + ned = construct_nd(cell) ned.to_fiat() for i, mat in ned.matrices[2][0].items(): print(i) diff --git a/test/test_polynomial_space.py b/test/test_polynomial_space.py index e7dad66..baf8120 100644 --- a/test/test_polynomial_space.py +++ b/test/test_polynomial_space.py @@ -41,6 +41,7 @@ def test_restriction(): res_on_set = restricted.to_ON_polynomial_set(cell) P3_on_set = P3.to_ON_polynomial_set(cell) + assert res_on_set.get_num_members() < P3_on_set.get_num_members() not_restricted = P3.restrict(0, 3) @@ -48,6 +49,16 @@ def test_restriction(): assert not_restricted.mindegree == 0 +def test_square_space(): + cell = polygon(3) + q2 = PolynomialSpace(3, 1) + + q2_on_set = q2.to_ON_polynomial_set(cell) + P3_on_set = P3.to_ON_polynomial_set(cell) + + assert q2_on_set.get_num_members() < P3_on_set.get_num_members() + + @pytest.mark.parametrize("deg", [1, 2, 3, 4]) def test_complete_space(deg): cell = polygon(3) diff --git a/test/test_tensor_prod.py b/test/test_tensor_prod.py index 7c52720..15e15a7 100644 --- a/test/test_tensor_prod.py +++ b/test/test_tensor_prod.py @@ -2,7 +2,7 @@ import numpy as np from fuse import * from firedrake import * -from test_2d_examples_docs import construct_cg1, construct_dg1 +from test_2d_examples_docs import construct_cg1, construct_dg1, construct_dg1_integral # from test_convert_to_fiat import create_cg1 @@ -32,27 +32,34 @@ def mass_solve(U): assemble(L) solve(a == L, out) assert np.allclose(out.dat.data, f.dat.data, rtol=1e-5) + return out.dat.data -@pytest.mark.parametrize("generator, code, deg", [(construct_cg1, "CG", 1), (construct_dg1, "DG", 1)]) -def test_tensor_product_ext_mesh(generator, code, deg): +@pytest.mark.parametrize("generator1, generator2, code1, code2, deg1, deg2", + [(construct_cg1, construct_cg1, "CG", "CG", 1, 1), + (construct_dg1, construct_dg1, "DG", "DG", 1, 1), + (construct_dg1, construct_cg1, "DG", "CG", 1, 1), + (construct_dg1_integral, construct_cg1, "DG", "CG", 1, 1)]) +def test_ext_mesh(generator1, generator2, code1, code2, deg1, deg2): m = UnitIntervalMesh(2) mesh = ExtrudedMesh(m, 2) # manual method of creating tensor product elements - horiz_elt = FiniteElement(code, as_cell("interval"), deg) - vert_elt = FiniteElement(code, as_cell("interval"), deg) + horiz_elt = FiniteElement(code1, as_cell("interval"), deg1) + vert_elt = FiniteElement(code2, as_cell("interval"), deg2) elt = TensorProductElement(horiz_elt, vert_elt) U = FunctionSpace(mesh, elt) - mass_solve(U) + res1 = mass_solve(U) # fuseonic way of creating tensor product elements - A = generator() - B = generator() + A = generator1() + B = generator2() elem = tensor_product(A, B) U = FunctionSpace(mesh, elem.to_ufl()) - mass_solve(U) + res2 = mass_solve(U) + + assert np.allclose(res1, res2) def test_helmholtz(): @@ -117,3 +124,72 @@ def test_quad_mesh_helmholtz(): conv = np.log2(res[:-1] / res[1:]) print("convergence order:", conv) assert (np.array(conv) > 1.8).all() + + +@pytest.mark.parametrize(["A", "B", "res"], [(Point(0), line(), False), + (line(), line(), True), + (polygon(3), line(), False),]) +def test_flattening(A, B, res): + tensor_cell = TensorProductPoint(A, B) + if not res: + with pytest.raises(AssertionError): + tensor_cell.flatten() + else: + cell = tensor_cell.flatten() + cell.construct_fuse_rep() + + +# def test_trace_galerkin_projection(): +# mesh = UnitSquareMesh(10, 10, quadrilateral=True) + +# x, y = SpatialCoordinate(mesh) +# A = construct_cg1() +# B = construct_dg0_integral() +# elem = tensor_product(A, B) +# elem = elem.flatten() + +# # Define the Trace Space +# T = FunctionSpace(mesh, elem.to_ufl()) + +# # Define trial and test functions +# lambdar = TrialFunction(T) +# gammar = TestFunction(T) + +# # Define right hand side function + +# V = FunctionSpace(mesh, "CG", 1) +# f = Function(V) +# f.interpolate(cos(x*pi*2)*cos(y*pi*2)) + +# # Construct bilinear form +# a = inner(lambdar, gammar) * ds + inner(lambdar('+'), gammar('+')) * dS + +# # Construct linear form +# l = inner(f, gammar) * ds + inner(f('+'), gammar('+')) * dS + +# # Compute the solution +# t = Function(T) +# solve(a == l, t, solver_parameters={'ksp_rtol': 1e-14}) + +# # Compute error in trace norm +# trace_error = sqrt(assemble(FacetArea(mesh)*inner((t - f)('+'), (t - f)('+')) * dS)) + +# assert trace_error < 1e-13 + +# def test_hdiv(): +# np.set_printoptions(linewidth=90, precision=4, suppress=True) +# m = UnitIntervalMesh(2) +# mesh = ExtrudedMesh(m, 2) +# CG_1 = FiniteElement("CG", interval, 1) +# DG_0 = FiniteElement("DG", interval, 0) +# P1P0 = TensorProductElement(CG_1, DG_0) +# RT_horiz = HDivElement(P1P0) +# P0P1 = TensorProductElement(DG_0, CG_1) +# RT_vert = HDivElement(P0P1) +# elt = RT_horiz + RT_vert +# V = FunctionSpace(mesh, elt) +# tabulation = V.finat_element.fiat_equivalent.tabulate(0, [(0, 0), (1, 0)]) +# for ent, arr in tabulation.items(): +# print(ent) +# for comp in arr: +# print(comp[0], comp[1])