From 7bfbb08d213c265f5384b15d1eb7eca5c81647f0 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Wed, 20 Mar 2024 11:43:48 -0400 Subject: [PATCH 1/5] Implementation of Hu-Zhang elements --- tsfc/finatinterface.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index fb4e9cef..3d298488 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -48,6 +48,7 @@ "Hellan-Herrmann-Johnson": finat.HellanHerrmannJohnson, "Nonconforming Arnold-Winther": finat.ArnoldWintherNC, "Conforming Arnold-Winther": finat.ArnoldWinther, + "Hu-Zhang": finat.HuZhang, "Hermite": finat.Hermite, "Kong-Mulder-Veldhuizen": finat.KongMulderVeldhuizen, "Argyris": finat.Argyris, From 988c85db4ff5764fb7f67c2c7170c3649d265f7f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 6 May 2024 10:46:24 +0100 Subject: [PATCH 2/5] Support covariant contravariant piola mapping for H(curl div) --- tsfc/ufl_utils.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tsfc/ufl_utils.py b/tsfc/ufl_utils.py index 0cec4dbf..bc5069b6 100644 --- a/tsfc/ufl_utils.py +++ b/tsfc/ufl_utils.py @@ -403,6 +403,10 @@ def apply_mapping(expression, element, domain): G(X) = det(J)^2 K g(x) K^T i.e. G_il(X)=(detJ)^2 K_ij g_jk K_lk + 'covariant contravariant piola' mapping for g: + + G(X) = det(J) J^T g(x) K^T i.e. G_il(X) = det(J) J_ji g_jk(x) K_lk + If 'contravariant piola' or 'covariant piola' (or their double variants) are applied to a matrix-valued function, the appropriate mappings are applied row-by-row. @@ -442,6 +446,13 @@ def apply_mapping(expression, element, domain): *k, i, j, m, n = indices(len(expression.ufl_shape) + 2) kmn = (*k, m, n) rexpression = as_tensor(detJ**2 * K[i, m] * expression[kmn] * K[j, n], (*k, i, j)) + elif mapping == "covariant contravariant piola": + J = Jacobian(mesh) + K = JacobianInverse(mesh) + detJ = JacobianDeterminant(mesh) + *k, i, j, m, n = indices(len(expression.ufl_shape) + 2) + kmn = (*k, m, n) + rexpression = as_tensor(detJ * J[m, i] * expression[kmn] * K[j, n], (*k, i, j)) elif mapping == "symmetries": # This tells us how to get from the pieces of the reference # space expression to the physical space one. From 8387a75443a697a63155a1ece42cb8999f331541 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 6 May 2024 11:09:57 +0100 Subject: [PATCH 3/5] Add GLS H(curl div) element --- tsfc/finatinterface.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index 320b3e32..0e6e741e 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -61,6 +61,7 @@ "Nedelec 2nd kind H(curl)": finat.NedelecSecondKind, "Raviart-Thomas": finat.RaviartThomas, "Regge": finat.Regge, + "Gopalakrishnan-Lederer-Schoberl": finat.GopalakrishnanLedererSchoberl, "BDMCE": finat.BrezziDouglasMariniCubeEdge, "BDMCF": finat.BrezziDouglasMariniCubeFace, # These require special treatment below From d7f064c828399131d2b7d168ef0b650e18d175fd Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 24 Oct 2024 17:58:26 +0100 Subject: [PATCH 4/5] Refactor variants --- tsfc/finatinterface.py | 75 ++++++++++++++++++++++-------------------- 1 file changed, 39 insertions(+), 36 deletions(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index f9f461fc..2e524735 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -20,7 +20,7 @@ # along with FFC. If not, see . import weakref -from functools import partial, singledispatch +from functools import singledispatch import FIAT import finat @@ -124,6 +124,23 @@ def convert(element, **kwargs): raise ValueError("Unsupported element type %s" % type(element)) +cg_interval_variants = { + "fdm": finat.FDMLagrange, + "fdm_ipdg": finat.FDMLagrange, + "fdm_quadrature": finat.FDMQuadrature, + "fdm_broken": finat.FDMBrokenH1, + "fdm_hermite": finat.FDMHermite, +} + + +dg_interval_variants = { + "fdm": finat.FDMDiscontinuousLagrange, + "fdm_quadrature": finat.FDMDiscontinuousLagrange, + "fdm_ipdg": lambda *args: finat.DiscontinuousElement(finat.FDMLagrange(*args)), + "fdm_broken": finat.FDMBrokenL2, +} + + # Base finite elements first @convert.register(finat.ufl.FiniteElement) def convert_finiteelement(element, **kwargs): @@ -152,30 +169,19 @@ def convert_finiteelement(element, **kwargs): finat_elem, deps = _create_element(element, **kwargs) return finat.FlattenedDimensions(finat_elem), deps + kw = {} kind = element.variant() if kind is None: kind = 'spectral' # default variant - is_interval = element.cell.cellname() == 'interval' - if element.family() in {"Raviart-Thomas", "Nedelec 1st kind H(curl)", - "Brezzi-Douglas-Marini", "Nedelec 2nd kind H(curl)", - "Argyris"}: - lmbda = partial(lmbda, variant=element.variant()) - elif element.family() == "Lagrange": + if element.family() == "Lagrange": if kind == 'spectral': lmbda = finat.GaussLobattoLegendre - elif kind.startswith('integral'): - lmbda = partial(finat.IntegratedLegendre, variant=kind) - elif kind in ['fdm', 'fdm_ipdg'] and is_interval: - lmbda = finat.FDMLagrange - elif kind == 'fdm_quadrature' and is_interval: - lmbda = finat.FDMQuadrature - elif kind == 'fdm_broken' and is_interval: - lmbda = finat.FDMBrokenH1 - elif kind == 'fdm_hermite' and is_interval: - lmbda = finat.FDMHermite - elif kind in ['demkowicz', 'fdm']: - lmbda = partial(finat.IntegratedLegendre, variant=kind) + elif element.cell.cellname() == "interval" and kind in cg_interval_variants: + lmbda = cg_interval_variants[kind] + elif kind.startswith('integral') or kind in ['demkowicz', 'fdm']: + lmbda = finat.IntegratedLegendre + kw["variant"] = kind elif kind in ['mgd', 'feec', 'qb', 'mse']: degree = element.degree() shift_axes = kwargs["shift_axes"] @@ -184,20 +190,17 @@ def convert_finiteelement(element, **kwargs): return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction), deps else: # Let FIAT handle the general case - lmbda = partial(finat.Lagrange, variant=kind) + lmbda = finat.Lagrange + kw["variant"] = kind + elif element.family() in ["Discontinuous Lagrange", "Discontinuous Lagrange L2"]: if kind == 'spectral': lmbda = finat.GaussLegendre - elif kind.startswith('integral'): - lmbda = partial(finat.Legendre, variant=kind) - elif kind in ['fdm', 'fdm_quadrature'] and is_interval: - lmbda = finat.FDMDiscontinuousLagrange - elif kind == 'fdm_ipdg' and is_interval: - lmbda = lambda *args: finat.DiscontinuousElement(finat.FDMLagrange(*args)) - elif kind in 'fdm_broken' and is_interval: - lmbda = finat.FDMBrokenL2 - elif kind in ['demkowicz', 'fdm']: - lmbda = partial(finat.Legendre, variant=kind) + elif element.cell.cellname() == "interval" and kind in dg_interval_variants: + lmbda = dg_interval_variants[kind] + elif kind.startswith('integral') or kind in ['demkowicz', 'fdm']: + lmbda = finat.Legendre + kw["variant"] = kind elif kind in ['mgd', 'feec', 'qb', 'mse']: degree = element.degree() shift_axes = kwargs["shift_axes"] @@ -206,13 +209,13 @@ def convert_finiteelement(element, **kwargs): return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction, continuous=False), deps else: # Let FIAT handle the general case - lmbda = partial(finat.DiscontinuousLagrange, variant=kind) - elif element.family() == ["DPC", "DPC L2", "S"]: - dim = element.cell.geometric_dimension() - if dim > 1: - element = element.reconstruct(cell=ufl.cell.hypercube(dim)) + lmbda = finat.DiscontinuousLagrange + kw["variant"] = kind + + elif element.variant() is not None: + kw["variant"] = element.variant() - return lmbda(cell, element.degree()), set() + return lmbda(cell, element.degree(), **kw), set() # Element modifiers and compound element types From 893bebbea1703e0b686e12a84917d3822473882a Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 30 Oct 2024 15:35:29 +0000 Subject: [PATCH 5/5] GLS 1st and 2nd kinds --- tsfc/finatinterface.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index a8863d9d..36d8f6d6 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -75,7 +75,8 @@ "Nedelec 2nd kind H(curl)": finat.NedelecSecondKind, "Raviart-Thomas": finat.RaviartThomas, "Regge": finat.Regge, - "Gopalakrishnan-Lederer-Schoberl": finat.GopalakrishnanLedererSchoberl, + "Gopalakrishnan-Lederer-Schoberl 1st kind": finat.GopalakrishnanLedererSchoberlFirstKind, + "Gopalakrishnan-Lederer-Schoberl 2nd kind": finat.GopalakrishnanLedererSchoberlSecondKind, "BDMCE": finat.BrezziDouglasMariniCubeEdge, "BDMCF": finat.BrezziDouglasMariniCubeFace, # These require special treatment below