From 306ee6acaf71eb6aa9256507eece62c7850484ee Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 27 Feb 2020 10:39:26 -0600 Subject: [PATCH 01/38] Add trimmed BDM for squares --- FIAT/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index f5c1677ce..96760042c 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -11,6 +11,7 @@ from FIAT.bell import Bell from FIAT.argyris import QuinticArgyris from FIAT.brezzi_douglas_marini import BrezziDouglasMarini +from FIAT.trimmed_brezzi_douglas_marini import TrimmedBrezziDouglasMariniCubeEdge, TrimmedBrezziDouglasMariniCubeFace from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor From e38963a2da9e7dcc4a90e0c3e312e808d47eed10 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 27 Feb 2020 10:42:52 -0600 Subject: [PATCH 02/38] Added trimmed BDM on square --- FIAT/trimmed_brezzi_douglas_marini.py | 296 ++++++++++++++++++++++++++ 1 file changed, 296 insertions(+) create mode 100644 FIAT/trimmed_brezzi_douglas_marini.py diff --git a/FIAT/trimmed_brezzi_douglas_marini.py b/FIAT/trimmed_brezzi_douglas_marini.py new file mode 100644 index 000000000..046f1c587 --- /dev/null +++ b/FIAT/trimmed_brezzi_douglas_marini.py @@ -0,0 +1,296 @@ +# Copyright (C) 2019 Cyrus Cheng (Imperial College London) +# +# This file is part of FIAT. +# +# FIAT is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# FIAT is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with FIAT. If not, see . +# +# Modified by David A. Ham (david.ham@imperial.ac.uk), 2019 + +from sympy import symbols, legendre, Array, diff +import numpy as np +from FIAT.finite_element import FiniteElement +from FIAT.dual_set import make_entity_closure_ids +from FIAT.polynomial_set import mis +from FIAT.reference_element import compute_unflattening_map, flatten_reference_cube + +x, y, z = symbols('x y z') +variables = (x, y, z) +leg = legendre + + +def triangular_number(n): + return int((n+1)*n/2) + + +class TrimmedBrezziDouglasMariniCube(FiniteElement): + def __init__(self, ref_el, degree, mapping): + if degree < 1: + raise Exception("Trimmed BDMce_k elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + raise Exception("Trimmed BDMce_k elements only valid for dimension 2") + + flat_topology = flat_el.get_topology() + entity_ids = {} + cur = 0 + + for top_dim, entities in flat_topology.items(): + entity_ids[top_dim] = {} + for entity in entities: + entity_ids[top_dim][entity] = [] + for j in sorted(flat_topology[1]): + entity_ids[1][j] = list(range(cur, cur + degree)) #assign entity ids to everything in the first dimension. + cur = cur + degree + + if(degree >= 2): + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) #Assigns an entity IDs to the face. I need to figure out the pattern here for trimmed serendipity. + #else: + # entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) #Assigns an entity IDs to the face. I need to figure out the pattern here for trimmed serendipity. + + #print(entity_ids) + cur += 2*triangular_number(degree - 2) + degree + #print("cur is equal to") + #print(cur) + + formdegree = 1 + + entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) + + super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, + mapping=mapping) + + topology = ref_el.get_topology() + unflattening_map = compute_unflattening_map(topology) + unflattened_entity_ids = {} + unflattened_entity_closure_ids = {} + + for dim, entities in sorted(topology.items()): + unflattened_entity_ids[dim] = {} + unflattened_entity_closure_ids[dim] = {} + for dim, entities in sorted(flat_topology.items()): + for entity in entities: + unflat_dim, unflat_entity = unflattening_map[(dim, entity)] + unflattened_entity_ids[unflat_dim][unflat_entity] = entity_ids[dim][entity] + unflattened_entity_closure_ids[unflat_dim][unflat_entity] = entity_closure_ids[dim][entity] + self.entity_ids = unflattened_entity_ids + self.entity_closure_ids = unflattened_entity_closure_ids + self._degree = degree + self.flat_el = flat_el + + def degree(self): + return self._degree + + def get_nodal_basis(self): + raise NotImplementedError("get_nodal_basis not implemented for bdmce") + + def get_dual_set(self): + raise NotImplementedError("get_dual_set is not implemented for bdmce") + + def get_coeffs(self): + raise NotImplementedError("get_coeffs not implemented for bdmce") + + def tabulate(self, order, points, entity=None): + + if entity is None: + entity = (self.ref_el.get_dimension(), 0) + + entity_dim, entity_id = entity + transform = self.ref_el.get_entity_transform(entity_dim, entity_id) + points = list(map(transform, points)) + + phivals = {} + + for o in range(order+1): + alphas = mis(2, o) + for alpha in alphas: + try: + polynomials = self.basis[alpha] + except KeyError: + polynomials = diff(self.basis[(0, 0)], *zip(variables, alpha)) + self.basis[alpha] = polynomials + T = np.zeros((len(polynomials[:, 0]), 2, len(points))) + for i in range(len(points)): + subs = {v: points[i][k] for k, v in enumerate(variables[:2])} + for j, f in enumerate(polynomials[:, 0]): + T[j, 0, i] = f.evalf(subs=subs) + for j, f in enumerate(polynomials[:, 1]): + T[j, 1, i] = f.evalf(subs=subs) + phivals[alpha] = T + + return phivals + + def entity_dofs(self): + """Return the map of topological entities to degrees of + freedom for the finite element.""" + return self.entity_ids + + def entity_closure_dofs(self): + """Return the map of topological entities to degrees of + freedom on the closure of those entities for the finite element.""" + return self.entity_closure_ids + + def value_shape(self): + return (2,) + + def dmats(self): + raise NotImplementedError + + def get_num_members(self, arg): + raise NotImplementedError + + def space_dimension(self): + return int(len(self.basis[(0, 0)])/2) + + +#Splitting the E Lambda function into two seperate functions for E Lambda and E tilde Lambda. +#Correlating with Andrew's paper, leg(j, x_mid) should be a polynomial x^i, leg(j, y_mid) should be y^i, +#dy[0] should represent y-1, dy[1] should represent y+1 (and similar for the dx and x+/- 1). +#Still not sure why we use for loops in only the EL tuple but not the ELTilde tuple. + +def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): + EL = tuple([(-leg(j, x_mid)*dy[0], 0) for j in range(deg)]+ + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]+ + [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)]+ + [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)]) + return EL + +def e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid): + ELTilde = tuple([(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]+ + [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]+ + [(-leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[0])]+ + [(leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[1])]) + return ELTilde + +def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + EL = e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid) + ELTilde = e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid) + + result = EL + ELTilde + return result + +def determine_f_lambda_portions_2d(deg): + if (deg < 2): + DegsOfIteration = [] + else: + DegsOfIteration = [] + for i in range(2, deg): + DegsOfIteration += [i] + + return DegsOfIteration + +def f_lambda_1_2d_pieces(current_deg, dx, dy, x_mid, y_mid): + if (current_deg == 2): + FLpiece = [(leg(0, x_mid) * leg(0, y_mid) * dy[0] * dy[1], 0)] + FLpiece += [(0, leg(0, x_mid) * leg(0, y_mid) * dx[0] * dx[1])] + else: + target_power = current_deg - 2 + FLpiece = tuple([]) + for j in range(0, target_power + 1): + k = target_power - j + FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dy[0] * dy[1], 0)]) + FLpiece += tuple([(0, leg(j, x_mid) * leg(k, y_mid) * dx[0] * dx[1])]) + return FLpiece + +def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): + DegsOfIteration = determine_f_lambda_portions_2d(deg) + FL = [] + for i in DegsOfIteration: + FL += f_lambda_1_2d_pieces(i, dx, dy, x_mid, y_mid) + return tuple(FL) + + + +def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + FL = [] + for k in range(2, deg+1): + for j in range(k-1): + FL += [(0, leg(j, x_mid)*leg(k-2-j, y_mid)*dx[0]*dx[1])] + FL += [(leg(k-2-j, x_mid)*leg(j, y_mid)*dy[0]*dy[1], 0)] + return tuple(FL) + +def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): + FLTilde = tuple([]) + FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) + FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) + for k in range (1, deg - 1): + FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) + + return tuple(FLTilde) + +def trimmed_f_lambda(deg, dx, dy, x_mid, y_mid): + FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) + FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) + result = FL + FLT + + return result + + + +class TrimmedBrezziDouglasMariniCubeEdge(TrimmedBrezziDouglasMariniCube): + def __init__(self, ref_el, degree): + if degree < 1: + raise Exception("Trimmed BDMce_k elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + raise Exception("Trimmed BDMce_k elements only valid for dimension 2") + + verts = flat_el.get_vertices() + + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) + dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + if degree >= 2: + FL = trimmed_f_lambda(degree, dx, dy, x_mid, y_mid) + else: + FL = () + + bdmce_list = EL + FL + self.basis = {(0, 0): Array(bdmce_list)} + super(TrimmedBrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") + + +class TrimmedBrezziDouglasMariniCubeFace(TrimmedBrezziDouglasMariniCube): + def __init__(self, ref_el, degree): + if degree < 1: + raise Exception("Trimmed BDMcf_k elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + raise Exception("Trimmed BDMcf_k elements only valid for dimension 2") + + verts = flat_el.get_vertices() + + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) + dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + if degree >= 2: + FL = trimmed_f_lambda(degree, dx, dy, x_mid, y_mid) + else: + FL = () + bdmcf_list = EL + FL + bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] + self.basis = {(0, 0): Array(bdmcf_list)} + super(TrimmedBrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + From e47be550b985e84928c50e0c1608105bce921f12 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 27 Feb 2020 11:05:08 -0600 Subject: [PATCH 03/38] Rename --- ...med_brezzi_douglas_marini.py => Sminus.py} | 30 +++++++++---------- FIAT/__init__.py | 2 +- 2 files changed, 16 insertions(+), 16 deletions(-) rename FIAT/{trimmed_brezzi_douglas_marini.py => Sminus.py} (90%) diff --git a/FIAT/trimmed_brezzi_douglas_marini.py b/FIAT/Sminus.py similarity index 90% rename from FIAT/trimmed_brezzi_douglas_marini.py rename to FIAT/Sminus.py index 046f1c587..c661a1601 100644 --- a/FIAT/trimmed_brezzi_douglas_marini.py +++ b/FIAT/Sminus.py @@ -33,15 +33,15 @@ def triangular_number(n): return int((n+1)*n/2) -class TrimmedBrezziDouglasMariniCube(FiniteElement): +class TrimmedSerendipityCube(FiniteElement): def __init__(self, ref_el, degree, mapping): if degree < 1: - raise Exception("Trimmed BDMce_k elements only valid for k >= 1") + raise Exception("Trimmed serendipity elements only valid for k >= 1") flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - raise Exception("Trimmed BDMce_k elements only valid for dimension 2") + raise Exception("Trimmed serendipity elements only valid for dimension 2") flat_topology = flat_el.get_topology() entity_ids = {} @@ -69,7 +69,7 @@ def __init__(self, ref_el, degree, mapping): entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, + super(SerendipityCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, mapping=mapping) topology = ref_el.get_topology() @@ -94,13 +94,13 @@ def degree(self): return self._degree def get_nodal_basis(self): - raise NotImplementedError("get_nodal_basis not implemented for bdmce") + raise NotImplementedError("get_nodal_basis not implemented for trimmed serendipity") def get_dual_set(self): - raise NotImplementedError("get_dual_set is not implemented for bdmce") + raise NotImplementedError("get_dual_set is not implemented for trimmed serendipity") def get_coeffs(self): - raise NotImplementedError("get_coeffs not implemented for bdmce") + raise NotImplementedError("get_coeffs not implemented for trimmed serendipity") def tabulate(self, order, points, entity=None): @@ -239,15 +239,15 @@ def trimmed_f_lambda(deg, dx, dy, x_mid, y_mid): -class TrimmedBrezziDouglasMariniCubeEdge(TrimmedBrezziDouglasMariniCube): +class TrimmedSerendipityCubeEdge(TrimmedSerendipityCube): def __init__(self, ref_el, degree): if degree < 1: - raise Exception("Trimmed BDMce_k elements only valid for k >= 1") + raise Exception("Trimmed Serendipity_k edge elements only valid for k >= 1") flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - raise Exception("Trimmed BDMce_k elements only valid for dimension 2") + raise Exception("Trimmed Serendipity_k edge elements only valid for dimension 2") verts = flat_el.get_vertices() @@ -264,18 +264,18 @@ def __init__(self, ref_el, degree): bdmce_list = EL + FL self.basis = {(0, 0): Array(bdmce_list)} - super(TrimmedBrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") + super(TrimmedSerendipityCubeEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") -class TrimmedBrezziDouglasMariniCubeFace(TrimmedBrezziDouglasMariniCube): +class TrimmedSerendipityCubeFace(TrimmedSerendipityCube): def __init__(self, ref_el, degree): if degree < 1: - raise Exception("Trimmed BDMcf_k elements only valid for k >= 1") + raise Exception("Trimmed serendipity face elements only valid for k >= 1") flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - raise Exception("Trimmed BDMcf_k elements only valid for dimension 2") + raise Exception("Trimmed serendipity face elements only valid for dimension 2") verts = flat_el.get_vertices() @@ -292,5 +292,5 @@ def __init__(self, ref_el, degree): bdmcf_list = EL + FL bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} - super(TrimmedBrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super(TrimmedSerendipityCubeFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 96760042c..d0b3abdc0 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -11,7 +11,7 @@ from FIAT.bell import Bell from FIAT.argyris import QuinticArgyris from FIAT.brezzi_douglas_marini import BrezziDouglasMarini -from FIAT.trimmed_brezzi_douglas_marini import TrimmedBrezziDouglasMariniCubeEdge, TrimmedBrezziDouglasMariniCubeFace +from FIAT.Sminus import TrimmedSerendipityCubeEdge, TrimmedSerendipityCubeFace from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor From ef7bf9db3e80d2da21866c1111892525e99b2a17 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 27 Feb 2020 11:08:15 -0600 Subject: [PATCH 04/38] Rename again --- FIAT/Sminus.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index c661a1601..d8fbfe762 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -33,7 +33,7 @@ def triangular_number(n): return int((n+1)*n/2) -class TrimmedSerendipityCube(FiniteElement): +class TrimmedSerendipity(FiniteElement): def __init__(self, ref_el, degree, mapping): if degree < 1: raise Exception("Trimmed serendipity elements only valid for k >= 1") @@ -69,7 +69,7 @@ def __init__(self, ref_el, degree, mapping): entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - super(SerendipityCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, + super(TrimmedSerendipity, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, mapping=mapping) topology = ref_el.get_topology() @@ -239,7 +239,7 @@ def trimmed_f_lambda(deg, dx, dy, x_mid, y_mid): -class TrimmedSerendipityCubeEdge(TrimmedSerendipityCube): +class TrimmedSerendipityEdge(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: raise Exception("Trimmed Serendipity_k edge elements only valid for k >= 1") @@ -264,10 +264,10 @@ def __init__(self, ref_el, degree): bdmce_list = EL + FL self.basis = {(0, 0): Array(bdmce_list)} - super(TrimmedSerendipityCubeEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") + super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") -class TrimmedSerendipityCubeFace(TrimmedSerendipityCube): +class TrimmedSerendipityFace(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: raise Exception("Trimmed serendipity face elements only valid for k >= 1") @@ -292,5 +292,5 @@ def __init__(self, ref_el, degree): bdmcf_list = EL + FL bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} - super(TrimmedSerendipityCubeFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") From 96423643efcbae7da02745e5c77f3a73d5390018 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 27 Feb 2020 12:43:46 -0600 Subject: [PATCH 05/38] Fix change of name in __init__.py --- FIAT/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index d0b3abdc0..d4f6bc28e 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -11,7 +11,7 @@ from FIAT.bell import Bell from FIAT.argyris import QuinticArgyris from FIAT.brezzi_douglas_marini import BrezziDouglasMarini -from FIAT.Sminus import TrimmedSerendipityCubeEdge, TrimmedSerendipityCubeFace +from FIAT.Sminus import TrimmedSerendipityEdge, TrimmedSerendipityFace from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor From f8841d16e9acbd3932d18ef2f97fb93e1d1f8c13 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 27 Feb 2020 14:19:46 -0600 Subject: [PATCH 06/38] Fix basis function ordering. --- FIAT/Sminus.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index d8fbfe762..3ae671334 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -56,14 +56,9 @@ def __init__(self, ref_el, degree, mapping): cur = cur + degree if(degree >= 2): - entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) #Assigns an entity IDs to the face. I need to figure out the pattern here for trimmed serendipity. - #else: - # entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) #Assigns an entity IDs to the face. I need to figure out the pattern here for trimmed serendipity. - - #print(entity_ids) + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + cur += 2*triangular_number(degree - 2) + degree - #print("cur is equal to") - #print(cur) formdegree = 1 @@ -161,10 +156,12 @@ def space_dimension(self): #Still not sure why we use for loops in only the EL tuple but not the ELTilde tuple. def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): - EL = tuple([(-leg(j, x_mid)*dy[0], 0) for j in range(deg)]+ - [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]+ - [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)]+ - [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)]) + EL = tuple( + [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)]+ + [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)]+ + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)]+ + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]) + return EL def e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid): @@ -261,7 +258,7 @@ def __init__(self, ref_el, degree): FL = trimmed_f_lambda(degree, dx, dy, x_mid, y_mid) else: FL = () - + bdmce_list = EL + FL self.basis = {(0, 0): Array(bdmce_list)} super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") From ab80b4e951a7f56e9b15fb68bec73ec8bd4143db Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 27 Feb 2020 14:25:42 -0600 Subject: [PATCH 07/38] Feeling less flake8-y --- FIAT/Sminus.py | 59 ++++++++++++++++++++++++++++-------------------- FIAT/__init__.py | 3 ++- 2 files changed, 37 insertions(+), 25 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 3ae671334..42b03a32b 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -52,20 +52,23 @@ def __init__(self, ref_el, degree, mapping): for entity in entities: entity_ids[top_dim][entity] = [] for j in sorted(flat_topology[1]): - entity_ids[1][j] = list(range(cur, cur + degree)) #assign entity ids to everything in the first dimension. + entity_ids[1][j] = list(range(cur, cur + degree)) # assign entity ids to everything in the first dimension. cur = cur + degree if(degree >= 2): entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) - + cur += 2*triangular_number(degree - 2) + degree formdegree = 1 entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - super(TrimmedSerendipity, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, - mapping=mapping) + super(TrimmedSerendipity, self).__init__(ref_el=ref_el, + dual=None, + order=degree, + formdegree=formdegree, + mapping=mapping) topology = ref_el.get_topology() unflattening_map = compute_unflattening_map(topology) @@ -150,34 +153,41 @@ def space_dimension(self): return int(len(self.basis[(0, 0)])/2) -#Splitting the E Lambda function into two seperate functions for E Lambda and E tilde Lambda. -#Correlating with Andrew's paper, leg(j, x_mid) should be a polynomial x^i, leg(j, y_mid) should be y^i, -#dy[0] should represent y-1, dy[1] should represent y+1 (and similar for the dx and x+/- 1). -#Still not sure why we use for loops in only the EL tuple but not the ELTilde tuple. +# Splitting the E Lambda function into two seperate functions for E Lambda and E tilde Lambda. +# Correlating with Andrew's paper, leg(j, x_mid) should be a polynomial x^i, leg(j, y_mid) should be y^i, +# dy[0] should represent y-1, dy[1] should represent y+1 (and similar for the dx and x+/- 1). +# Still not sure why we use for loops in only the EL tuple but not the ELTilde tuple. def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): EL = tuple( - [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)]+ - [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)]+ - [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)]+ - [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]) + [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)] + + [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)] + + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]) return EL + def e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid): - ELTilde = tuple([(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]+ - [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]+ - [(-leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[0])]+ - [(leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[1])]) + ELTilde = tuple([(-leg(deg, x_mid) * dy[0], + -leg(deg-1, x_mid) * dx[0] * dx[1] / (deg+1))] + + [(-leg(deg, x_mid) * dy[1], + leg(deg-1, x_mid) * dx[0] * dx[1] / (deg+1))] + + [(-leg(deg-1, y_mid) * dy[0] * dy[1] / (deg+1), + -leg(deg, y_mid) * dx[0])] + + [(leg(deg-1, y_mid) * dy[0] * dy[1] / (deg+1), + -leg(deg, y_mid) * dx[1])]) return ELTilde + def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): EL = e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid) ELTilde = e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid) - + result = EL + ELTilde return result + def determine_f_lambda_portions_2d(deg): if (deg < 2): DegsOfIteration = [] @@ -185,13 +195,14 @@ def determine_f_lambda_portions_2d(deg): DegsOfIteration = [] for i in range(2, deg): DegsOfIteration += [i] - + return DegsOfIteration + def f_lambda_1_2d_pieces(current_deg, dx, dy, x_mid, y_mid): if (current_deg == 2): - FLpiece = [(leg(0, x_mid) * leg(0, y_mid) * dy[0] * dy[1], 0)] - FLpiece += [(0, leg(0, x_mid) * leg(0, y_mid) * dx[0] * dx[1])] + FLpiece = [(leg(0, x_mid) * leg(0, y_mid) * dy[0] * dy[1], 0)] + FLpiece += [(0, leg(0, x_mid) * leg(0, y_mid) * dx[0] * dx[1])] else: target_power = current_deg - 2 FLpiece = tuple([]) @@ -201,6 +212,7 @@ def f_lambda_1_2d_pieces(current_deg, dx, dy, x_mid, y_mid): FLpiece += tuple([(0, leg(j, x_mid) * leg(k, y_mid) * dx[0] * dx[1])]) return FLpiece + def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): DegsOfIteration = determine_f_lambda_portions_2d(deg) FL = [] @@ -209,7 +221,6 @@ def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): return tuple(FL) - def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): FL = [] for k in range(2, deg+1): @@ -218,15 +229,17 @@ def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): FL += [(leg(k-2-j, x_mid)*leg(j, y_mid)*dy[0]*dy[1], 0)] return tuple(FL) + def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): FLTilde = tuple([]) FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) - for k in range (1, deg - 1): + for k in range(1, deg - 1): FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) return tuple(FLTilde) + def trimmed_f_lambda(deg, dx, dy, x_mid, y_mid): FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) @@ -235,7 +248,6 @@ def trimmed_f_lambda(deg, dx, dy, x_mid, y_mid): return result - class TrimmedSerendipityEdge(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: @@ -290,4 +302,3 @@ def __init__(self, ref_el, degree): bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") - diff --git a/FIAT/__init__.py b/FIAT/__init__.py index d4f6bc28e..0cd27f747 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -11,7 +11,8 @@ from FIAT.bell import Bell from FIAT.argyris import QuinticArgyris from FIAT.brezzi_douglas_marini import BrezziDouglasMarini -from FIAT.Sminus import TrimmedSerendipityEdge, TrimmedSerendipityFace +from FIAT.Sminus import TrimmedSerendipityEdge # noqa: F401 +from FIAT.Sminus import TrimmedSerendipityFace # noqa: F401 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor From deece33b7a02d784c2208b3efef2128cb57b5495 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Fri, 13 Mar 2020 11:07:05 -0700 Subject: [PATCH 08/38] Working on implementing the 3d version of 1-forms. --- FIAT/Sminus.py | 171 +++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 159 insertions(+), 12 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 42b03a32b..be4a86fb6 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -221,15 +221,6 @@ def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): return tuple(FL) -def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - FL = [] - for k in range(2, deg+1): - for j in range(k-1): - FL += [(0, leg(j, x_mid)*leg(k-2-j, y_mid)*dx[0]*dx[1])] - FL += [(leg(k-2-j, x_mid)*leg(j, y_mid)*dy[0]*dy[1], 0)] - return tuple(FL) - - def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): FLTilde = tuple([]) FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) @@ -240,7 +231,7 @@ def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): return tuple(FLTilde) -def trimmed_f_lambda(deg, dx, dy, x_mid, y_mid): +def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) result = FL + FLT @@ -248,6 +239,162 @@ def trimmed_f_lambda(deg, dx, dy, x_mid, y_mid): return result +def e_lambda_3d_piece(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): + ELpiece = tuple([(leg(current_deg, x_mid) * dy[0] * dz[0], 0, 0)] + + [(leg(current_deg, x_mid)) * dy[1] * dz[0], 0, 0)] + + [(leg(current_deg, x_mid) * dy[1] * dz[1], 0, 0)] + + [(leg(current_deg, x_mid) * dy[0] * dz[1], 0, 0)] + + [(0, leg(current_deg, y_mid) * dx[0] * dz[1], 0)] + + [(0, leg(current_deg, y_mid) * dx[1] * dz[1], 0)] + + [(0, leg(current_deg, y_mid) * dx[1] * dz[0], 0)] + + [(0, leg(current_deg, y_mid) * dx[0] * dz[0], 0)] + + [(0, 0, leg(current_deg, z_mid) * dx[0] * dy[1])] + + [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[1])] + + [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[0])] + + [(0, 0, leg(current_deg, z_mid) * dx[0] * dy[0])]) + return ELpiece + +def e_lambda_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): + EL = [] + for i in range(max_deg): + Elpiece = e_lambda_3d_piece(i, dx, dy, dz, x_mid, y_mid, z_mid) + EL += Elpiece + return EL + +def determine_f_lambda_portions_3d(deg): + if (deg < 2): + DegsOfIteration = [] + else: + DegsOfIteration = [] + for i in range(2, deg): + DegsOfIteration += [i] + + return DegsOfIteration + +def f_lambda_1_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): + if (current_deg == 2): + j = 0 + k = 0 + FLpiece = [(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)] + FLpiece += [(leg(j, x_mid) * leg(k, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)] + FLpiece += [(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)] + FLpiece += [(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)] + FLpiece += [(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)] + FLpiece += [(0, leg(j, y_mid) * leg(k, x_mid) * dz[1] * dx[0] * dx[1], 0)] + FLpiece += [(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)] + FLpiece += [(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)] + FLpiece += [(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[0] * dx[0] * dx[1])] + FLpiece += [(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[1] * dx[0] * dx[1])] + FLpiece += [(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[0] * dy[0] * dy[1])] + FLpiece += [(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[1] * dy[0] * dy[1])] + else: + target_power = current_deg - 2 + FLpiece = tuple([]) + for j in range(0, target_power + 1): + k = target_power - j + FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) + FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) + FLpiece += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) + FLpiece += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) + FLpiece += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)]) + FLpiece += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[1] * dx[0] * dx[1], 0)]) + FLpiece += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) + FLpiece += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)]) + FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[0] * dx[0] * dx[1])]) + FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[1] * dx[0] * dx[1])]) + FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[0] * dy[0] * dy[1])]) + FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[1] * dy[0] * dy[1])]) + return FLpiece + + +def f_lambda_1_3d_tilde(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): + FLtilde = tuple([]) + FLtilde += tuple([(leg(max_deg - 2, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) + FLtilde += tuple([(leg(max_deg - 2, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) + FLtilde += tuple([(leg(max_deg - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) + FLtilde += tuple([(leg(max_deg - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) + FLtilde += tuple([(0, leg(max_deg - 2, x_mid) * dz[0] * dx[0] * dx[1], 0)]) + FLtilde += tuple([(0, leg(max_deg - 2, x_mid) * dz[1] * dx[0] * dx[1], 0)]) + FLtilde += tuple([(0, leg(max_deg - 2, z_mid) * dx[0] * dz[0] * dz[1], 0)]) + FLtilde += tuple([(0, leg(max_deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) + FLtilde += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[0] * dx[0] * dx[1])]) + FLtilde += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[1] * dx[0] * dx[1])]) + FLtilde += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[0] * dy[0] * dy[1])]) + FLtilde += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) + for j in range(1, max_deg - 1): + for k in range(0, 2): + FLtilde += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, y_mid) * dz[k] * dy[0] * dy[1], + -leg(j-1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[k] * dx[0] * dx[1], 0)]) + FLtilde += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[k] * dz[0] * dz[1], 0, + -leg(j-1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[k] * dx[0] * dx[1])]) + FLtilde += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[k] * dz[0] * dz[1], + -leg(j-1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[k] * dy[0] * dy[1])]) + return FLtilde + + +def f_lambda_1_3d_trim(deg, dx, dy, dz, x_mid, y_mid, z_mid): + DegsOfIteration = determine_f_lambda_portions_3d(deg) + FL = [] + for i in DegsOfIteration: + FL += f_lambda_1_3d_pieces(i, dx, dy, dz, x_mid, y_mid, z_mid) + return tuple(FL) + + +def trimmed_f_lambda_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): + FL = f_lambda_1_3d_trim(deg, dx, dy, dz, x_mid, y_mid, z_mid) + FLT = f_lambda_1_3d_tilde(deg, dx, dy, dz, x_mid, y_mid, z_mid) + result = FL + FLT + + return result + + +def determine_I_lambda_1_portions_3d(deg): + if (deg < 4): + DegsOfIteration = [] + else: + Degs + DegsOfIteration = tuple([]) + for x in range(0, deg - 3): + for y in range(0, deg - 3 - x): + for z in range(0, deg - 3 - x - y): + Degs += tuple([(x, y, z)]) + for degs in Degs: + if(degs[0] + degs[1] + degs[2] == deg - 4): + DegsOfIteration += tuple([degs]) + return DegsOfIteration + + +def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): + DegsOfIteration = determine_I_lambda_1_portions_3d(deg) + IL = tuple([]) + for Degs in DegsOfIteration: + IL += tuple([(leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * + dy[0] * dy[1] * dz[0] * dz[1], 0, 0)]) + IL += tuple([(0, leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * + dx[0] * dx[1] * dz[0] * dz[1], 0)]) + IL += tuple([(0, 0, leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * + dy[0] * dy[1] * dy[0] * dy[1])]) + + +def I_lambda_1_tilde_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): + ILtilde = tuple([]) + ILtilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)]) + ILtilde += tuple([(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)]) + ILtilde += tuple([(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)]) + ILtilde += tuple([(0, leg(deg - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)]) + ILtilde += tuple([(0, 0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + ILtilde += tuple([(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + for j in range(1, deg - 3): + ILtilde += tuple([(leg(j, x_mid) * leg(deg - j - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], + -leg(j - 1, x_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)]) + ILtilde += tuple([(leg(j, x_mid) * leg(deg - j - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, + -leg(j - 1, x_mid) * leg(deg - j - 3, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + if(deg > 5): + ILtilde += tuple([(0, leg(j, y_mid) * leg(deg - j - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], + -leg(j - 1, y_mid) * leg(deg - j - 3, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + return ILtilde + + class TrimmedSerendipityEdge(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: @@ -267,7 +414,7 @@ def __init__(self, ref_el, degree): EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: - FL = trimmed_f_lambda(degree, dx, dy, x_mid, y_mid) + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) else: FL = () @@ -295,7 +442,7 @@ def __init__(self, ref_el, degree): EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: - FL = trimmed_f_lambda(degree, dx, dy, x_mid, y_mid) + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) else: FL = () bdmcf_list = EL + FL From c782b9c84149eb787646ef728bc80af0f5f83172 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Fri, 20 Mar 2020 09:35:06 -0700 Subject: [PATCH 09/38] Continuing work on updating Sminus to work in 3d --- FIAT/Sminus.py | 108 ++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 84 insertions(+), 24 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index be4a86fb6..72bf2a77b 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -41,7 +41,8 @@ def __init__(self, ref_el, degree, mapping): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - raise Exception("Trimmed serendipity elements only valid for dimension 2") + if dim != 3: + raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3") flat_topology = flat_el.get_topology() entity_ids = {} @@ -51,14 +52,32 @@ def __init__(self, ref_el, degree, mapping): entity_ids[top_dim] = {} for entity in entities: entity_ids[top_dim][entity] = [] - for j in sorted(flat_topology[1]): - entity_ids[1][j] = list(range(cur, cur + degree)) # assign entity ids to everything in the first dimension. - cur = cur + degree - - if(degree >= 2): - entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) - - cur += 2*triangular_number(degree - 2) + degree + if dim == 2: + for j in sorted(flat_topology[1]): + entity_ids[1][j] = list(range(cur, cur + degree)) + cur = cur + degree + + if(degree >= 2): + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + + cur += 2*triangular_number(degree - 2) + degree + else: + #3-d case. + for j in sorted(flat_topology[1]): + entity_ids[1][j] = list(range(cur, cur + degree)) + cur = cur + degree + + if (degree >= 2): + for j in sorted(flat_topology[2]): + entity_ids[2][j] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + + if (degree >= 4): + if (degree == 4): + entity_ids[3][0] = list(range(cur, cur + 6)) + elif (degree == 5): + entity_ids[3][0] = list(range(cur, cur + 8)) + else: + entity_ids[3][0] = list(range(cur, cur + 6 + (degree - 4) * 3)) formdegree = 1 @@ -239,9 +258,9 @@ def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): return result -def e_lambda_3d_piece(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): +def e_lambda_1_3d_piece(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): ELpiece = tuple([(leg(current_deg, x_mid) * dy[0] * dz[0], 0, 0)] + - [(leg(current_deg, x_mid)) * dy[1] * dz[0], 0, 0)] + + [(leg(current_deg, x_mid) * dy[1] * dz[0], 0, 0)] + [(leg(current_deg, x_mid) * dy[1] * dz[1], 0, 0)] + [(leg(current_deg, x_mid) * dy[0] * dz[1], 0, 0)] + [(0, leg(current_deg, y_mid) * dx[0] * dz[1], 0)] + @@ -254,14 +273,14 @@ def e_lambda_3d_piece(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): [(0, 0, leg(current_deg, z_mid) * dx[0] * dy[0])]) return ELpiece -def e_lambda_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): +def e_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): EL = [] for i in range(max_deg): - Elpiece = e_lambda_3d_piece(i, dx, dy, dz, x_mid, y_mid, z_mid) + Elpiece = e_lambda_1_3d_piece(i, dx, dy, dz, x_mid, y_mid, z_mid) EL += Elpiece return EL -def determine_f_lambda_portions_3d(deg): +def determine_f_lambda_1_portions_3d(deg): if (deg < 2): DegsOfIteration = [] else: @@ -403,7 +422,8 @@ def __init__(self, ref_el, degree): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - raise Exception("Trimmed Serendipity_k edge elements only valid for dimension 2") + if dim != 3: + raise Exception("Trimmed Serendipity_k edge elements only valid for dimensions 2 and 3") verts = flat_el.get_vertices() @@ -411,14 +431,33 @@ def __init__(self, ref_el, degree): dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) x_mid = 2*x-(verts[-1][0] + verts[0][0]) y_mid = 2*y-(verts[-1][1] + verts[0][1]) - - EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + try: + dz = ((verts[-1][2] - z)/(verts[-1][2] - verts[0][2]), (z - verts[0][2])/(verts[-1][2] - verts[0][2])) + z_mid = 2*z-(verts[-1][2] + verts[0][2]) + except IndexError: + dz = None + z_mid = None + + if dim == 2: + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + else: + EL = e_lambda_1_3d_trimmed(degree, dx, dy, dz, x_mid, y_mid, z_mid) if degree >= 2: - FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + if dim == 2: + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + else: + FL = trimmed_f_lambda_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) else: FL = () + if dim == 3: + if degree >= 4: + IL = I_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + I_lambda_1_tilde_3d(degree, dx, dy, + dz, x_mid, + y_mid, z_mid) + else: + IL = () - bdmce_list = EL + FL + bdmce_list = EL + FL + IL self.basis = {(0, 0): Array(bdmce_list)} super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") @@ -431,7 +470,8 @@ def __init__(self, ref_el, degree): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - raise Exception("Trimmed serendipity face elements only valid for dimension 2") + if dim != 3: + raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") verts = flat_el.get_vertices() @@ -439,13 +479,33 @@ def __init__(self, ref_el, degree): dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) x_mid = 2*x-(verts[-1][0] + verts[0][0]) y_mid = 2*y-(verts[-1][1] + verts[0][1]) - - EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + try: + dz = ((verts[-1][2] - z)/(verts[-1][2] - verts[0][2]), (z - verts[0][2])/(verts[-1][2] - verts[0][2])) + z_mid = 2*z-(verts[-1][2] + verts[0][2]) + except IndexError: + dz = None + z_mid = None + + if (dim == 2): + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + else: + EL = e_lambda_1_3d_trimmed(degree, dx, dy, dz, x_mid, y_mid, z_mid) if degree >= 2: - FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + if dim == 2: + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + else: + FL = trimmed_f_lambda_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) else: FL = () - bdmcf_list = EL + FL + if dim == 3: + if degree >= 4: + IL = I_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + I_lambda_1_tilde_3d(degree, dx, dy, + dz, x_mid, + y_mid, z_mid) + else: + IL = () + + bdmcf_list = EL + FL + IL bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") From 6bb477c02969f546d0b97ef38f2c6c8e9482fb1e Mon Sep 17 00:00:00 2001 From: Justincrum Date: Wed, 1 Apr 2020 14:20:09 -0700 Subject: [PATCH 10/38] Continuing work on Sminus.py --- FIAT/Sminus.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 72bf2a77b..70c600b12 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -274,7 +274,7 @@ def e_lambda_1_3d_piece(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): return ELpiece def e_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): - EL = [] + EL = tuple([]) for i in range(max_deg): Elpiece = e_lambda_1_3d_piece(i, dx, dy, dz, x_mid, y_mid, z_mid) EL += Elpiece @@ -352,7 +352,7 @@ def f_lambda_1_3d_tilde(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): def f_lambda_1_3d_trim(deg, dx, dy, dz, x_mid, y_mid, z_mid): - DegsOfIteration = determine_f_lambda_portions_3d(deg) + DegsOfIteration = determine_f_lambda_1_portions_3d(deg) FL = [] for i in DegsOfIteration: FL += f_lambda_1_3d_pieces(i, dx, dy, dz, x_mid, y_mid, z_mid) @@ -371,7 +371,7 @@ def determine_I_lambda_1_portions_3d(deg): if (deg < 4): DegsOfIteration = [] else: - Degs + Degs = tuple([]) DegsOfIteration = tuple([]) for x in range(0, deg - 3): for y in range(0, deg - 3 - x): @@ -393,6 +393,7 @@ def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): dx[0] * dx[1] * dz[0] * dz[1], 0)]) IL += tuple([(0, 0, leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * dy[0] * dy[1] * dy[0] * dy[1])]) + return IL def I_lambda_1_tilde_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): @@ -456,7 +457,6 @@ def __init__(self, ref_el, degree): y_mid, z_mid) else: IL = () - bdmce_list = EL + FL + IL self.basis = {(0, 0): Array(bdmce_list)} super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") @@ -505,6 +505,9 @@ def __init__(self, ref_el, degree): else: IL = () + print(len(EL)) + print(len(FL)) + print(len(IL)) bdmcf_list = EL + FL + IL bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} From 2ad80e2548c09bbb6ccab515f4af4a473675b1f5 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Fri, 10 Apr 2020 08:31:58 -0700 Subject: [PATCH 11/38] Sminus.py continues to bug me. --- FIAT/Sminus.py | 80 ++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 58 insertions(+), 22 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 70c600b12..588afdacd 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -32,6 +32,15 @@ def triangular_number(n): return int((n+1)*n/2) +def choose_ijk_total(degree): + top = 1 + for i in range(1, 2 + degree + 1): + top = i * top + bottom = 1 + for i in range(1, degree + 1): + bottom = i * bottom + return int(top /(2 * bottom)) + class TrimmedSerendipity(FiniteElement): def __init__(self, ref_el, degree, mapping): @@ -61,23 +70,37 @@ def __init__(self, ref_el, degree, mapping): entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) cur += 2*triangular_number(degree - 2) + degree + else: #3-d case. + entity_ids[3] = {} for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) cur = cur + degree - if (degree >= 2): - for j in sorted(flat_topology[2]): - entity_ids[2][j] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + if (degree >= 2 ): + if (degree < 4): + for j in sorted(flat_topology[2]): + entity_ids[2][j] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + cur = cur + 2*triangular_number(degree - 2) + degree + else: + for j in sorted(flat_topology[2]): + entity_ids[2][j] = list(range(cur, cur + 3 * degree - 4)) + cur = cur + 3*degree - 4 + if (degree >= 4): if (degree == 4): entity_ids[3][0] = list(range(cur, cur + 6)) elif (degree == 5): - entity_ids[3][0] = list(range(cur, cur + 8)) + entity_ids[3][0] = list(range(cur, cur + 11)) else: - entity_ids[3][0] = list(range(cur, cur + 6 + (degree - 4) * 3)) + interior_ids = 0 + for i in range(0, degree - 4): + interior_ids += 3 * choose_ijk_total(i) + entity_ids[3][0] = list(range(cur, cur + 6 + (degree - 4) * 3 + interior_ids)) + else: + entity_ids[3][0] = list(range(cur, cur)) formdegree = 1 @@ -258,19 +281,29 @@ def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): return result +""" def e_lambda_1_3d_part_one(deg, dx, dy, dz, x_mid, y_mid, z_mid): + EL = tuple( + [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)] + + [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)] + + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]) + + return EL """ + + def e_lambda_1_3d_piece(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): - ELpiece = tuple([(leg(current_deg, x_mid) * dy[0] * dz[0], 0, 0)] + - [(leg(current_deg, x_mid) * dy[1] * dz[0], 0, 0)] + - [(leg(current_deg, x_mid) * dy[1] * dz[1], 0, 0)] + - [(leg(current_deg, x_mid) * dy[0] * dz[1], 0, 0)] + + ELpiece = tuple([(0, 0, leg(current_deg, z_mid) * dx[0] * dy[1])] + + [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[1])] + + [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[0])] + + [(0, 0, leg(current_deg, z_mid) * dx[0] * dy[0])] + [(0, leg(current_deg, y_mid) * dx[0] * dz[1], 0)] + [(0, leg(current_deg, y_mid) * dx[1] * dz[1], 0)] + [(0, leg(current_deg, y_mid) * dx[1] * dz[0], 0)] + [(0, leg(current_deg, y_mid) * dx[0] * dz[0], 0)] + - [(0, 0, leg(current_deg, z_mid) * dx[0] * dy[1])] + - [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[1])] + - [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[0])] + - [(0, 0, leg(current_deg, z_mid) * dx[0] * dy[0])]) + [(leg(current_deg, x_mid) * dy[0] * dz[0], 0, 0)] + + [(leg(current_deg, x_mid) * dy[1] * dz[0], 0, 0)] + + [(leg(current_deg, x_mid) * dy[1] * dz[1], 0, 0)] + + [(leg(current_deg, x_mid) * dy[0] * dz[1], 0, 0)]) return ELpiece def e_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): @@ -415,6 +448,7 @@ def I_lambda_1_tilde_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): return ILtilde +#This is always 1-forms regardless of 2 or 3 dimensions. class TrimmedSerendipityEdge(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: @@ -455,13 +489,18 @@ def __init__(self, ref_el, degree): IL = I_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + I_lambda_1_tilde_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) - else: - IL = () - bdmce_list = EL + FL + IL - self.basis = {(0, 0): Array(bdmce_list)} + else: + IL = () + + Sminus_list = EL + FL + if dim == 3: + Sminus_list = Sminus_list + IL + + self.basis = {(0, 0): Array(Sminus_list)} super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") +#This is 1 forms in 2d (rotated) and 2 forms in 3d (not just a rotation) class TrimmedSerendipityFace(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: @@ -502,12 +541,9 @@ def __init__(self, ref_el, degree): IL = I_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + I_lambda_1_tilde_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) - else: - IL = () + else: + IL = () - print(len(EL)) - print(len(FL)) - print(len(IL)) bdmcf_list = EL + FL + IL bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} From 7965603f4e71bd73947d6d1f2e7dd4ac16a37a87 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Mon, 13 Apr 2020 11:31:58 -0700 Subject: [PATCH 12/38] Sminus.py more changes. --- FIAT/Sminus.py | 23 +++++------------------ 1 file changed, 5 insertions(+), 18 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 588afdacd..992b8b2f3 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -525,26 +525,13 @@ def __init__(self, ref_el, degree): dz = None z_mid = None - if (dim == 2): - EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) - else: - EL = e_lambda_1_3d_trimmed(degree, dx, dy, dz, x_mid, y_mid, z_mid) + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: - if dim == 2: - FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) - else: - FL = trimmed_f_lambda_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) else: FL = () - if dim == 3: - if degree >= 4: - IL = I_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + I_lambda_1_tilde_3d(degree, dx, dy, - dz, x_mid, - y_mid, z_mid) - else: - IL = () - bdmcf_list = EL + FL + IL - bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] - self.basis = {(0, 0): Array(bdmcf_list)} + Sminus_list = EL + FL + Sminus_list = [[-a[1], a[0]] for a in Sminus_list] + self.basis = {(0, 0): Array(Sminus_list)} super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") From da802537613284c35e1882bdcf66551bd40593f8 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Tue, 14 Apr 2020 11:17:18 -0700 Subject: [PATCH 13/38] Continuing to work on 3d codes for Sminus elements, including the plumbing for SminusDiv elements. --- FIAT/Sminus.py | 12 +- FIAT/SminusDiv.py | 361 ++++++++++++++++++++++++++++++++++++++++++++++ FIAT/__init__.py | 1 + 3 files changed, 365 insertions(+), 9 deletions(-) create mode 100644 FIAT/SminusDiv.py diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 992b8b2f3..ea5c2a7d5 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -496,11 +496,12 @@ def __init__(self, ref_el, degree): if dim == 3: Sminus_list = Sminus_list + IL + print(dim) + self.basis = {(0, 0): Array(Sminus_list)} super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") -#This is 1 forms in 2d (rotated) and 2 forms in 3d (not just a rotation) class TrimmedSerendipityFace(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: @@ -509,8 +510,7 @@ def __init__(self, ref_el, degree): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - if dim != 3: - raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") + raise Exception("Trimmed serendipity face elements only valid for dimensions 2") verts = flat_el.get_vertices() @@ -518,12 +518,6 @@ def __init__(self, ref_el, degree): dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) x_mid = 2*x-(verts[-1][0] + verts[0][0]) y_mid = 2*y-(verts[-1][1] + verts[0][1]) - try: - dz = ((verts[-1][2] - z)/(verts[-1][2] - verts[0][2]), (z - verts[0][2])/(verts[-1][2] - verts[0][2])) - z_mid = 2*z-(verts[-1][2] + verts[0][2]) - except IndexError: - dz = None - z_mid = None EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py new file mode 100644 index 000000000..74a275c22 --- /dev/null +++ b/FIAT/SminusDiv.py @@ -0,0 +1,361 @@ +#Working on 3d trimmed serendipity 2 forms now. +from sympy import symbols, legendre, Array, diff +import numpy as np +from FIAT.finite_element import FiniteElement +from FIAT.dual_set import make_entity_closure_ids +from FIAT.polynomial_set import mis +from FIAT.reference_element import compute_unflattening_map, flatten_reference_cube + +x, y, z = symbols('x y z') +variables = (x, y, z) +leg = legendre + + +def triangular_number(n): + return int((n+1)*n/2) + +def choose_ijk_total(degree): + top = 1 + for i in range(1, 2 + degree + 1): + top = i * top + bottom = 1 + for i in range(1, degree + 1): + bottom = i * bottom + return int(top /(2 * bottom)) + + +class TrimmedSerendipity(FiniteElement): + def __init__(self, ref_el, degree, mapping): + if degree < 1: + raise Exception("Trimmed serendipity elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 3: + if dim !=2: + raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3") + + flat_topology = flat_el.get_topology() + entity_ids = {} + cur = 0 + + for top_dim, entities in flat_topology.items(): + entity_ids[top_dim] = {} + for entity in entities: + entity_ids[top_dim][entity] = [] + + #3-d case. + if dim == 3: + entity_ids[3] = {} + for j in sorted(flat_topology[2]): + entity_ids[2][j] = list(range(cur, cur + triangular_number(degree))) + cur = cur + triangular_number(degree) + interior_ids = 0 + for k in range(2, degree): + interior_ids = interior_ids + 3 * choose_ijk_total(k - 2) + if (degree > 1): + interior_tilde_ids = 3 + for k in range(1, degree - 1): + interior_tilde_ids = interior_tilde_ids + 3 + if (degree > 3): + interior_tilde_ids += choose_ijk_total(degree - 2) - (2 * degree - 1) + if degree == 1: + interior_tilde_ids = 0 + + entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) + cur = cur + interior_ids + interior_tilde_ids + else: + for j in sorted(flat_topology[1]): + entity_ids[1][j] = list(range(cur, cur + degree)) + cur = cur + degree + + if(degree >= 2): + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + + cur += 2*triangular_number(degree - 2) + degree + + formdegree = dim - 1 + + entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) + + super(TrimmedSerendipity, self).__init__(ref_el=ref_el, + dual=None, + order=degree, + formdegree=formdegree, + mapping=mapping) + + topology = ref_el.get_topology() + unflattening_map = compute_unflattening_map(topology) + unflattened_entity_ids = {} + unflattened_entity_closure_ids = {} + + for dim, entities in sorted(topology.items()): + unflattened_entity_ids[dim] = {} + unflattened_entity_closure_ids[dim] = {} + for dim, entities in sorted(flat_topology.items()): + for entity in entities: + unflat_dim, unflat_entity = unflattening_map[(dim, entity)] + unflattened_entity_ids[unflat_dim][unflat_entity] = entity_ids[dim][entity] + unflattened_entity_closure_ids[unflat_dim][unflat_entity] = entity_closure_ids[dim][entity] + self.entity_ids = unflattened_entity_ids + self.entity_closure_ids = unflattened_entity_closure_ids + self._degree = degree + self.flat_el = flat_el + + def degree(self): + return self._degree + + def get_nodal_basis(self): + raise NotImplementedError("get_nodal_basis not implemented for trimmed serendipity") + + def get_dual_set(self): + raise NotImplementedError("get_dual_set is not implemented for trimmed serendipity") + + def get_coeffs(self): + raise NotImplementedError("get_coeffs not implemented for trimmed serendipity") + + def tabulate(self, order, points, entity=None): + + if entity is None: + entity = (self.ref_el.get_dimension(), 0) + + entity_dim, entity_id = entity + transform = self.ref_el.get_entity_transform(entity_dim, entity_id) + points = list(map(transform, points)) + + phivals = {} + + for o in range(order+1): + alphas = mis(2, o) + for alpha in alphas: + try: + polynomials = self.basis[alpha] + except KeyError: + polynomials = diff(self.basis[(0, 0)], *zip(variables, alpha)) + self.basis[alpha] = polynomials + T = np.zeros((len(polynomials[:, 0]), 2, len(points))) + for i in range(len(points)): + subs = {v: points[i][k] for k, v in enumerate(variables[:2])} + for j, f in enumerate(polynomials[:, 0]): + T[j, 0, i] = f.evalf(subs=subs) + for j, f in enumerate(polynomials[:, 1]): + T[j, 1, i] = f.evalf(subs=subs) + phivals[alpha] = T + + return phivals + + def entity_dofs(self): + """Return the map of topological entities to degrees of + freedom for the finite element.""" + return self.entity_ids + + def entity_closure_dofs(self): + """Return the map of topological entities to degrees of + freedom on the closure of those entities for the finite element.""" + return self.entity_closure_ids + + def value_shape(self): + return (2,) + + def dmats(self): + raise NotImplementedError + + def get_num_members(self, arg): + raise NotImplementedError + + def space_dimension(self): + return int(len(self.basis[(0, 0)])/2) + + +class TrimmedSerendipityDiv(TrimmedSerendipity): + def __init__(self, ref_el, degree): + if degree < 1: + raise Exception("Trimmed serendipity face elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + if dim !=3: + raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") + + verts = flat_el.get_vertices() + + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) + dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + try: + dz = ((verts[-1][2] - z)/(verts[-1][2] - verts[0][2]), (z - verts[0][2])/(verts[-1][2] - verts[0][2])) + z_mid = 2*z-(verts[-1][2] + verts[0][2]) + except IndexError: + dz = None + z_mid = None + + if dim == 3: + FL = f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + if (degree > 1): + IL = I_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + else: + IL = () + + Sminus_list = FL + IL + self.basis = {(0, 0): Array(Sminus_list)} + super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + + else: + ##Put all 2 dimensional stuff here. + if degree < 1: + raise Exception("Trimmed serendipity face elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + verts = flat_el.get_vertices() + + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + if degree >= 2: + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + else: + FL = () + + Sminus_list = EL + FL + Sminus_list = [[-a[1], a[0]] for a in Sminus_list] + self.basis = {(0, 0): Array(Sminus_list)} + super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + + +def f_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): + FLpiece = tuple([]) + for j in range(0, current_deg + 1): + FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(current_deg - j, y_mid) * c) for c in dz] + + [(0, leg(j, x_mid) * leg(current_deg - j, z_mid) * b, 0) for b in dy] + + [(leg(j, y_mid) * leg(current_deg - j, z_mid) * a, 0, 0) for a in dx]) + return FLpiece + + +def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): + FL = tuple([]) + for j in range(0, degree): + FL += f_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) + return FL + + +def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): + assert current_deg > 1, 'invalid for i = 1' + ILpiece = tuple([]) + for j in range(0, deg -1): + for k in range(0, deg - 1 - j): + ILpiece += tuple([(0, 0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + dz[0] * dz[1])]+ + [(0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + dy[0] * dy[1] ,0)] + + [(-leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * + dx[1], 0, 0)]) + return ILpiece + + +def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): + assert degree > 1, 'invalid for i = 1' + IL_tilde = tuple([(-leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + + [(0, -leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + + [(0, 0, -leg(degree - 2, z_mid) * dz[0] * dz[1])]) + IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], -leg(degree - j - 1, x_mid) * + leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + + [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * + leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)] + + [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * d[y0] * dy[1], -leg(degree - j - 1, y_mid) * + leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) + for k in range(1, degree - 2): + for l in range(1, degree - 2 - k): + j = degree - 2 - k - l + IL_tilde += tuple([(-leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], + leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], + leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) + return IL_tilde + + +def I_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): + IL = tuple([]) + for j in range(2, degree): + IL += I_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) + IL += I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid) + return IL + +#Everything for 2-d should work already. +def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): + EL = tuple( + [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)] + + [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)] + + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]) + + return EL + + +def e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid): + ELTilde = tuple([(-leg(deg, x_mid) * dy[0], + -leg(deg-1, x_mid) * dx[0] * dx[1] / (deg+1))] + + [(-leg(deg, x_mid) * dy[1], + leg(deg-1, x_mid) * dx[0] * dx[1] / (deg+1))] + + [(-leg(deg-1, y_mid) * dy[0] * dy[1] / (deg+1), + -leg(deg, y_mid) * dx[0])] + + [(leg(deg-1, y_mid) * dy[0] * dy[1] / (deg+1), + -leg(deg, y_mid) * dx[1])]) + return ELTilde + + +def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + EL = e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid) + ELTilde = e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid) + + result = EL + ELTilde + return result + + +def determine_f_lambda_portions_2d(deg): + if (deg < 2): + DegsOfIteration = [] + else: + DegsOfIteration = [] + for i in range(2, deg): + DegsOfIteration += [i] + + return DegsOfIteration + + +def f_lambda_1_2d_pieces(current_deg, dx, dy, x_mid, y_mid): + if (current_deg == 2): + FLpiece = [(leg(0, x_mid) * leg(0, y_mid) * dy[0] * dy[1], 0)] + FLpiece += [(0, leg(0, x_mid) * leg(0, y_mid) * dx[0] * dx[1])] + else: + target_power = current_deg - 2 + FLpiece = tuple([]) + for j in range(0, target_power + 1): + k = target_power - j + FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dy[0] * dy[1], 0)]) + FLpiece += tuple([(0, leg(j, x_mid) * leg(k, y_mid) * dx[0] * dx[1])]) + return FLpiece + + +def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): + DegsOfIteration = determine_f_lambda_portions_2d(deg) + FL = [] + for i in DegsOfIteration: + FL += f_lambda_1_2d_pieces(i, dx, dy, x_mid, y_mid) + return tuple(FL) + + +def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): + FLTilde = tuple([]) + FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) + FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) + for k in range(1, deg - 1): + FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) + + return tuple(FLTilde) + + +def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): + FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) + FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) + result = FL + FLT + + return result \ No newline at end of file diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 0cd27f747..6ac974573 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -13,6 +13,7 @@ from FIAT.brezzi_douglas_marini import BrezziDouglasMarini from FIAT.Sminus import TrimmedSerendipityEdge # noqa: F401 from FIAT.Sminus import TrimmedSerendipityFace # noqa: F401 +from FIAT.SminusDiv import TrimmedSerendipityDiv #noqa: F401 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor From 2d4d38b962535ece4c3ec2e3ea235ce7419863dd Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Mon, 20 Apr 2020 14:53:27 -0500 Subject: [PATCH 14/38] Work on fixing some dimension stuff --- FIAT/SminusDiv.py | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 74a275c22..d9ef29221 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -31,6 +31,7 @@ def __init__(self, ref_el, degree, mapping): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() + self.fdim = dim if dim != 3: if dim !=2: raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3") @@ -124,22 +125,24 @@ def tabulate(self, order, points, entity=None): points = list(map(transform, points)) phivals = {} - + for o in range(order+1): - alphas = mis(2, o) + alphas = mis(self.fdim, o) for alpha in alphas: try: polynomials = self.basis[alpha] except KeyError: - polynomials = diff(self.basis[(0, 0)], *zip(variables, alpha)) + zr = tuple([0] * self.fdim) + polynomials = diff(self.basis[zr], *zip(variables, alpha)) self.basis[alpha] = polynomials - T = np.zeros((len(polynomials[:, 0]), 2, len(points))) + T = np.zeros((len(polynomials[:, 0]), self.fdim, len(points))) for i in range(len(points)): - subs = {v: points[i][k] for k, v in enumerate(variables[:2])} - for j, f in enumerate(polynomials[:, 0]): - T[j, 0, i] = f.evalf(subs=subs) - for j, f in enumerate(polynomials[:, 1]): - T[j, 1, i] = f.evalf(subs=subs) + subs = {v: points[i][k] for k, v in enumerate(variables[:self.fdim])} + for ell in range(self.fdim): + for j, f in enumerate(polynomials[:, ell]): + print(subs) + print(f.subs(subs)) + T[j, ell, i] = f.evalf(subs=subs) phivals[alpha] = T return phivals @@ -164,7 +167,7 @@ def get_num_members(self, arg): raise NotImplementedError def space_dimension(self): - return int(len(self.basis[(0, 0)])/2) + return int(len(self.basis[tuple([0] * self.fdim)])/2) class TrimmedSerendipityDiv(TrimmedSerendipity): @@ -199,7 +202,7 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = FL + IL - self.basis = {(0, 0): Array(Sminus_list)} + self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") else: @@ -358,4 +361,4 @@ def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) result = FL + FLT - return result \ No newline at end of file + return result From e08c005d0511a8c37fbc12de663f02ee4643de7e Mon Sep 17 00:00:00 2001 From: Justincrum Date: Mon, 20 Apr 2020 14:05:06 -0700 Subject: [PATCH 15/38] more changes. --- FIAT/SminusDiv.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 74a275c22..53e1489a6 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -125,8 +125,10 @@ def tabulate(self, order, points, entity=None): phivals = {} + print(entity_dim) + for o in range(order+1): - alphas = mis(2, o) + alphas = mis(self.fdim, o) for alpha in alphas: try: polynomials = self.basis[alpha] From b4d2e8563a1d794d4e416ba531c6a62b5ca3ba2e Mon Sep 17 00:00:00 2001 From: Justincrum Date: Mon, 20 Apr 2020 14:28:15 -0700 Subject: [PATCH 16/38] Fixed the sizing issues in SminusDiv.py to account for the different dimensions being used by the elements. --- FIAT/SminusDiv.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index a3671a447..3b58a6901 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -140,8 +140,6 @@ def tabulate(self, order, points, entity=None): subs = {v: points[i][k] for k, v in enumerate(variables[:self.fdim])} for ell in range(self.fdim): for j, f in enumerate(polynomials[:, ell]): - print(subs) - print(f.subs(subs)) T[j, ell, i] = f.evalf(subs=subs) phivals[alpha] = T @@ -158,7 +156,7 @@ def entity_closure_dofs(self): return self.entity_closure_ids def value_shape(self): - return (2,) + return (self.fdim,) def dmats(self): raise NotImplementedError @@ -167,7 +165,7 @@ def get_num_members(self, arg): raise NotImplementedError def space_dimension(self): - return int(len(self.basis[tuple([0] * self.fdim)])/2) + return int(len(self.basis[tuple([0] * self.fdim)])/self.fdim) class TrimmedSerendipityDiv(TrimmedSerendipity): @@ -203,6 +201,7 @@ def __init__(self, ref_el, degree): Sminus_list = FL + IL self.basis = {(0, 0, 0): Array(Sminus_list)} + print(len(Array(Sminus_list))) super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") else: From 95952c989b637faeba229470858401af35094ffe Mon Sep 17 00:00:00 2001 From: Justincrum Date: Mon, 27 Apr 2020 13:19:23 -0700 Subject: [PATCH 17/38] Updating face association for SminusDiv elements. --- FIAT/SminusDiv.py | 48 ++++++++++++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 13 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 3b58a6901..6010553fa 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -74,6 +74,9 @@ def __init__(self, ref_el, degree, mapping): entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) cur += 2*triangular_number(degree - 2) + degree + print("entity_ids = :", entity_ids) + print("topology of mesh:", flat_topology) + print("vertex locations:", flat_el.vertices) formdegree = dim - 1 @@ -200,8 +203,8 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = FL + IL + print(Sminus_list) self.basis = {(0, 0, 0): Array(Sminus_list)} - print(len(Array(Sminus_list))) super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") else: @@ -224,19 +227,38 @@ def __init__(self, ref_el, degree): super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") -def f_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): +""" def f_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): FLpiece = tuple([]) for j in range(0, current_deg + 1): FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(current_deg - j, y_mid) * c) for c in dz] + [(0, leg(j, x_mid) * leg(current_deg - j, z_mid) * b, 0) for b in dy] + [(leg(j, y_mid) * leg(current_deg - j, z_mid) * a, 0, 0) for a in dx]) + return FLpiece """ + + +def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): + FLpiece = tuple([]) + #Dx's first. + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(0, leg(j, x_mid) * leg(i - j, z_mid) * dy[k], 0)]) + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) return FLpiece def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): FL = tuple([]) - for j in range(0, degree): - FL += f_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) + #for j in range(0, degree): + # FL += f_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) + FL += f_lambda_2_3d_pieces(degree, dx, dy, dz, x_mid, y_mid, z_mid) return FL @@ -245,30 +267,30 @@ def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): ILpiece = tuple([]) for j in range(0, deg -1): for k in range(0, deg - 1 - j): - ILpiece += tuple([(0, 0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + ILpiece += tuple([(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dz[0] * dz[1])]+ - [(0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + [(0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dy[0] * dy[1] ,0)] + - [(-leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * + [(leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * dx[1], 0, 0)]) return ILpiece def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): assert degree > 1, 'invalid for i = 1' - IL_tilde = tuple([(-leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + - [(0, -leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + - [(0, 0, -leg(degree - 2, z_mid) * dz[0] * dz[1])]) - IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], -leg(degree - j - 1, x_mid) * + IL_tilde = tuple([(leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + + [(0, leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + + [(0, 0, leg(degree - 2, z_mid) * dz[0] * dz[1])]) + IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], leg(degree - j - 1, x_mid) * leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)] + - [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * d[y0] * dy[1], -leg(degree - j - 1, y_mid) * + [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * d[y0] * dy[1], leg(degree - j - 1, y_mid) * leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) for k in range(1, degree - 2): for l in range(1, degree - 2 - k): j = degree - 2 - k - l - IL_tilde += tuple([(-leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], + IL_tilde += tuple([(leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) return IL_tilde From 155a63393ff50bc265a5c21505fb92255ef25d85 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Thu, 30 Apr 2020 14:05:07 -0700 Subject: [PATCH 18/38] Changes to fix the association of entities and basis functions. This should be working correctly now for both 2 and 3 dimensions. --- FIAT/Sminus.py | 250 +++++++++++++++++++++++----------------------- FIAT/SminusDiv.py | 2 +- 2 files changed, 128 insertions(+), 124 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index ea5c2a7d5..c587587b7 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -49,6 +49,7 @@ def __init__(self, ref_el, degree, mapping): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() + self.fdim = dim if dim != 2: if dim != 3: raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3") @@ -103,6 +104,9 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur)) formdegree = 1 + print("entities", entity_ids) + print("topology", flat_topology) + print("vertices", flat_el.vertices) entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) @@ -154,24 +158,24 @@ def tabulate(self, order, points, entity=None): phivals = {} for o in range(order+1): - alphas = mis(2, o) + alphas = mis(self.fdim, o) for alpha in alphas: try: polynomials = self.basis[alpha] except KeyError: - polynomials = diff(self.basis[(0, 0)], *zip(variables, alpha)) + zr = tuple([0] * self.fdim) + polynomials = diff(self.basis[zr], *zip(variables, alpha)) self.basis[alpha] = polynomials - T = np.zeros((len(polynomials[:, 0]), 2, len(points))) + T = np.zeros((len(polynomials[:, 0]), self.fdim, len(points))) for i in range(len(points)): - subs = {v: points[i][k] for k, v in enumerate(variables[:2])} - for j, f in enumerate(polynomials[:, 0]): - T[j, 0, i] = f.evalf(subs=subs) - for j, f in enumerate(polynomials[:, 1]): - T[j, 1, i] = f.evalf(subs=subs) + subs = {v: points[i][k] for k, v in enumerate(variables[:self.fdim])} + for ell in range(self.fdim): + for j, f in enumerate(polynomials[:, ell]): + T[j, ell, i] = f.evalf(subs=subs) phivals[alpha] = T - return phivals + def entity_dofs(self): """Return the map of topological entities to degrees of freedom for the finite element.""" @@ -183,7 +187,7 @@ def entity_closure_dofs(self): return self.entity_closure_ids def value_shape(self): - return (2,) + return (self.fdim,) def dmats(self): raise NotImplementedError @@ -192,7 +196,7 @@ def get_num_members(self, arg): raise NotImplementedError def space_dimension(self): - return int(len(self.basis[(0, 0)])/2) + return int(len(self.basis[tuple([0] * self.fdim)])/self.fdim) # Splitting the E Lambda function into two seperate functions for E Lambda and E tilde Lambda. @@ -281,123 +285,120 @@ def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): return result -""" def e_lambda_1_3d_part_one(deg, dx, dy, dz, x_mid, y_mid, z_mid): - EL = tuple( - [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)] + - [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)] + - [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + - [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]) - - return EL """ - - -def e_lambda_1_3d_piece(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): - ELpiece = tuple([(0, 0, leg(current_deg, z_mid) * dx[0] * dy[1])] + - [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[1])] + - [(0, 0, leg(current_deg, z_mid) * dx[1] * dy[0])] + - [(0, 0, leg(current_deg, z_mid) * dx[0] * dy[0])] + - [(0, leg(current_deg, y_mid) * dx[0] * dz[1], 0)] + - [(0, leg(current_deg, y_mid) * dx[1] * dz[1], 0)] + - [(0, leg(current_deg, y_mid) * dx[1] * dz[0], 0)] + - [(0, leg(current_deg, y_mid) * dx[0] * dz[0], 0)] + - [(leg(current_deg, x_mid) * dy[0] * dz[0], 0, 0)] + - [(leg(current_deg, x_mid) * dy[1] * dz[0], 0, 0)] + - [(leg(current_deg, x_mid) * dy[1] * dz[1], 0, 0)] + - [(leg(current_deg, x_mid) * dy[0] * dz[1], 0, 0)]) - return ELpiece - def e_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): EL = tuple([]) - for i in range(max_deg): - Elpiece = e_lambda_1_3d_piece(i, dx, dy, dz, x_mid, y_mid, z_mid) - EL += Elpiece + ##assignment to edge x=y=0 + for i in range (0, max_deg): + EL += tuple([(0, 0, leg(i, z_mid) * dx[0] * dy[0])]) + ##assignment to edge x=0, y=1 + for i in range(0, max_deg): + EL += tuple([(0, 0, leg(i, z_mid) * dy[1] * dx[0])]) + ##assignment to edge x = 1, y = 0 + for i in range(0, max_deg): + EL += tuple([(0, 0, leg(i, z_mid) * dx[1] * dy[0])]) + ##assignment to edge x = 1, y = 1 + for i in range(0, max_deg): + EL += tuple([(0, 0, leg(i, z_mid) * dx[1] * dy[1])]) + ##assignment to edge x = 0, z = 0 + for i in range(0, max_deg): + EL += tuple([(0, leg(i, y_mid) * dx[0] * dz[0], 0)]) + ##assignment to edge x = 0, z = 1 + for i in range(0, max_deg): + EL += tuple([(0, leg(i, y_mid) * dx[0] * dz[1], 0)]) + ##assignment to edge x = 1, z = 0 + for i in range(0, max_deg): + EL += tuple([(0, leg(i, y_mid) * dx[1] * dz[0], 0)]) + ##assignment to edge x = 1, z = 1 + for i in range(0, max_deg): + EL += tuple([(0, leg(i, y_mid) * dx[1] * dz[1], 0)]) + ##assignment to edge y = 0, z = 0 + for i in range(0, max_deg): + EL += tuple([(leg(i, x_mid) * dy[0] * dz[0], 0, 0)]) + ##assignment to edge y = 0, z = 1 + for i in range(0, max_deg): + EL += tuple([(leg(i, x_mid) * dy[0] * dz[1], 0, 0)]) + ##assignment to edge y = 1, z = 0 + for i in range(0, max_deg): + EL += tuple([(leg(i, x_mid) * dy[1] * dz[0], 0, 0)]) + ##assignment to edge y = 1, z = 1 + for i in range(0, max_deg): + EL += tuple([(leg(i, x_mid) * dy[1] * dz[1], 0, 0)]) return EL -def determine_f_lambda_1_portions_3d(deg): - if (deg < 2): - DegsOfIteration = [] - else: - DegsOfIteration = [] - for i in range(2, deg): - DegsOfIteration += [i] - - return DegsOfIteration - -def f_lambda_1_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): - if (current_deg == 2): - j = 0 - k = 0 - FLpiece = [(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)] - FLpiece += [(leg(j, x_mid) * leg(k, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)] - FLpiece += [(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)] - FLpiece += [(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)] - FLpiece += [(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)] - FLpiece += [(0, leg(j, y_mid) * leg(k, x_mid) * dz[1] * dx[0] * dx[1], 0)] - FLpiece += [(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)] - FLpiece += [(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)] - FLpiece += [(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[0] * dx[0] * dx[1])] - FLpiece += [(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[1] * dx[0] * dx[1])] - FLpiece += [(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[0] * dy[0] * dy[1])] - FLpiece += [(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[1] * dy[0] * dy[1])] - else: - target_power = current_deg - 2 - FLpiece = tuple([]) - for j in range(0, target_power + 1): - k = target_power - j - FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) - FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) - FLpiece += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) - FLpiece += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) - FLpiece += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)]) - FLpiece += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[1] * dx[0] * dx[1], 0)]) - FLpiece += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) - FLpiece += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)]) - FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[0] * dx[0] * dx[1])]) - FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[1] * dx[0] * dx[1])]) - FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[0] * dy[0] * dy[1])]) - FLpiece += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[1] * dy[0] * dy[1])]) - return FLpiece +def f_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): + FL = tuple([]) + ##Assignment to face x = 0, Ftilde + FL += tuple([(0, leg(max_deg - 2, z_mid) * dx[0] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[0] * dy[0] * dy[1])]) + for j in range(1, max_deg - 1): + FL += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[0] * dz[0] * dz[1], + -leg(j - 1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[0] * dy[0] * dy[1])]) + ##Assignment to face x = 0, F + for j in range(1, max_deg - 1): + k = max_deg - j - 2 + FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[0] * dy[0] * dy[1])]) + ##Assignment to face x = 1, Ftilde + FL += tuple([(0, leg(max_deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) + for j in range(1, max_deg - 1): + FL += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[1] * dz[0] * dz[1], + -leg(j - 1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[1] * dy[0] * dy[1])]) + ##Assignment to face x = 1, F + for j in range(1, max_deg - 1): + k = max_deg - j - 2 + FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[1] * dy[0] * dy[1])]) -def f_lambda_1_3d_tilde(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): - FLtilde = tuple([]) - FLtilde += tuple([(leg(max_deg - 2, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) - FLtilde += tuple([(leg(max_deg - 2, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) - FLtilde += tuple([(leg(max_deg - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) - FLtilde += tuple([(leg(max_deg - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) - FLtilde += tuple([(0, leg(max_deg - 2, x_mid) * dz[0] * dx[0] * dx[1], 0)]) - FLtilde += tuple([(0, leg(max_deg - 2, x_mid) * dz[1] * dx[0] * dx[1], 0)]) - FLtilde += tuple([(0, leg(max_deg - 2, z_mid) * dx[0] * dz[0] * dz[1], 0)]) - FLtilde += tuple([(0, leg(max_deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) - FLtilde += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[0] * dx[0] * dx[1])]) - FLtilde += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[1] * dx[0] * dx[1])]) - FLtilde += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[0] * dy[0] * dy[1])]) - FLtilde += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) + ##Assignment to face y = 0, Ftilde + FL += tuple([(leg(max_deg - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[0] * dx[0] * dx[1])]) for j in range(1, max_deg - 1): - for k in range(0, 2): - FLtilde += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, y_mid) * dz[k] * dy[0] * dy[1], - -leg(j-1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[k] * dx[0] * dx[1], 0)]) - FLtilde += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[k] * dz[0] * dz[1], 0, - -leg(j-1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[k] * dx[0] * dx[1])]) - FLtilde += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[k] * dz[0] * dz[1], - -leg(j-1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[k] * dy[0] * dy[1])]) - return FLtilde - - -def f_lambda_1_3d_trim(deg, dx, dy, dz, x_mid, y_mid, z_mid): - DegsOfIteration = determine_f_lambda_1_portions_3d(deg) - FL = [] - for i in DegsOfIteration: - FL += f_lambda_1_3d_pieces(i, dx, dy, dz, x_mid, y_mid, z_mid) - return tuple(FL) + FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, + -leg(j - 1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[0] * dx[0] * dx[1])]) + ##Assignment to face y = 0, F + for j in range(1, max_deg - 1): + k = max_deg - j - 2 + FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[0] * dx[0] * dx[1])]) + ##Assignment to face y = 1, Ftilde + FL += tuple([(leg(max_deg - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[1] * dx[0] * dx[1])]) + for j in range(1, max_deg - 1): + FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, + -leg(j - 1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[1] * dx[0] * dx[1])]) + ##Assignment to face y = 1, F + for j in range(1, max_deg - 1): + k = max_deg - j - 2 + FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[1] * dx[0] * dx[1])]) -def trimmed_f_lambda_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): - FL = f_lambda_1_3d_trim(deg, dx, dy, dz, x_mid, y_mid, z_mid) - FLT = f_lambda_1_3d_tilde(deg, dx, dy, dz, x_mid, y_mid, z_mid) - result = FL + FLT + ##Assignment to face z = 0, Ftilde + FL += tuple([(leg(max_deg - 2, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(max_deg - 2, x_mid) * dz[0] * dx[0] * dx[1], 0)]) + for j in range(1, max_deg - 1): + FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, y_mid) * dz[0] * dy[0] * dy[1], + -leg(j - 1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[0] * dx[0] * dx[1], 0)]) + ##Assignment to face z = 0, F + for j in range(1, max_deg - 1): + k = max_deg - j - 2 + FL += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)]) - return result + ##Assignment to face z = 1, Ftilde + FL += tuple([(leg(max_deg - 2, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(max_deg - 2, x_mid) * dz[1] * dx[0] * dx[1], 0)]) + for j in range(1, max_deg - 1): + FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, y_mid) * dz[1] * dy[0] * dy[1], + -leg(j - 1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[1] * dx[0] * dx[1], 0)]) + ##Assignment to face z = 1, F + for j in range(1, max_deg - 1): + k = max_deg - j - 2 + FL += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[1] * dx[0] * dx[1], 0)]) + return FL def determine_I_lambda_1_portions_3d(deg): @@ -481,7 +482,7 @@ def __init__(self, ref_el, degree): if dim == 2: FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) else: - FL = trimmed_f_lambda_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + FL = f_lambda_1_3d_trimmed(degree, dx, dy, dz, x_mid, y_mid, z_mid) else: FL = () if dim == 3: @@ -495,10 +496,13 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL if dim == 3: Sminus_list = Sminus_list + IL + print("EL", EL) + print("FL", FL) - print(dim) - - self.basis = {(0, 0): Array(Sminus_list)} + if dim == 2: + self.basis = {(0, 0): Array(Sminus_list)} + else: + self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") @@ -528,4 +532,4 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL Sminus_list = [[-a[1], a[0]] for a in Sminus_list] self.basis = {(0, 0): Array(Sminus_list)} - super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") \ No newline at end of file diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 6010553fa..f198445f6 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -77,7 +77,7 @@ def __init__(self, ref_el, degree, mapping): print("entity_ids = :", entity_ids) print("topology of mesh:", flat_topology) print("vertex locations:", flat_el.vertices) - + formdegree = dim - 1 entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) From e796ef5fe1bf73c2b61c5159007c6a8d38b20383 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Thu, 30 Apr 2020 14:07:59 -0700 Subject: [PATCH 19/38] Removing extraneous print messages. --- FIAT/Sminus.py | 5 ----- FIAT/SminusDiv.py | 4 ---- 2 files changed, 9 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index c587587b7..728886914 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -104,9 +104,6 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur)) formdegree = 1 - print("entities", entity_ids) - print("topology", flat_topology) - print("vertices", flat_el.vertices) entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) @@ -496,8 +493,6 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL if dim == 3: Sminus_list = Sminus_list + IL - print("EL", EL) - print("FL", FL) if dim == 2: self.basis = {(0, 0): Array(Sminus_list)} diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index f198445f6..d48fa0c52 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -74,9 +74,6 @@ def __init__(self, ref_el, degree, mapping): entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) cur += 2*triangular_number(degree - 2) + degree - print("entity_ids = :", entity_ids) - print("topology of mesh:", flat_topology) - print("vertex locations:", flat_el.vertices) formdegree = dim - 1 @@ -203,7 +200,6 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = FL + IL - print(Sminus_list) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") From 8c899b7f6066f4941e718c9e8ff9429081eb9de2 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Tue, 8 Sep 2020 13:25:33 -0700 Subject: [PATCH 20/38] More changes for SminusCurl. --- FIAT/SminusCurl.py | 588 +++++++++++++++++++++++++++++++++++++++++++++ FIAT/SminusDiv.py | 6 +- FIAT/__init__.py | 1 + 3 files changed, 592 insertions(+), 3 deletions(-) create mode 100644 FIAT/SminusCurl.py diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py new file mode 100644 index 000000000..66ef97aad --- /dev/null +++ b/FIAT/SminusCurl.py @@ -0,0 +1,588 @@ +#Working on 3d trimmed serendipity 2 forms now. +from sympy import symbols, legendre, Array, diff +import numpy as np +from FIAT.finite_element import FiniteElement +from FIAT.dual_set import make_entity_closure_ids +from FIAT.polynomial_set import mis +from FIAT.reference_element import compute_unflattening_map, flatten_reference_cube + +x, y, z = symbols('x y z') +variables = (x, y, z) +leg = legendre + + +def triangular_number(n): + return int((n+1)*n/2) + +def choose_ijk_total(degree): + top = 1 + for i in range(1, 2 + degree + 1): + top = i * top + bottom = 1 + for i in range(1, degree + 1): + bottom = i * bottom + return int(top /(2 * bottom)) + + +class TrimmedSerendipity(FiniteElement): + def __init__(self, ref_el, degree, mapping): + if degree < 1: + raise Exception("Trimmed serendipity elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + self.fdim = dim + if dim != 3: + if dim !=2: + raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3") + + flat_topology = flat_el.get_topology() + entity_ids = {} + cur = 0 + + for top_dim, entities in flat_topology.items(): + entity_ids[top_dim] = {} + for entity in entities: + entity_ids[top_dim][entity] = [] + + #3-d case. + if dim == 3: + entity_ids[3] = {} + for l in sorted(flat_topology[1]): + entity_ids[1][l] = list(range(cur, cur + degree)) + cur = cur + degree + if (degree > 1): + face_ids = degree + for k in range(2, degree): + face_ids += 2 * (k - 1) + for j in sorted(flat_topology[2]): + entity_ids[2][j] = list(range(cur, cur + face_ids)) + cur += face_ids + + interior_ids = 0 + for k in range(4, degree): + interior_ids = interior_ids + 3 * choose_ijk_total(k - 4) + if (degree > 3): + if (degree == 4): + interior_tilde_ids = 6 + elif (degree == 5): + interior_tilde_ids = 8 + else: + interior_tilde_ids = 6 + 3 * (degree - 4) + else: + interior_tilde_ids = 0 + + entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) + cur = cur + interior_ids + interior_tilde_ids + else: + for j in sorted(flat_topology[1]): + entity_ids[1][j] = list(range(cur, cur + degree)) + cur = cur + degree + + if(degree >= 2): + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + + cur += 2*triangular_number(degree - 2) + degree + formdegree = 1 + + entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) + + super(TrimmedSerendipity, self).__init__(ref_el=ref_el, + dual=None, + order=degree, + formdegree=formdegree, + mapping=mapping) + + topology = ref_el.get_topology() + unflattening_map = compute_unflattening_map(topology) + unflattened_entity_ids = {} + unflattened_entity_closure_ids = {} + + for dim, entities in sorted(topology.items()): + unflattened_entity_ids[dim] = {} + unflattened_entity_closure_ids[dim] = {} + for dim, entities in sorted(flat_topology.items()): + for entity in entities: + unflat_dim, unflat_entity = unflattening_map[(dim, entity)] + unflattened_entity_ids[unflat_dim][unflat_entity] = entity_ids[dim][entity] + unflattened_entity_closure_ids[unflat_dim][unflat_entity] = entity_closure_ids[dim][entity] + self.entity_ids = unflattened_entity_ids + self.entity_closure_ids = unflattened_entity_closure_ids + self._degree = degree + self.flat_el = flat_el + + def degree(self): + return self._degree + + def get_nodal_basis(self): + raise NotImplementedError("get_nodal_basis not implemented for trimmed serendipity") + + def get_dual_set(self): + raise NotImplementedError("get_dual_set is not implemented for trimmed serendipity") + + def get_coeffs(self): + raise NotImplementedError("get_coeffs not implemented for trimmed serendipity") + + def tabulate(self, order, points, entity=None): + + if entity is None: + entity = (self.ref_el.get_dimension(), 0) + + entity_dim, entity_id = entity + transform = self.ref_el.get_entity_transform(entity_dim, entity_id) + points = list(map(transform, points)) + + phivals = {} + + for o in range(order+1): + alphas = mis(self.fdim, o) + for alpha in alphas: + try: + polynomials = self.basis[alpha] + except KeyError: + zr = tuple([0] * self.fdim) + polynomials = diff(self.basis[zr], *zip(variables, alpha)) + self.basis[alpha] = polynomials + T = np.zeros((len(polynomials[:, 0]), self.fdim, len(points))) + for i in range(len(points)): + subs = {v: points[i][k] for k, v in enumerate(variables[:self.fdim])} + for ell in range(self.fdim): + for j, f in enumerate(polynomials[:, ell]): + T[j, ell, i] = f.evalf(subs=subs) + phivals[alpha] = T + + return phivals + + def entity_dofs(self): + """Return the map of topological entities to degrees of + freedom for the finite element.""" + return self.entity_ids + + def entity_closure_dofs(self): + """Return the map of topological entities to degrees of + freedom on the closure of those entities for the finite element.""" + return self.entity_closure_ids + + def value_shape(self): + return (self.fdim,) + + def dmats(self): + raise NotImplementedError + + def get_num_members(self, arg): + raise NotImplementedError + + def space_dimension(self): + return int(len(self.basis[tuple([0] * self.fdim)])/self.fdim) + + +class TrimmedSerendipityCurl(TrimmedSerendipity): + #3d Curl element needs to be fixed entirely. This should correspond to 1-forms in 3d. + def __init__(self, ref_el, degree): + if degree < 1: + raise Exception("Trimmed serendipity face elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + if dim !=3: + raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") + + verts = flat_el.get_vertices() + + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) + dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + try: + dz = ((verts[-1][2] - z)/(verts[-1][2] - verts[0][2]), (z - verts[0][2])/(verts[-1][2] - verts[0][2])) + z_mid = 2*z-(verts[-1][2] + verts[0][2]) + except IndexError: + dz = None + z_mid = None + + if dim == 3: + if degree < 1: + raise Exception ("Trimmed Serendipity fce elements only valid for k >= 1") + EL = e_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + if (degree > 1): + FL = f_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + else: + FL = () + if (degree > 3): + IL = I_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) + else: + IL = () + + Sminus_list = EL + FL + IL + self.basis = {(0, 0, 0): Array(Sminus_list)} + super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + + else: + ##Put all 2 dimensional stuff here. + ##This part might be good, need to test it. + if degree < 1: + raise Exception("Trimmed serendipity face elements only valid for k >= 1") + + #flat_el = flatten_reference_cube(ref_el) + #verts = flat_el.get_vertices() + + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) + if degree >= 2: + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + else: + FL = () + + Sminus_list = EL + FL + self.basis = {(0, 0): Array(Sminus_list)} + super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + +#def e_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): +# EL = () +# for j in range(0, deg): +# for k in range(0, 2): +# for l in range(0, 2): +# EL += tuple([(leg(j, x_mid) * dy[k] * dz[l], 0, 0)] + +# [(0, leg(j, y_mid) * dx[k] * dz[l], 0)] + +# [(0, 0, leg(j, z_mid) * dx[k] * dy[l])]) +# return EL + + +def e_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): + EL = () + #IDs associated to edge 0. + for j in range(0, deg): + EL += tuple([(0, 0, leg(j, z_mid) * dx[0] * dy[0])]) + #IDs associated to edge 1. + for j in range(0, deg): + EL += tuple([(0, 0, leg(j, z_mid) * dx[0] * dy[1])]) + #IDs associated to edge 2. + for j in range(0, deg): + EL += tuple([(0, 0, leg(j, z_mid) * dx[1] * dy[0])]) + #IDs associated to edge 3. + for j in range(0, deg): + EL += tuple([(0, 0, leg(j, z_mid) * dx[1] * dy[1])]) + #IDs associated edge 4. + for j in range(0, deg): + EL += tuple([(0, leg(j, y_mid) * dx[0] * dz[0], 0)]) + #IDs associated to edge 5. + for j in range(0, deg): + EL += tuple([(0, leg(j, y_mid) * dx[0] * dz[1], 0)]) + #IDs associated to edge 6. + for j in range(0, deg): + EL += tuple([(0, leg(j, y_mid) * dx[1] * dz[0], 0)]) + #IDs associated to edge 7. + for j in range(0, deg): + EL += tuple([(0, leg(j, y_mid) * dx[1] * dz[1], 0)]) + #IDs associated to edge 8. + for j in range(0, deg): + EL += tuple([(leg(j, x_mid) * dz[0] * dy[0], 0, 0)]) + #IDs associated to edge 9. + for j in range(0, deg): + EL += tuple([(leg(j, x_mid) * dz[1] * dy[0], 0, 0)]) + #IDs associated to edge 10. + for j in range(0, deg): + EL += tuple([(leg(j, x_mid) * dz[0] * dy[1], 0, 0)]) + #IDs associated to edge 11. + for j in range(0, deg): + EL += tuple([(leg(j, x_mid) * dz[1] * dy[1], 0, 0)]) + return EL + + +#def f_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): +# FLpiece = () +# for l in range(0, 2): +# for j in range(0, deg - 2): +# k = deg - 2 - j +# FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[l] * dy[0] * dy[1], 0, 0)] + +# [(leg(j, x_mid) * leg(k, z_mid) * dy[l] * dz[0] * dz[1], 0, 0)] + +# [(0, leg(j, y_mid) * leg(k, x_mid) * dz[l] * dx[0] * dx[1], 0)] + +# [(0, leg(j, y_mid) * leg(k, z_mid) * dx[l] * dz[0] * dz[1], 0)] + +# [(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[l] * dx[0] * dx[1])] + +# [(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[l] * dy[0] * dy[1])]) +# return FLpiece + + +#def f_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): +# FTilde = () +# for l in range(0, 2): +# FTilde += tuple([(leg(deg - 2, y_mid) * dz[l] * dy[0] * dy[1], 0, 0)] + +# [(leg(deg - 2, z_mid) * dy[l] * dz[0] * dz[1], 0, 0)] + +# [(0, leg(deg - 2, x_mid) * dz[l] * dx[0] * dx[1], 0)] + +# [(0, leg(deg - 2, z_mid) * dx[l] * dz[0] * dz[1], 0)] + +# [(0, 0, leg(deg - 2, x_mid) * dy[l] * dx[0] * dx[1])] + +# [(0, 0, leg(deg - 2, y_mid) * dx[l] * dy[0] * dy[1])]) +# for j in range(1, deg - 1): +# for l in range(0, 2): +# FTilde += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[l] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg - j - 1, y_mid) * dz[l] * dx[0] * dx[1], 0)] + +# [(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[l] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[l] * dy[0] * dy[1])]) +# return FTilde + + +def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): + FL = () + #IDs associated to face 0. + #FTilde first. + FL += tuple([(0, leg(deg - 2, z_mid) * dx[0] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(deg - 2, y_mid) * dx[0] * dy[0] * dy[1])]) + for j in range(1, deg - 2): + FL += tuple([(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[0] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[0] * dy[0] * dy[1])]) + #F next + for i in range(2, deg): + for j in range(0, i - 1): + k = i - 2 - j + FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dx[0] * dy[0] * dy[1])]) + #IDs associated to face 1. + FL += tuple([(0, leg(deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) + for j in range(1, deg - 2): + FL += tuple([(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[1] * dy[0] * dy[1])]) + #F next. + for i in range(2, deg): + for j in range(0, i - 1): + k = i - 2 - j + FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dx[1] * dy[0] * dy[1])]) + #IDs associated to face 2. + #FTilde first. + FL += tuple([(leg(deg - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(deg - 2, x_mid) * dy[0] * dx[0] * dx[1])]) + for j in range(1, deg - 2): + FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 1, z_mid) * dy[0] * dx[0] * dx[1])]) + #F next. + for i in range(2, deg): + for j in range(0, i - 1): + k = i - 2 - j + FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dy[0] * dx[0] * dx[1])]) + #IDs associated to face 3. + FL += tuple([(leg(deg - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(deg - 2, x_mid) * dy[1] * dx[0] * dx[1])]) + for j in range(1, deg - 2): + FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 1, z_mid) * dy[1] * dx[0] * dx[1])]) + #F next. + for i in range(2, deg): + for j in range(0, i - 1): + k = i - 2 - j + FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) + FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dy[1] * dx[0] * dx[1])]) + #IDs associated to face 4. + FL += tuple([(leg(deg - 2, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(deg - 2, x_mid) * dz[0] * dx[0] * dx[1], 0)]) + for j in range(1, deg - 2): + FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[0] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg - j - 1, y_mid) * dz[0] * dx[0] * dx[1], 0)]) + for i in range(2, deg): + for j in range(0, i - 1): + k = i - 2 - j + FL += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)]) + #IDs associated to face 5. + FL += tuple([(leg(deg - 2, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(deg - 2, x_mid) * dz[1] * dx[0] * dx[1], 0)]) + for j in range(1, deg - 2): + FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[1] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg -j - 1, y_mid) * dz[1] * dx[0] * dx[1], 0)]) + for i in range(2, deg): + for j in range(0, i - 1): + k = i - 2 - j + FL += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) + FL += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[1] * dx[0] * dx[1], 0)]) + return FL + + + +#def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): +# F = () +# for i in range(2, deg): +# F += f_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid) +# F += f_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid) +# return F + + +def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): + I = () + for j in range(0, deg - 3): + for k in range(0, deg - 3 - j): + l = deg - 4 - j - k + I += tuple([(leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(0, leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + return I + + +def I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): + ITilde = () + ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(0, leg(deg - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(0, 0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dy[0] * dy[1])] + + [(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + for j in range(1, deg - 3): + ITilde += tuple([(leg(j, x_mid) * leg(deg - j - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], -leg(j - 1, x_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(leg(j, x_mid) * leg(deg - j - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 3, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + if deg > 5: + ITilde += [(0, leg(j, y_mid) * leg(deg - j - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1])] + return ITilde + + +def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): + I = () + for i in range(4, deg): + I += I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid) + I += I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid) + return I + + +""" +def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): + FLpiece = tuple([]) + #Dx's first. + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(0, leg(j, x_mid) * leg(i - j, z_mid) * dy[k], 0)]) + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) + return FLpiece + + +def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): + FL = tuple([]) + #for j in range(0, degree): + # FL += f_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) + FL += f_lambda_2_3d_pieces(degree, dx, dy, dz, x_mid, y_mid, z_mid) + return FL + + +def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): + assert current_deg > 1, invalid for i = 1 + ILpiece = tuple([]) + for j in range(0, current_deg -1): + for k in range(0, current_deg - 1 - j): + ILpiece += tuple([(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + dz[0] * dz[1])]+ + [(0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + dy[0] * dy[1] ,0)] + + [(leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * + dx[1], 0, 0)]) + return ILpiece + + + + +def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): + assert degree > 1, invalid for i = 1 + IL_tilde = tuple([(leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + + [(0, leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + + [(0, 0, leg(degree - 2, z_mid) * dz[0] * dz[1])]) + IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], leg(degree - j - 1, x_mid) * + leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + + [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * + leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)] + + [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * dy[0] * dy[1], leg(degree - j - 1, y_mid) * + leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) + for k in range(1, degree - 2): + for l in range(1, degree - 2 - k): + j = degree - 2 - k - l + IL_tilde += tuple([(leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], + leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], + leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) + return IL_tilde + + +def I_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): + IL = tuple([]) + for j in range(2, degree): + IL += I_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) + IL += I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid) + return IL + +""" + +#Everything for 2-d should work already. +def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): + EL = tuple( + [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)] + + [(0, -leg(j, y_mid) * dx[1]) for j in range(deg)] + + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)]) + + return EL + + +def e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid): + ELTilde = tuple([(-leg(deg, x_mid) * dy[0], + -leg(deg-1, x_mid) * dx[0] * dx[1] / (deg+1))] + + [(-leg(deg, x_mid) * dy[1], + leg(deg-1, x_mid) * dx[0] * dx[1] / (deg+1))] + + [(-leg(deg-1, y_mid) * dy[0] * dy[1] / (deg+1), + -leg(deg, y_mid) * dx[0])] + + [(leg(deg-1, y_mid) * dy[0] * dy[1] / (deg+1), + -leg(deg, y_mid) * dx[1])]) + return ELTilde + + +def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + EL = e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid) + ELTilde = e_lambda_tilde_1_2d_part_two(deg, dx, dy, x_mid, y_mid) + + result = EL + ELTilde + return result + + +def determine_f_lambda_portions_2d(deg): + if (deg < 2): + DegsOfIteration = [] + else: + DegsOfIteration = [] + for i in range(2, deg): + DegsOfIteration += [i] + + return DegsOfIteration + + +def f_lambda_1_2d_pieces(current_deg, dx, dy, x_mid, y_mid): + if (current_deg == 2): + FLpiece = [(leg(0, x_mid) * leg(0, y_mid) * dy[0] * dy[1], 0)] + FLpiece += [(0, leg(0, x_mid) * leg(0, y_mid) * dx[0] * dx[1])] + else: + target_power = current_deg - 2 + FLpiece = tuple([]) + for j in range(0, target_power + 1): + k = target_power - j + FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dy[0] * dy[1], 0)]) + FLpiece += tuple([(0, leg(j, x_mid) * leg(k, y_mid) * dx[0] * dx[1])]) + return FLpiece + + +def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): + DegsOfIteration = determine_f_lambda_portions_2d(deg) + FL = [] + for i in DegsOfIteration: + FL += f_lambda_1_2d_pieces(i, dx, dy, x_mid, y_mid) + return tuple(FL) + + +def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): + FLTilde = tuple([]) + FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) + FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) + for k in range(1, deg - 1): + FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) + + return tuple(FLTilde) + + +def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): + FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) + FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) + result = FL + FLT + + return result diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index d48fa0c52..4776bd3dc 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -261,8 +261,8 @@ def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): assert current_deg > 1, 'invalid for i = 1' ILpiece = tuple([]) - for j in range(0, deg -1): - for k in range(0, deg - 1 - j): + for j in range(0, current_deg -1): + for k in range(0, current_deg - 1 - j): ILpiece += tuple([(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dz[0] * dz[1])]+ [(0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * @@ -281,7 +281,7 @@ def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)] + - [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * d[y0] * dy[1], leg(degree - j - 1, y_mid) * + [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * dy[0] * dy[1], leg(degree - j - 1, y_mid) * leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) for k in range(1, degree - 2): for l in range(1, degree - 2 - k): diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 6ac974573..74c439f1f 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -14,6 +14,7 @@ from FIAT.Sminus import TrimmedSerendipityEdge # noqa: F401 from FIAT.Sminus import TrimmedSerendipityFace # noqa: F401 from FIAT.SminusDiv import TrimmedSerendipityDiv #noqa: F401 +from FIAT.SminusCurl import TrimmedSerendipityCurl #noqa: F401 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor From 42e8899ba147ba738c8bb149ff387ae3d02a3405 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Mon, 21 Sep 2020 15:03:19 -0700 Subject: [PATCH 21/38] Fixed a small Curl issue. --- FIAT/SminusCurl.py | 1 - 1 file changed, 1 deletion(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 66ef97aad..875942541 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -10,7 +10,6 @@ variables = (x, y, z) leg = legendre - def triangular_number(n): return int((n+1)*n/2) From 3738c000ed130d2352631fd60043408eaefd985a Mon Sep 17 00:00:00 2001 From: Justincrum Date: Fri, 25 Sep 2020 16:32:24 -0700 Subject: [PATCH 22/38] Working on testing Sminus Div in 3-d. --- FIAT/SminusDiv.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 4776bd3dc..ea1f09e44 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -39,7 +39,6 @@ def __init__(self, ref_el, degree, mapping): flat_topology = flat_el.get_topology() entity_ids = {} cur = 0 - for top_dim, entities in flat_topology.items(): entity_ids[top_dim] = {} for entity in entities: @@ -180,7 +179,6 @@ def __init__(self, ref_el, degree): raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") verts = flat_el.get_vertices() - dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) x_mid = 2*x-(verts[-1][0] + verts[0][0]) @@ -198,7 +196,6 @@ def __init__(self, ref_el, degree): IL = I_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) else: IL = () - Sminus_list = FL + IL self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -263,20 +260,20 @@ def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): ILpiece = tuple([]) for j in range(0, current_deg -1): for k in range(0, current_deg - 1 - j): - ILpiece += tuple([(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + ILpiece += tuple([(0, 0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dz[0] * dz[1])]+ - [(0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + [(0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dy[0] * dy[1] ,0)] + - [(leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * + [(-leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * dx[1], 0, 0)]) return ILpiece def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): assert degree > 1, 'invalid for i = 1' - IL_tilde = tuple([(leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + - [(0, leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + - [(0, 0, leg(degree - 2, z_mid) * dz[0] * dz[1])]) + IL_tilde = tuple([(-leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + + [(0, -leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + + [(0, 0, -leg(degree - 2, z_mid) * dz[0] * dz[1])]) IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], leg(degree - j - 1, x_mid) * leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * @@ -286,9 +283,9 @@ def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): for k in range(1, degree - 2): for l in range(1, degree - 2 - k): j = degree - 2 - k - l - IL_tilde += tuple([(leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], - leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], - leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) + IL_tilde += tuple([(-leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], + -leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], + -leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) return IL_tilde From d7b68b13b0ada8e34e86c0cdfbe5109c94148b4b Mon Sep 17 00:00:00 2001 From: Justincrum Date: Tue, 6 Oct 2020 13:09:09 -0700 Subject: [PATCH 23/38] Update for making sure basis functions are working properly now at all orders. Still need to test order 5 and 6 to make sure those work. --- FIAT/SminusDiv.py | 58 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 48 insertions(+), 10 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index ea1f09e44..4a5e1d8a6 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -57,7 +57,9 @@ def __init__(self, ref_el, degree, mapping): interior_tilde_ids = 3 for k in range(1, degree - 1): interior_tilde_ids = interior_tilde_ids + 3 - if (degree > 3): + if (degree == 4): + interior_tilde_ids += choose_ijk_total(degree - 2) - (2 * degree - 2) + if (degree > 4): interior_tilde_ids += choose_ijk_total(degree - 2) - (2 * degree - 1) if degree == 1: interior_tilde_ids = 0 @@ -197,6 +199,8 @@ def __init__(self, ref_el, degree): else: IL = () Sminus_list = FL + IL + #print(FL) + print(IL) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -213,29 +217,38 @@ def __init__(self, ref_el, degree): FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) else: FL = () - + Sminus_list = EL + FL Sminus_list = [[-a[1], a[0]] for a in Sminus_list] self.basis = {(0, 0): Array(Sminus_list)} super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") -""" def f_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): +""" def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): FLpiece = tuple([]) - for j in range(0, current_deg + 1): - FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(current_deg - j, y_mid) * c) for c in dz] + - [(0, leg(j, x_mid) * leg(current_deg - j, z_mid) * b, 0) for b in dy] + - [(leg(j, y_mid) * leg(current_deg - j, z_mid) * a, 0, 0) for a in dx]) + #Dx's first. + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(0, leg(j, x_mid) * leg(i - j, z_mid) * dy[k], 0)]) + for k in range(0, 2): + for i in range(0, deg): + for j in range(0, i + 1): + FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) return FLpiece """ -def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): +""" def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): FLpiece = tuple([]) #Dx's first. for k in range(0, 2): for i in range(0, deg): for j in range(0, i + 1): - FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) + FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) for k in range(0, 2): for i in range(0, deg): for j in range(0, i + 1): @@ -243,7 +256,7 @@ def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): for k in range(0, 2): for i in range(0, deg): for j in range(0, i + 1): - FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) + FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) return FLpiece @@ -252,6 +265,16 @@ def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): #for j in range(0, degree): # FL += f_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) FL += f_lambda_2_3d_pieces(degree, dx, dy, dz, x_mid, y_mid, z_mid) + return FL """ + +def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): + + FL = tuple([(-leg(j, y_mid) * leg(k, z_mid) * a, 0, 0) + for a in dx for k in range(0, degree) for j in range(0, degree - k)] + + [(0, leg(j, x_mid) * leg(k, z_mid) * b, 0) + for b in dy for k in range(0, degree) for j in range(0, degree - k)] + + [(0, 0, -leg(j, x_mid) * leg(k, y_mid) * c) + for c in dz for k in range(0, degree) for j in range(0, degree - k)]) return FL @@ -268,6 +291,20 @@ def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): dx[1], 0, 0)]) return ILpiece +""" def i_lambda_2_3d_normal(degree, dx, dy, dz, x_mid, y_mid, z_mid): + + IL = tuple([(-leg(l-2-j, x_mid) * leg(j-k, y_mid) * leg(k, z_mid) * + dx[0] * dx[1], 0, 0) + for l in range(2, degree) for j in range(l-1) for k in range(j+1)] + + [(0,-leg(l-2-j, x_mid) * leg(j-k, y_mid) * leg(k, z_mid) * + dy[0] * dy[1], 0, 0) + for l in range(2, degree) for j in range(l-1) for k in range(j+1)] + + [(0, 0, -leg(l-2-j, x_mid) * leg(j-k, y_mid) * leg(k, z_mid) * + dz[0] * dz[1]) + for l in range(2, degree) for j in range(l-1) for k in range(j+1)]) + + return IL """ + def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): assert degree > 1, 'invalid for i = 1' @@ -293,6 +330,7 @@ def I_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): IL = tuple([]) for j in range(2, degree): IL += I_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) + #IL += i_lambda_2_3d_normal(degree, dx, dy, dz, x_mid, y_mid, z_mid) IL += I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid) return IL From 6e0e0a1c1c8a33e0f01e4bfa9a51608312529e86 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Thu, 8 Oct 2020 13:29:22 -0700 Subject: [PATCH 24/38] More bug testing in 3d. --- FIAT/SminusDiv.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 4a5e1d8a6..6228dea9b 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -58,14 +58,15 @@ def __init__(self, ref_el, degree, mapping): for k in range(1, degree - 1): interior_tilde_ids = interior_tilde_ids + 3 if (degree == 4): - interior_tilde_ids += choose_ijk_total(degree - 2) - (2 * degree - 2) + interior_tilde_ids += choose_ijk_total(degree - 2) - (degree - 1) - (degree - 1) + 1 if (degree > 4): - interior_tilde_ids += choose_ijk_total(degree - 2) - (2 * degree - 1) + #interior_tilde_ids += choose_ijk_total(degree - 2) - (2 * degree - 1) + interior_tilde_ids += choose_ijk_total(degree - 2) - (degree - 1) - (degree - 1) + 1 if degree == 1: interior_tilde_ids = 0 - entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids + print(cur) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -199,8 +200,6 @@ def __init__(self, ref_el, degree): else: IL = () Sminus_list = FL + IL - #print(FL) - print(IL) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -217,7 +216,6 @@ def __init__(self, ref_el, degree): FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) else: FL = () - Sminus_list = EL + FL Sminus_list = [[-a[1], a[0]] for a in Sminus_list] self.basis = {(0, 0): Array(Sminus_list)} @@ -318,10 +316,10 @@ def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * dy[0] * dy[1], leg(degree - j - 1, y_mid) * leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) for k in range(1, degree - 2): - for l in range(1, degree - 2 - k): + for l in range(1, degree - 1 - k): j = degree - 2 - k - l IL_tilde += tuple([(-leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], - -leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], + leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], -leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) return IL_tilde From 1dcaf644300703dac7b6cc3646cd1b2a3ece1d3a Mon Sep 17 00:00:00 2001 From: Justincrum Date: Thu, 8 Oct 2020 15:53:59 -0700 Subject: [PATCH 25/38] Bug testing in 2d. --- FIAT/SminusDiv.py | 36 ++++++++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 6228dea9b..7413aceb8 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -73,10 +73,10 @@ def __init__(self, ref_el, degree, mapping): cur = cur + degree if(degree >= 2): - entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) - - cur += 2*triangular_number(degree - 2) + degree + #entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2))) + cur += 2*triangular_number(degree - 2) formdegree = dim - 1 entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) @@ -208,12 +208,13 @@ def __init__(self, ref_el, degree): if degree < 1: raise Exception("Trimmed serendipity face elements only valid for k >= 1") - flat_el = flatten_reference_cube(ref_el) - verts = flat_el.get_vertices() + #flat_el = flatten_reference_cube(ref_el) + #verts = flat_el.get_vertices() EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: - FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + #FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + FL = F_lambda_1_2d(degree, dx, dy, x_mid, y_mid) else: FL = () Sminus_list = EL + FL @@ -398,17 +399,32 @@ def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): FLTilde = tuple([]) - FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) - FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) + #FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) + #FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) + FLTilde += tuple([(x_mid, 0)]) + FLTilde += tuple([(0, x_mid)]) for k in range(1, deg - 1): - FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) + #FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) + FLTilde += tuple([(x_mid, x_mid)]) + return tuple(FLTilde) +def F_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + FL = [] + for k in range(2, deg): + for j in range(k-1): + FL += [(0, leg(j, x_mid)*leg(k-2-j, y_mid)*dx[0]*dx[1])] + FL += [(leg(k-2-j, x_mid)*leg(j, y_mid)*dy[0]*dy[1], 0)] + + return tuple(FL) + + def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) result = FL + FLT - return result + return FL +# return result From 8076688afdc57b694ac6a464b1b1d6c88c16883b Mon Sep 17 00:00:00 2001 From: Justincrum Date: Thu, 8 Oct 2020 17:02:18 -0700 Subject: [PATCH 26/38] More bug testing. --- FIAT/SminusDiv.py | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 7413aceb8..f0b325f0d 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -66,17 +66,16 @@ def __init__(self, ref_el, degree, mapping): interior_tilde_ids = 0 entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids - print(cur) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) cur = cur + degree if(degree >= 2): - #entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) - entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2))) + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) + #entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2))) - cur += 2*triangular_number(degree - 2) + cur += 2*triangular_number(degree - 2) + degree formdegree = dim - 1 entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) @@ -213,8 +212,8 @@ def __init__(self, ref_el, degree): EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: - #FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) - FL = F_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) + #FL = F_lambda_1_2d(degree, dx, dy, x_mid, y_mid) else: FL = () Sminus_list = EL + FL @@ -399,14 +398,10 @@ def f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid): def f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid): FLTilde = tuple([]) - #FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) - #FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) - FLTilde += tuple([(x_mid, 0)]) - FLTilde += tuple([(0, x_mid)]) + FLTilde += tuple([(leg(deg - 2, y_mid)*dy[0]*dy[1], 0)]) + FLTilde += tuple([(0, leg(deg - 2, x_mid)*dx[0]*dx[1])]) for k in range(1, deg - 1): - #FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) - FLTilde += tuple([(x_mid, x_mid)]) - + FLTilde += tuple([(leg(k, x_mid) * leg(deg - k - 2, y_mid) * dy[0] * dy[1], -leg(k - 1, x_mid) * leg(deg - k - 1, y_mid) * dx[0] * dx[1])]) return tuple(FLTilde) @@ -422,9 +417,10 @@ def F_lambda_1_2d(deg, dx, dy, x_mid, y_mid): def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): - FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) + #FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) + FL = F_lambda_1_2d(deg, dx, dy, x_mid, y_mid) FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) result = FL + FLT - - return FL -# return result + print(FLT) + # return FL + return result From 25b2af0a6ff446a22e6925c1bf99064168d6d9d0 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Wed, 21 Oct 2020 11:46:19 -0700 Subject: [PATCH 27/38] Getting rid of some print statements for debugging. --- FIAT/SminusDiv.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index f0b325f0d..dcfd8a722 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -306,9 +306,9 @@ def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): assert degree > 1, 'invalid for i = 1' - IL_tilde = tuple([(-leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + - [(0, -leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + - [(0, 0, -leg(degree - 2, z_mid) * dz[0] * dz[1])]) + IL_tilde = tuple([(0, 0, leg(degree - 2, z_mid) * dz[0] * dz[1])] + + [(0, leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + + [(leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)]) IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], leg(degree - j - 1, x_mid) * leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * @@ -421,6 +421,5 @@ def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): FL = F_lambda_1_2d(deg, dx, dy, x_mid, y_mid) FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) result = FL + FLT - print(FLT) # return FL return result From 575e380d9e569a49f9abc3c8c97a51c991bedd4d Mon Sep 17 00:00:00 2001 From: Justincrum Date: Thu, 22 Oct 2020 11:24:33 -0700 Subject: [PATCH 28/38] Fixing some basis functions for SminusCurl. --- FIAT/SminusCurl.py | 98 +++++----------------------------------------- 1 file changed, 10 insertions(+), 88 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 875942541..3c04f7df0 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -73,6 +73,8 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids + print("cur = ") + print(cur) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -214,6 +216,7 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL + print(len(FL)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -236,16 +239,6 @@ def __init__(self, ref_el, degree): self.basis = {(0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") -#def e_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): -# EL = () -# for j in range(0, deg): -# for k in range(0, 2): -# for l in range(0, 2): -# EL += tuple([(leg(j, x_mid) * dy[k] * dz[l], 0, 0)] + -# [(0, leg(j, y_mid) * dx[k] * dz[l], 0)] + -# [(0, 0, leg(j, z_mid) * dx[k] * dy[l])]) -# return EL - def e_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): EL = () @@ -324,7 +317,7 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): #FTilde first. FL += tuple([(0, leg(deg - 2, z_mid) * dx[0] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(deg - 2, y_mid) * dx[0] * dy[0] * dy[1])]) - for j in range(1, deg - 2): + for j in range(1, deg - 1): FL += tuple([(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[0] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[0] * dy[0] * dy[1])]) #F next for i in range(2, deg): @@ -332,10 +325,11 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): k = i - 2 - j FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dx[0] * dy[0] * dy[1])]) + print(FL) #IDs associated to face 1. FL += tuple([(0, leg(deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) - for j in range(1, deg - 2): + for j in range(1, deg - 1): FL += tuple([(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[1] * dy[0] * dy[1])]) #F next. for i in range(2, deg): @@ -347,7 +341,7 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): #FTilde first. FL += tuple([(leg(deg - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(deg - 2, x_mid) * dy[0] * dx[0] * dx[1])]) - for j in range(1, deg - 2): + for j in range(1, deg - 1): FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 1, z_mid) * dy[0] * dx[0] * dx[1])]) #F next. for i in range(2, deg): @@ -358,7 +352,7 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): #IDs associated to face 3. FL += tuple([(leg(deg - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(deg - 2, x_mid) * dy[1] * dx[0] * dx[1])]) - for j in range(1, deg - 2): + for j in range(1, deg - 1): FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 1, z_mid) * dy[1] * dx[0] * dx[1])]) #F next. for i in range(2, deg): @@ -369,7 +363,7 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): #IDs associated to face 4. FL += tuple([(leg(deg - 2, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(deg - 2, x_mid) * dz[0] * dx[0] * dx[1], 0)]) - for j in range(1, deg - 2): + for j in range(1, deg - 1): FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[0] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg - j - 1, y_mid) * dz[0] * dx[0] * dx[1], 0)]) for i in range(2, deg): for j in range(0, i - 1): @@ -379,7 +373,7 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): #IDs associated to face 5. FL += tuple([(leg(deg - 2, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(deg - 2, x_mid) * dz[1] * dx[0] * dx[1], 0)]) - for j in range(1, deg - 2): + for j in range(1, deg - 1): FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[1] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg -j - 1, y_mid) * dz[1] * dx[0] * dx[1], 0)]) for i in range(2, deg): for j in range(0, i - 1): @@ -433,78 +427,6 @@ def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): return I -""" -def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): - FLpiece = tuple([]) - #Dx's first. - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(0, leg(j, x_mid) * leg(i - j, z_mid) * dy[k], 0)]) - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) - return FLpiece - - -def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): - FL = tuple([]) - #for j in range(0, degree): - # FL += f_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) - FL += f_lambda_2_3d_pieces(degree, dx, dy, dz, x_mid, y_mid, z_mid) - return FL - - -def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): - assert current_deg > 1, invalid for i = 1 - ILpiece = tuple([]) - for j in range(0, current_deg -1): - for k in range(0, current_deg - 1 - j): - ILpiece += tuple([(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * - dz[0] * dz[1])]+ - [(0, leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * - dy[0] * dy[1] ,0)] + - [(leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * - dx[1], 0, 0)]) - return ILpiece - - - - -def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): - assert degree > 1, invalid for i = 1 - IL_tilde = tuple([(leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)] + - [(0, leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + - [(0, 0, leg(degree - 2, z_mid) * dz[0] * dz[1])]) - IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], leg(degree - j - 1, x_mid) * - leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + - [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * - leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)] + - [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * dy[0] * dy[1], leg(degree - j - 1, y_mid) * - leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) - for k in range(1, degree - 2): - for l in range(1, degree - 2 - k): - j = degree - 2 - k - l - IL_tilde += tuple([(leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], - leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], - leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) - return IL_tilde - - -def I_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): - IL = tuple([]) - for j in range(2, degree): - IL += I_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) - IL += I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid) - return IL - -""" - #Everything for 2-d should work already. def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): EL = tuple( From c7e5ca5933980b35ae87034ad5fdc4badc76ff94 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Thu, 22 Oct 2020 11:33:39 -0700 Subject: [PATCH 29/38] Getting rid of some prints. --- FIAT/SminusCurl.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 3c04f7df0..b8d64df42 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -216,7 +216,7 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL - print(len(FL)) + print(len(Sminus_list)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -325,7 +325,6 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): k = i - 2 - j FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dx[0] * dy[0] * dy[1])]) - print(FL) #IDs associated to face 1. FL += tuple([(0, leg(deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) From fbeca61213dcbeb9dd5c2795b12b90002ade3fc0 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Mon, 26 Oct 2020 09:31:09 -0700 Subject: [PATCH 30/38] Fixing SminusCurl.py number of basis functions for higher degrees. --- FIAT/SminusCurl.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index b8d64df42..ac679d5ef 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -57,7 +57,6 @@ def __init__(self, ref_el, degree, mapping): for j in sorted(flat_topology[2]): entity_ids[2][j] = list(range(cur, cur + face_ids)) cur += face_ids - interior_ids = 0 for k in range(4, degree): interior_ids = interior_ids + 3 * choose_ijk_total(k - 4) @@ -73,8 +72,8 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids - print("cur = ") - print(cur) + print("InteriorIDs =", interior_ids +interior_tilde_ids) + print("Cur =", cur) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -216,7 +215,7 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL - print(len(Sminus_list)) + print(len(IL)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -393,9 +392,10 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): I = () - for j in range(0, deg - 3): - for k in range(0, deg - 3 - j): - l = deg - 4 - j - k + D = deg - 1 + for j in range(0, D - 3): + for k in range(0, D - 3 - j): + l = D - 4 - j - k I += tuple([(leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + [(0, leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + [(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) From c1f477d26d1568abeb1efc39a0523d3f9996fcbb Mon Sep 17 00:00:00 2001 From: Justincrum Date: Fri, 30 Oct 2020 14:30:44 -0700 Subject: [PATCH 31/38] Fixing some basis functions in SminusCurl from getting double counted at degree 4. --- FIAT/SminusCurl.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index ac679d5ef..fa09d1b9a 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -62,7 +62,7 @@ def __init__(self, ref_el, degree, mapping): interior_ids = interior_ids + 3 * choose_ijk_total(k - 4) if (degree > 3): if (degree == 4): - interior_tilde_ids = 6 + interior_tilde_ids = 3 elif (degree == 5): interior_tilde_ids = 8 else: @@ -215,7 +215,8 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL - print(len(IL)) + print("Interior IDs = ", len(IL)) + print("IDs =", len(Sminus_list)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -404,12 +405,17 @@ def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): ITilde = () - ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + - [(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + - [(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + - [(0, leg(deg - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + - [(0, 0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dy[0] * dy[1])] + - [(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + if(deg==4): + ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + if(deg > 4): + ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(0, leg(deg - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(0, 0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dy[0] * dy[1])] + + [(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) for j in range(1, deg - 3): ITilde += tuple([(leg(j, x_mid) * leg(deg - j - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], -leg(j - 1, x_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + [(leg(j, x_mid) * leg(deg - j - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 3, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) @@ -422,7 +428,9 @@ def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): I = () for i in range(4, deg): I += I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid) + print("Len of I normal", len(I)) I += I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid) + print("Len of total I", len(I)) return I From e82c3fee36c2d21aabd6ca69ea0871385a4c7225 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Fri, 30 Oct 2020 14:32:15 -0700 Subject: [PATCH 32/38] Got rid of debug print statements. --- FIAT/SminusCurl.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index fa09d1b9a..6c2eb56a9 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -72,8 +72,6 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids - print("InteriorIDs =", interior_ids +interior_tilde_ids) - print("Cur =", cur) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -215,8 +213,6 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL - print("Interior IDs = ", len(IL)) - print("IDs =", len(Sminus_list)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -428,9 +424,7 @@ def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): I = () for i in range(4, deg): I += I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid) - print("Len of I normal", len(I)) I += I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid) - print("Len of total I", len(I)) return I From 7b4950ae2889c5909525e50581aa4955c30f1895 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Mon, 2 Nov 2020 16:16:40 -0700 Subject: [PATCH 33/38] Fixing order 6+ --- FIAT/SminusCurl.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 6c2eb56a9..1c7626849 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -50,6 +50,7 @@ def __init__(self, ref_el, degree, mapping): for l in sorted(flat_topology[1]): entity_ids[1][l] = list(range(cur, cur + degree)) cur = cur + degree + print("Total number of edge IDs available", cur) if (degree > 1): face_ids = degree for k in range(2, degree): @@ -72,6 +73,10 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids + print("Total spots for basis functions", cur) + print("The spots for interior basis functions", interior_ids) + print("The spots for interior tilde basis functions", interior_tilde_ids) + print("The total spots for interior functions", interior_ids + interior_tilde_ids) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -213,6 +218,8 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL + print("The total number of edge functions", len(EL)) + print("The number of interior basis functions", len(IL)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -389,10 +396,9 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): I = () - D = deg - 1 - for j in range(0, D - 3): - for k in range(0, D - 3 - j): - l = D - 4 - j - k + for j in range(0, deg - 3): + for k in range(0, deg - 3 - j): + l = deg - 4 - j - k I += tuple([(leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + [(0, leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + [(0, 0, leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) @@ -402,9 +408,9 @@ def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): ITilde = () if(deg==4): - ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + ITilde += tuple([(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + [(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + - [(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + [(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) if(deg > 4): ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + [(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + @@ -416,15 +422,16 @@ def I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): ITilde += tuple([(leg(j, x_mid) * leg(deg - j - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], -leg(j - 1, x_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + [(leg(j, x_mid) * leg(deg - j - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 3, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) if deg > 5: - ITilde += [(0, leg(j, y_mid) * leg(deg - j - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1])] + ITilde += tuple([(0, leg(j, y_mid) * leg(deg - j - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1])]) return ITilde def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): I = () for i in range(4, deg): - I += I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid) + I += I_lambda_1_3d_pieces(i, dx, dy, dz, x_mid, y_mid, z_mid) I += I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid) + print("Total interior functions", len(I)) return I From dd67cde8438a34e45eb41cef61f39f28151573c7 Mon Sep 17 00:00:00 2001 From: Justincrum Date: Wed, 4 Nov 2020 13:46:33 -0700 Subject: [PATCH 34/38] Trying to fix order 4. --- FIAT/SminusCurl.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 1c7626849..880574fda 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -58,6 +58,7 @@ def __init__(self, ref_el, degree, mapping): for j in sorted(flat_topology[2]): entity_ids[2][j] = list(range(cur, cur + face_ids)) cur += face_ids + print("Total number of edge + face IDs available", cur) interior_ids = 0 for k in range(4, degree): interior_ids = interior_ids + 3 * choose_ijk_total(k - 4) @@ -408,9 +409,9 @@ def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): ITilde = () if(deg==4): - ITilde += tuple([(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + - [(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + - [(0, 0, leg(deg - 4, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) + ITilde += tuple([(dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(0, dx[0] * dx[1] * dz[0] * dz[1], 0)] + + [(0, 0, dx[0] * dx[1] * dy[0] * dy[1])]) if(deg > 4): ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + [(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + From fead6e5d876209b98afa8a57793f7bc095e0900e Mon Sep 17 00:00:00 2001 From: Justincrum Date: Fri, 13 Nov 2020 10:46:08 -0700 Subject: [PATCH 35/38] Adding in the missing basis functions for SminusCurl at order 6. --- FIAT/SminusCurl.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 880574fda..06fb93f10 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -50,7 +50,6 @@ def __init__(self, ref_el, degree, mapping): for l in sorted(flat_topology[1]): entity_ids[1][l] = list(range(cur, cur + degree)) cur = cur + degree - print("Total number of edge IDs available", cur) if (degree > 1): face_ids = degree for k in range(2, degree): @@ -58,7 +57,6 @@ def __init__(self, ref_el, degree, mapping): for j in sorted(flat_topology[2]): entity_ids[2][j] = list(range(cur, cur + face_ids)) cur += face_ids - print("Total number of edge + face IDs available", cur) interior_ids = 0 for k in range(4, degree): interior_ids = interior_ids + 3 * choose_ijk_total(k - 4) @@ -67,17 +65,15 @@ def __init__(self, ref_el, degree, mapping): interior_tilde_ids = 3 elif (degree == 5): interior_tilde_ids = 8 + elif (degree == 6): + interior_tilde_ids = 6 + 3 * (degree - 3) else: - interior_tilde_ids = 6 + 3 * (degree - 4) + interior_tilde_ids = 6 + 3 * (degree - 4) else: interior_tilde_ids = 0 entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids - print("Total spots for basis functions", cur) - print("The spots for interior basis functions", interior_ids) - print("The spots for interior tilde basis functions", interior_tilde_ids) - print("The total spots for interior functions", interior_ids + interior_tilde_ids) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -219,8 +215,6 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL - print("The total number of edge functions", len(EL)) - print("The number of interior basis functions", len(IL)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") @@ -424,6 +418,10 @@ def I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): [(leg(j, x_mid) * leg(deg - j - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 3, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) if deg > 5: ITilde += tuple([(0, leg(j, y_mid) * leg(deg - j - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 3, y_mid) * dx[0] * dx[1] * dz[0] * dz[1])]) + if (deg == 6): + ITilde += tuple([(leg(1, y_mid) * leg(1, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)]) + ITilde += tuple([(0, leg(1, x_mid) * leg(1, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)]) + ITilde += tuple([(0, 0, leg(1, x_mid) * leg(1, y_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) return ITilde @@ -432,7 +430,6 @@ def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): for i in range(4, deg): I += I_lambda_1_3d_pieces(i, dx, dy, dz, x_mid, y_mid, z_mid) I += I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid) - print("Total interior functions", len(I)) return I From 4f457bcbb588165081f1439d242c90bf78e62e70 Mon Sep 17 00:00:00 2001 From: Justin Crum Date: Wed, 16 Dec 2020 09:51:11 -0700 Subject: [PATCH 36/38] Some change --- FIAT/SminusCurl.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 06fb93f10..2d37cf708 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -74,6 +74,7 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids + #print("cur =", cur) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -215,6 +216,7 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL + #print("Number of basis functions", len(Sminus_list)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") From 26feeed26906464c0a04b572c89795a71c7fdbcb Mon Sep 17 00:00:00 2001 From: Justin Crum Date: Tue, 30 Mar 2021 14:18:33 -0700 Subject: [PATCH 37/38] Getting rid of print statements. --- FIAT/SminusCurl.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 2d37cf708..cbf46dc3b 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -74,7 +74,6 @@ def __init__(self, ref_el, degree, mapping): entity_ids[3][0] = list(range(cur, cur + interior_ids + interior_tilde_ids)) cur = cur + interior_ids + interior_tilde_ids - #print("cur =", cur) else: for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) @@ -190,6 +189,7 @@ def __init__(self, ref_el, degree): raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") verts = flat_el.get_vertices() + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) @@ -216,7 +216,6 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = EL + FL + IL - #print("Number of basis functions", len(Sminus_list)) self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") From d3cc392060823e70fa29a52eba804e8b3c45b2d9 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 23 Jun 2022 23:17:57 -0500 Subject: [PATCH 38/38] flake8 cleanup --- FIAT/Sminus.py | 103 +++++++++++++++++++------------------ FIAT/SminusCurl.py | 124 +++++++++++++++------------------------------ FIAT/SminusDiv.py | 105 ++++++++------------------------------ FIAT/__init__.py | 4 +- 4 files changed, 113 insertions(+), 223 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 728886914..9373e0a02 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -1,3 +1,4 @@ + # Copyright (C) 2019 Cyrus Cheng (Imperial College London) # # This file is part of FIAT. @@ -32,6 +33,7 @@ def triangular_number(n): return int((n+1)*n/2) + def choose_ijk_total(degree): top = 1 for i in range(1, 2 + degree + 1): @@ -39,7 +41,7 @@ def choose_ijk_total(degree): bottom = 1 for i in range(1, degree + 1): bottom = i * bottom - return int(top /(2 * bottom)) + return int(top / (2 * bottom)) class TrimmedSerendipity(FiniteElement): @@ -69,17 +71,17 @@ def __init__(self, ref_el, degree, mapping): if(degree >= 2): entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) - + cur += 2*triangular_number(degree - 2) + degree else: - #3-d case. + # 3-d case. entity_ids[3] = {} for j in sorted(flat_topology[1]): entity_ids[1][j] = list(range(cur, cur + degree)) cur = cur + degree - - if (degree >= 2 ): + + if (degree >= 2): if (degree < 4): for j in sorted(flat_topology[2]): entity_ids[2][j] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) @@ -88,8 +90,6 @@ def __init__(self, ref_el, degree, mapping): for j in sorted(flat_topology[2]): entity_ids[2][j] = list(range(cur, cur + 3 * degree - 4)) cur = cur + 3*degree - 4 - - if (degree >= 4): if (degree == 4): entity_ids[3][0] = list(range(cur, cur + 6)) @@ -172,7 +172,6 @@ def tabulate(self, order, points, entity=None): phivals[alpha] = T return phivals - def entity_dofs(self): """Return the map of topological entities to degrees of freedom for the finite element.""" @@ -284,40 +283,40 @@ def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): def e_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): EL = tuple([]) - ##assignment to edge x=y=0 - for i in range (0, max_deg): + # assignment to edge x=y=0 + for i in range(0, max_deg): EL += tuple([(0, 0, leg(i, z_mid) * dx[0] * dy[0])]) - ##assignment to edge x=0, y=1 + # assignment to edge x=0, y=1 for i in range(0, max_deg): EL += tuple([(0, 0, leg(i, z_mid) * dy[1] * dx[0])]) - ##assignment to edge x = 1, y = 0 + # assignment to edge x = 1, y = 0 for i in range(0, max_deg): EL += tuple([(0, 0, leg(i, z_mid) * dx[1] * dy[0])]) - ##assignment to edge x = 1, y = 1 + # assignment to edge x = 1, y = 1 for i in range(0, max_deg): EL += tuple([(0, 0, leg(i, z_mid) * dx[1] * dy[1])]) - ##assignment to edge x = 0, z = 0 + # assignment to edge x = 0, z = 0 for i in range(0, max_deg): EL += tuple([(0, leg(i, y_mid) * dx[0] * dz[0], 0)]) - ##assignment to edge x = 0, z = 1 + # assignment to edge x = 0, z = 1 for i in range(0, max_deg): EL += tuple([(0, leg(i, y_mid) * dx[0] * dz[1], 0)]) - ##assignment to edge x = 1, z = 0 + # assignment to edge x = 1, z = 0 for i in range(0, max_deg): EL += tuple([(0, leg(i, y_mid) * dx[1] * dz[0], 0)]) - ##assignment to edge x = 1, z = 1 + # assignment to edge x = 1, z = 1 for i in range(0, max_deg): EL += tuple([(0, leg(i, y_mid) * dx[1] * dz[1], 0)]) - ##assignment to edge y = 0, z = 0 + # assignment to edge y = 0, z = 0 for i in range(0, max_deg): EL += tuple([(leg(i, x_mid) * dy[0] * dz[0], 0, 0)]) - ##assignment to edge y = 0, z = 1 + # assignment to edge y = 0, z = 1 for i in range(0, max_deg): EL += tuple([(leg(i, x_mid) * dy[0] * dz[1], 0, 0)]) - ##assignment to edge y = 1, z = 0 + # assignment to edge y = 1, z = 0 for i in range(0, max_deg): EL += tuple([(leg(i, x_mid) * dy[1] * dz[0], 0, 0)]) - ##assignment to edge y = 1, z = 1 + # assignment to edge y = 1, z = 1 for i in range(0, max_deg): EL += tuple([(leg(i, x_mid) * dy[1] * dz[1], 0, 0)]) return EL @@ -325,72 +324,72 @@ def e_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): def f_lambda_1_3d_trimmed(max_deg, dx, dy, dz, x_mid, y_mid, z_mid): FL = tuple([]) - ##Assignment to face x = 0, Ftilde + # Assignment to face x = 0, Ftilde FL += tuple([(0, leg(max_deg - 2, z_mid) * dx[0] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[0] * dy[0] * dy[1])]) for j in range(1, max_deg - 1): - FL += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[0] * dz[0] * dz[1], - -leg(j - 1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[0] * dy[0] * dy[1])]) - ##Assignment to face x = 0, F + FL += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[0] * dz[0] * dz[1], + -leg(j - 1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[0] * dy[0] * dy[1])]) + # Assignment to face x = 0, F for j in range(1, max_deg - 1): k = max_deg - j - 2 FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[0] * dy[0] * dy[1])]) - ##Assignment to face x = 1, Ftilde + # Assignment to face x = 1, Ftilde FL += tuple([(0, leg(max_deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(max_deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) for j in range(1, max_deg - 1): - FL += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[1] * dz[0] * dz[1], - -leg(j - 1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[1] * dy[0] * dy[1])]) - ##Assignment to face x = 1, F + FL += tuple([(0, leg(j, y_mid) * leg(max_deg - j - 2, z_mid) * dx[1] * dz[0] * dz[1], + -leg(j - 1, y_mid) * leg(max_deg - j - 1, z_mid) * dx[1] * dy[0] * dy[1])]) + # Assignment to face x = 1, F for j in range(1, max_deg - 1): k = max_deg - j - 2 FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[1] * dy[0] * dy[1])]) - ##Assignment to face y = 0, Ftilde + # Assignment to face y = 0, Ftilde FL += tuple([(leg(max_deg - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[0] * dx[0] * dx[1])]) for j in range(1, max_deg - 1): - FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, - -leg(j - 1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[0] * dx[0] * dx[1])]) - ##Assignment to face y = 0, F + FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, + -leg(j - 1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[0] * dx[0] * dx[1])]) + # Assignment to face y = 0, F for j in range(1, max_deg - 1): k = max_deg - j - 2 FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[0] * dx[0] * dx[1])]) - ##Assignment to face y = 1, Ftilde + # Assignment to face y = 1, Ftilde FL += tuple([(leg(max_deg - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(max_deg - 2, x_mid) * dy[1] * dx[0] * dx[1])]) for j in range(1, max_deg - 1): - FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, - -leg(j - 1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[1] * dx[0] * dx[1])]) - ##Assignment to face y = 1, F + FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, + -leg(j - 1, x_mid) * leg(max_deg - j - 1, z_mid) * dy[1] * dx[0] * dx[1])]) + # Assignment to face y = 1, F for j in range(1, max_deg - 1): k = max_deg - j - 2 FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[1] * dx[0] * dx[1])]) - ##Assignment to face z = 0, Ftilde + # Assignment to face z = 0, Ftilde FL += tuple([(leg(max_deg - 2, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(max_deg - 2, x_mid) * dz[0] * dx[0] * dx[1], 0)]) for j in range(1, max_deg - 1): FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, y_mid) * dz[0] * dy[0] * dy[1], - -leg(j - 1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[0] * dx[0] * dx[1], 0)]) - ##Assignment to face z = 0, F + -leg(j - 1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[0] * dx[0] * dx[1], 0)]) + # Assignment to face z = 0, F for j in range(1, max_deg - 1): k = max_deg - j - 2 FL += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)]) - ##Assignment to face z = 1, Ftilde + # Assignment to face z = 1, Ftilde FL += tuple([(leg(max_deg - 2, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(max_deg - 2, x_mid) * dz[1] * dx[0] * dx[1], 0)]) for j in range(1, max_deg - 1): FL += tuple([(leg(j, x_mid) * leg(max_deg - j - 2, y_mid) * dz[1] * dy[0] * dy[1], - -leg(j - 1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[1] * dx[0] * dx[1], 0)]) - ##Assignment to face z = 1, F + -leg(j - 1, x_mid) * leg(max_deg - j - 1, y_mid) * dz[1] * dx[0] * dx[1], 0)]) + # Assignment to face z = 1, F for j in range(1, max_deg - 1): k = max_deg - j - 2 FL += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) @@ -418,11 +417,11 @@ def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): DegsOfIteration = determine_I_lambda_1_portions_3d(deg) IL = tuple([]) for Degs in DegsOfIteration: - IL += tuple([(leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * + IL += tuple([(leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)]) - IL += tuple([(0, leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * + IL += tuple([(0, leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)]) - IL += tuple([(0, 0, leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * + IL += tuple([(0, 0, leg(Degs[0], x_mid) * leg(Degs[1], y_mid) * leg(Degs[2], z_mid) * dy[0] * dy[1] * dy[0] * dy[1])]) return IL @@ -443,10 +442,10 @@ def I_lambda_1_tilde_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): if(deg > 5): ILtilde += tuple([(0, leg(j, y_mid) * leg(deg - j - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 3, z_mid) * dx[0] * dx[1] * dy[0] * dy[1])]) - return ILtilde + return ILtilde -#This is always 1-forms regardless of 2 or 3 dimensions. +# This is always 1-forms regardless of 2 or 3 dimensions. class TrimmedSerendipityEdge(TrimmedSerendipity): def __init__(self, ref_el, degree): if degree < 1: @@ -489,11 +488,11 @@ def __init__(self, ref_el, degree): y_mid, z_mid) else: IL = () - + Sminus_list = EL + FL if dim == 3: Sminus_list = Sminus_list + IL - + if dim == 2: self.basis = {(0, 0): Array(Sminus_list)} else: @@ -517,7 +516,7 @@ def __init__(self, ref_el, degree): dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) x_mid = 2*x-(verts[-1][0] + verts[0][0]) y_mid = 2*y-(verts[-1][1] + verts[0][1]) - + EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) @@ -527,4 +526,4 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL Sminus_list = [[-a[1], a[0]] for a in Sminus_list] self.basis = {(0, 0): Array(Sminus_list)} - super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") \ No newline at end of file + super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index cbf46dc3b..a03f1529c 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -1,4 +1,3 @@ -#Working on 3d trimmed serendipity 2 forms now. from sympy import symbols, legendre, Array, diff import numpy as np from FIAT.finite_element import FiniteElement @@ -10,9 +9,11 @@ variables = (x, y, z) leg = legendre + def triangular_number(n): return int((n+1)*n/2) + def choose_ijk_total(degree): top = 1 for i in range(1, 2 + degree + 1): @@ -20,7 +21,7 @@ def choose_ijk_total(degree): bottom = 1 for i in range(1, degree + 1): bottom = i * bottom - return int(top /(2 * bottom)) + return int(top / (2 * bottom)) class TrimmedSerendipity(FiniteElement): @@ -32,7 +33,7 @@ def __init__(self, ref_el, degree, mapping): dim = flat_el.get_spatial_dimension() self.fdim = dim if dim != 3: - if dim !=2: + if dim != 2: raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3") flat_topology = flat_el.get_topology() @@ -44,7 +45,7 @@ def __init__(self, ref_el, degree, mapping): for entity in entities: entity_ids[top_dim][entity] = [] - #3-d case. + # 3-d case. if dim == 3: entity_ids[3] = {} for l in sorted(flat_topology[1]): @@ -68,7 +69,7 @@ def __init__(self, ref_el, degree, mapping): elif (degree == 6): interior_tilde_ids = 6 + 3 * (degree - 3) else: - interior_tilde_ids = 6 + 3 * (degree - 4) + interior_tilde_ids = 6 + 3 * (degree - 4) else: interior_tilde_ids = 0 @@ -81,7 +82,7 @@ def __init__(self, ref_el, degree, mapping): if(degree >= 2): entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) - + cur += 2*triangular_number(degree - 2) + degree formdegree = 1 @@ -177,7 +178,6 @@ def space_dimension(self): class TrimmedSerendipityCurl(TrimmedSerendipity): - #3d Curl element needs to be fixed entirely. This should correspond to 1-forms in 3d. def __init__(self, ref_el, degree): if degree < 1: raise Exception("Trimmed serendipity face elements only valid for k >= 1") @@ -185,11 +185,10 @@ def __init__(self, ref_el, degree): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - if dim !=3: + if dim != 3: raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") verts = flat_el.get_vertices() - dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) @@ -204,7 +203,7 @@ def __init__(self, ref_el, degree): if dim == 3: if degree < 1: - raise Exception ("Trimmed Serendipity fce elements only valid for k >= 1") + raise Exception("Trimmed Serendipity fce elements only valid for k >= 1") EL = e_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) if (degree > 1): FL = f_lambda_1_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid) @@ -218,16 +217,12 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL + IL self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") - + else: - ##Put all 2 dimensional stuff here. - ##This part might be good, need to test it. + # Put all 2 dimensional stuff here. if degree < 1: raise Exception("Trimmed serendipity face elements only valid for k >= 1") - #flat_el = flatten_reference_cube(ref_el) - #verts = flat_el.get_vertices() - EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) @@ -241,124 +236,94 @@ def __init__(self, ref_el, degree): def e_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): EL = () - #IDs associated to edge 0. + # IDs associated to edge 0. for j in range(0, deg): EL += tuple([(0, 0, leg(j, z_mid) * dx[0] * dy[0])]) - #IDs associated to edge 1. + # IDs associated to edge 1. for j in range(0, deg): EL += tuple([(0, 0, leg(j, z_mid) * dx[0] * dy[1])]) - #IDs associated to edge 2. + # IDs associated to edge 2. for j in range(0, deg): EL += tuple([(0, 0, leg(j, z_mid) * dx[1] * dy[0])]) - #IDs associated to edge 3. + # IDs associated to edge 3. for j in range(0, deg): EL += tuple([(0, 0, leg(j, z_mid) * dx[1] * dy[1])]) - #IDs associated edge 4. + # IDs associated edge 4. for j in range(0, deg): EL += tuple([(0, leg(j, y_mid) * dx[0] * dz[0], 0)]) - #IDs associated to edge 5. + # IDs associated to edge 5. for j in range(0, deg): EL += tuple([(0, leg(j, y_mid) * dx[0] * dz[1], 0)]) - #IDs associated to edge 6. + # IDs associated to edge 6. for j in range(0, deg): EL += tuple([(0, leg(j, y_mid) * dx[1] * dz[0], 0)]) - #IDs associated to edge 7. + # IDs associated to edge 7. for j in range(0, deg): EL += tuple([(0, leg(j, y_mid) * dx[1] * dz[1], 0)]) - #IDs associated to edge 8. + # IDs associated to edge 8. for j in range(0, deg): EL += tuple([(leg(j, x_mid) * dz[0] * dy[0], 0, 0)]) - #IDs associated to edge 9. + # IDs associated to edge 9. for j in range(0, deg): EL += tuple([(leg(j, x_mid) * dz[1] * dy[0], 0, 0)]) - #IDs associated to edge 10. + # IDs associated to edge 10. for j in range(0, deg): EL += tuple([(leg(j, x_mid) * dz[0] * dy[1], 0, 0)]) - #IDs associated to edge 11. + # IDs associated to edge 11. for j in range(0, deg): EL += tuple([(leg(j, x_mid) * dz[1] * dy[1], 0, 0)]) return EL -#def f_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): -# FLpiece = () -# for l in range(0, 2): -# for j in range(0, deg - 2): -# k = deg - 2 - j -# FLpiece += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[l] * dy[0] * dy[1], 0, 0)] + -# [(leg(j, x_mid) * leg(k, z_mid) * dy[l] * dz[0] * dz[1], 0, 0)] + -# [(0, leg(j, y_mid) * leg(k, x_mid) * dz[l] * dx[0] * dx[1], 0)] + -# [(0, leg(j, y_mid) * leg(k, z_mid) * dx[l] * dz[0] * dz[1], 0)] + -# [(0, 0, leg(j, z_mid) * leg(k, y_mid) * dy[l] * dx[0] * dx[1])] + -# [(0, 0, leg(j, z_mid) * leg(k, x_mid) * dx[l] * dy[0] * dy[1])]) -# return FLpiece - - -#def f_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): -# FTilde = () -# for l in range(0, 2): -# FTilde += tuple([(leg(deg - 2, y_mid) * dz[l] * dy[0] * dy[1], 0, 0)] + -# [(leg(deg - 2, z_mid) * dy[l] * dz[0] * dz[1], 0, 0)] + -# [(0, leg(deg - 2, x_mid) * dz[l] * dx[0] * dx[1], 0)] + -# [(0, leg(deg - 2, z_mid) * dx[l] * dz[0] * dz[1], 0)] + -# [(0, 0, leg(deg - 2, x_mid) * dy[l] * dx[0] * dx[1])] + -# [(0, 0, leg(deg - 2, y_mid) * dx[l] * dy[0] * dy[1])]) -# for j in range(1, deg - 1): -# for l in range(0, 2): -# FTilde += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[l] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg - j - 1, y_mid) * dz[l] * dx[0] * dx[1], 0)] + -# [(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[l] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[l] * dy[0] * dy[1])]) -# return FTilde - - def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): FL = () - #IDs associated to face 0. - #FTilde first. + # IDs associated to face 0. + # FTilde first. FL += tuple([(0, leg(deg - 2, z_mid) * dx[0] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(deg - 2, y_mid) * dx[0] * dy[0] * dy[1])]) for j in range(1, deg - 1): FL += tuple([(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[0] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[0] * dy[0] * dy[1])]) - #F next + # F next for i in range(2, deg): for j in range(0, i - 1): k = i - 2 - j FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[0] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dx[0] * dy[0] * dy[1])]) - #IDs associated to face 1. + # IDs associated to face 1. FL += tuple([(0, leg(deg - 2, z_mid) * dx[1] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(deg - 2, y_mid) * dx[1] * dy[0] * dy[1])]) for j in range(1, deg - 1): FL += tuple([(0, leg(j, y_mid) * leg(deg - j - 2, z_mid) * dx[1] * dz[0] * dz[1], -leg(j - 1, y_mid) * leg(deg - j - 1, z_mid) * dx[1] * dy[0] * dy[1])]) - #F next. + # F next. for i in range(2, deg): for j in range(0, i - 1): k = i - 2 - j FL += tuple([(0, leg(j, y_mid) * leg(k, z_mid) * dx[1] * dz[0] * dz[1], 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, y_mid) * dx[1] * dy[0] * dy[1])]) - #IDs associated to face 2. - #FTilde first. + # IDs associated to face 2. + # FTilde first. FL += tuple([(leg(deg - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(deg - 2, x_mid) * dy[0] * dx[0] * dx[1])]) for j in range(1, deg - 1): FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, z_mid) * dy[0] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 1, z_mid) * dy[0] * dx[0] * dx[1])]) - #F next. + # F next. for i in range(2, deg): for j in range(0, i - 1): k = i - 2 - j FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[0] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dy[0] * dx[0] * dx[1])]) - #IDs associated to face 3. + # IDs associated to face 3. FL += tuple([(leg(deg - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(deg - 2, x_mid) * dy[1] * dx[0] * dx[1])]) for j in range(1, deg - 1): FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, z_mid) * dy[1] * dz[0] * dz[1], 0, -leg(j - 1, x_mid) * leg(deg - j - 1, z_mid) * dy[1] * dx[0] * dx[1])]) - #F next. + # F next. for i in range(2, deg): for j in range(0, i - 1): k = i - 2 - j FL += tuple([(leg(j, x_mid) * leg(k, z_mid) * dy[1] * dz[0] * dz[1], 0, 0)]) FL += tuple([(0, 0, leg(j, z_mid) * leg(k, x_mid) * dy[1] * dx[0] * dx[1])]) - #IDs associated to face 4. + # IDs associated to face 4. FL += tuple([(leg(deg - 2, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(deg - 2, x_mid) * dz[0] * dx[0] * dx[1], 0)]) for j in range(1, deg - 1): @@ -368,11 +333,11 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): k = i - 2 - j FL += tuple([(leg(j, x_mid) * leg(k, y_mid) * dz[0] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(j, y_mid) * leg(k, x_mid) * dz[0] * dx[0] * dx[1], 0)]) - #IDs associated to face 5. + # IDs associated to face 5. FL += tuple([(leg(deg - 2, y_mid) * dz[1] * dy[0] * dy[1], 0, 0)]) FL += tuple([(0, leg(deg - 2, x_mid) * dz[1] * dx[0] * dx[1], 0)]) for j in range(1, deg - 1): - FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[1] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg -j - 1, y_mid) * dz[1] * dx[0] * dx[1], 0)]) + FL += tuple([(leg(j, x_mid) * leg(deg - j - 2, y_mid) * dz[1] * dy[0] * dy[1], -leg(j - 1, x_mid) * leg(deg - j - 1, y_mid) * dz[1] * dx[0] * dx[1], 0)]) for i in range(2, deg): for j in range(0, i - 1): k = i - 2 - j @@ -381,15 +346,6 @@ def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): return FL - -#def f_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): -# F = () -# for i in range(2, deg): -# F += f_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid) -# F += f_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid) -# return F - - def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): I = () for j in range(0, deg - 3): @@ -403,13 +359,13 @@ def I_lambda_1_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_tilde_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): ITilde = () - if(deg==4): + if(deg == 4): ITilde += tuple([(dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + [(0, dx[0] * dx[1] * dz[0] * dz[1], 0)] + - [(0, 0, dx[0] * dx[1] * dy[0] * dy[1])]) + [(0, 0, dx[0] * dx[1] * dy[0] * dy[1])]) if(deg > 4): ITilde += tuple([(leg(deg - 4, y_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + - [(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + + [(leg(deg - 4, z_mid) * dy[0] * dy[1] * dz[0] * dz[1], 0, 0)] + [(0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + [(0, leg(deg - 4, z_mid) * dx[0] * dx[1] * dz[0] * dz[1], 0)] + [(0, 0, leg(deg - 4, x_mid) * dx[0] * dx[1] * dy[0] * dy[1])] + @@ -434,7 +390,7 @@ def I_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): return I -#Everything for 2-d should work already. +# Everything for 2-d should work already. def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): EL = tuple( [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)] + diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index dcfd8a722..83461cb1e 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -1,4 +1,3 @@ -#Working on 3d trimmed serendipity 2 forms now. from sympy import symbols, legendre, Array, diff import numpy as np from FIAT.finite_element import FiniteElement @@ -14,6 +13,7 @@ def triangular_number(n): return int((n+1)*n/2) + def choose_ijk_total(degree): top = 1 for i in range(1, 2 + degree + 1): @@ -21,7 +21,7 @@ def choose_ijk_total(degree): bottom = 1 for i in range(1, degree + 1): bottom = i * bottom - return int(top /(2 * bottom)) + return int(top / (2 * bottom)) class TrimmedSerendipity(FiniteElement): @@ -33,7 +33,7 @@ def __init__(self, ref_el, degree, mapping): dim = flat_el.get_spatial_dimension() self.fdim = dim if dim != 3: - if dim !=2: + if dim != 2: raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3") flat_topology = flat_el.get_topology() @@ -44,7 +44,7 @@ def __init__(self, ref_el, degree, mapping): for entity in entities: entity_ids[top_dim][entity] = [] - #3-d case. + # 3-d case. if dim == 3: entity_ids[3] = {} for j in sorted(flat_topology[2]): @@ -60,7 +60,6 @@ def __init__(self, ref_el, degree, mapping): if (degree == 4): interior_tilde_ids += choose_ijk_total(degree - 2) - (degree - 1) - (degree - 1) + 1 if (degree > 4): - #interior_tilde_ids += choose_ijk_total(degree - 2) - (2 * degree - 1) interior_tilde_ids += choose_ijk_total(degree - 2) - (degree - 1) - (degree - 1) + 1 if degree == 1: interior_tilde_ids = 0 @@ -73,7 +72,6 @@ def __init__(self, ref_el, degree, mapping): if(degree >= 2): entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2) + degree)) - #entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 2))) cur += 2*triangular_number(degree - 2) + degree formdegree = dim - 1 @@ -177,7 +175,7 @@ def __init__(self, ref_el, degree): flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - if dim !=3: + if dim != 3: raise Exception("Trimmed serendipity face elements only valid for dimensions 2 and 3") verts = flat_el.get_vertices() @@ -201,19 +199,15 @@ def __init__(self, ref_el, degree): Sminus_list = FL + IL self.basis = {(0, 0, 0): Array(Sminus_list)} super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") - + else: - ##Put all 2 dimensional stuff here. + # Put all 2 dimensional stuff here. if degree < 1: raise Exception("Trimmed serendipity face elements only valid for k >= 1") - #flat_el = flatten_reference_cube(ref_el) - #verts = flat_el.get_vertices() - EL = e_lambda_1_2d_part_one(degree, dx, dy, x_mid, y_mid) if degree >= 2: FL = trimmed_f_lambda_2d(degree, dx, dy, x_mid, y_mid) - #FL = F_lambda_1_2d(degree, dx, dy, x_mid, y_mid) else: FL = () Sminus_list = EL + FL @@ -222,54 +216,11 @@ def __init__(self, ref_el, degree): super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") -""" def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): - FLpiece = tuple([]) - #Dx's first. - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(0, leg(j, x_mid) * leg(i - j, z_mid) * dy[k], 0)]) - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) - return FLpiece """ - - -""" def f_lambda_2_3d_pieces(deg, dx, dy, dz, x_mid, y_mid, z_mid): - FLpiece = tuple([]) - #Dx's first. - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(0, 0, leg(j, x_mid) * leg(i - j, y_mid) * dz[k])]) - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(0, leg(j, x_mid) * leg(i - j, z_mid) * dy[k], 0)]) - for k in range(0, 2): - for i in range(0, deg): - for j in range(0, i + 1): - FLpiece += tuple([(leg(j, y_mid) * leg(i - j, z_mid) * dx[k], 0, 0)]) - return FLpiece - - -def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): - FL = tuple([]) - #for j in range(0, degree): - # FL += f_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) - FL += f_lambda_2_3d_pieces(degree, dx, dy, dz, x_mid, y_mid, z_mid) - return FL """ - def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): FL = tuple([(-leg(j, y_mid) * leg(k, z_mid) * a, 0, 0) for a in dx for k in range(0, degree) for j in range(0, degree - k)] + - [(0, leg(j, x_mid) * leg(k, z_mid) * b, 0) + [(0, leg(j, x_mid) * leg(k, z_mid) * b, 0) for b in dy for k in range(0, degree) for j in range(0, degree - k)] + [(0, 0, -leg(j, x_mid) * leg(k, y_mid) * c) for c in dz for k in range(0, degree) for j in range(0, degree - k)]) @@ -279,30 +230,16 @@ def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): def I_lambda_2_3d_pieces(current_deg, dx, dy, dz, x_mid, y_mid, z_mid): assert current_deg > 1, 'invalid for i = 1' ILpiece = tuple([]) - for j in range(0, current_deg -1): + for j in range(0, current_deg - 1): for k in range(0, current_deg - 1 - j): ILpiece += tuple([(0, 0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * - dz[0] * dz[1])]+ - [(0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * - dy[0] * dy[1] ,0)] + - [(-leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * - dx[1], 0, 0)]) + dz[0] * dz[1])] + + [(0, -leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * + dy[0] * dy[1], 0)] + + [(-leg(j, x_mid) * leg(k, y_mid) * leg(current_deg - 2 - j - k, z_mid) * dx[0] * + dx[1], 0, 0)]) return ILpiece -""" def i_lambda_2_3d_normal(degree, dx, dy, dz, x_mid, y_mid, z_mid): - - IL = tuple([(-leg(l-2-j, x_mid) * leg(j-k, y_mid) * leg(k, z_mid) * - dx[0] * dx[1], 0, 0) - for l in range(2, degree) for j in range(l-1) for k in range(j+1)] + - [(0,-leg(l-2-j, x_mid) * leg(j-k, y_mid) * leg(k, z_mid) * - dy[0] * dy[1], 0, 0) - for l in range(2, degree) for j in range(l-1) for k in range(j+1)] + - [(0, 0, -leg(l-2-j, x_mid) * leg(j-k, y_mid) * leg(k, z_mid) * - dz[0] * dz[1]) - for l in range(2, degree) for j in range(l-1) for k in range(j+1)]) - - return IL """ - def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): assert degree > 1, 'invalid for i = 1' @@ -310,15 +247,15 @@ def I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid): [(0, leg(degree - 2, y_mid) * dy[0] * dy[1], 0)] + [(leg(degree - 2, x_mid) * dx[0] * dx[1], 0, 0)]) IL_tilde += tuple([(leg(degree - j - 2, x_mid) * leg(j, y_mid) * dx[0] * dx[1], leg(degree - j - 1, x_mid) * - leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + + leg(j - 1, y_mid) * dy[0] * dy[1], 0) for j in range(1, degree - 1)] + [(leg(degree - j - 2, x_mid) * leg(j, z_mid) * dx[0] * dx[1], 0, leg(degree - j - 1, x_mid) * - leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)] + + leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)] + [(0, leg(degree - j - 2, y_mid) * leg(j, z_mid) * dy[0] * dy[1], leg(degree - j - 1, y_mid) * - leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) + leg(j - 1, z_mid) * dz[0] * dz[1]) for j in range(1, degree - 1)]) for k in range(1, degree - 2): for l in range(1, degree - 1 - k): j = degree - 2 - k - l - IL_tilde += tuple([(-leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], + IL_tilde += tuple([(-leg(j, x_mid) * leg(k, y_mid) * leg(l, z_mid) * dx[0] * dx[1], leg(j + 1, x_mid) * leg(k - 1, y_mid) * leg(l, z_mid) * dy[0] * dy[1], -leg(j + 1, x_mid) * leg(k, y_mid) * leg(l - 1, z_mid) * dz[0] * dz[1])]) return IL_tilde @@ -328,11 +265,11 @@ def I_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): IL = tuple([]) for j in range(2, degree): IL += I_lambda_2_3d_pieces(j, dx, dy, dz, x_mid, y_mid, z_mid) - #IL += i_lambda_2_3d_normal(degree, dx, dy, dz, x_mid, y_mid, z_mid) IL += I_lambda_2_3d_tilde(degree, dx, dy, dz, x_mid, y_mid, z_mid) return IL -#Everything for 2-d should work already. + +# Everything for 2-d should work already. def e_lambda_1_2d_part_one(deg, dx, dy, x_mid, y_mid): EL = tuple( [(0, -leg(j, y_mid) * dx[0]) for j in range(deg)] + @@ -417,9 +354,7 @@ def F_lambda_1_2d(deg, dx, dy, x_mid, y_mid): def trimmed_f_lambda_2d(deg, dx, dy, x_mid, y_mid): - #FL = f_lambda_1_2d_trim(deg, dx, dy, x_mid, y_mid) FL = F_lambda_1_2d(deg, dx, dy, x_mid, y_mid) FLT = f_lambda_1_2d_tilde(deg, dx, dy, x_mid, y_mid) result = FL + FLT - # return FL return result diff --git a/FIAT/__init__.py b/FIAT/__init__.py index b2bbd9a83..38c156143 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -13,8 +13,8 @@ from FIAT.brezzi_douglas_marini import BrezziDouglasMarini from FIAT.Sminus import TrimmedSerendipityEdge # noqa: F401 from FIAT.Sminus import TrimmedSerendipityFace # noqa: F401 -from FIAT.SminusDiv import TrimmedSerendipityDiv #noqa: F401 -from FIAT.SminusCurl import TrimmedSerendipityCurl #noqa: F401 +from FIAT.SminusDiv import TrimmedSerendipityDiv # noqa: F401 +from FIAT.SminusCurl import TrimmedSerendipityCurl # noqa: F401 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor