From 54de90707dd1e046fd67f43c93d3a62fd0919a29 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 10 Aug 2023 11:18:49 +0100 Subject: [PATCH 01/14] Work on removing switch/case from integrals/ids --- ffcx/codegeneration/form.py | 124 +++++++++++++++++++++------ ffcx/codegeneration/form_template.py | 9 +- ffcx/codegeneration/ufcx.h | 5 ++ 3 files changed, 112 insertions(+), 26 deletions(-) diff --git a/ffcx/codegeneration/form.py b/ffcx/codegeneration/form.py index a2bac06e5..5922dac5d 100644 --- a/ffcx/codegeneration/form.py +++ b/ffcx/codegeneration/form.py @@ -25,7 +25,7 @@ def generator(ir, options): d = {} d["factory_name"] = ir.name d["name_from_uflfile"] = ir.name_from_uflfile - d["signature"] = f"\"{ir.signature}\"" + d["signature"] = f'"{ir.signature}"' d["rank"] = ir.rank d["num_coefficients"] = ir.num_coefficients d["num_constants"] = ir.num_constants @@ -39,8 +39,11 @@ def generator(ir, options): if len(ir.original_coefficient_position) > 0: d["original_coefficient_position_init"] = L.ArrayDecl( - "int", f"original_coefficient_position_{ir.name}", - values=ir.original_coefficient_position, sizes=len(ir.original_coefficient_position)) + "int", + f"original_coefficient_position_{ir.name}", + values=ir.original_coefficient_position, + sizes=len(ir.original_coefficient_position), + ) d["original_coefficient_position"] = f"original_coefficient_position_{ir.name}" else: d["original_coefficient_position_init"] = "" @@ -49,7 +52,7 @@ def generator(ir, options): cnames = ir.coefficient_names assert ir.num_coefficients == len(cnames) names = L.Symbol("names") - if (len(cnames) == 0): + if len(cnames) == 0: code = [L.Return(L.Null())] else: code = [L.ArrayDecl("static const char*", names, len(cnames), cnames)] @@ -67,37 +70,95 @@ def generator(ir, options): if len(ir.finite_elements) > 0: d["finite_elements"] = f"finite_elements_{ir.name}" - d["finite_elements_init"] = L.ArrayDecl("ufcx_finite_element*", f"finite_elements_{ir.name}", values=[ - L.AddressOf(L.Symbol(el)) for el in ir.finite_elements], - sizes=len(ir.finite_elements)) + d["finite_elements_init"] = L.ArrayDecl( + "ufcx_finite_element*", + f"finite_elements_{ir.name}", + values=[L.AddressOf(L.Symbol(el)) for el in ir.finite_elements], + sizes=len(ir.finite_elements), + ) else: d["finite_elements"] = L.Null() d["finite_elements_init"] = "" if len(ir.dofmaps) > 0: d["dofmaps"] = f"dofmaps_{ir.name}" - d["dofmaps_init"] = L.ArrayDecl("ufcx_dofmap*", f"dofmaps_{ir.name}", values=[ - L.AddressOf(L.Symbol(dofmap)) for dofmap in ir.dofmaps], sizes=len(ir.dofmaps)) + d["dofmaps_init"] = L.ArrayDecl( + "ufcx_dofmap*", + f"dofmaps_{ir.name}", + values=[L.AddressOf(L.Symbol(dofmap)) for dofmap in ir.dofmaps], + sizes=len(ir.dofmaps), + ) else: d["dofmaps"] = L.Null() d["dofmaps_init"] = "" + integrals = [] + integral_ids = [] + integral_offsets = [0] + for itg_type in ("cell", "interior_facet", "exterior_facet"): + integrals += [L.AddressOf(L.Symbol(itg)) for itg in ir.integral_names[itg_type]] + integral_ids += ir.subdomain_ids[itg_type] + integral_offsets.append(len(integrals)) + + d["form_integrals_init"] = L.ArrayDecl( + "static ufcx_integral*", + f"integrals_{ir.name}", + values=integrals, + sizes=len(integrals), + ) + + d["form_integral_ids_init"] = L.ArrayDecl( + "int", + f"integral_ids_{ir.name}", + values=integral_ids, + sizes=len(integral_ids), + ) + + d["form_integral_offsets_init"] = L.ArrayDecl( + "int", + f"integral_offsets_{ir.name}", + values=integral_offsets, + sizes=len(integral_offsets), + ) + code = [] cases = [] code_ids = [] cases_ids = [] for itg_type in ("cell", "interior_facet", "exterior_facet"): if len(ir.integral_names[itg_type]) > 0: - code += [L.ArrayDecl( - "static ufcx_integral*", f"integrals_{itg_type}_{ir.name}", - values=[L.AddressOf(L.Symbol(itg)) for itg in ir.integral_names[itg_type]], - sizes=len(ir.integral_names[itg_type]))] - cases.append((L.Symbol(itg_type), L.Return(L.Symbol(f"integrals_{itg_type}_{ir.name}")))) - - code_ids += [L.ArrayDecl( - "static int", f"integral_ids_{itg_type}_{ir.name}", - values=ir.subdomain_ids[itg_type], sizes=len(ir.subdomain_ids[itg_type]))] - cases_ids.append((L.Symbol(itg_type), L.Return(L.Symbol(f"integral_ids_{itg_type}_{ir.name}")))) + code += [ + L.ArrayDecl( + "static ufcx_integral*", + f"integrals_{itg_type}_{ir.name}", + values=[ + L.AddressOf(L.Symbol(itg)) + for itg in ir.integral_names[itg_type] + ], + sizes=len(ir.integral_names[itg_type]), + ) + ] + cases.append( + ( + L.Symbol(itg_type), + L.Return(L.Symbol(f"integrals_{itg_type}_{ir.name}")), + ) + ) + + code_ids += [ + L.ArrayDecl( + "static int", + f"integral_ids_{itg_type}_{ir.name}", + values=ir.subdomain_ids[itg_type], + sizes=len(ir.subdomain_ids[itg_type]), + ) + ] + cases_ids.append( + ( + L.Symbol(itg_type), + L.Return(L.Symbol(f"integral_ids_{itg_type}_{ir.name}")), + ) + ) code += [L.Switch("integral_type", cases, default=L.Return(L.Null()))] code_ids += [L.Switch("integral_type", cases_ids, default=L.Return(L.Null()))] @@ -110,12 +171,19 @@ def generator(ir, options): # FIXME: Should be handled differently, revise how # ufcx_function_space is generated - for (name, (element, dofmap, cmap_family, cmap_degree, cmap_celltype, cmap_variant)) in ir.function_spaces.items(): + for name, ( + element, + dofmap, + cmap_family, + cmap_degree, + cmap_celltype, + cmap_variant, + ) in ir.function_spaces.items(): code += [f"static ufcx_function_space functionspace_{name} ="] code += ["{"] code += [f".finite_element = &{element},"] code += [f".dofmap = &{dofmap},"] - code += [f".geometry_family = \"{cmap_family}\","] + code += [f'.geometry_family = "{cmap_family}",'] code += [f".geometry_degree = {cmap_degree},"] code += [f".geometry_basix_cell = {int(cmap_celltype)},"] code += [f".geometry_basix_variant = {int(cmap_variant)}"] @@ -133,14 +201,20 @@ def generator(ir, options): # Check that no keys are redundant or have been missed from string import Formatter - fields = [fname for _, fname, _, _ in Formatter().parse(form_template.factory) if fname] - assert set(fields) == set(d.keys()), "Mismatch between keys in template and in formatting dict" + + fields = [ + fname for _, fname, _, _ in Formatter().parse(form_template.factory) if fname + ] + assert set(fields) == set( + d.keys() + ), "Mismatch between keys in template and in formatting dict" # Format implementation code implementation = form_template.factory.format_map(d) # Format declaration - declaration = form_template.declaration.format(factory_name=d["factory_name"], - name_from_uflfile=d["name_from_uflfile"]) + declaration = form_template.declaration.format( + factory_name=d["factory_name"], name_from_uflfile=d["name_from_uflfile"] + ) return declaration, implementation diff --git a/ffcx/codegeneration/form_template.py b/ffcx/codegeneration/form_template.py index 749d3a812..ce80894dd 100644 --- a/ffcx/codegeneration/form_template.py +++ b/ffcx/codegeneration/form_template.py @@ -24,6 +24,9 @@ {original_coefficient_position_init} {dofmaps_init} {finite_elements_init} +{form_integral_offsets_init} +{form_integrals_init} +{form_integral_ids_init} // Return a list of the coefficient names. const char** coefficient_name_{factory_name}(void) @@ -70,7 +73,11 @@ .integral_ids = integral_ids_{factory_name}, .num_integrals = num_integrals_{factory_name}, - .integrals = integrals_{factory_name} + .integrals = integrals_{factory_name}, + + .form_integrals = form_integrals_{factory_name}, + .form_integral_ids = form_integral_ids_{factory_name}, + .form_integral_offsets = form_integral_offsets_{factory_name} }}; // Alias name diff --git a/ffcx/codegeneration/ufcx.h b/ffcx/codegeneration/ufcx.h index 8e84112d7..9b678aedc 100644 --- a/ffcx/codegeneration/ufcx.h +++ b/ffcx/codegeneration/ufcx.h @@ -462,6 +462,11 @@ extern "C" /// Get an integral on sub domain subdomain_id ufcx_integral** (*integrals)(ufcx_integral_type); + /// New interface, static data, no function calls + ufcx_integral* form_integrals; + int* form_integral_ids; + int* form_integral_offsets; + } ufcx_form; // FIXME: Formalise a UFCX 'function space' From 8b7fef2ee034ea48cdb668feaaa88cd511c33e7b Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 10 Aug 2023 11:26:31 +0100 Subject: [PATCH 02/14] Fix typos --- ffcx/codegeneration/form.py | 6 +++--- ffcx/codegeneration/ufcx.h | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ffcx/codegeneration/form.py b/ffcx/codegeneration/form.py index 5922dac5d..4b9a7a9b3 100644 --- a/ffcx/codegeneration/form.py +++ b/ffcx/codegeneration/form.py @@ -102,21 +102,21 @@ def generator(ir, options): d["form_integrals_init"] = L.ArrayDecl( "static ufcx_integral*", - f"integrals_{ir.name}", + f"form_integrals_{ir.name}", values=integrals, sizes=len(integrals), ) d["form_integral_ids_init"] = L.ArrayDecl( "int", - f"integral_ids_{ir.name}", + f"form_integral_ids_{ir.name}", values=integral_ids, sizes=len(integral_ids), ) d["form_integral_offsets_init"] = L.ArrayDecl( "int", - f"integral_offsets_{ir.name}", + f"form_integral_offsets_{ir.name}", values=integral_offsets, sizes=len(integral_offsets), ) diff --git a/ffcx/codegeneration/ufcx.h b/ffcx/codegeneration/ufcx.h index 9b678aedc..ae48e4716 100644 --- a/ffcx/codegeneration/ufcx.h +++ b/ffcx/codegeneration/ufcx.h @@ -463,7 +463,7 @@ extern "C" ufcx_integral** (*integrals)(ufcx_integral_type); /// New interface, static data, no function calls - ufcx_integral* form_integrals; + ufcx_integral** form_integrals; int* form_integral_ids; int* form_integral_offsets; From 76131845655dec916386a64825606de0f3895421 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 10 Aug 2023 12:11:00 +0100 Subject: [PATCH 03/14] Fix for zero integrals case --- ffcx/codegeneration/form.py | 33 +++++++++++++++++----------- ffcx/codegeneration/form_template.py | 4 ++-- 2 files changed, 22 insertions(+), 15 deletions(-) diff --git a/ffcx/codegeneration/form.py b/ffcx/codegeneration/form.py index 4b9a7a9b3..962666f24 100644 --- a/ffcx/codegeneration/form.py +++ b/ffcx/codegeneration/form.py @@ -100,19 +100,26 @@ def generator(ir, options): integral_ids += ir.subdomain_ids[itg_type] integral_offsets.append(len(integrals)) - d["form_integrals_init"] = L.ArrayDecl( - "static ufcx_integral*", - f"form_integrals_{ir.name}", - values=integrals, - sizes=len(integrals), - ) - - d["form_integral_ids_init"] = L.ArrayDecl( - "int", - f"form_integral_ids_{ir.name}", - values=integral_ids, - sizes=len(integral_ids), - ) + if len(integrals) > 0: + d["form_integrals_init"] = L.ArrayDecl( + "static ufcx_integral*", + f"form_integrals_{ir.name}", + values=integrals, + sizes=len(integrals), + ) + d["form_integrals"] = f"form_integrals_{ir.name}" + d["form_integral_ids_init"] = L.ArrayDecl( + "int", + f"form_integral_ids_{ir.name}", + values=integral_ids, + sizes=len(integral_ids), + ) + d["form_integral_ids"] = f"form_integral_ids_{ir.name}" + else: + d["form_integrals_init"] = "" + d["form_integals"] = "NULL" + d["form_integrals_ids_init"] = "" + d["form_integral_ids"] = "NULL" d["form_integral_offsets_init"] = L.ArrayDecl( "int", diff --git a/ffcx/codegeneration/form_template.py b/ffcx/codegeneration/form_template.py index ce80894dd..5c0b1d22c 100644 --- a/ffcx/codegeneration/form_template.py +++ b/ffcx/codegeneration/form_template.py @@ -75,8 +75,8 @@ .integrals = integrals_{factory_name}, - .form_integrals = form_integrals_{factory_name}, - .form_integral_ids = form_integral_ids_{factory_name}, + .form_integrals = {form_integrals}, + .form_integral_ids = {form_integral_ids}, .form_integral_offsets = form_integral_offsets_{factory_name} }}; From ba2cab9c9cbdd38fa036f0e752b918039ee4c608 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 10 Aug 2023 12:38:47 +0100 Subject: [PATCH 04/14] Typos --- ffcx/codegeneration/form.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ffcx/codegeneration/form.py b/ffcx/codegeneration/form.py index 962666f24..6de4c46b9 100644 --- a/ffcx/codegeneration/form.py +++ b/ffcx/codegeneration/form.py @@ -117,8 +117,8 @@ def generator(ir, options): d["form_integral_ids"] = f"form_integral_ids_{ir.name}" else: d["form_integrals_init"] = "" - d["form_integals"] = "NULL" - d["form_integrals_ids_init"] = "" + d["form_integrals"] = "NULL" + d["form_integral_ids_init"] = "" d["form_integral_ids"] = "NULL" d["form_integral_offsets_init"] = L.ArrayDecl( From 10d357da6b96e02106869e42a54ee2c9b9e96c16 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 10 Aug 2023 14:01:08 +0100 Subject: [PATCH 05/14] Improve documentation --- ffcx/codegeneration/ufcx.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/ffcx/codegeneration/ufcx.h b/ffcx/codegeneration/ufcx.h index ae48e4716..170961b33 100644 --- a/ffcx/codegeneration/ufcx.h +++ b/ffcx/codegeneration/ufcx.h @@ -462,9 +462,13 @@ extern "C" /// Get an integral on sub domain subdomain_id ufcx_integral** (*integrals)(ufcx_integral_type); - /// New interface, static data, no function calls + /// List of cell, interior facet and exterior facet integrals ufcx_integral** form_integrals; + + /// IDs for each integral in form_integrals list int* form_integral_ids; + + /// Offsets for cell, interior facet and exterior facet integrals in form_integrals list int* form_integral_offsets; } ufcx_form; From 814b8565b5303a27ca71fb8a2b182cd94f5d6195 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Fri, 11 Aug 2023 07:29:47 +0100 Subject: [PATCH 06/14] Remove from form.py --- ffcx/codegeneration/dofmap_template.py | 18 ------------------ ffcx/codegeneration/ufcx.h | 12 ------------ 2 files changed, 30 deletions(-) diff --git a/ffcx/codegeneration/dofmap_template.py b/ffcx/codegeneration/dofmap_template.py index 0088d9a5f..abe5563f7 100644 --- a/ffcx/codegeneration/dofmap_template.py +++ b/ffcx/codegeneration/dofmap_template.py @@ -12,20 +12,6 @@ {sub_dofmaps_initialization} -void tabulate_entity_dofs_{factory_name}(int* restrict dofs, int d, int i) -{{ -{tabulate_entity_dofs} -}} - -void tabulate_entity_closure_dofs_{factory_name}(int* restrict dofs, int d, int i) -{{ -{tabulate_entity_closure_dofs} -}} - -{num_entity_dofs_init} - -{num_entity_closure_dofs_init} - {entity_dofs_init} {entity_dof_offsets_init} @@ -44,10 +30,6 @@ .entity_dof_offsets = {entity_dof_offsets}, .entity_closure_dofs = {entity_closure_dofs}, .entity_closure_dof_offsets = {entity_closure_dof_offsets}, - .num_entity_dofs = {num_entity_dofs}, - .tabulate_entity_dofs = tabulate_entity_dofs_{factory_name}, - .num_entity_closure_dofs = {num_entity_closure_dofs}, - .tabulate_entity_closure_dofs = tabulate_entity_closure_dofs_{factory_name}, .num_sub_dofmaps = {num_sub_dofmaps}, .sub_dofmaps = {sub_dofmaps} }}; diff --git a/ffcx/codegeneration/ufcx.h b/ffcx/codegeneration/ufcx.h index 170961b33..e70206d0b 100644 --- a/ffcx/codegeneration/ufcx.h +++ b/ffcx/codegeneration/ufcx.h @@ -230,18 +230,6 @@ extern "C" /// Offset for closure dofs of each entity in entity_closure_dofs int *entity_closure_dof_offsets; - /// Number of dofs associated with each cell entity of dimension d - int *num_entity_dofs; - - /// Tabulate the local-to-local mapping of dofs on entity (d, i) - void (*tabulate_entity_dofs)(int* restrict dofs, int d, int i); - - /// Number of dofs associated with the closure of each cell entity of dimension d - int *num_entity_closure_dofs; - - /// Tabulate the local-to-local mapping of dofs on the closure of entity (d, i) - void (*tabulate_entity_closure_dofs)(int* restrict dofs, int d, int i); - /// Number of sub dofmaps (for a mixed element) int num_sub_dofmaps; From fb9d2f03d0a1075f7ccb1add3bc9f4e36799a731 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Fri, 11 Aug 2023 07:51:24 +0100 Subject: [PATCH 07/14] Fix order --- ffcx/codegeneration/form.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ffcx/codegeneration/form.py b/ffcx/codegeneration/form.py index 6de4c46b9..14a806de1 100644 --- a/ffcx/codegeneration/form.py +++ b/ffcx/codegeneration/form.py @@ -95,7 +95,7 @@ def generator(ir, options): integrals = [] integral_ids = [] integral_offsets = [0] - for itg_type in ("cell", "interior_facet", "exterior_facet"): + for itg_type in ("cell", "exterior_facet", "interior_facet"): integrals += [L.AddressOf(L.Symbol(itg)) for itg in ir.integral_names[itg_type]] integral_ids += ir.subdomain_ids[itg_type] integral_offsets.append(len(integrals)) From 1f8a46425e639cf225e84b65a7391526348968ce Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Fri, 11 Aug 2023 09:10:41 +0100 Subject: [PATCH 08/14] Remove for integrals --- ffcx/codegeneration/form_template.py | 20 -------------------- ffcx/codegeneration/ufcx.h | 24 +++++++++++++++--------- 2 files changed, 15 insertions(+), 29 deletions(-) diff --git a/ffcx/codegeneration/form_template.py b/ffcx/codegeneration/form_template.py index 5c0b1d22c..02ad31d6f 100644 --- a/ffcx/codegeneration/form_template.py +++ b/ffcx/codegeneration/form_template.py @@ -40,21 +40,6 @@ {constant_name_map} }} -int* integral_ids_{factory_name}(ufcx_integral_type integral_type) -{{ -{integral_ids} -}} - -int num_integrals_{factory_name}(ufcx_integral_type integral_type) -{{ -{num_integrals} -}} - -ufcx_integral** integrals_{factory_name}(ufcx_integral_type integral_type) -{{ -{integrals} -}} - ufcx_form {factory_name} = {{ @@ -70,11 +55,6 @@ .finite_elements = {finite_elements}, .dofmaps = {dofmaps}, - .integral_ids = integral_ids_{factory_name}, - .num_integrals = num_integrals_{factory_name}, - - .integrals = integrals_{factory_name}, - .form_integrals = {form_integrals}, .form_integral_ids = {form_integral_ids}, .form_integral_offsets = form_integral_offsets_{factory_name} diff --git a/ffcx/codegeneration/ufcx.h b/ffcx/codegeneration/ufcx.h index e70206d0b..f3e3d1e00 100644 --- a/ffcx/codegeneration/ufcx.h +++ b/ffcx/codegeneration/ufcx.h @@ -200,6 +200,9 @@ extern "C" /// The highest degree of a polynomial in the element int highest_degree; + + /// The polyset type of the element + int polyset_type; } ufcx_basix_custom_finite_element; typedef struct ufcx_dofmap @@ -230,6 +233,18 @@ extern "C" /// Offset for closure dofs of each entity in entity_closure_dofs int *entity_closure_dof_offsets; + /// Number of dofs associated with each cell entity of dimension d + int *num_entity_dofs; + + /// Tabulate the local-to-local mapping of dofs on entity (d, i) + void (*tabulate_entity_dofs)(int* restrict dofs, int d, int i); + + /// Number of dofs associated with the closure of each cell entity of dimension d + int *num_entity_closure_dofs; + + /// Tabulate the local-to-local mapping of dofs on the closure of entity (d, i) + void (*tabulate_entity_closure_dofs)(int* restrict dofs, int d, int i); + /// Number of sub dofmaps (for a mixed element) int num_sub_dofmaps; @@ -441,15 +456,6 @@ extern "C" /// Coefficient number j=i-r if r+j <= i < r+n ufcx_dofmap** dofmaps; - /// All ids for integrals - int* (*integral_ids)(ufcx_integral_type); - - /// Number of integrals - int (*num_integrals)(ufcx_integral_type); - - /// Get an integral on sub domain subdomain_id - ufcx_integral** (*integrals)(ufcx_integral_type); - /// List of cell, interior facet and exterior facet integrals ufcx_integral** form_integrals; From 1fdecdcee08a7b396cecb663ae8173344fca76fb Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Fri, 11 Aug 2023 09:12:47 +0100 Subject: [PATCH 09/14] Remove all remaining cases --- ffcx/codegeneration/dofmap.py | 66 ----------------------------------- ffcx/codegeneration/form.py | 52 --------------------------- 2 files changed, 118 deletions(-) diff --git a/ffcx/codegeneration/dofmap.py b/ffcx/codegeneration/dofmap.py index 3ca6e3174..d6556e2be 100644 --- a/ffcx/codegeneration/dofmap.py +++ b/ffcx/codegeneration/dofmap.py @@ -15,47 +15,6 @@ logger = logging.getLogger("ffcx") -def tabulate_entity_dofs( - L, - entity_dofs: typing.List[typing.List[typing.List[int]]], - num_dofs_per_entity: typing.List[int], -): - # Output argument array - dofs = L.Symbol("dofs") - - # Input arguments - d = L.Symbol("d") - i = L.Symbol("i") - - # TODO: Removed check for (d <= tdim + 1) - tdim = len(num_dofs_per_entity) - 1 - - # Generate cases for each dimension: - all_cases = [] - for dim in range(tdim + 1): - # Ignore if no entities for this dimension - if num_dofs_per_entity[dim] == 0: - continue - - # Generate cases for each mesh entity - cases = [] - for entity in range(len(entity_dofs[dim])): - casebody = [] - for j, dof in enumerate(entity_dofs[dim][entity]): - casebody += [L.Assign(dofs[j], dof)] - cases.append((entity, L.StatementList(casebody))) - - # Generate inner switch - # TODO: Removed check for (i <= num_entities-1) - inner_switch = L.Switch(i, cases, autoscope=False) - all_cases.append((dim, inner_switch)) - - if all_cases: - return L.Switch(d, all_cases, autoscope=False) - else: - return L.NoOp() - - def generator(ir, options): """Generate UFC code for a dofmap.""" logger.info("Generating code for dofmap:") @@ -73,23 +32,6 @@ def generator(ir, options): import ffcx.codegeneration.C.cnodes as L - num_entity_dofs = ir.num_entity_dofs + [0, 0, 0, 0] - num_entity_dofs = num_entity_dofs[:4] - d["num_entity_dofs"] = f"num_entity_dofs_{ir.name}" - d["num_entity_dofs_init"] = L.ArrayDecl( - "int", f"num_entity_dofs_{ir.name}", values=num_entity_dofs, sizes=4 - ) - - num_entity_closure_dofs = ir.num_entity_closure_dofs + [0, 0, 0, 0] - num_entity_closure_dofs = num_entity_closure_dofs[:4] - d["num_entity_closure_dofs"] = f"num_entity_closure_dofs_{ir.name}" - d["num_entity_closure_dofs_init"] = L.ArrayDecl( - "int", - f"num_entity_closure_dofs_{ir.name}", - values=num_entity_closure_dofs, - sizes=4, - ) - flattened_entity_dofs = [] entity_dof_offsets = [0] for dim in ir.entity_dofs: @@ -137,14 +79,6 @@ def generator(ir, options): d["block_size"] = ir.block_size - # Functions - d["tabulate_entity_dofs"] = tabulate_entity_dofs( - L, ir.entity_dofs, ir.num_entity_dofs - ) - d["tabulate_entity_closure_dofs"] = tabulate_entity_dofs( - L, ir.entity_closure_dofs, ir.num_entity_closure_dofs - ) - if len(ir.sub_dofmaps) > 0: d["sub_dofmaps_initialization"] = L.ArrayDecl( "ufcx_dofmap*", diff --git a/ffcx/codegeneration/form.py b/ffcx/codegeneration/form.py index 6de4c46b9..4e0d1a73c 100644 --- a/ffcx/codegeneration/form.py +++ b/ffcx/codegeneration/form.py @@ -30,13 +30,6 @@ def generator(ir, options): d["num_coefficients"] = ir.num_coefficients d["num_constants"] = ir.num_constants - code = [] - cases = [] - for itg_type in ("cell", "interior_facet", "exterior_facet"): - cases += [(L.Symbol(itg_type), L.Return(len(ir.subdomain_ids[itg_type])))] - code += [L.Switch("integral_type", cases, default=L.Return(0))] - d["num_integrals"] = L.StatementList(code) - if len(ir.original_coefficient_position) > 0: d["original_coefficient_position_init"] = L.ArrayDecl( "int", @@ -128,51 +121,6 @@ def generator(ir, options): sizes=len(integral_offsets), ) - code = [] - cases = [] - code_ids = [] - cases_ids = [] - for itg_type in ("cell", "interior_facet", "exterior_facet"): - if len(ir.integral_names[itg_type]) > 0: - code += [ - L.ArrayDecl( - "static ufcx_integral*", - f"integrals_{itg_type}_{ir.name}", - values=[ - L.AddressOf(L.Symbol(itg)) - for itg in ir.integral_names[itg_type] - ], - sizes=len(ir.integral_names[itg_type]), - ) - ] - cases.append( - ( - L.Symbol(itg_type), - L.Return(L.Symbol(f"integrals_{itg_type}_{ir.name}")), - ) - ) - - code_ids += [ - L.ArrayDecl( - "static int", - f"integral_ids_{itg_type}_{ir.name}", - values=ir.subdomain_ids[itg_type], - sizes=len(ir.subdomain_ids[itg_type]), - ) - ] - cases_ids.append( - ( - L.Symbol(itg_type), - L.Return(L.Symbol(f"integral_ids_{itg_type}_{ir.name}")), - ) - ) - - code += [L.Switch("integral_type", cases, default=L.Return(L.Null()))] - code_ids += [L.Switch("integral_type", cases_ids, default=L.Return(L.Null()))] - d["integrals"] = L.StatementList(code) - - d["integral_ids"] = L.StatementList(code_ids) - code = [] function_name = L.Symbol("function_name") From 6f23e402f8bff2fd7210550d1b1712a1d3fd8dba Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Fri, 11 Aug 2023 09:15:18 +0100 Subject: [PATCH 10/14] Remove from Cnodes --- ffcx/codegeneration/C/cnodes.py | 274 ++++++++++++++++---------------- ffcx/codegeneration/dofmap.py | 1 - 2 files changed, 141 insertions(+), 134 deletions(-) diff --git a/ffcx/codegeneration/C/cnodes.py b/ffcx/codegeneration/C/cnodes.py index 860945923..f8b927763 100644 --- a/ffcx/codegeneration/C/cnodes.py +++ b/ffcx/codegeneration/C/cnodes.py @@ -10,8 +10,7 @@ import numpy as np from ffcx.codegeneration.C.format_lines import Indented, format_indented_lines -from ffcx.codegeneration.C.format_value import (format_float, format_int, - format_value) +from ffcx.codegeneration.C.format_value import format_float, format_int, format_value from ffcx.codegeneration.C.precedence import PRECEDENCE logger = logging.getLogger("ffcx") @@ -32,18 +31,21 @@ def is_zero_cexpr(cexpr): - return ((isinstance(cexpr, LiteralFloat) and cexpr.value == 0.0) - or (isinstance(cexpr, LiteralInt) and cexpr.value == 0)) + return (isinstance(cexpr, LiteralFloat) and cexpr.value == 0.0) or ( + isinstance(cexpr, LiteralInt) and cexpr.value == 0 + ) def is_one_cexpr(cexpr): - return ((isinstance(cexpr, LiteralFloat) and cexpr.value == 1.0) - or (isinstance(cexpr, LiteralInt) and cexpr.value == 1)) + return (isinstance(cexpr, LiteralFloat) and cexpr.value == 1.0) or ( + isinstance(cexpr, LiteralInt) and cexpr.value == 1 + ) def is_negative_one_cexpr(cexpr): - return ((isinstance(cexpr, LiteralFloat) and cexpr.value == -1.0) - or (isinstance(cexpr, LiteralInt) and cexpr.value == -1)) + return (isinstance(cexpr, LiteralFloat) and cexpr.value == -1.0) or ( + isinstance(cexpr, LiteralInt) and cexpr.value == -1 + ) def float_product(factors): @@ -58,6 +60,8 @@ def float_product(factors): if is_zero_cexpr(f): return f return Product(factors) + + # CNode core @@ -77,6 +81,7 @@ def __eq__(self, other): def __ne__(self, other): return not self.__eq__(other) + # CExpr base classes @@ -261,7 +266,7 @@ def __eq__(self, other): class LiteralFloat(CExprLiteral): """A floating point literal value.""" - __slots__ = ("value", ) + __slots__ = ("value",) precedence = PRECEDENCE.LITERAL def __init__(self, value): @@ -289,7 +294,7 @@ def flops(self): class LiteralInt(CExprLiteral): """An integer literal value.""" - __slots__ = ("value", ) + __slots__ = ("value",) precedence = PRECEDENCE.LITERAL def __init__(self, value): @@ -323,11 +328,11 @@ def __hash__(self): class LiteralBool(CExprLiteral): """A boolean literal value.""" - __slots__ = ("value", ) + __slots__ = ("value",) precedence = PRECEDENCE.LITERAL def __init__(self, value): - assert isinstance(value, (bool, )) + assert isinstance(value, (bool,)) self.value = value def ce_format(self, precision=None): @@ -345,16 +350,16 @@ def __bool__(self): class LiteralString(CExprLiteral): """A boolean literal value.""" - __slots__ = ("value", ) + __slots__ = ("value",) precedence = PRECEDENCE.LITERAL def __init__(self, value): - assert isinstance(value, (str, )) + assert isinstance(value, (str,)) assert '"' not in value self.value = value def ce_format(self, precision=None): - return '"%s"' % (self.value, ) + return '"%s"' % (self.value,) def __eq__(self, other): return isinstance(other, LiteralString) and self.value == other.value @@ -363,7 +368,7 @@ def __eq__(self, other): class Symbol(CExprTerminal): """A named symbol.""" - __slots__ = ("name", ) + __slots__ = ("name",) precedence = PRECEDENCE.SYMBOL def __init__(self, name): @@ -385,10 +390,11 @@ def __hash__(self): # CExprOperator base classes + class UnaryOp(CExprOperator): """Base class for unary operators.""" - __slots__ = ("arg", ) + __slots__ = ("arg",) def __init__(self, arg): self.arg = as_cexpr(arg) @@ -408,7 +414,7 @@ class PrefixUnaryOp(UnaryOp): def ce_format(self, precision=None): arg = self.arg.ce_format(precision) if self.arg.precedence >= self.precedence: - arg = '(' + arg + ')' + arg = "(" + arg + ")" return self.op + arg def __eq__(self, other): @@ -423,7 +429,7 @@ class PostfixUnaryOp(UnaryOp): def ce_format(self, precision=None): arg = self.arg.ce_format(precision) if self.arg.precedence >= self.precedence: - arg = '(' + arg + ')' + arg = "(" + arg + ")" return arg + self.op def __eq__(self, other): @@ -444,15 +450,19 @@ def ce_format(self, precision=None): # Apply parentheses if self.lhs.precedence >= self.precedence: - lhs = '(' + lhs + ')' + lhs = "(" + lhs + ")" if self.rhs.precedence >= self.precedence: - rhs = '(' + rhs + ')' + rhs = "(" + rhs + ")" # Return combined string return lhs + (" " + self.op + " ") + rhs def __eq__(self, other): - return (isinstance(other, type(self)) and self.lhs == other.lhs and self.rhs == other.rhs) + return ( + isinstance(other, type(self)) + and self.lhs == other.lhs + and self.rhs == other.rhs + ) def __hash__(self): return hash(self.ce_format()) @@ -464,7 +474,7 @@ def flops(self): class NaryOp(CExprOperator): """Base class for special n-ary operators.""" - __slots__ = ("args", ) + __slots__ = ("args",) def __init__(self, args): self.args = [as_cexpr(arg) for arg in args] @@ -476,7 +486,7 @@ def ce_format(self, precision=None): # Apply parentheses for i in range(len(args)): if self.args[i].precedence >= self.precedence: - args[i] = '(' + args[i] + ')' + args[i] = "(" + args[i] + ")" # Return combined string op = " " + self.op + " " @@ -486,8 +496,11 @@ def ce_format(self, precision=None): return s def __eq__(self, other): - return (isinstance(other, type(self)) and len(self.args) == len(other.args) - and all(a == b for a, b in zip(self.args, other.args))) + return ( + isinstance(other, type(self)) + and len(self.args) == len(other.args) + and all(a == b for a, b in zip(self.args, other.args)) + ) def flops(self): flops = len(self.args) - 1 @@ -716,6 +729,7 @@ class AssignDiv(AssignOp): __slots__ = () op = "/=" + # CExpr operators @@ -763,7 +777,7 @@ def __init__(self, array, dummy=None, dims=None, strides=None, offset=None): def __getitem__(self, indices): if not isinstance(indices, (list, tuple)): - indices = (indices, ) + indices = (indices,) n = len(indices) if n == 0: # Handle scalar case, allowing dims=() and indices=() for A[0] @@ -773,7 +787,7 @@ def __getitem__(self, indices): else: i, s = (indices[0], self.strides[0]) literal_one = LiteralInt(1) - flat = (i if s == literal_one else s * i) + flat = i if s == literal_one else s * i if self.offset is not None: flat = self.offset + flat for i, s in zip(indices[1:n], self.strides[1:n]): @@ -796,11 +810,11 @@ def __init__(self, array, indices): elif isinstance(array, ArrayDecl): self.array = array.symbol else: - raise ValueError("Unexpected array type %s." % (type(array).__name__, )) + raise ValueError("Unexpected array type %s." % (type(array).__name__,)) # Allow expressions or literals as indices if not isinstance(indices, (list, tuple)): - indices = (indices, ) + indices = (indices,) self.indices = tuple(as_cexpr_or_string_symbol(i) for i in indices) # Early error checking for negative array dimensions @@ -812,8 +826,10 @@ def __init__(self, array, indices): if len(self.indices) != len(array.sizes): raise ValueError("Invalid number of indices.") ints = (int, LiteralInt) - if any((isinstance(i, ints) and isinstance(d, ints) and int(i) >= int(d)) - for i, d in zip(self.indices, array.sizes)): + if any( + (isinstance(i, ints) and isinstance(d, ints) and int(i) >= int(d)) + for i, d in zip(self.indices, array.sizes) + ): raise ValueError("Index value >= array dimension.") def __getitem__(self, indices): @@ -821,7 +837,7 @@ def __getitem__(self, indices): if isinstance(indices, list): indices = tuple(indices) elif not isinstance(indices, tuple): - indices = (indices, ) + indices = (indices,) return ArrayAccess(self.array, self.indices + indices) def ce_format(self, precision=None): @@ -831,8 +847,11 @@ def ce_format(self, precision=None): return s def __eq__(self, other): - return (isinstance(other, type(self)) and self.array == other.array - and self.indices == other.indices) + return ( + isinstance(other, type(self)) + and self.array == other.array + and self.indices == other.indices + ) def __hash__(self): return hash(self.ce_format()) @@ -858,18 +877,22 @@ def ce_format(self, precision=None): # Apply parentheses if self.condition.precedence >= self.precedence: - c = '(' + c + ')' + c = "(" + c + ")" if self.true.precedence >= self.precedence: - t = '(' + t + ')' + t = "(" + t + ")" if self.false.precedence >= self.precedence: - f = '(' + f + ')' + f = "(" + f + ")" # Return combined string return c + " ? " + t + " : " + f def __eq__(self, other): - return (isinstance(other, type(self)) and self.condition == other.condition - and self.true == other.true and self.false == other.false) + return ( + isinstance(other, type(self)) + and self.condition == other.condition + and self.true == other.true + and self.false == other.false + ) def flops(self): raise NotImplementedError("Flop count is not implemented for conditionals") @@ -887,7 +910,7 @@ def __init__(self, function, arguments=None): if arguments is None: arguments = () elif not isinstance(arguments, (tuple, list)): - arguments = (arguments, ) + arguments = (arguments,) self.arguments = [as_cexpr(arg) for arg in arguments] def ce_format(self, precision=None): @@ -895,8 +918,11 @@ def ce_format(self, precision=None): return self.function.ce_format(precision) + "(" + args + ")" def __eq__(self, other): - return (isinstance(other, type(self)) and self.function == other.function - and self.arguments == other.arguments) + return ( + isinstance(other, type(self)) + and self.function == other.function + and self.arguments == other.arguments + ) def flops(self): return 1 @@ -934,7 +960,7 @@ def as_cexpr(node): elif isinstance(node, numbers.Real): return LiteralFloat(node) elif isinstance(node, str): - raise RuntimeError("Got string for CExpr, this is ambiguous: %s" % (node, )) + raise RuntimeError("Got string for CExpr, this is ambiguous: %s" % (node,)) else: raise RuntimeError("Unexpected CExpr type %s:\n%s" % (type(node), str(node))) @@ -1005,7 +1031,9 @@ class CStatement(CNode): def cs_format(self, precision=None): """Return S: string | list(S) | Indented(S).""" - raise NotImplementedError("Missing implementation of cs_format() in CStatement.") + raise NotImplementedError( + "Missing implementation of cs_format() in CStatement." + ) def __str__(self): try: @@ -1025,7 +1053,7 @@ def flops(self): class VerbatimStatement(CStatement): """Wraps a source code string to be pasted verbatim into the source code.""" - __slots__ = ("codestring", ) + __slots__ = ("codestring",) is_scoped = False def __init__(self, codestring): @@ -1036,13 +1064,13 @@ def cs_format(self, precision=None): return self.codestring def __eq__(self, other): - return (isinstance(other, type(self)) and self.codestring == other.codestring) + return isinstance(other, type(self)) and self.codestring == other.codestring class Statement(CStatement): """Make an expression into a statement.""" - __slots__ = ("expr", ) + __slots__ = ("expr",) is_scoped = False def __init__(self, expr): @@ -1052,7 +1080,7 @@ def cs_format(self, precision=None): return self.expr.ce_format(precision) + ";" def __eq__(self, other): - return (isinstance(other, type(self)) and self.expr == other.expr) + return isinstance(other, type(self)) and self.expr == other.expr def flops(self): # print(self.expr.rhs.flops()) @@ -1062,7 +1090,7 @@ def flops(self): class StatementList(CStatement): """A simple sequence of statements. No new scopes are introduced.""" - __slots__ = ("statements", ) + __slots__ = ("statements",) def __init__(self, statements): self.statements = [as_cstatement(st) for st in statements] @@ -1075,7 +1103,7 @@ def cs_format(self, precision=None): return [st.cs_format(precision) for st in self.statements] def __eq__(self, other): - return (isinstance(other, type(self)) and self.statements == other.statements) + return isinstance(other, type(self)) and self.statements == other.statements def flops(self): flops = 0 @@ -1116,7 +1144,7 @@ def flops(self): class Return(CStatement): - __slots__ = ("value", ) + __slots__ = ("value",) is_scoped = True def __init__(self, value=None): @@ -1129,10 +1157,10 @@ def cs_format(self, precision=None): if self.value is None: return "return;" else: - return "return %s;" % (self.value.ce_format(precision), ) + return "return %s;" % (self.value.ce_format(precision),) def __eq__(self, other): - return (isinstance(other, type(self)) and self.value == other.value) + return isinstance(other, type(self)) and self.value == other.value def flops(self): return 0 @@ -1141,7 +1169,7 @@ def flops(self): class Comment(CStatement): """Line comment(s) used for annotating the generated code with human readable remarks.""" - __slots__ = ("comment", ) + __slots__ = ("comment",) is_scoped = True def __init__(self, comment): @@ -1153,7 +1181,7 @@ def cs_format(self, precision=None): return ["// " + line.strip() for line in lines] def __eq__(self, other): - return (isinstance(other, type(self)) and self.comment == other.comment) + return isinstance(other, type(self)) and self.comment == other.comment def flops(self): return 0 @@ -1179,7 +1207,7 @@ def commented_code_list(code, comments): class Pragma(CStatement): """Pragma comments used for compiler-specific annotations.""" - __slots__ = ("comment", ) + __slots__ = ("comment",) is_scoped = True def __init__(self, comment): @@ -1191,7 +1219,7 @@ def cs_format(self, precision=None): return "#pragma " + self.comment def __eq__(self, other): - return (isinstance(other, type(self)) and self.comment == other.comment) + return isinstance(other, type(self)) and self.comment == other.comment def flops(self): return 0 @@ -1207,7 +1235,6 @@ class VariableDecl(CStatement): is_scoped = False def __init__(self, typename, symbol, value=None): - # No type system yet, just using strings assert isinstance(typename, str) self.typename = typename @@ -1226,8 +1253,12 @@ def cs_format(self, precision=None): return code + ";" def __eq__(self, other): - return (isinstance(other, type(self)) and self.typename == other.typename - and self.symbol == other.symbol and self.value == other.value) + return ( + isinstance(other, type(self)) + and self.typename == other.typename + and self.symbol == other.symbol + and self.value == other.value + ) def flops(self): if self.value is not None: @@ -1305,13 +1336,18 @@ def formatter(x, p): r = len(sizes) assert r > 0 if r == 1: - return [build_1d_initializer_list(values, formatter, padlen=padlen, precision=precision)] + return [ + build_1d_initializer_list( + values, formatter, padlen=padlen, precision=precision + ) + ] else: # Render all sublists parts = [] for val in values: sublist = build_initializer_lists( - val, sizes[1:], level + 1, formatter, padlen=padlen, precision=precision) + val, sizes[1:], level + 1, formatter, padlen=padlen, precision=precision + ) parts.append(sublist) # Add comma after last line in each part except the last one for part in parts[:-1]: @@ -1357,7 +1393,7 @@ def __init__(self, typename, symbol, sizes=None, values=None, padlen=0): self.symbol = as_symbol(symbol) if isinstance(sizes, int): - sizes = (sizes, ) + sizes = (sizes,) self.sizes = tuple(sizes) # NB! No type checking, assuming nested lists of literal values. Not applying as_cexpr. @@ -1370,13 +1406,15 @@ def __init__(self, typename, symbol, sizes=None, values=None, padlen=0): def cs_format(self, precision=None): if not all(self.sizes): - raise RuntimeError(f"Detected an array {self.symbol} dimension of zero. This is not valid in C.") + raise RuntimeError( + f"Detected an array {self.symbol} dimension of zero. This is not valid in C." + ) # Pad innermost array dimension sizes = pad_innermost_dim(self.sizes, self.padlen) # Add brackets - brackets = ''.join("[%d]" % n for n in sizes) + brackets = "".join("[%d]" % n for n in sizes) # Join declaration decl = self.typename + " " + self.symbol.name + brackets @@ -1398,13 +1436,21 @@ def cs_format(self, precision=None): elif self.values.dtype.kind == "i": formatter = format_int elif self.values.dtype == np.bool_: + def format_bool(x, precision=None): return "true" if x is True else "false" + formatter = format_bool else: formatter = format_value initializer_lists = build_initializer_lists( - self.values, self.sizes, 0, formatter, padlen=self.padlen, precision=precision) + self.values, + self.sizes, + 0, + formatter, + padlen=self.padlen, + precision=precision, + ) if len(initializer_lists) == 1: return decl + " = " + initializer_lists[0] + ";" else: @@ -1413,8 +1459,9 @@ def format_bool(x, precision=None): def __eq__(self, other): attributes = ("typename", "symbol", "sizes", "padlen", "values") - return (isinstance(other, type(self)) - and all(getattr(self, name) == getattr(self, name) for name in attributes)) + return isinstance(other, type(self)) and all( + getattr(self, name) == getattr(self, name) for name in attributes + ) def flops(self): return 0 @@ -1424,7 +1471,7 @@ def flops(self): class Scope(CStatement): - __slots__ = ("body", ) + __slots__ = ("body",) is_scoped = True def __init__(self, body): @@ -1434,7 +1481,7 @@ def cs_format(self, precision=None): return ("{", Indented(self.body.cs_format(precision)), "}") def __eq__(self, other): - return (isinstance(other, type(self)) and self.body == other.body) + return isinstance(other, type(self)) and self.body == other.body def flops(self): return 0 @@ -1444,7 +1491,7 @@ def _is_simple_if_body(body): if isinstance(body, StatementList): if len(body.statements) > 1: return False - body, = body.statements + (body,) = body.statements return isinstance(body, (Return, AssignOp, Break, Continue)) @@ -1465,8 +1512,11 @@ def cs_format(self, precision=None): return (statement, "{", body_fmt, "}") def __eq__(self, other): - return (isinstance(other, type(self)) and self.condition == other.condition - and self.body == other.body) + return ( + isinstance(other, type(self)) + and self.condition == other.condition + and self.body == other.body + ) class ElseIf(CStatement): @@ -1486,12 +1536,15 @@ def cs_format(self, precision=None): return (statement, "{", body_fmt, "}") def __eq__(self, other): - return (isinstance(other, type(self)) and self.condition == other.condition - and self.body == other.body) + return ( + isinstance(other, type(self)) + and self.condition == other.condition + and self.body == other.body + ) class Else(CStatement): - __slots__ = ("body", ) + __slots__ = ("body",) is_scoped = True def __init__(self, body): @@ -1506,7 +1559,7 @@ def cs_format(self, precision=None): return (statement, "{", body_fmt, "}") def __eq__(self, other): - return (isinstance(other, type(self)) and self.body == other.body) + return isinstance(other, type(self)) and self.body == other.body def is_simple_inner_loop(code): @@ -1517,56 +1570,6 @@ def is_simple_inner_loop(code): return False -class Switch(CStatement): - __slots__ = ("arg", "cases", "default", "autobreak", "autoscope") - is_scoped = True - - def __init__(self, arg, cases, default=None, autobreak=True, autoscope=True): - self.arg = as_cexpr_or_string_symbol(arg) - self.cases = [(as_cexpr(value), as_cstatement(body)) for value, body in cases] - if default is not None: - default = as_cstatement(default) - defcase = [(None, default)] - else: - defcase = [] - self.default = default - # If this is a switch where every case returns, scopes or breaks are never needed - if all(isinstance(case[1], Return) for case in self.cases + defcase): - autobreak = False - autoscope = False - if all(case[1].is_scoped for case in self.cases + defcase): - autoscope = False - assert autobreak in (True, False) - assert autoscope in (True, False) - self.autobreak = autobreak - self.autoscope = autoscope - - def cs_format(self, precision=None): - cases = [] - for case in self.cases: - caseheader = "case " + case[0].ce_format(precision) + ":" - casebody = case[1].cs_format(precision) - if self.autoscope: - casebody = ("{", Indented(casebody), "}") - if self.autobreak: - casebody = (casebody, "break;") - cases.extend([caseheader, Indented(casebody)]) - - if self.default is not None: - caseheader = "default:" - casebody = self.default.cs_format(precision) - if self.autoscope: - casebody = ("{", Indented(casebody), "}") - cases.extend([caseheader, Indented(casebody)]) - - return ("switch (" + self.arg.ce_format(precision) + ")", "{", cases, "}") - - def __eq__(self, other): - attributes = ("arg", "cases", "default", "autobreak", "autoscope") - return (isinstance(other, type(self)) - and all(getattr(self, name) == getattr(self, name) for name in attributes)) - - class ForRange(CStatement): """Slightly higher-level for loop assuming incrementing an index over a range.""" @@ -1603,8 +1606,9 @@ def cs_format(self, precision=None): def __eq__(self, other): attributes = ("index", "begin", "end", "body", "index_type") - return (isinstance(other, type(self)) - and all(getattr(self, name) == getattr(self, name) for name in attributes)) + return isinstance(other, type(self)) and all( + getattr(self, name) == getattr(self, name) for name in attributes + ) def flops(self): return (self.end.value - self.begin.value) * self.body.flops() @@ -1626,8 +1630,10 @@ def as_cstatement(node): # Special case for using assignment expressions as statements return Statement(node) else: - raise RuntimeError("Trying to create a statement of CExprOperator type %s:\n%s" % - (type(node), str(node))) + raise RuntimeError( + "Trying to create a statement of CExprOperator type %s:\n%s" + % (type(node), str(node)) + ) elif isinstance(node, list): # Convenience case for list of statements if len(node) == 1: @@ -1639,4 +1645,6 @@ def as_cstatement(node): # Backdoor for flexibility in code generation to allow verbatim pasted statements return VerbatimStatement(node) else: - raise RuntimeError("Unexpected CStatement type %s:\n%s" % (type(node), str(node))) + raise RuntimeError( + "Unexpected CStatement type %s:\n%s" % (type(node), str(node)) + ) diff --git a/ffcx/codegeneration/dofmap.py b/ffcx/codegeneration/dofmap.py index d6556e2be..de3d21dc4 100644 --- a/ffcx/codegeneration/dofmap.py +++ b/ffcx/codegeneration/dofmap.py @@ -8,7 +8,6 @@ # old implementation in FFC import logging -import typing import ffcx.codegeneration.dofmap_template as ufcx_dofmap From 4618c33a592f1c8fd6b8fb83ee5ec1eb42088691 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Fri, 11 Aug 2023 11:03:28 +0100 Subject: [PATCH 11/14] Fix tests --- test/test_add_mode.py | 95 ++-- test/test_blocked_elements.py | 65 +-- test/test_jit_forms.py | 892 +++++++++++++++++++++++----------- 3 files changed, 683 insertions(+), 369 deletions(-) diff --git a/test/test_add_mode.py b/test/test_add_mode.py index 7ee4aba31..778b81a00 100644 --- a/test/test_add_mode.py +++ b/test/test_add_mode.py @@ -13,21 +13,17 @@ from ffcx.naming import cdtype_to_numpy, scalar_to_value_type -@pytest.mark.parametrize("mode", - [ - "double", - "float", - "long double", - "double _Complex", - "float _Complex" - ]) +@pytest.mark.parametrize( + "mode", ["double", "float", "long double", "double _Complex", "float _Complex"] +) def test_additive_facet_integral(mode, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) a = ufl.inner(u, v) * ufl.ds forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -35,11 +31,14 @@ def test_additive_facet_integral(mode, compile_args): ffi = module.ffi form0 = compiled_forms[0] - assert form0.num_integrals(module.lib.exterior_facet) == 1 - ids = form0.integral_ids(module.lib.exterior_facet) - assert ids[0] == -1 + integral_offsets = form0.form_integral_offsets + ex = module.lib.exterior_facet + num_integrals = integral_offsets[ex + 1] - integral_offsets[ex] + assert num_integrals == 1 + integral_id = form0.form_integral_ids[integral_offsets[ex]] + assert integral_id == -1 - default_integral = form0.integrals(module.lib.exterior_facet)[0] + default_integral = form0.form_integrals[integral_offsets[ex]] np_type = cdtype_to_numpy(mode) A = np.zeros((3, 3), dtype=np_type) @@ -50,32 +49,38 @@ def test_additive_facet_integral(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([0.0, 2.0, 0.0, - np.sqrt(3.0), -1.0, 0.0, - -np.sqrt(3.0), -1.0, 0.0], dtype=np_gtype) + coords = np.array( + [0.0, 2.0, 0.0, np.sqrt(3.0), -1.0, 0.0, -np.sqrt(3.0), -1.0, 0.0], + dtype=np_gtype, + ) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") for i in range(3): facets[0] = i - kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), - ffi.cast('int *', facets.ctypes.data), - ffi.cast('uint8_t *', perm.ctypes.data)) + kernel( + ffi.cast("{type} *".format(type=mode), A.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.cast("int *", facets.ctypes.data), + ffi.cast("uint8_t *", perm.ctypes.data), + ) assert np.isclose(A.sum(), np.sqrt(12) * (i + 1)) -@pytest.mark.parametrize("mode", ["double", "float", "long double", "double _Complex", "float _Complex"]) +@pytest.mark.parametrize( + "mode", ["double", "float", "long double", "double _Complex", "float _Complex"] +) def test_additive_cell_integral(mode, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -83,11 +88,14 @@ def test_additive_cell_integral(mode, compile_args): ffi = module.ffi form0 = compiled_forms[0] - assert form0.num_integrals(module.lib.cell) == 1 - ids = form0.integral_ids(module.lib.cell) - assert ids[0] == -1 + offsets = form0.form_integral_offsets + cell = module.lib.cell + num_integrals = offsets[cell + 1] - offsets[cell] + assert num_integrals == 1 + integral_id = form0.form_integral_ids[offsets[cell]] + assert integral_id == -1 - default_integral = form0.integrals(0)[0] + default_integral = form0.form_integrals[offsets[cell]] np_type = cdtype_to_numpy(mode) A = np.zeros((3, 3), dtype=np_type) @@ -96,22 +104,31 @@ def test_additive_cell_integral(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([0.0, 2.0, 0.0, - np.sqrt(3.0), -1.0, 0.0, - -np.sqrt(3.0), -1.0, 0.0], dtype=np_gtype) + coords = np.array( + [0.0, 2.0, 0.0, np.sqrt(3.0), -1.0, 0.0, -np.sqrt(3.0), -1.0, 0.0], + dtype=np_gtype, + ) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("{type} *".format(type=mode), A.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) A0 = np.array(A) for i in range(3): - kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("{type} *".format(type=mode), A.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.all(np.isclose(A, (i + 2) * A0)) diff --git a/test/test_blocked_elements.py b/test/test_blocked_elements.py index fd4db7b48..3388a8638 100644 --- a/test/test_blocked_elements.py +++ b/test/test_blocked_elements.py @@ -15,7 +15,8 @@ def test_finite_element(compile_args): ufl_element = basix.ufl.element("Lagrange", "triangle", 1) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args) + [ufl_element], cffi_extra_compile_args=compile_args + ) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 2 @@ -32,22 +33,18 @@ def test_finite_element(compile_args): assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_element_support_dofs == 3 - assert ufcx_dofmap.num_entity_dofs[0] == 1 - assert ufcx_dofmap.num_entity_dofs[1] == 0 - assert ufcx_dofmap.num_entity_dofs[2] == 0 - assert ufcx_dofmap.num_entity_dofs[3] == 0 + off = np.array([ufcx_dofmap.entity_dof_offsets[i] for i in range(8)]) + assert np.all(np.diff(off) == [1, 1, 1, 0, 0, 0, 0]) for v in range(3): - vals = np.zeros(1, dtype=np.int32) - vals_ptr = module.ffi.cast("int *", module.ffi.from_buffer(vals)) - ufcx_dofmap.tabulate_entity_dofs(vals_ptr, 0, v) - assert vals[0] == v + assert ufcx_dofmap.entity_dofs[v] == v assert ufcx_dofmap.num_sub_dofmaps == 0 def test_vector_element(compile_args): ufl_element = basix.ufl.element("Lagrange", "triangle", 1, rank=1) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args) + [ufl_element], cffi_extra_compile_args=compile_args + ) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 2 @@ -66,22 +63,19 @@ def test_vector_element(compile_args): assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_element_support_dofs == 3 - assert ufcx_dofmap.num_entity_dofs[0] == 1 - assert ufcx_dofmap.num_entity_dofs[1] == 0 - assert ufcx_dofmap.num_entity_dofs[2] == 0 - assert ufcx_dofmap.num_entity_dofs[3] == 0 + off = np.array([ufcx_dofmap.entity_dof_offsets[i] for i in range(8)]) + assert np.all(np.diff(off) == [1, 1, 1, 0, 0, 0, 0]) + for v in range(3): - vals = np.zeros(1, dtype=np.int32) - vals_ptr = module.ffi.cast("int *", module.ffi.from_buffer(vals)) - ufcx_dofmap.tabulate_entity_dofs(vals_ptr, 0, v) - assert vals[0] == v + assert ufcx_dofmap.entity_dofs[v] == v assert ufcx_dofmap.num_sub_dofmaps == 2 def test_tensor_element(compile_args): ufl_element = basix.ufl.element("Lagrange", "triangle", 1, rank=2) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args) + [ufl_element], cffi_extra_compile_args=compile_args + ) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 2 @@ -102,22 +96,21 @@ def test_tensor_element(compile_args): assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_element_support_dofs == 3 - assert ufcx_dofmap.num_entity_dofs[0] == 1 - assert ufcx_dofmap.num_entity_dofs[1] == 0 - assert ufcx_dofmap.num_entity_dofs[2] == 0 - assert ufcx_dofmap.num_entity_dofs[3] == 0 + off = np.array([ufcx_dofmap.entity_dof_offsets[i] for i in range(8)]) + assert np.all(np.diff(off) == [1, 1, 1, 0, 0, 0, 0]) + for v in range(3): - vals = np.zeros(1, dtype=np.int32) - vals_ptr = module.ffi.cast("int *", module.ffi.from_buffer(vals)) - ufcx_dofmap.tabulate_entity_dofs(vals_ptr, 0, v) - assert vals[0] == v + assert ufcx_dofmap.entity_dofs[v] == v assert ufcx_dofmap.num_sub_dofmaps == 4 def test_vector_quadrature_element(compile_args): - ufl_element = ufl.VectorElement(ufl.FiniteElement("Quadrature", "tetrahedron", degree=2, quad_scheme="default")) + ufl_element = ufl.VectorElement( + ufl.FiniteElement("Quadrature", "tetrahedron", degree=2, quad_scheme="default") + ) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args) + [ufl_element], cffi_extra_compile_args=compile_args + ) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 3 @@ -136,14 +129,10 @@ def test_vector_quadrature_element(compile_args): assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_global_support_dofs == 0 assert ufcx_dofmap.num_element_support_dofs == 4 - assert ufcx_dofmap.num_entity_dofs[0] == 0 - assert ufcx_dofmap.num_entity_dofs[1] == 0 - assert ufcx_dofmap.num_entity_dofs[2] == 0 - assert ufcx_dofmap.num_entity_dofs[3] == 4 - - vals = np.zeros(4, dtype=np.int32) - vals_ptr = module.ffi.cast("int *", module.ffi.from_buffer(vals)) - ufcx_dofmap.tabulate_entity_dofs(vals_ptr, 3, 0) - assert (vals == [0, 1, 2, 3]).all() + off = np.array([ufcx_dofmap.entity_dof_offsets[i] for i in range(16)]) + assert np.all(np.diff(off) == [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4]) + + for i in range(4): + assert ufcx_dofmap.entity_dofs[i] == i assert ufcx_dofmap.num_sub_dofmaps == 3 diff --git a/test/test_jit_forms.py b/test/test_jit_forms.py index b6fe9518b..8f67ec4ba 100644 --- a/test/test_jit_forms.py +++ b/test/test_jit_forms.py @@ -15,14 +15,29 @@ from ffcx.naming import cdtype_to_numpy, scalar_to_value_type -@pytest.mark.parametrize("mode,expected_result", [ - ("double", np.array([[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], dtype=np.float64)), - ("double _Complex", - np.array( - [[1.0 + 0j, -0.5 + 0j, -0.5 + 0j], [-0.5 + 0j, 0.5 + 0j, 0.0 + 0j], - [-0.5 + 0j, 0.0 + 0j, 0.5 + 0j]], - dtype=np.complex128)), -]) +@pytest.mark.parametrize( + "mode,expected_result", + [ + ( + "double", + np.array( + [[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], + dtype=np.float64, + ), + ), + ( + "double _Complex", + np.array( + [ + [1.0 + 0j, -0.5 + 0j, -0.5 + 0j], + [-0.5 + 0j, 0.5 + 0j, 0.0 + 0j], + [-0.5 + 0j, 0.0 + 0j, 0.5 + 0j], + ], + dtype=np.complex128, + ), + ), + ], +) def test_laplace_bilinear_form_2d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) kappa = ufl.Constant(ufl.triangle, shape=(2, 2)) @@ -31,7 +46,8 @@ def test_laplace_bilinear_form_2d(mode, expected_result, compile_args): a = ufl.tr(kappa) * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -39,11 +55,13 @@ def test_laplace_bilinear_form_2d(mode, expected_result, compile_args): ffi = module.ffi form0 = compiled_forms[0] - assert form0.num_integrals(module.lib.cell) == 1 - ids = form0.integral_ids(module.lib.cell) - assert ids[0] == -1 + offsets = form0.form_integral_offsets + cell = module.lib.cell + assert offsets[cell + 1] - offsets[cell] == 1 + integral_id = form0.form_integral_ids[offsets[cell]] + assert integral_id == -1 - default_integral = form0.integrals(module.lib.cell)[0] + default_integral = form0.form_integrals[offsets[cell]] np_type = cdtype_to_numpy(mode) A = np.zeros((3, 3), dtype=np_type) @@ -54,47 +72,84 @@ def test_laplace_bilinear_form_2d(mode, expected_result, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([[0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np_gtype) + coords = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype + ) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("{type} *".format(type=mode), A.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.allclose(A, np.trace(kappa_value) * expected_result) -@pytest.mark.parametrize("mode,expected_result", [ - ("float", - np.array( - [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], - dtype=np.float32)), - ("long double", - np.array( - [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], - dtype=np.longdouble)), - ("double", - np.array( - [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], - dtype=np.float64)), - ("double _Complex", - np.array( - [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], - dtype=np.complex128)), - ("float _Complex", - np.array( - [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], - dtype=np.complex64)), -]) +@pytest.mark.parametrize( + "mode,expected_result", + [ + ( + "float", + np.array( + [ + [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], + ], + dtype=np.float32, + ), + ), + ( + "long double", + np.array( + [ + [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], + ], + dtype=np.longdouble, + ), + ), + ( + "double", + np.array( + [ + [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], + ], + dtype=np.float64, + ), + ), + ( + "double _Complex", + np.array( + [ + [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], + ], + dtype=np.complex128, + ), + ), + ( + "float _Complex", + np.array( + [ + [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], + ], + dtype=np.complex64, + ), + ), + ], +) def test_mass_bilinear_form_2d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) @@ -102,13 +157,14 @@ def test_mass_bilinear_form_2d(mode, expected_result, compile_args): L = ufl.conj(v) * ufl.dx forms = [a, L] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) - form0 = compiled_forms[0].integrals(module.lib.cell)[0] - form1 = compiled_forms[1].integrals(module.lib.cell)[0] + form0 = compiled_forms[0].form_integrals[0] + form1 = compiled_forms[1].form_integrals[0] np_type = cdtype_to_numpy(mode) A = np.zeros((3, 3), dtype=np_type) @@ -119,34 +175,64 @@ def test_mass_bilinear_form_2d(mode, expected_result, compile_args): np_gtype = cdtype_to_numpy(geom_type) ffi = module.ffi - coords = np.array([[0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np_gtype) + coords = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype + ) - kernel0 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form0, f"tabulate_tensor_{np_type}")) - kernel0(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel0 = ffi.cast( + f"ufcx_tabulate_tensor_{np_type} *", + getattr(form0, f"tabulate_tensor_{np_type}"), + ) + kernel0( + ffi.cast("{type} *".format(type=mode), A.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) b = np.zeros(3, dtype=np_type) - kernel1 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form1, f"tabulate_tensor_{np_type}")) - kernel1(ffi.cast('{type} *'.format(type=mode), b.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel1 = ffi.cast( + f"ufcx_tabulate_tensor_{np_type} *", + getattr(form1, f"tabulate_tensor_{np_type}"), + ) + kernel1( + ffi.cast("{type} *".format(type=mode), b.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.allclose(A, expected_result) assert np.allclose(b, 1.0 / 6.0) -@pytest.mark.parametrize("mode,expected_result", [ - ("double", np.array([[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], dtype=np.float64) - - (1.0 / 24.0) * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.float64)), - ("double _Complex", - np.array([[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], dtype=np.complex128) - - (1.0j / 24.0) * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.complex128)), -]) +@pytest.mark.parametrize( + "mode,expected_result", + [ + ( + "double", + np.array( + [[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], + dtype=np.float64, + ) + - (1.0 / 24.0) + * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.float64), + ), + ( + "double _Complex", + np.array( + [[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], + dtype=np.complex128, + ) + - (1.0j / 24.0) + * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.complex128), + ), + ], +) def test_helmholtz_form_2d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) @@ -160,12 +246,13 @@ def test_helmholtz_form_2d(mode, expected_result, compile_args): a = (ufl.inner(ufl.grad(u), ufl.grad(v)) - ufl.inner(k * u, v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) - form0 = compiled_forms[0].integrals(module.lib.cell)[0] + form0 = compiled_forms[0].form_integrals[0] np_type = cdtype_to_numpy(mode) A = np.zeros((3, 3), dtype=np_type) @@ -176,44 +263,65 @@ def test_helmholtz_form_2d(mode, expected_result, compile_args): np_gtype = cdtype_to_numpy(geom_type) ffi = module.ffi - coords = np.array([[0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np_gtype) + coords = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype + ) kernel = getattr(form0, f"tabulate_tensor_{np_type}") - kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("{type} *".format(type=mode), A.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.allclose(A, expected_result) -@pytest.mark.parametrize("mode,expected_result", [ - ("double", np.array([[0.5, -1 / 6, -1 / 6, -1 / 6], - [-1 / 6, 1 / 6, 0.0, 0.0], - [-1 / 6, 0.0, 1 / 6, 0.0], - [-1 / 6, 0.0, 0.0, 1 / 6]], dtype=np.float64)), - ("double _Complex", - np.array( - [[0.5 + 0j, -1 / 6 + 0j, -1 / 6 + 0j, -1 / 6 + 0j], - [-1 / 6 + 0j, 1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j], - [-1 / 6 + 0j, 0.0 + 0j, 1 / 6 + 0j, 0.0 + 0j], - [-1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j, 1 / 6 + 0j]], - dtype=np.complex128)), -]) +@pytest.mark.parametrize( + "mode,expected_result", + [ + ( + "double", + np.array( + [ + [0.5, -1 / 6, -1 / 6, -1 / 6], + [-1 / 6, 1 / 6, 0.0, 0.0], + [-1 / 6, 0.0, 1 / 6, 0.0], + [-1 / 6, 0.0, 0.0, 1 / 6], + ], + dtype=np.float64, + ), + ), + ( + "double _Complex", + np.array( + [ + [0.5 + 0j, -1 / 6 + 0j, -1 / 6 + 0j, -1 / 6 + 0j], + [-1 / 6 + 0j, 1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j], + [-1 / 6 + 0j, 0.0 + 0j, 1 / 6 + 0j, 0.0 + 0j], + [-1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j, 1 / 6 + 0j], + ], + dtype=np.complex128, + ), + ), + ], +) def test_laplace_bilinear_form_3d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "tetrahedron", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) - form0 = compiled_forms[0].integrals(module.lib.cell)[0] + form0 = compiled_forms[0].form_integrals[0] np_type = cdtype_to_numpy(mode) A = np.zeros((4, 4), dtype=np_type) @@ -224,16 +332,19 @@ def test_laplace_bilinear_form_3d(mode, expected_result, compile_args): np_gtype = cdtype_to_numpy(geom_type) ffi = module.ffi - coords = np.array([0.0, 0.0, 0.0, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0], dtype=np_gtype) + coords = np.array( + [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], dtype=np_gtype + ) kernel = getattr(form0, f"tabulate_tensor_{np_type}") - kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("{type} *".format(type=mode), A.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.allclose(A, expected_result) @@ -244,31 +355,36 @@ def test_form_coefficient(compile_args): g = ufl.Coefficient(element) a = g * ufl.inner(u, v) * ufl.dx forms = [a] - compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms(forms, cffi_extra_compile_args=compile_args) + compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( + forms, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) - form0 = compiled_forms[0].integrals(module.lib.cell)[0] + form0 = compiled_forms[0].form_integrals[0] A = np.zeros((3, 3), dtype=np.float64) w = np.array([1.0, 1.0, 1.0], dtype=np.float64) c = np.array([], dtype=np.float64) perm = np.array([0], dtype=np.uint8) ffi = module.ffi - coords = np.array([[0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np.float64) + coords = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np.float64 + ) kernel = getattr(form0, "tabulate_tensor_float64") - kernel(ffi.cast('double *', A.ctypes.data), - ffi.cast('double *', w.ctypes.data), - ffi.cast('double *', c.ctypes.data), - ffi.cast('double *', coords.ctypes.data), ffi.NULL, - ffi.cast('uint8_t *', perm.ctypes.data)) + kernel( + ffi.cast("double *", A.ctypes.data), + ffi.cast("double *", w.ctypes.data), + ffi.cast("double *", c.ctypes.data), + ffi.cast("double *", coords.ctypes.data), + ffi.NULL, + ffi.cast("uint8_t *", perm.ctypes.data), + ) A_analytic = np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.float64) / 24.0 - A_diff = (A - A_analytic) + A_diff = A - A_analytic assert np.isclose(A_diff.max(), 0.0) assert np.isclose(A_diff.min(), 0.0) @@ -282,27 +398,53 @@ def test_subdomains(compile_args): a3 = ufl.inner(u, v) * ufl.ds(210) + ufl.inner(u, v) * ufl.ds(0) forms = [a0, a1, a2, a3] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': 'double'}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": "double"}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) + cell = module.lib.cell form0 = compiled_forms[0] - ids = form0.integral_ids(module.lib.cell) + ids = [ + form0.form_integral_ids[j] + for j in range( + form0.form_integral_offsets[cell], form0.form_integral_offsets[cell + 1] + ) + ] assert ids[0] == -1 and ids[1] == 2 form1 = compiled_forms[1] - ids = form1.integral_ids(module.lib.cell) + ids = [ + form1.form_integral_ids[j] + for j in range( + form1.form_integral_offsets[cell], form1.form_integral_offsets[cell + 1] + ) + ] assert ids[0] == -1 and ids[1] == 2 form2 = compiled_forms[2] - ids = form2.integral_ids(module.lib.cell) + ids = [ + form2.form_integral_ids[j] + for j in range( + form2.form_integral_offsets[cell], form2.form_integral_offsets[cell + 1] + ) + ] assert ids[0] == 1 and ids[1] == 2 form3 = compiled_forms[3] - assert form3.num_integrals(module.lib.cell) == 0 + assert ( + form3.form_integral_offsets[cell + 1] - form3.form_integral_offsets[cell] == 0 + ) - ids = form3.integral_ids(module.lib.exterior_facet) + ext_facet = module.lib.exterior_facet + ids = [ + form3.form_integral_ids[j] + for j in range( + form3.form_integral_offsets[ext_facet], + form3.form_integral_offsets[ext_facet + 1], + ) + ] assert ids[0] == 0 and ids[1] == 210 @@ -313,7 +455,8 @@ def test_interior_facet_integral(mode, compile_args): a0 = ufl.inner(ufl.jump(ufl.grad(u)), ufl.jump(ufl.grad(v))) * ufl.dS forms = [a0] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -325,7 +468,7 @@ def test_interior_facet_integral(mode, compile_args): ffi = module.ffi np_type = cdtype_to_numpy(mode) - integral0 = form0.integrals(module.lib.interior_facet)[0] + integral0 = form0.form_integrals[0] A = np.zeros((6, 6), dtype=np_type) w = np.array([], dtype=np_type) c = np.array([], dtype=np.float64) @@ -336,19 +479,23 @@ def test_interior_facet_integral(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([[0.0, 0.0, 0.0, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0], - [1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 1.0, 1.0, 0.0]], dtype=np_gtype) + coords = np.array( + [ + [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0], + ], + dtype=np_gtype, + ) kernel = getattr(integral0, f"tabulate_tensor_{np_type}") - kernel(ffi.cast(f'{mode} *', A.ctypes.data), - ffi.cast(f'{mode} *', w.ctypes.data), - ffi.cast(f'{mode} *', c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.cast('int *', facets.ctypes.data), - ffi.cast('uint8_t *', perms.ctypes.data)) + kernel( + ffi.cast(f"{mode} *", A.ctypes.data), + ffi.cast(f"{mode} *", w.ctypes.data), + ffi.cast(f"{mode} *", c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.cast("int *", facets.ctypes.data), + ffi.cast("uint8_t *", perms.ctypes.data), + ) @pytest.mark.parametrize("mode", ["double", "double _Complex"]) @@ -356,8 +503,9 @@ def test_conditional(mode, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) x = ufl.SpatialCoordinate(ufl.triangle) - condition = ufl.Or(ufl.ge(ufl.real(x[0] + x[1]), 0.1), - ufl.ge(ufl.real(x[1] + x[1]**2), 0.1)) + condition = ufl.Or( + ufl.ge(ufl.real(x[0] + x[1]), 0.1), ufl.ge(ufl.real(x[1] + x[1] ** 2), 0.1) + ) c1 = ufl.conditional(condition, 2.0, 1.0) a = c1 * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx @@ -368,10 +516,11 @@ def test_conditional(mode, compile_args): forms = [a, b] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) - form0 = compiled_forms[0].integrals(module.lib.cell)[0] - form1 = compiled_forms[1].integrals(module.lib.cell)[0] + form0 = compiled_forms[0].form_integrals[0] + form1 = compiled_forms[1].form_integrals[0] ffi = module.ffi np_type = cdtype_to_numpy(mode) @@ -383,15 +532,22 @@ def test_conditional(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([[0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np_gtype) + coords = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype + ) - kernel0 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form0, f"tabulate_tensor_{np_type}")) - kernel0(ffi.cast('{type} *'.format(type=mode), A1.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w1.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel0 = ffi.cast( + f"ufcx_tabulate_tensor_{np_type} *", + getattr(form0, f"tabulate_tensor_{np_type}"), + ) + kernel0( + ffi.cast("{type} *".format(type=mode), A1.ctypes.data), + ffi.cast("{type} *".format(type=mode), w1.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) expected_result = np.array([[2, -1, -1], [-1, 1, 0], [-1, 0, 1]], dtype=np_type) assert np.allclose(A1, expected_result) @@ -399,11 +555,18 @@ def test_conditional(mode, compile_args): A2 = np.zeros(3, dtype=np_type) w2 = np.array([1.0, 1.0, 1.0], dtype=np_type) - kernel1 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form1, f"tabulate_tensor_{np_type}")) - kernel1(ffi.cast('{type} *'.format(type=mode), A2.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w2.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel1 = ffi.cast( + f"ufcx_tabulate_tensor_{np_type} *", + getattr(form1, f"tabulate_tensor_{np_type}"), + ) + kernel1( + ffi.cast("{type} *".format(type=mode), A2.ctypes.data), + ffi.cast("{type} *".format(type=mode), w2.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) expected_result = np.ones(3, dtype=np_type) assert np.allclose(A2, expected_result) @@ -419,29 +582,44 @@ def test_custom_quadrature(compile_args): points = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.5, 0.5], [0.0, 0.5], [0.5, 0.0]] weights = [1 / 12] * 6 - a = u * v * ufl.dx(metadata={"quadrature_rule": "custom", - "quadrature_points": points, "quadrature_weights": weights}) + a = ( + u + * v + * ufl.dx( + metadata={ + "quadrature_rule": "custom", + "quadrature_points": points, + "quadrature_weights": weights, + } + ) + ) forms = [a] - compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms(forms, cffi_extra_compile_args=compile_args) + compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( + forms, cffi_extra_compile_args=compile_args + ) ffi = module.ffi form = compiled_forms[0] - default_integral = form.integrals(module.lib.cell)[0] + default_integral = form.form_integrals[0] A = np.zeros((6, 6), dtype=np.float64) w = np.array([], dtype=np.float64) c = np.array([], dtype=np.float64) - coords = np.array([[0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np.float64) + coords = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np.float64 + ) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel(ffi.cast("double *", A.ctypes.data), - ffi.cast("double *", w.ctypes.data), - ffi.cast("double *", c.ctypes.data), - ffi.cast("double *", coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("double *", A.ctypes.data), + ffi.cast("double *", w.ctypes.data), + ffi.cast("double *", c.ctypes.data), + ffi.cast("double *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) # Check that A is diagonal assert np.count_nonzero(A - np.diag(np.diagonal(A))) == 0 @@ -453,12 +631,19 @@ def test_curl_curl(compile_args): a = ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx forms = [a] - compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms(forms, cffi_extra_compile_args=compile_args) + compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( + forms, cffi_extra_compile_args=compile_args + ) -def lagrange_triangle_symbolic(order, corners=[(1, 0), (2, 0), (0, 1)], fun=lambda i: i): +def lagrange_triangle_symbolic( + order, corners=[(1, 0), (2, 0), (0, 1)], fun=lambda i: i +): from sympy import S - poly_basis = [x**i * y**j for i in range(order + 1) for j in range(order + 1 - i)] + + poly_basis = [ + x**i * y**j for i in range(order + 1) for j in range(order + 1 - i) + ] # vertices eval_points = [S(c) for c in corners] # edges @@ -468,24 +653,39 @@ def lagrange_triangle_symbolic(order, corners=[(1, 0), (2, 0), (0, 1)], fun=lamb if order > 3: raise NotImplementedError elif order == 3: - eval_points += [tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) - for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2]] + eval_points += [ + tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) + for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2] + ] else: - eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) - for i in range(1, order)] + eval_points += [ + tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) + for i in range(1, order) + ] # face for f in [(0, 1, 2)]: p0 = corners[f[0]] p1 = corners[f[1]] p2 = corners[f[2]] - eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) - + sympy.Rational((c - a) * j, order) for a, b, c in zip(p0, p1, p2)) - for i in range(1, order) for j in range(1, order - i)] - - dual_mat = [[f.subs(x, p[0]).subs(y, p[1]) for p in eval_points] for f in poly_basis] + eval_points += [ + tuple( + S(a) + + sympy.Rational((b - a) * i, order) + + sympy.Rational((c - a) * j, order) + for a, b, c in zip(p0, p1, p2) + ) + for i in range(1, order) + for j in range(1, order - i) + ] + + dual_mat = [ + [f.subs(x, p[0]).subs(y, p[1]) for p in eval_points] for f in poly_basis + ] dual_mat = sympy.Matrix(dual_mat) mat = dual_mat.inv() - functions = [sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows)] + functions = [ + sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows) + ] results = [] for f in functions: integrand = fun(f) @@ -494,10 +694,14 @@ def lagrange_triangle_symbolic(order, corners=[(1, 0), (2, 0), (0, 1)], fun=lamb @pytest.mark.parametrize("mode", ["double"]) -@pytest.mark.parametrize("sym_fun,ufl_fun", [ - (lambda i: i, lambda i: i), - (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), - (lambda i: i.diff(y), lambda i: ufl.grad(i)[1])]) +@pytest.mark.parametrize( + "sym_fun,ufl_fun", + [ + (lambda i: i, lambda i: i), + (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), + (lambda i: i.diff(y), lambda i: ufl.grad(i)[1]), + ], +) @pytest.mark.parametrize("order", [1, 2, 3]) def test_lagrange_triangle(compile_args, order, mode, sym_fun, ufl_fun): sym = lagrange_triangle_symbolic(order, fun=sym_fun) @@ -507,13 +711,14 @@ def test_lagrange_triangle(compile_args, order, mode, sym_fun, ufl_fun): a = ufl_fun(v) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) ffi = module.ffi form0 = compiled_forms[0] - assert form0.num_integrals(module.lib.cell) == 1 - default_integral = form0.integrals(module.lib.cell)[0] + assert form0.form_integral_offsets[module.lib.cell + 1] == 1 + default_integral = form0.form_integrals[0] np_type = cdtype_to_numpy(mode) b = np.zeros((order + 2) * (order + 1) // 2, dtype=np_type) @@ -522,25 +727,35 @@ def test_lagrange_triangle(compile_args, order, mode, sym_fun, ufl_fun): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([[1.0, 0.0, 0.0], - [2.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np_gtype) + coords = np.array( + [[1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype + ) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel(ffi.cast('{type} *'.format(type=mode), b.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.NULL, - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("{type} *".format(type=mode), b.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.NULL, + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) # Check that the result is the same as for sympy assert np.allclose(b, [float(i) for i in sym]) -def lagrange_tetrahedron_symbolic(order, corners=[(1, 0, 0), (2, 0, 0), (0, 1, 0), (0, 0, 1)], fun=lambda i: i): +def lagrange_tetrahedron_symbolic( + order, corners=[(1, 0, 0), (2, 0, 0), (0, 1, 0), (0, 0, 1)], fun=lambda i: i +): from sympy import S + poly_basis = [ - x**i * y**j * z**k for i in range(order + 1) for j in range(order + 1 - i) - for k in range(order + 1 - i - j)] + x**i * y**j * z**k + for i in range(order + 1) + for j in range(order + 1 - i) + for k in range(order + 1 - i - j) + ] # vertices eval_points = [S(c) for c in corners] # edges @@ -550,45 +765,78 @@ def lagrange_tetrahedron_symbolic(order, corners=[(1, 0, 0), (2, 0, 0), (0, 1, 0 if order > 3: raise NotImplementedError elif order == 3: - eval_points += [tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) - for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2]] + eval_points += [ + tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) + for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2] + ] else: - eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) - for i in range(1, order)] + eval_points += [ + tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) + for i in range(1, order) + ] # face for f in [(1, 2, 3), (0, 2, 3), (0, 1, 3), (0, 1, 2)]: p0 = corners[f[0]] p1 = corners[f[1]] p2 = corners[f[2]] - eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) - + sympy.Rational((c - a) * j, order) for a, b, c in zip(p0, p1, p2)) - for i in range(1, order) for j in range(1, order - i)] + eval_points += [ + tuple( + S(a) + + sympy.Rational((b - a) * i, order) + + sympy.Rational((c - a) * j, order) + for a, b, c in zip(p0, p1, p2) + ) + for i in range(1, order) + for j in range(1, order - i) + ] # interior for v in [(0, 1, 2, 3)]: p0 = corners[v[0]] p1 = corners[v[1]] p2 = corners[v[2]] p3 = corners[v[3]] - eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) + sympy.Rational((c - a) * j, order) - + sympy.Rational((d - a) * k, order) for a, b, c, d in zip(p0, p1, p2, p3)) - for i in range(1, order) for j in range(1, order - i) for k in range(1, order - i - j)] - - dual_mat = [[f.subs(x, p[0]).subs(y, p[1]).subs(z, p[2]) for p in eval_points] for f in poly_basis] + eval_points += [ + tuple( + S(a) + + sympy.Rational((b - a) * i, order) + + sympy.Rational((c - a) * j, order) + + sympy.Rational((d - a) * k, order) + for a, b, c, d in zip(p0, p1, p2, p3) + ) + for i in range(1, order) + for j in range(1, order - i) + for k in range(1, order - i - j) + ] + + dual_mat = [ + [f.subs(x, p[0]).subs(y, p[1]).subs(z, p[2]) for p in eval_points] + for f in poly_basis + ] dual_mat = sympy.Matrix(dual_mat) mat = dual_mat.inv() - functions = [sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows)] + functions = [ + sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows) + ] results = [] for f in functions: integrand = fun(f) - results.append(integrand.integrate((x, 1 - y - z, 2 - 2 * y - 2 * z), (y, 0, 1 - z), (z, 0, 1))) + results.append( + integrand.integrate( + (x, 1 - y - z, 2 - 2 * y - 2 * z), (y, 0, 1 - z), (z, 0, 1) + ) + ) return results @pytest.mark.parametrize("mode", ["double"]) -@pytest.mark.parametrize("sym_fun,ufl_fun", [ - (lambda i: i, lambda i: i), - (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), - (lambda i: i.diff(y), lambda i: ufl.grad(i)[1])]) +@pytest.mark.parametrize( + "sym_fun,ufl_fun", + [ + (lambda i: i, lambda i: i), + (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), + (lambda i: i.diff(y), lambda i: ufl.grad(i)[1]), + ], +) @pytest.mark.parametrize("order", [1, 2, 3]) def test_lagrange_tetrahedron(compile_args, order, mode, sym_fun, ufl_fun): sym = lagrange_tetrahedron_symbolic(order, fun=sym_fun) @@ -598,14 +846,15 @@ def test_lagrange_tetrahedron(compile_args, order, mode, sym_fun, ufl_fun): a = ufl_fun(v) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) ffi = module.ffi form0 = compiled_forms[0] - assert form0.num_integrals(module.lib.cell) == 1 + assert form0.form_integral_offsets[module.lib.cell + 1] == 1 - default_integral = form0.integrals(module.lib.cell)[0] + default_integral = form0.form_integrals[0] np_type = cdtype_to_numpy(mode) b = np.zeros((order + 3) * (order + 2) * (order + 1) // 6, dtype=np_type) @@ -614,16 +863,19 @@ def test_lagrange_tetrahedron(compile_args, order, mode, sym_fun, ufl_fun): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([1.0, 0.0, 0.0, - 2.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0], dtype=np_gtype) + coords = np.array( + [1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], dtype=np_gtype + ) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel(ffi.cast('{type} *'.format(type=mode), b.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w.ctypes.data), - ffi.NULL, - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("{type} *".format(type=mode), b.ctypes.data), + ffi.cast("{type} *".format(type=mode), w.ctypes.data), + ffi.NULL, + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) # Check that the result is the same as for sympy assert np.allclose(b, [float(i) for i in sym]) @@ -636,26 +888,48 @@ def test_prism(compile_args): L = v * ufl.dx forms = [L] compiled_forms, module, _ = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': 'double'}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": "double"}, cffi_extra_compile_args=compile_args + ) ffi = module.ffi form0 = compiled_forms[0] - assert form0.num_integrals(module.lib.cell) == 1 + assert form0.form_integral_offsets[module.lib.cell + 1] == 1 - default_integral = form0.integrals(module.lib.cell)[0] + default_integral = form0.form_integrals[0] b = np.zeros(6, dtype=np.float64) - coords = np.array([1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, - 1.0, 0.0, 1.0, - 0.0, 1.0, 1.0, - 0.0, 0.0, 1.0], dtype=np.float64) + coords = np.array( + [ + 1.0, + 0.0, + 0.0, + 0.0, + 1.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 1.0, + 0.0, + 0.0, + 1.0, + ], + dtype=np.float64, + ) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel(ffi.cast('double *', b.ctypes.data), - ffi.NULL, - ffi.NULL, - ffi.cast('double *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("double *", b.ctypes.data), + ffi.NULL, + ffi.NULL, + ffi.cast("double *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.isclose(sum(b), 0.5) @@ -673,10 +947,11 @@ def test_complex_operations(compile_args): forms = [J1, J2] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) + forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args + ) - form0 = compiled_forms[0].integrals(module.lib.cell)[0] - form1 = compiled_forms[1].integrals(module.lib.cell)[0] + form0 = compiled_forms[0].form_integrals[0] + form1 = compiled_forms[1].form_integrals[0] ffi = module.ffi np_type = cdtype_to_numpy(mode) @@ -686,27 +961,48 @@ def test_complex_operations(compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array([[0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0]], dtype=np_gtype) + coords = np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype + ) J_1 = np.zeros((1), dtype=np_type) - kernel0 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form0, f"tabulate_tensor_{np_type}")) - kernel0(ffi.cast('{type} *'.format(type=mode), J_1.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w1.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) - - expected_result = np.array([0.5 * np.real(w1[0]) * np.imag(w1[1]) - * (np.real(w1[0]) - 1j * np.imag(w1[0]))], dtype=np_type) + kernel0 = ffi.cast( + f"ufcx_tabulate_tensor_{np_type} *", + getattr(form0, f"tabulate_tensor_{np_type}"), + ) + kernel0( + ffi.cast("{type} *".format(type=mode), J_1.ctypes.data), + ffi.cast("{type} *".format(type=mode), w1.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) + + expected_result = np.array( + [ + 0.5 + * np.real(w1[0]) + * np.imag(w1[1]) + * (np.real(w1[0]) - 1j * np.imag(w1[0])) + ], + dtype=np_type, + ) assert np.allclose(J_1, expected_result) J_2 = np.zeros((1), dtype=np_type) - kernel1 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form1, f"tabulate_tensor_{np_type}")) - kernel1(ffi.cast('{type} *'.format(type=mode), J_2.ctypes.data), - ffi.cast('{type} *'.format(type=mode), w1.ctypes.data), - ffi.cast('{type} *'.format(type=mode), c.ctypes.data), - ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel1 = ffi.cast( + f"ufcx_tabulate_tensor_{np_type} *", + getattr(form1, f"tabulate_tensor_{np_type}"), + ) + kernel1( + ffi.cast("{type} *".format(type=mode), J_2.ctypes.data), + ffi.cast("{type} *".format(type=mode), w1.ctypes.data), + ffi.cast("{type} *".format(type=mode), c.ctypes.data), + ffi.cast(f"{geom_type} *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.allclose(J_2, expected_result) @@ -726,7 +1022,8 @@ def test_invalid_function_name(compile_args): try: compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args) + forms, cffi_extra_compile_args=compile_args + ) except ValueError: pass except Exception: @@ -737,36 +1034,37 @@ def test_invalid_function_name(compile_args): def test_interval_vertex_quadrature(compile_args): - c_el = basix.ufl.element("Lagrange", "interval", 1, rank=1) mesh = ufl.Mesh(c_el) x = ufl.SpatialCoordinate(mesh) - dx = ufl.Measure( - "dx", metadata={"quadrature_rule": "vertex"}) + dx = ufl.Measure("dx", metadata={"quadrature_rule": "vertex"}) b = x[0] * dx forms = [b] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args) + forms, cffi_extra_compile_args=compile_args + ) ffi = module.ffi form0 = compiled_forms[0] - assert form0.num_integrals(module.lib.cell) == 1 + assert form0.form_integral_offsets[module.lib.cell + 1] == 1 - default_integral = form0.integrals(module.lib.cell)[0] + default_integral = form0.form_integrals[0] J = np.zeros(1, dtype=np.float64) a = np.pi b = np.exp(1) - coords = np.array([a, 0.0, 0.0, - - b, 0.0, 0.0], dtype=np.float64) + coords = np.array([a, 0.0, 0.0, b, 0.0, 0.0], dtype=np.float64) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel(ffi.cast('double *', J.ctypes.data), - ffi.NULL, - ffi.NULL, - ffi.cast('double *', coords.ctypes.data), ffi.NULL, ffi.NULL) + kernel( + ffi.cast("double *", J.ctypes.data), + ffi.NULL, + ffi.NULL, + ffi.cast("double *", coords.ctypes.data), + ffi.NULL, + ffi.NULL, + ) assert np.isclose(J[0], (0.5 * a + 0.5 * b) * np.abs(b - a)) @@ -778,9 +1076,8 @@ def test_facet_vertex_quadrature(compile_args): mesh = ufl.Mesh(c_el) x = ufl.SpatialCoordinate(mesh) - ds = ufl.Measure( - "ds", metadata={"quadrature_rule": "vertex"}) - expr = (x[0] + ufl.cos(x[1])) + ds = ufl.Measure("ds", metadata={"quadrature_rule": "vertex"}) + expr = x[0] + ufl.cos(x[1]) b1 = expr * ds ds_c = ufl.Measure( "ds", @@ -788,40 +1085,47 @@ def test_facet_vertex_quadrature(compile_args): "quadrature_rule": "custom", "quadrature_points": np.array([[0.0], [1.0]]), "quadrature_weights": np.array([1.0 / 2.0, 1.0 / 2.0]), - } + }, ) b2 = expr * ds_c forms = [b1, b2] compiled_forms, module, _ = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args) + forms, cffi_extra_compile_args=compile_args + ) ffi = module.ffi assert len(compiled_forms) == 2 solutions = [] for form in compiled_forms: - assert form.num_integrals(module.lib.exterior_facet) == 1 + assert form.form_integral_offsets[module.lib.exterior_facet + 1] == 1 - default_integral = form.integrals(module.lib.exterior_facet)[0] + default_integral = form.form_integrals[0] J = np.zeros(1, dtype=np.float64) a = np.pi b = np.exp(1) - coords = np.array([a, 0.1, 0.0, - a + b, 0.0, 0.0, - a, a, 0., - a + 2 * b, a, 0.], dtype=np.float64) + coords = np.array( + [a, 0.1, 0.0, a + b, 0.0, 0.0, a, a, 0.0, a + 2 * b, a, 0.0], + dtype=np.float64, + ) # First facet is between vertex 0 and 1 in coords facets = np.array([0], dtype=np.intc) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel(ffi.cast('double *', J.ctypes.data), - ffi.NULL, - ffi.NULL, - ffi.cast('double *', coords.ctypes.data), - ffi.cast('int *', facets.ctypes.data), - ffi.NULL) + kernel( + ffi.cast("double *", J.ctypes.data), + ffi.NULL, + ffi.NULL, + ffi.cast("double *", coords.ctypes.data), + ffi.cast("int *", facets.ctypes.data), + ffi.NULL, + ) solutions.append(J[0]) # Test against exact result - assert np.isclose(J[0], (0.5 * (a + np.cos(0.1)) + 0.5 * (a + b + np.cos(0))) * np.sqrt(b**2 + 0.1**2)) + assert np.isclose( + J[0], + (0.5 * (a + np.cos(0.1)) + 0.5 * (a + b + np.cos(0))) + * np.sqrt(b**2 + 0.1**2), + ) # Compare custom quadrature with vertex quadrature assert np.isclose(solutions[0], solutions[1]) @@ -843,14 +1147,15 @@ def test_manifold_derivatives(compile_args): u = ufl.Coefficient(V) d = 5.3 - f_ex = d * order * (order - 1) * x[1]**(order - 2) + f_ex = d * order * (order - 1) * x[1] ** (order - 2) expr = u.dx(1).dx(1) - f_ex J = expr * expr * dx compiled_forms, module, _ = ffcx.codegeneration.jit.compile_forms( - [J], cffi_extra_compile_args=compile_args) + [J], cffi_extra_compile_args=compile_args + ) - default_integral = compiled_forms[0].integrals(module.lib.cell)[0] + default_integral = compiled_forms[0].form_integrals[0] scale = 2.5 coords = np.array([0.0, 0.0, 0.0, 0.0, scale, 0.0], dtype=np.float64) dof_coords = el.element.points.reshape(-1) @@ -863,10 +1168,13 @@ def test_manifold_derivatives(compile_args): ffi = module.ffi J = np.zeros(1, dtype=np.float64) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel(ffi.cast('double *', J.ctypes.data), - ffi.cast('double *', w.ctypes.data), - ffi.cast('double *', c.ctypes.data), - ffi.cast('double *', coords.ctypes.data), ffi.NULL, - ffi.cast('uint8_t *', perm.ctypes.data)) + kernel( + ffi.cast("double *", J.ctypes.data), + ffi.cast("double *", w.ctypes.data), + ffi.cast("double *", c.ctypes.data), + ffi.cast("double *", coords.ctypes.data), + ffi.NULL, + ffi.cast("uint8_t *", perm.ctypes.data), + ) assert np.isclose(J[0], 0.0) From d0014fb2e6c16045fa7c698e115e4450390a63c4 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Mon, 14 Aug 2023 17:55:05 +0100 Subject: [PATCH 12/14] Remove break/continue --- ffcx/codegeneration/C/cnodes.py | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/ffcx/codegeneration/C/cnodes.py b/ffcx/codegeneration/C/cnodes.py index f8b927763..4820b64f1 100644 --- a/ffcx/codegeneration/C/cnodes.py +++ b/ffcx/codegeneration/C/cnodes.py @@ -1115,34 +1115,6 @@ def flops(self): # Simple statements -class Break(CStatement): - __slots__ = () - is_scoped = True - - def cs_format(self, precision=None): - return "break;" - - def __eq__(self, other): - return isinstance(other, type(self)) - - def flops(self): - return 0 - - -class Continue(CStatement): - __slots__ = () - is_scoped = True - - def cs_format(self, precision=None): - return "continue;" - - def __eq__(self, other): - return isinstance(other, type(self)) - - def flops(self): - return 0 - - class Return(CStatement): __slots__ = ("value",) is_scoped = True From 58ed1c274dbf64b2f7e1cf67b8695b54c47b4071 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Mon, 14 Aug 2023 18:00:01 +0100 Subject: [PATCH 13/14] fix --- ffcx/codegeneration/C/cnodes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ffcx/codegeneration/C/cnodes.py b/ffcx/codegeneration/C/cnodes.py index 4820b64f1..3d6eb7666 100644 --- a/ffcx/codegeneration/C/cnodes.py +++ b/ffcx/codegeneration/C/cnodes.py @@ -1464,7 +1464,7 @@ def _is_simple_if_body(body): if len(body.statements) > 1: return False (body,) = body.statements - return isinstance(body, (Return, AssignOp, Break, Continue)) + return isinstance(body, (Return, AssignOp)) class If(CStatement): From e265aa9c24bc325a57c155777f70b1e984d7dac9 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Tue, 15 Aug 2023 12:17:45 +0100 Subject: [PATCH 14/14] Reset files in tests --- ffcx/codegeneration/form.py | 1 + test/test_add_mode.py | 78 ++-- test/test_blocked_elements.py | 17 +- test/test_jit_forms.py | 851 +++++++++++----------------------- 4 files changed, 316 insertions(+), 631 deletions(-) diff --git a/ffcx/codegeneration/form.py b/ffcx/codegeneration/form.py index 0014aa38e..2e95edadd 100644 --- a/ffcx/codegeneration/form.py +++ b/ffcx/codegeneration/form.py @@ -88,6 +88,7 @@ def generator(ir, options): integrals = [] integral_ids = [] integral_offsets = [0] + # Note: the order of this list is defined by the enum ufcx_integral_type in ufcx.h for itg_type in ("cell", "exterior_facet", "interior_facet"): integrals += [L.AddressOf(L.Symbol(itg)) for itg in ir.integral_names[itg_type]] integral_ids += ir.subdomain_ids[itg_type] diff --git a/test/test_add_mode.py b/test/test_add_mode.py index 778b81a00..87d159a78 100644 --- a/test/test_add_mode.py +++ b/test/test_add_mode.py @@ -13,17 +13,21 @@ from ffcx.naming import cdtype_to_numpy, scalar_to_value_type -@pytest.mark.parametrize( - "mode", ["double", "float", "long double", "double _Complex", "float _Complex"] -) +@pytest.mark.parametrize("mode", + [ + "double", + "float", + "long double", + "double _Complex", + "float _Complex" + ]) def test_additive_facet_integral(mode, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) a = ufl.inner(u, v) * ufl.ds forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -33,8 +37,7 @@ def test_additive_facet_integral(mode, compile_args): integral_offsets = form0.form_integral_offsets ex = module.lib.exterior_facet - num_integrals = integral_offsets[ex + 1] - integral_offsets[ex] - assert num_integrals == 1 + assert integral_offsets[ex + 1] - integral_offsets[ex] == 1 integral_id = form0.form_integral_ids[integral_offsets[ex]] assert integral_id == -1 @@ -49,38 +52,32 @@ def test_additive_facet_integral(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [0.0, 2.0, 0.0, np.sqrt(3.0), -1.0, 0.0, -np.sqrt(3.0), -1.0, 0.0], - dtype=np_gtype, - ) + coords = np.array([0.0, 2.0, 0.0, + np.sqrt(3.0), -1.0, 0.0, + -np.sqrt(3.0), -1.0, 0.0], dtype=np_gtype) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") for i in range(3): facets[0] = i - kernel( - ffi.cast("{type} *".format(type=mode), A.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.cast("int *", facets.ctypes.data), - ffi.cast("uint8_t *", perm.ctypes.data), - ) + kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), + ffi.cast('int *', facets.ctypes.data), + ffi.cast('uint8_t *', perm.ctypes.data)) assert np.isclose(A.sum(), np.sqrt(12) * (i + 1)) -@pytest.mark.parametrize( - "mode", ["double", "float", "long double", "double _Complex", "float _Complex"] -) +@pytest.mark.parametrize("mode", ["double", "float", "long double", "double _Complex", "float _Complex"]) def test_additive_cell_integral(mode, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -88,8 +85,8 @@ def test_additive_cell_integral(mode, compile_args): ffi = module.ffi form0 = compiled_forms[0] - offsets = form0.form_integral_offsets cell = module.lib.cell + offsets = form0.form_integral_offsets num_integrals = offsets[cell + 1] - offsets[cell] assert num_integrals == 1 integral_id = form0.form_integral_ids[offsets[cell]] @@ -104,31 +101,22 @@ def test_additive_cell_integral(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [0.0, 2.0, 0.0, np.sqrt(3.0), -1.0, 0.0, -np.sqrt(3.0), -1.0, 0.0], - dtype=np_gtype, - ) + coords = np.array([0.0, 2.0, 0.0, + np.sqrt(3.0), -1.0, 0.0, + -np.sqrt(3.0), -1.0, 0.0], dtype=np_gtype) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel( - ffi.cast("{type} *".format(type=mode), A.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) A0 = np.array(A) for i in range(3): - kernel( - ffi.cast("{type} *".format(type=mode), A.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.all(np.isclose(A, (i + 2) * A0)) diff --git a/test/test_blocked_elements.py b/test/test_blocked_elements.py index 3388a8638..c67660a54 100644 --- a/test/test_blocked_elements.py +++ b/test/test_blocked_elements.py @@ -15,8 +15,7 @@ def test_finite_element(compile_args): ufl_element = basix.ufl.element("Lagrange", "triangle", 1) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args - ) + [ufl_element], cffi_extra_compile_args=compile_args) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 2 @@ -35,6 +34,7 @@ def test_finite_element(compile_args): assert ufcx_dofmap.num_element_support_dofs == 3 off = np.array([ufcx_dofmap.entity_dof_offsets[i] for i in range(8)]) assert np.all(np.diff(off) == [1, 1, 1, 0, 0, 0, 0]) + for v in range(3): assert ufcx_dofmap.entity_dofs[v] == v assert ufcx_dofmap.num_sub_dofmaps == 0 @@ -43,8 +43,7 @@ def test_finite_element(compile_args): def test_vector_element(compile_args): ufl_element = basix.ufl.element("Lagrange", "triangle", 1, rank=1) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args - ) + [ufl_element], cffi_extra_compile_args=compile_args) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 2 @@ -74,8 +73,7 @@ def test_vector_element(compile_args): def test_tensor_element(compile_args): ufl_element = basix.ufl.element("Lagrange", "triangle", 1, rank=2) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args - ) + [ufl_element], cffi_extra_compile_args=compile_args) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 2 @@ -105,12 +103,9 @@ def test_tensor_element(compile_args): def test_vector_quadrature_element(compile_args): - ufl_element = ufl.VectorElement( - ufl.FiniteElement("Quadrature", "tetrahedron", degree=2, quad_scheme="default") - ) + ufl_element = ufl.VectorElement(ufl.FiniteElement("Quadrature", "tetrahedron", degree=2, quad_scheme="default")) jit_compiled_elements, module, code = ffcx.codegeneration.jit.compile_elements( - [ufl_element], cffi_extra_compile_args=compile_args - ) + [ufl_element], cffi_extra_compile_args=compile_args) ufcx_element, ufcx_dofmap = jit_compiled_elements[0] assert ufcx_element.topological_dimension == 3 diff --git a/test/test_jit_forms.py b/test/test_jit_forms.py index 8f67ec4ba..caa085d35 100644 --- a/test/test_jit_forms.py +++ b/test/test_jit_forms.py @@ -15,29 +15,14 @@ from ffcx.naming import cdtype_to_numpy, scalar_to_value_type -@pytest.mark.parametrize( - "mode,expected_result", - [ - ( - "double", - np.array( - [[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], - dtype=np.float64, - ), - ), - ( - "double _Complex", - np.array( - [ - [1.0 + 0j, -0.5 + 0j, -0.5 + 0j], - [-0.5 + 0j, 0.5 + 0j, 0.0 + 0j], - [-0.5 + 0j, 0.0 + 0j, 0.5 + 0j], - ], - dtype=np.complex128, - ), - ), - ], -) +@pytest.mark.parametrize("mode,expected_result", [ + ("double", np.array([[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], dtype=np.float64)), + ("double _Complex", + np.array( + [[1.0 + 0j, -0.5 + 0j, -0.5 + 0j], [-0.5 + 0j, 0.5 + 0j, 0.0 + 0j], + [-0.5 + 0j, 0.0 + 0j, 0.5 + 0j]], + dtype=np.complex128)), +]) def test_laplace_bilinear_form_2d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) kappa = ufl.Constant(ufl.triangle, shape=(2, 2)) @@ -46,8 +31,7 @@ def test_laplace_bilinear_form_2d(mode, expected_result, compile_args): a = ufl.tr(kappa) * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -72,84 +56,47 @@ def test_laplace_bilinear_form_2d(mode, expected_result, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype - ) + coords = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np_gtype) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel( - ffi.cast("{type} *".format(type=mode), A.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.allclose(A, np.trace(kappa_value) * expected_result) -@pytest.mark.parametrize( - "mode,expected_result", - [ - ( - "float", - np.array( - [ - [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], - ], - dtype=np.float32, - ), - ), - ( - "long double", - np.array( - [ - [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], - ], - dtype=np.longdouble, - ), - ), - ( - "double", - np.array( - [ - [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], - ], - dtype=np.float64, - ), - ), - ( - "double _Complex", - np.array( - [ - [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], - ], - dtype=np.complex128, - ), - ), - ( - "float _Complex", - np.array( - [ - [1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], - [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0], - ], - dtype=np.complex64, - ), - ), - ], -) +@pytest.mark.parametrize("mode,expected_result", [ + ("float", + np.array( + [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], + dtype=np.float32)), + ("long double", + np.array( + [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], + dtype=np.longdouble)), + ("double", + np.array( + [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], + dtype=np.float64)), + ("double _Complex", + np.array( + [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], + dtype=np.complex128)), + ("float _Complex", + np.array( + [[1.0 / 12.0, 1.0 / 24.0, 1.0 / 24.0], [1.0 / 24.0, 1.0 / 12.0, 1.0 / 24.0], + [1.0 / 24.0, 1.0 / 24.0, 1.0 / 12.0]], + dtype=np.complex64)), +]) def test_mass_bilinear_form_2d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) @@ -157,8 +104,7 @@ def test_mass_bilinear_form_2d(mode, expected_result, compile_args): L = ufl.conj(v) * ufl.dx forms = [a, L] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -175,64 +121,34 @@ def test_mass_bilinear_form_2d(mode, expected_result, compile_args): np_gtype = cdtype_to_numpy(geom_type) ffi = module.ffi - coords = np.array( - [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype - ) + coords = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np_gtype) - kernel0 = ffi.cast( - f"ufcx_tabulate_tensor_{np_type} *", - getattr(form0, f"tabulate_tensor_{np_type}"), - ) - kernel0( - ffi.cast("{type} *".format(type=mode), A.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel0 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form0, f"tabulate_tensor_{np_type}")) + kernel0(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) b = np.zeros(3, dtype=np_type) - kernel1 = ffi.cast( - f"ufcx_tabulate_tensor_{np_type} *", - getattr(form1, f"tabulate_tensor_{np_type}"), - ) - kernel1( - ffi.cast("{type} *".format(type=mode), b.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel1 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form1, f"tabulate_tensor_{np_type}")) + kernel1(ffi.cast('{type} *'.format(type=mode), b.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.allclose(A, expected_result) assert np.allclose(b, 1.0 / 6.0) -@pytest.mark.parametrize( - "mode,expected_result", - [ - ( - "double", - np.array( - [[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], - dtype=np.float64, - ) - - (1.0 / 24.0) - * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.float64), - ), - ( - "double _Complex", - np.array( - [[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], - dtype=np.complex128, - ) - - (1.0j / 24.0) - * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.complex128), - ), - ], -) +@pytest.mark.parametrize("mode,expected_result", [ + ("double", np.array([[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], dtype=np.float64) + - (1.0 / 24.0) * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.float64)), + ("double _Complex", + np.array([[1.0, -0.5, -0.5], [-0.5, 0.5, 0.0], [-0.5, 0.0, 0.5]], dtype=np.complex128) + - (1.0j / 24.0) * np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.complex128)), +]) def test_helmholtz_form_2d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) @@ -246,8 +162,7 @@ def test_helmholtz_form_2d(mode, expected_result, compile_args): a = (ufl.inner(ufl.grad(u), ufl.grad(v)) - ufl.inner(k * u, v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -263,60 +178,39 @@ def test_helmholtz_form_2d(mode, expected_result, compile_args): np_gtype = cdtype_to_numpy(geom_type) ffi = module.ffi - coords = np.array( - [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype - ) + coords = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np_gtype) kernel = getattr(form0, f"tabulate_tensor_{np_type}") - kernel( - ffi.cast("{type} *".format(type=mode), A.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.allclose(A, expected_result) -@pytest.mark.parametrize( - "mode,expected_result", - [ - ( - "double", - np.array( - [ - [0.5, -1 / 6, -1 / 6, -1 / 6], - [-1 / 6, 1 / 6, 0.0, 0.0], - [-1 / 6, 0.0, 1 / 6, 0.0], - [-1 / 6, 0.0, 0.0, 1 / 6], - ], - dtype=np.float64, - ), - ), - ( - "double _Complex", - np.array( - [ - [0.5 + 0j, -1 / 6 + 0j, -1 / 6 + 0j, -1 / 6 + 0j], - [-1 / 6 + 0j, 1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j], - [-1 / 6 + 0j, 0.0 + 0j, 1 / 6 + 0j, 0.0 + 0j], - [-1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j, 1 / 6 + 0j], - ], - dtype=np.complex128, - ), - ), - ], -) +@pytest.mark.parametrize("mode,expected_result", [ + ("double", np.array([[0.5, -1 / 6, -1 / 6, -1 / 6], + [-1 / 6, 1 / 6, 0.0, 0.0], + [-1 / 6, 0.0, 1 / 6, 0.0], + [-1 / 6, 0.0, 0.0, 1 / 6]], dtype=np.float64)), + ("double _Complex", + np.array( + [[0.5 + 0j, -1 / 6 + 0j, -1 / 6 + 0j, -1 / 6 + 0j], + [-1 / 6 + 0j, 1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j], + [-1 / 6 + 0j, 0.0 + 0j, 1 / 6 + 0j, 0.0 + 0j], + [-1 / 6 + 0j, 0.0 + 0j, 0.0 + 0j, 1 / 6 + 0j]], + dtype=np.complex128)), +]) def test_laplace_bilinear_form_3d(mode, expected_result, compile_args): element = basix.ufl.element("Lagrange", "tetrahedron", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -332,19 +226,16 @@ def test_laplace_bilinear_form_3d(mode, expected_result, compile_args): np_gtype = cdtype_to_numpy(geom_type) ffi = module.ffi - coords = np.array( - [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], dtype=np_gtype - ) + coords = np.array([0.0, 0.0, 0.0, + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0], dtype=np_gtype) kernel = getattr(form0, f"tabulate_tensor_{np_type}") - kernel( - ffi.cast("{type} *".format(type=mode), A.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('{type} *'.format(type=mode), A.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.allclose(A, expected_result) @@ -355,9 +246,7 @@ def test_form_coefficient(compile_args): g = ufl.Coefficient(element) a = g * ufl.inner(u, v) * ufl.dx forms = [a] - compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args - ) + compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms(forms, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -369,22 +258,19 @@ def test_form_coefficient(compile_args): perm = np.array([0], dtype=np.uint8) ffi = module.ffi - coords = np.array( - [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np.float64 - ) + coords = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np.float64) kernel = getattr(form0, "tabulate_tensor_float64") - kernel( - ffi.cast("double *", A.ctypes.data), - ffi.cast("double *", w.ctypes.data), - ffi.cast("double *", c.ctypes.data), - ffi.cast("double *", coords.ctypes.data), - ffi.NULL, - ffi.cast("uint8_t *", perm.ctypes.data), - ) + kernel(ffi.cast('double *', A.ctypes.data), + ffi.cast('double *', w.ctypes.data), + ffi.cast('double *', c.ctypes.data), + ffi.cast('double *', coords.ctypes.data), ffi.NULL, + ffi.cast('uint8_t *', perm.ctypes.data)) A_analytic = np.array([[2, 1, 1], [1, 2, 1], [1, 1, 2]], dtype=np.float64) / 24.0 - A_diff = A - A_analytic + A_diff = (A - A_analytic) assert np.isclose(A_diff.max(), 0.0) assert np.isclose(A_diff.min(), 0.0) @@ -398,53 +284,32 @@ def test_subdomains(compile_args): a3 = ufl.inner(u, v) * ufl.ds(210) + ufl.inner(u, v) * ufl.ds(0) forms = [a0, a1, a2, a3] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": "double"}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': 'double'}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) - cell = module.lib.cell form0 = compiled_forms[0] - ids = [ - form0.form_integral_ids[j] - for j in range( - form0.form_integral_offsets[cell], form0.form_integral_offsets[cell + 1] - ) - ] + offsets = form0.form_integral_offsets + cell = module.lib.cell + ids = [form0.form_integral_ids[j] for j in range(offsets[cell], offsets[cell + 1])] assert ids[0] == -1 and ids[1] == 2 form1 = compiled_forms[1] - ids = [ - form1.form_integral_ids[j] - for j in range( - form1.form_integral_offsets[cell], form1.form_integral_offsets[cell + 1] - ) - ] + offsets = form1.form_integral_offsets + ids = [form1.form_integral_ids[j] for j in range(offsets[cell], offsets[cell + 1])] assert ids[0] == -1 and ids[1] == 2 form2 = compiled_forms[2] - ids = [ - form2.form_integral_ids[j] - for j in range( - form2.form_integral_offsets[cell], form2.form_integral_offsets[cell + 1] - ) - ] + offsets = form2.form_integral_offsets + ids = [form2.form_integral_ids[j] for j in range(offsets[cell], offsets[cell + 1])] assert ids[0] == 1 and ids[1] == 2 form3 = compiled_forms[3] - assert ( - form3.form_integral_offsets[cell + 1] - form3.form_integral_offsets[cell] == 0 - ) - - ext_facet = module.lib.exterior_facet - ids = [ - form3.form_integral_ids[j] - for j in range( - form3.form_integral_offsets[ext_facet], - form3.form_integral_offsets[ext_facet + 1], - ) - ] + offsets = form3.form_integral_offsets + assert offsets[cell + 1] - offsets[cell] == 0 + exf = module.lib.exterior_facet + ids = [form3.form_integral_ids[j] for j in range(offsets[exf], offsets[exf + 1])] assert ids[0] == 0 and ids[1] == 210 @@ -455,8 +320,7 @@ def test_interior_facet_integral(mode, compile_args): a0 = ufl.inner(ufl.jump(ufl.grad(u)), ufl.jump(ufl.grad(v))) * ufl.dS forms = [a0] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) for f, compiled_f in zip(forms, compiled_forms): assert compiled_f.rank == len(f.arguments()) @@ -479,23 +343,19 @@ def test_interior_facet_integral(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [ - [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0], - [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0], - ], - dtype=np_gtype, - ) + coords = np.array([[0.0, 0.0, 0.0, + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0], + [1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 1.0, 1.0, 0.0]], dtype=np_gtype) kernel = getattr(integral0, f"tabulate_tensor_{np_type}") - kernel( - ffi.cast(f"{mode} *", A.ctypes.data), - ffi.cast(f"{mode} *", w.ctypes.data), - ffi.cast(f"{mode} *", c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.cast("int *", facets.ctypes.data), - ffi.cast("uint8_t *", perms.ctypes.data), - ) + kernel(ffi.cast(f'{mode} *', A.ctypes.data), + ffi.cast(f'{mode} *', w.ctypes.data), + ffi.cast(f'{mode} *', c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.cast('int *', facets.ctypes.data), + ffi.cast('uint8_t *', perms.ctypes.data)) @pytest.mark.parametrize("mode", ["double", "double _Complex"]) @@ -503,9 +363,8 @@ def test_conditional(mode, compile_args): element = basix.ufl.element("Lagrange", "triangle", 1) u, v = ufl.TrialFunction(element), ufl.TestFunction(element) x = ufl.SpatialCoordinate(ufl.triangle) - condition = ufl.Or( - ufl.ge(ufl.real(x[0] + x[1]), 0.1), ufl.ge(ufl.real(x[1] + x[1] ** 2), 0.1) - ) + condition = ufl.Or(ufl.ge(ufl.real(x[0] + x[1]), 0.1), + ufl.ge(ufl.real(x[1] + x[1]**2), 0.1)) c1 = ufl.conditional(condition, 2.0, 1.0) a = c1 * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx @@ -516,8 +375,7 @@ def test_conditional(mode, compile_args): forms = [a, b] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) form0 = compiled_forms[0].form_integrals[0] form1 = compiled_forms[1].form_integrals[0] @@ -532,22 +390,15 @@ def test_conditional(mode, compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype - ) + coords = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np_gtype) - kernel0 = ffi.cast( - f"ufcx_tabulate_tensor_{np_type} *", - getattr(form0, f"tabulate_tensor_{np_type}"), - ) - kernel0( - ffi.cast("{type} *".format(type=mode), A1.ctypes.data), - ffi.cast("{type} *".format(type=mode), w1.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel0 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form0, f"tabulate_tensor_{np_type}")) + kernel0(ffi.cast('{type} *'.format(type=mode), A1.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w1.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) expected_result = np.array([[2, -1, -1], [-1, 1, 0], [-1, 0, 1]], dtype=np_type) assert np.allclose(A1, expected_result) @@ -555,18 +406,11 @@ def test_conditional(mode, compile_args): A2 = np.zeros(3, dtype=np_type) w2 = np.array([1.0, 1.0, 1.0], dtype=np_type) - kernel1 = ffi.cast( - f"ufcx_tabulate_tensor_{np_type} *", - getattr(form1, f"tabulate_tensor_{np_type}"), - ) - kernel1( - ffi.cast("{type} *".format(type=mode), A2.ctypes.data), - ffi.cast("{type} *".format(type=mode), w2.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel1 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form1, f"tabulate_tensor_{np_type}")) + kernel1(ffi.cast('{type} *'.format(type=mode), A2.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w2.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) expected_result = np.ones(3, dtype=np_type) assert np.allclose(A2, expected_result) @@ -582,22 +426,11 @@ def test_custom_quadrature(compile_args): points = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.5, 0.5], [0.0, 0.5], [0.5, 0.0]] weights = [1 / 12] * 6 - a = ( - u - * v - * ufl.dx( - metadata={ - "quadrature_rule": "custom", - "quadrature_points": points, - "quadrature_weights": weights, - } - ) - ) + a = u * v * ufl.dx(metadata={"quadrature_rule": "custom", + "quadrature_points": points, "quadrature_weights": weights}) forms = [a] - compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args - ) + compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms(forms, cffi_extra_compile_args=compile_args) ffi = module.ffi form = compiled_forms[0] @@ -607,19 +440,15 @@ def test_custom_quadrature(compile_args): w = np.array([], dtype=np.float64) c = np.array([], dtype=np.float64) - coords = np.array( - [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np.float64 - ) + coords = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np.float64) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel( - ffi.cast("double *", A.ctypes.data), - ffi.cast("double *", w.ctypes.data), - ffi.cast("double *", c.ctypes.data), - ffi.cast("double *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast("double *", A.ctypes.data), + ffi.cast("double *", w.ctypes.data), + ffi.cast("double *", c.ctypes.data), + ffi.cast("double *", coords.ctypes.data), ffi.NULL, ffi.NULL) # Check that A is diagonal assert np.count_nonzero(A - np.diag(np.diagonal(A))) == 0 @@ -631,19 +460,12 @@ def test_curl_curl(compile_args): a = ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx forms = [a] - compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args - ) + compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms(forms, cffi_extra_compile_args=compile_args) -def lagrange_triangle_symbolic( - order, corners=[(1, 0), (2, 0), (0, 1)], fun=lambda i: i -): +def lagrange_triangle_symbolic(order, corners=[(1, 0), (2, 0), (0, 1)], fun=lambda i: i): from sympy import S - - poly_basis = [ - x**i * y**j for i in range(order + 1) for j in range(order + 1 - i) - ] + poly_basis = [x**i * y**j for i in range(order + 1) for j in range(order + 1 - i)] # vertices eval_points = [S(c) for c in corners] # edges @@ -653,39 +475,24 @@ def lagrange_triangle_symbolic( if order > 3: raise NotImplementedError elif order == 3: - eval_points += [ - tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) - for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2] - ] + eval_points += [tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) + for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2]] else: - eval_points += [ - tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) - for i in range(1, order) - ] + eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) + for i in range(1, order)] # face for f in [(0, 1, 2)]: p0 = corners[f[0]] p1 = corners[f[1]] p2 = corners[f[2]] - eval_points += [ - tuple( - S(a) - + sympy.Rational((b - a) * i, order) - + sympy.Rational((c - a) * j, order) - for a, b, c in zip(p0, p1, p2) - ) - for i in range(1, order) - for j in range(1, order - i) - ] - - dual_mat = [ - [f.subs(x, p[0]).subs(y, p[1]) for p in eval_points] for f in poly_basis - ] + eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) + + sympy.Rational((c - a) * j, order) for a, b, c in zip(p0, p1, p2)) + for i in range(1, order) for j in range(1, order - i)] + + dual_mat = [[f.subs(x, p[0]).subs(y, p[1]) for p in eval_points] for f in poly_basis] dual_mat = sympy.Matrix(dual_mat) mat = dual_mat.inv() - functions = [ - sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows) - ] + functions = [sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows)] results = [] for f in functions: integrand = fun(f) @@ -694,14 +501,10 @@ def lagrange_triangle_symbolic( @pytest.mark.parametrize("mode", ["double"]) -@pytest.mark.parametrize( - "sym_fun,ufl_fun", - [ - (lambda i: i, lambda i: i), - (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), - (lambda i: i.diff(y), lambda i: ufl.grad(i)[1]), - ], -) +@pytest.mark.parametrize("sym_fun,ufl_fun", [ + (lambda i: i, lambda i: i), + (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), + (lambda i: i.diff(y), lambda i: ufl.grad(i)[1])]) @pytest.mark.parametrize("order", [1, 2, 3]) def test_lagrange_triangle(compile_args, order, mode, sym_fun, ufl_fun): sym = lagrange_triangle_symbolic(order, fun=sym_fun) @@ -711,8 +514,7 @@ def test_lagrange_triangle(compile_args, order, mode, sym_fun, ufl_fun): a = ufl_fun(v) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) ffi = module.ffi form0 = compiled_forms[0] @@ -727,35 +529,25 @@ def test_lagrange_triangle(compile_args, order, mode, sym_fun, ufl_fun): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [[1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype - ) + coords = np.array([[1.0, 0.0, 0.0], + [2.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np_gtype) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel( - ffi.cast("{type} *".format(type=mode), b.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.NULL, - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('{type} *'.format(type=mode), b.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.NULL, + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) # Check that the result is the same as for sympy assert np.allclose(b, [float(i) for i in sym]) -def lagrange_tetrahedron_symbolic( - order, corners=[(1, 0, 0), (2, 0, 0), (0, 1, 0), (0, 0, 1)], fun=lambda i: i -): +def lagrange_tetrahedron_symbolic(order, corners=[(1, 0, 0), (2, 0, 0), (0, 1, 0), (0, 0, 1)], fun=lambda i: i): from sympy import S - poly_basis = [ - x**i * y**j * z**k - for i in range(order + 1) - for j in range(order + 1 - i) - for k in range(order + 1 - i - j) - ] + x**i * y**j * z**k for i in range(order + 1) for j in range(order + 1 - i) + for k in range(order + 1 - i - j)] # vertices eval_points = [S(c) for c in corners] # edges @@ -765,78 +557,45 @@ def lagrange_tetrahedron_symbolic( if order > 3: raise NotImplementedError elif order == 3: - eval_points += [ - tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) - for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2] - ] + eval_points += [tuple(S(a) + (b - a) * i for a, b in zip(p0, p1)) + for i in [(1 - 1 / sympy.sqrt(5)) / 2, (1 + 1 / sympy.sqrt(5)) / 2]] else: - eval_points += [ - tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) - for i in range(1, order) - ] + eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) for a, b in zip(p0, p1)) + for i in range(1, order)] # face for f in [(1, 2, 3), (0, 2, 3), (0, 1, 3), (0, 1, 2)]: p0 = corners[f[0]] p1 = corners[f[1]] p2 = corners[f[2]] - eval_points += [ - tuple( - S(a) - + sympy.Rational((b - a) * i, order) - + sympy.Rational((c - a) * j, order) - for a, b, c in zip(p0, p1, p2) - ) - for i in range(1, order) - for j in range(1, order - i) - ] + eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) + + sympy.Rational((c - a) * j, order) for a, b, c in zip(p0, p1, p2)) + for i in range(1, order) for j in range(1, order - i)] # interior for v in [(0, 1, 2, 3)]: p0 = corners[v[0]] p1 = corners[v[1]] p2 = corners[v[2]] p3 = corners[v[3]] - eval_points += [ - tuple( - S(a) - + sympy.Rational((b - a) * i, order) - + sympy.Rational((c - a) * j, order) - + sympy.Rational((d - a) * k, order) - for a, b, c, d in zip(p0, p1, p2, p3) - ) - for i in range(1, order) - for j in range(1, order - i) - for k in range(1, order - i - j) - ] - - dual_mat = [ - [f.subs(x, p[0]).subs(y, p[1]).subs(z, p[2]) for p in eval_points] - for f in poly_basis - ] + eval_points += [tuple(S(a) + sympy.Rational((b - a) * i, order) + sympy.Rational((c - a) * j, order) + + sympy.Rational((d - a) * k, order) for a, b, c, d in zip(p0, p1, p2, p3)) + for i in range(1, order) for j in range(1, order - i) for k in range(1, order - i - j)] + + dual_mat = [[f.subs(x, p[0]).subs(y, p[1]).subs(z, p[2]) for p in eval_points] for f in poly_basis] dual_mat = sympy.Matrix(dual_mat) mat = dual_mat.inv() - functions = [ - sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows) - ] + functions = [sum(i * j for i, j in zip(mat.row(k), poly_basis)) for k in range(mat.rows)] results = [] for f in functions: integrand = fun(f) - results.append( - integrand.integrate( - (x, 1 - y - z, 2 - 2 * y - 2 * z), (y, 0, 1 - z), (z, 0, 1) - ) - ) + results.append(integrand.integrate((x, 1 - y - z, 2 - 2 * y - 2 * z), (y, 0, 1 - z), (z, 0, 1))) return results @pytest.mark.parametrize("mode", ["double"]) -@pytest.mark.parametrize( - "sym_fun,ufl_fun", - [ - (lambda i: i, lambda i: i), - (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), - (lambda i: i.diff(y), lambda i: ufl.grad(i)[1]), - ], -) +@pytest.mark.parametrize("sym_fun,ufl_fun", [ + (lambda i: i, lambda i: i), + (lambda i: i.diff(x), lambda i: ufl.grad(i)[0]), + (lambda i: i.diff(y), lambda i: ufl.grad(i)[1])]) @pytest.mark.parametrize("order", [1, 2, 3]) def test_lagrange_tetrahedron(compile_args, order, mode, sym_fun, ufl_fun): sym = lagrange_tetrahedron_symbolic(order, fun=sym_fun) @@ -846,8 +605,7 @@ def test_lagrange_tetrahedron(compile_args, order, mode, sym_fun, ufl_fun): a = ufl_fun(v) * ufl.dx forms = [a] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) ffi = module.ffi form0 = compiled_forms[0] @@ -863,19 +621,16 @@ def test_lagrange_tetrahedron(compile_args, order, mode, sym_fun, ufl_fun): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], dtype=np_gtype - ) + coords = np.array([1.0, 0.0, 0.0, + 2.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0], dtype=np_gtype) kernel = getattr(default_integral, f"tabulate_tensor_{np_type}") - kernel( - ffi.cast("{type} *".format(type=mode), b.ctypes.data), - ffi.cast("{type} *".format(type=mode), w.ctypes.data), - ffi.NULL, - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('{type} *'.format(type=mode), b.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w.ctypes.data), + ffi.NULL, + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) # Check that the result is the same as for sympy assert np.allclose(b, [float(i) for i in sym]) @@ -888,8 +643,7 @@ def test_prism(compile_args): L = v * ufl.dx forms = [L] compiled_forms, module, _ = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": "double"}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': 'double'}, cffi_extra_compile_args=compile_args) ffi = module.ffi form0 = compiled_forms[0] @@ -897,39 +651,18 @@ def test_prism(compile_args): default_integral = form0.form_integrals[0] b = np.zeros(6, dtype=np.float64) - coords = np.array( - [ - 1.0, - 0.0, - 0.0, - 0.0, - 1.0, - 0.0, - 0.0, - 0.0, - 0.0, - 1.0, - 0.0, - 1.0, - 0.0, - 1.0, - 1.0, - 0.0, - 0.0, - 1.0, - ], - dtype=np.float64, - ) + coords = np.array([1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, + 1.0, 0.0, 1.0, + 0.0, 1.0, 1.0, + 0.0, 0.0, 1.0], dtype=np.float64) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel( - ffi.cast("double *", b.ctypes.data), - ffi.NULL, - ffi.NULL, - ffi.cast("double *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('double *', b.ctypes.data), + ffi.NULL, + ffi.NULL, + ffi.cast('double *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.isclose(sum(b), 0.5) @@ -947,8 +680,7 @@ def test_complex_operations(compile_args): forms = [J1, J2] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, options={"scalar_type": mode}, cffi_extra_compile_args=compile_args - ) + forms, options={'scalar_type': mode}, cffi_extra_compile_args=compile_args) form0 = compiled_forms[0].form_integrals[0] form1 = compiled_forms[1].form_integrals[0] @@ -961,48 +693,27 @@ def test_complex_operations(compile_args): geom_type = scalar_to_value_type(mode) np_gtype = cdtype_to_numpy(geom_type) - coords = np.array( - [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np_gtype - ) + coords = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np_gtype) J_1 = np.zeros((1), dtype=np_type) - kernel0 = ffi.cast( - f"ufcx_tabulate_tensor_{np_type} *", - getattr(form0, f"tabulate_tensor_{np_type}"), - ) - kernel0( - ffi.cast("{type} *".format(type=mode), J_1.ctypes.data), - ffi.cast("{type} *".format(type=mode), w1.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) - - expected_result = np.array( - [ - 0.5 - * np.real(w1[0]) - * np.imag(w1[1]) - * (np.real(w1[0]) - 1j * np.imag(w1[0])) - ], - dtype=np_type, - ) + kernel0 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form0, f"tabulate_tensor_{np_type}")) + kernel0(ffi.cast('{type} *'.format(type=mode), J_1.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w1.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) + + expected_result = np.array([0.5 * np.real(w1[0]) * np.imag(w1[1]) + * (np.real(w1[0]) - 1j * np.imag(w1[0]))], dtype=np_type) assert np.allclose(J_1, expected_result) J_2 = np.zeros((1), dtype=np_type) - kernel1 = ffi.cast( - f"ufcx_tabulate_tensor_{np_type} *", - getattr(form1, f"tabulate_tensor_{np_type}"), - ) - kernel1( - ffi.cast("{type} *".format(type=mode), J_2.ctypes.data), - ffi.cast("{type} *".format(type=mode), w1.ctypes.data), - ffi.cast("{type} *".format(type=mode), c.ctypes.data), - ffi.cast(f"{geom_type} *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel1 = ffi.cast(f"ufcx_tabulate_tensor_{np_type} *", getattr(form1, f"tabulate_tensor_{np_type}")) + kernel1(ffi.cast('{type} *'.format(type=mode), J_2.ctypes.data), + ffi.cast('{type} *'.format(type=mode), w1.ctypes.data), + ffi.cast('{type} *'.format(type=mode), c.ctypes.data), + ffi.cast(f'{geom_type} *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.allclose(J_2, expected_result) @@ -1022,8 +733,7 @@ def test_invalid_function_name(compile_args): try: compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args - ) + forms, cffi_extra_compile_args=compile_args) except ValueError: pass except Exception: @@ -1034,17 +744,18 @@ def test_invalid_function_name(compile_args): def test_interval_vertex_quadrature(compile_args): + c_el = basix.ufl.element("Lagrange", "interval", 1, rank=1) mesh = ufl.Mesh(c_el) x = ufl.SpatialCoordinate(mesh) - dx = ufl.Measure("dx", metadata={"quadrature_rule": "vertex"}) + dx = ufl.Measure( + "dx", metadata={"quadrature_rule": "vertex"}) b = x[0] * dx forms = [b] compiled_forms, module, code = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args - ) + forms, cffi_extra_compile_args=compile_args) ffi = module.ffi form0 = compiled_forms[0] @@ -1054,17 +765,15 @@ def test_interval_vertex_quadrature(compile_args): J = np.zeros(1, dtype=np.float64) a = np.pi b = np.exp(1) - coords = np.array([a, 0.0, 0.0, b, 0.0, 0.0], dtype=np.float64) + coords = np.array([a, 0.0, 0.0, + + b, 0.0, 0.0], dtype=np.float64) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel( - ffi.cast("double *", J.ctypes.data), - ffi.NULL, - ffi.NULL, - ffi.cast("double *", coords.ctypes.data), - ffi.NULL, - ffi.NULL, - ) + kernel(ffi.cast('double *', J.ctypes.data), + ffi.NULL, + ffi.NULL, + ffi.cast('double *', coords.ctypes.data), ffi.NULL, ffi.NULL) assert np.isclose(J[0], (0.5 * a + 0.5 * b) * np.abs(b - a)) @@ -1076,8 +785,9 @@ def test_facet_vertex_quadrature(compile_args): mesh = ufl.Mesh(c_el) x = ufl.SpatialCoordinate(mesh) - ds = ufl.Measure("ds", metadata={"quadrature_rule": "vertex"}) - expr = x[0] + ufl.cos(x[1]) + ds = ufl.Measure( + "ds", metadata={"quadrature_rule": "vertex"}) + expr = (x[0] + ufl.cos(x[1])) b1 = expr * ds ds_c = ufl.Measure( "ds", @@ -1085,47 +795,42 @@ def test_facet_vertex_quadrature(compile_args): "quadrature_rule": "custom", "quadrature_points": np.array([[0.0], [1.0]]), "quadrature_weights": np.array([1.0 / 2.0, 1.0 / 2.0]), - }, + } ) b2 = expr * ds_c forms = [b1, b2] compiled_forms, module, _ = ffcx.codegeneration.jit.compile_forms( - forms, cffi_extra_compile_args=compile_args - ) + forms, cffi_extra_compile_args=compile_args) ffi = module.ffi assert len(compiled_forms) == 2 solutions = [] for form in compiled_forms: - assert form.form_integral_offsets[module.lib.exterior_facet + 1] == 1 + offsets = form.form_integral_offsets + exf = module.lib.exterior_facet + assert offsets[exf + 1] - offsets[exf] == 1 - default_integral = form.form_integrals[0] + default_integral = form.form_integrals[offsets[exf]] J = np.zeros(1, dtype=np.float64) a = np.pi b = np.exp(1) - coords = np.array( - [a, 0.1, 0.0, a + b, 0.0, 0.0, a, a, 0.0, a + 2 * b, a, 0.0], - dtype=np.float64, - ) + coords = np.array([a, 0.1, 0.0, + a + b, 0.0, 0.0, + a, a, 0., + a + 2 * b, a, 0.], dtype=np.float64) # First facet is between vertex 0 and 1 in coords facets = np.array([0], dtype=np.intc) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel( - ffi.cast("double *", J.ctypes.data), - ffi.NULL, - ffi.NULL, - ffi.cast("double *", coords.ctypes.data), - ffi.cast("int *", facets.ctypes.data), - ffi.NULL, - ) + kernel(ffi.cast('double *', J.ctypes.data), + ffi.NULL, + ffi.NULL, + ffi.cast('double *', coords.ctypes.data), + ffi.cast('int *', facets.ctypes.data), + ffi.NULL) solutions.append(J[0]) # Test against exact result - assert np.isclose( - J[0], - (0.5 * (a + np.cos(0.1)) + 0.5 * (a + b + np.cos(0))) - * np.sqrt(b**2 + 0.1**2), - ) + assert np.isclose(J[0], (0.5 * (a + np.cos(0.1)) + 0.5 * (a + b + np.cos(0))) * np.sqrt(b**2 + 0.1**2)) # Compare custom quadrature with vertex quadrature assert np.isclose(solutions[0], solutions[1]) @@ -1147,13 +852,12 @@ def test_manifold_derivatives(compile_args): u = ufl.Coefficient(V) d = 5.3 - f_ex = d * order * (order - 1) * x[1] ** (order - 2) + f_ex = d * order * (order - 1) * x[1]**(order - 2) expr = u.dx(1).dx(1) - f_ex J = expr * expr * dx compiled_forms, module, _ = ffcx.codegeneration.jit.compile_forms( - [J], cffi_extra_compile_args=compile_args - ) + [J], cffi_extra_compile_args=compile_args) default_integral = compiled_forms[0].form_integrals[0] scale = 2.5 @@ -1168,13 +872,10 @@ def test_manifold_derivatives(compile_args): ffi = module.ffi J = np.zeros(1, dtype=np.float64) kernel = getattr(default_integral, "tabulate_tensor_float64") - kernel( - ffi.cast("double *", J.ctypes.data), - ffi.cast("double *", w.ctypes.data), - ffi.cast("double *", c.ctypes.data), - ffi.cast("double *", coords.ctypes.data), - ffi.NULL, - ffi.cast("uint8_t *", perm.ctypes.data), - ) + kernel(ffi.cast('double *', J.ctypes.data), + ffi.cast('double *', w.ctypes.data), + ffi.cast('double *', c.ctypes.data), + ffi.cast('double *', coords.ctypes.data), ffi.NULL, + ffi.cast('uint8_t *', perm.ctypes.data)) assert np.isclose(J[0], 0.0)