From 239ed126690c426384d9c5d2234020cb0eafb02c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 28 Mar 2024 18:23:15 -0500 Subject: [PATCH 01/13] Let FIAT handle general CG/DG variants --- tsfc/finatinterface.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index b0c813ed..b4d18fdb 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -144,9 +144,7 @@ def convert_finiteelement(element, **kwargs): kind = 'spectral' # default variant if element.family() == "Lagrange": - if kind == 'equispaced': - lmbda = finat.Lagrange - elif kind == 'spectral': + if kind == 'spectral': lmbda = finat.GaussLobattoLegendre elif kind == 'integral': lmbda = finat.IntegratedLegendre @@ -167,14 +165,13 @@ def convert_finiteelement(element, **kwargs): deps = {"shift_axes", "restriction"} return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction), deps else: - raise ValueError("Variant %r not supported on %s" % (kind, element.cell)) + # Let FIAT handle the general case + lmbda = finat.Lagrange elif element.family() in {"Raviart-Thomas", "Nedelec 1st kind H(curl)", "Brezzi-Douglas-Marini", "Nedelec 2nd kind H(curl)"}: lmbda = partial(lmbda, variant=element.variant()) elif element.family() in ["Discontinuous Lagrange", "Discontinuous Lagrange L2"]: - if kind == 'equispaced': - lmbda = finat.DiscontinuousLagrange - elif kind == 'spectral': + if kind == 'spectral': lmbda = finat.GaussLegendre elif kind == 'integral': lmbda = finat.Legendre @@ -191,7 +188,8 @@ def convert_finiteelement(element, **kwargs): deps = {"shift_axes", "restriction"} return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction, continuous=False), deps else: - raise ValueError("Variant %r not supported on %s" % (kind, element.cell)) + # Let FIAT handle the general case + lmbda = finat.DiscontinuousLagrange elif element.family() == ["DPC", "DPC L2"]: if element.cell.geometric_dimension() == 2: element = element.reconstruct(cell=ufl.cell.hypercube(2)) From 23ae55f8d5e7b8022aec7992a8f11c5ff2962f83 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 28 Mar 2024 18:47:47 -0500 Subject: [PATCH 02/13] pass the variant --- tsfc/finatinterface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index b4d18fdb..b0501356 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -166,7 +166,7 @@ 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 = finat.Lagrange + lmbda = partial(finat.Lagrange, variant=kind) elif element.family() in {"Raviart-Thomas", "Nedelec 1st kind H(curl)", "Brezzi-Douglas-Marini", "Nedelec 2nd kind H(curl)"}: lmbda = partial(lmbda, variant=element.variant()) @@ -189,7 +189,7 @@ 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 = finat.DiscontinuousLagrange + lmbda = partial(finat.DiscontinuousLagrange, variant=kind) elif element.family() == ["DPC", "DPC L2"]: if element.cell.geometric_dimension() == 2: element = element.reconstruct(cell=ufl.cell.hypercube(2)) From f3add68f65f4a2d47b4b6b9d74a2f524004191dc Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 2 Apr 2024 13:50:27 -0500 Subject: [PATCH 03/13] fix tests --- tests/test_create_fiat_element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_create_fiat_element.py b/tests/test_create_fiat_element.py index 28356b00..f8a7d6ef 100644 --- a/tests/test_create_fiat_element.py +++ b/tests/test_create_fiat_element.py @@ -1,7 +1,7 @@ import pytest import FIAT -from FIAT.discontinuous_lagrange import HigherOrderDiscontinuousLagrange as FIAT_DiscontinuousLagrange +from FIAT.discontinuous_lagrange import DiscontinuousLagrange as FIAT_DiscontinuousLagrange import ufl import finat.ufl From be216cc83a96754c231e18516bd8b05ff2bc187f Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Tue, 2 Apr 2024 22:50:11 -0500 Subject: [PATCH 04/13] WIP: enable macro-quadrature --- tsfc/kernel_interface/common.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 67f7dac8..6a1514be 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -1,29 +1,24 @@ import collections -import string import operator +import string from functools import reduce from itertools import chain, product +import gem +import gem.impero_utils as impero_utils import numpy -from numpy import asarray - -from ufl.utils.sequences import max_degree - from FIAT.reference_element import TensorProductCell - from finat.quadrature import AbstractQuadratureRule, make_quadrature - -import gem - from gem.node import traversal +from gem.optimise import constant_fold_zero +from gem.optimise import remove_componenttensors as prune from gem.utils import cached_property -import gem.impero_utils as impero_utils -from gem.optimise import remove_componenttensors as prune, constant_fold_zero - +from numpy import asarray from tsfc import fem, ufl_utils +from tsfc.finatinterface import as_fiat_cell, convert, create_element from tsfc.kernel_interface import KernelInterface -from tsfc.finatinterface import as_fiat_cell, create_element from tsfc.logging import logger +from ufl.utils.sequences import max_degree class KernelBuilderBase(KernelInterface): @@ -301,7 +296,8 @@ def set_quad_rule(params, cell, integral_type, functions): quadrature_degree = params["quadrature_degree"] except KeyError: quadrature_degree = params["estimated_polynomial_degree"] - function_degrees = [f.ufl_function_space().ufl_element().degree() for f in functions] + function_degrees = [f.ufl_function_space().ufl_element().degree() + for f in functions] if all((asarray(quadrature_degree) > 10 * asarray(degree)).all() for degree in function_degrees): logger.warning("Estimated quadrature degree %s more " @@ -314,9 +310,14 @@ def set_quad_rule(params, cell, integral_type, functions): quad_rule = params["quadrature_rule"] except KeyError: fiat_cell = as_fiat_cell(cell) + finat_elements = set(create_element(f.ufl_element()) for f in functions) + print(list(f.degree for f in finat_elements)) + + fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] + integration_dim, _ = lower_integral_type(fiat_cell, integral_type) - integration_cell = fiat_cell.construct_subelement(integration_dim) - quad_rule = make_quadrature(integration_cell, quadrature_degree) + + quad_rule = make_quadrature(fiat_cells, quadrature_degree, dim=integration_dim) params["quadrature_rule"] = quad_rule if not isinstance(quad_rule, AbstractQuadratureRule): From 27cb54fa4676253f9f9bc43f67f618b4c33f2bbf Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Tue, 2 Apr 2024 23:28:29 -0500 Subject: [PATCH 05/13] flake8 --- tsfc/kernel_interface/common.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 6a1514be..b455516c 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -15,7 +15,7 @@ from gem.utils import cached_property from numpy import asarray from tsfc import fem, ufl_utils -from tsfc.finatinterface import as_fiat_cell, convert, create_element +from tsfc.finatinterface import as_fiat_cell, create_element from tsfc.kernel_interface import KernelInterface from tsfc.logging import logger from ufl.utils.sequences import max_degree From b55e03a65ed5b9829528cf1f7ce93a61844b3512 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 3 Apr 2024 13:22:59 -0500 Subject: [PATCH 06/13] Handle quadrature for macro elements --- tsfc/kernel_interface/common.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index b455516c..2e097efc 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -304,20 +304,22 @@ def set_quad_rule(params, cell, integral_type, functions): "than tenfold greater than any " "argument/coefficient degree (max %s)", quadrature_degree, max_degree(function_degrees)) - if params.get("quadrature_rule") == "default": - del params["quadrature_rule"] - try: - quad_rule = params["quadrature_rule"] - except KeyError: + quad_rule = params.get("quadrature_rule", "default") + if isinstance(quad_rule, str): + scheme = quad_rule fiat_cell = as_fiat_cell(cell) - finat_elements = set(create_element(f.ufl_element()) for f in functions) - print(list(f.degree for f in finat_elements)) + integration_dim, _ = lower_integral_type(fiat_cell, integral_type) + finat_elements = set(create_element(f.ufl_element()) for f in functions) fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] + max_cell = max(fiat_cells) + if all(max_cell >= b for b in fiat_cells): + fiat_cell = max_cell + else: + raise ValueError("Can't find a maximal complex") - integration_dim, _ = lower_integral_type(fiat_cell, integral_type) - - quad_rule = make_quadrature(fiat_cells, quadrature_degree, dim=integration_dim) + fiat_cell = fiat_cell.construct_subcomplex(integration_dim) + quad_rule = make_quadrature(fiat_cell, quadrature_degree, scheme=scheme) params["quadrature_rule"] = quad_rule if not isinstance(quad_rule, AbstractQuadratureRule): From 89e9820fdafdfa56805779e10f3d01eb5fc24144 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 3 Apr 2024 20:18:25 -0500 Subject: [PATCH 07/13] do not create complex on Real (Constant) spaces --- tsfc/kernel_interface/common.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 2e097efc..73a76374 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -310,7 +310,8 @@ def set_quad_rule(params, cell, integral_type, functions): fiat_cell = as_fiat_cell(cell) integration_dim, _ = lower_integral_type(fiat_cell, integral_type) - finat_elements = set(create_element(f.ufl_element()) for f in functions) + finat_elements = set(create_element(f.ufl_element()) for f in functions + if f.ufl_element().family() != "Real") fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] max_cell = max(fiat_cells) if all(max_cell >= b for b in fiat_cells): From b43f74e30ca87d8089dd4721e490e3983c0c6436 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 3 Apr 2024 21:58:43 -0500 Subject: [PATCH 08/13] handle integral variants --- tsfc/finatinterface.py | 8 ++++---- tsfc/kernel_interface/common.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index b0501356..e570bbea 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -146,8 +146,8 @@ def convert_finiteelement(element, **kwargs): if element.family() == "Lagrange": if kind == 'spectral': lmbda = finat.GaussLobattoLegendre - elif kind == 'integral': - lmbda = finat.IntegratedLegendre + 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: @@ -173,8 +173,8 @@ def convert_finiteelement(element, **kwargs): elif element.family() in ["Discontinuous Lagrange", "Discontinuous Lagrange L2"]: if kind == 'spectral': lmbda = finat.GaussLegendre - elif kind == 'integral': - lmbda = finat.Legendre + 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: diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 73a76374..54aac54a 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -319,7 +319,7 @@ def set_quad_rule(params, cell, integral_type, functions): else: raise ValueError("Can't find a maximal complex") - fiat_cell = fiat_cell.construct_subcomplex(integration_dim) + integration_cell = fiat_cell.construct_subcomplex(integration_dim) quad_rule = make_quadrature(fiat_cell, quadrature_degree, scheme=scheme) params["quadrature_rule"] = quad_rule From c96e6e546cd3aebceea28660757f7b78ff29f7d3 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 4 Apr 2024 18:33:18 -0500 Subject: [PATCH 09/13] use finat's max complex logic --- tsfc/kernel_interface/common.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 54aac54a..4fe9f543 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -8,6 +8,7 @@ import gem.impero_utils as impero_utils import numpy from FIAT.reference_element import TensorProductCell +from finat.cell_tools import max_complex from finat.quadrature import AbstractQuadratureRule, make_quadrature from gem.node import traversal from gem.optimise import constant_fold_zero @@ -313,14 +314,10 @@ def set_quad_rule(params, cell, integral_type, functions): finat_elements = set(create_element(f.ufl_element()) for f in functions if f.ufl_element().family() != "Real") fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] - max_cell = max(fiat_cells) - if all(max_cell >= b for b in fiat_cells): - fiat_cell = max_cell - else: - raise ValueError("Can't find a maximal complex") + max_cell = max_complex(fiat_cells) integration_cell = fiat_cell.construct_subcomplex(integration_dim) - quad_rule = make_quadrature(fiat_cell, quadrature_degree, scheme=scheme) + quad_rule = make_quadrature(integration_cell, quadrature_degree, scheme=scheme) params["quadrature_rule"] = quad_rule if not isinstance(quad_rule, AbstractQuadratureRule): From fe988ef7f63f91b5d8f73d195831e1acb8200f2e Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 4 Apr 2024 18:36:01 -0500 Subject: [PATCH 10/13] fix an error, flake8 --- tsfc/kernel_interface/common.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 4fe9f543..4b0370e9 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -316,7 +316,7 @@ def set_quad_rule(params, cell, integral_type, functions): fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] max_cell = max_complex(fiat_cells) - integration_cell = fiat_cell.construct_subcomplex(integration_dim) + integration_cell = max_cell.construct_subcomplex(integration_dim) quad_rule = make_quadrature(integration_cell, quadrature_degree, scheme=scheme) params["quadrature_rule"] = quad_rule From cfef1d7a184c87f957a1f6a375186f6c2200eb08 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 4 Apr 2024 18:43:55 -0500 Subject: [PATCH 11/13] lower integral type on the final cell --- tsfc/kernel_interface/common.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 4fe9f543..18fd363f 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -309,13 +309,12 @@ def set_quad_rule(params, cell, integral_type, functions): if isinstance(quad_rule, str): scheme = quad_rule fiat_cell = as_fiat_cell(cell) - integration_dim, _ = lower_integral_type(fiat_cell, integral_type) - finat_elements = set(create_element(f.ufl_element()) for f in functions if f.ufl_element().family() != "Real") fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] - max_cell = max_complex(fiat_cells) + fiat_cell = max_complex(fiat_cells) + integration_dim, _ = lower_integral_type(fiat_cell, integral_type) integration_cell = fiat_cell.construct_subcomplex(integration_dim) quad_rule = make_quadrature(integration_cell, quadrature_degree, scheme=scheme) params["quadrature_rule"] = quad_rule From 436328129006f03737a2b56c096c1f3578d004f3 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 12 Apr 2024 11:35:08 -0500 Subject: [PATCH 12/13] register HCT --- tsfc/finatinterface.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index e570bbea..ac73a681 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -52,6 +52,7 @@ "Hermite": finat.Hermite, "Kong-Mulder-Veldhuizen": finat.KongMulderVeldhuizen, "Argyris": finat.Argyris, + "Hsieh-Clough-Tocher": finat.HsiehCloughTocher, "Mardal-Tai-Winther": finat.MardalTaiWinther, "Morley": finat.Morley, "Bell": finat.Bell, From 916e773117b669711f45cfb25e0d70cfae06f72c Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 19 Apr 2024 22:34:09 -0500 Subject: [PATCH 13/13] Add reduced element --- tsfc/finatinterface.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index ac73a681..320b3e32 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -19,14 +19,13 @@ # You should have received a copy of the GNU Lesser General Public License # along with FFC. If not, see . -from functools import singledispatch, partial import weakref +from functools import partial, singledispatch import FIAT import finat -import ufl import finat.ufl - +import ufl __all__ = ("as_fiat_cell", "create_base_element", "create_element", "supported_elements") @@ -53,6 +52,7 @@ "Kong-Mulder-Veldhuizen": finat.KongMulderVeldhuizen, "Argyris": finat.Argyris, "Hsieh-Clough-Tocher": finat.HsiehCloughTocher, + "Reduced-Hsieh-Clough-Tocher": finat.ReducedHsiehCloughTocher, "Mardal-Tai-Winther": finat.MardalTaiWinther, "Morley": finat.Morley, "Bell": finat.Bell,