mscroggs
commented
1 year ago Any chance this is fixed by https://github.com/FEniCS/ffcx/pull/622?
Closed jorgensd closed 1 year ago
mscroggs
commented
1 year ago Any chance this is fixed by https://github.com/FEniCS/ffcx/pull/622?
jorgensd
commented
1 year ago Any chance this is fixed by FEniCS/ffcx#622?
It does not:
In [1]: import basix.ufl
...: import dolfinx
...: import ufl
...: from mpi4py import MPI
...:
...: mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
...: V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
...: u, v = dolfinx.fem.Function(V), ufl.TestFunction(V)
...: q_el = basix.ufl.quadrature_element(mesh.topology.cell_name(),value_shape=(),
...: scheme="default", degree=1)
...: W = dolfinx.fem.functionspace(mesh, q_el)
...: c = dolfinx.fem.Function(W)
...: a = c * ufl.dot(u, v) * ufl.dx(domain=mesh,
...: metadata={"quadrature_degree": q_el.degree})
...: dolfinx.fem.form(a)
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
Cell In[1], line 15
12 c = dolfinx.fem.Function(W)
13 a = c * ufl.dot(u, v) * ufl.dx(domain=mesh,
14 metadata={"quadrature_degree": q_el.degree})
---> 15 dolfinx.fem.form(a)
File ~/shared/dolfinx/python/dolfinx/fem/forms.py:188, in form(form, dtype, form_compiler_options, jit_options)
185 return list(map(lambda sub_form: _create_form(sub_form), form))
186 return form
--> 188 return _create_form(form)
File ~/shared/dolfinx/python/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
180 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
181 return form argument"""
182 if isinstance(form, ufl.Form):
--> 183 return _form(form)
184 elif isinstance(form, collections.abc.Iterable):
185 return list(map(lambda sub_form: _create_form(sub_form), form))
File ~/shared/dolfinx/python/dolfinx/fem/forms.py:141, in form.<locals>._form(form)
139 if mesh is None:
140 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 141 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
142 form_compiler_options=form_compiler_options,
143 jit_options=jit_options)
145 # For each argument in form extract its function space
146 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File ~/shared/dolfinx/python/dolfinx/jit.py:56, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
51 @functools.wraps(local_jit)
52 def mpi_jit(comm, *args, **kwargs):
53
54 # Just call JIT compiler when running in serial
55 if comm.size == 1:
---> 56 return local_jit(*args, **kwargs)
58 # Default status (0 == ok, 1 == fail)
59 status = 0
File ~/shared/dolfinx/python/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
202 # Switch on type and compile, returning cffi object
203 if isinstance(ufl_object, ufl.Form):
--> 204 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
205 elif isinstance(ufl_object, ufl.AbstractFiniteElement):
206 r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)
File ~/shared/ffcx/ffcx/codegeneration/jit.py:199, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
197 except Exception:
198 pass
--> 199 raise e
201 obj, module = _load_objects(cache_dir, module_name, form_names)
202 return obj, module, (decl, impl)
File ~/shared/ffcx/ffcx/codegeneration/jit.py:190, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
187 for name in form_names:
188 decl += form_template.format(name=name)
--> 190 impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
191 cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
192 except Exception as e:
193 try:
194 # remove c file so that it will not timeout next time
File ~/shared/ffcx/ffcx/codegeneration/jit.py:260, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
256 import ffcx.compiler
258 # JIT uses module_name as prefix, which is needed to make names of all struct/function
259 # unique across modules
--> 260 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
262 ffibuilder = cffi.FFI()
263 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
264 extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)
File ~/shared/ffcx/ffcx/compiler.py:97, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise)
95 # Stage 1: analysis
96 cpu_time = time()
---> 97 analysis = analyze_ufl_objects(ufl_objects, options)
98 _print_timing(1, time() - cpu_time)
100 # Stage 2: intermediate representation
File ~/shared/ffcx/ffcx/analysis.py:86, in analyze_ufl_objects(ufl_objects, options)
83 else:
84 raise TypeError("UFL objects not recognised.")
---> 86 form_data = tuple(_analyze_form(form, options) for form in forms)
87 for data in form_data:
88 elements += data.unique_sub_elements
File ~/shared/ffcx/ffcx/analysis.py:86, in <genexpr>(.0)
83 else:
84 raise TypeError("UFL objects not recognised.")
---> 86 form_data = tuple(_analyze_form(form, options) for form in forms)
87 for data in form_data:
88 elements += data.unique_sub_elements
File ~/shared/ffcx/ffcx/analysis.py:160, in _analyze_form(form, options)
157 complex_mode = "_Complex" in options["scalar_type"]
159 # Compute form metadata
--> 160 form_data = ufl.algorithms.compute_form_data(
161 form,
162 do_apply_function_pullbacks=True,
163 do_apply_integral_scaling=True,
164 do_apply_geometry_lowering=True,
165 preserve_geometry_types=(ufl.classes.Jacobian,),
166 do_apply_restrictions=True,
167 do_append_everywhere_integrals=False, # do not add dx integrals to dx(i) in UFL
168 complex_mode=complex_mode)
170 # Determine unique quadrature degree and quadrature scheme
171 # per each integral data
172 for id, integral_data in enumerate(form_data.integral_data):
173 # Iterate through groups of integral data. There is one integral
174 # data for all integrals with same domain, itype, subdomain_id
(...)
178 # all integrals in this integral data group, i.e. must be the
179 # same for for the same (domain, itype, subdomain_id)
File ~/shared/ufl/ufl/algorithms/compute_form_data.py:281, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
277 # Estimate polynomial degree of integrands now, before applying
278 # any pullbacks and geometric lowering. Otherwise quad degrees
279 # blow up horrifically.
280 if do_estimate_degrees:
--> 281 form = attach_estimated_degrees(form)
283 if do_apply_function_pullbacks:
284 # Rewrite coefficients and arguments in terms of their
285 # reference cell values with Piola transforms and symmetry
(...)
288 # Domain. Current dolfin works if Expression has a
289 # cell but this should be changed to a mesh.
290 form = apply_function_pullbacks(form)
File ~/shared/ufl/ufl/algorithms/compute_form_data.py:207, in attach_estimated_degrees(form)
205 md = {}
206 md.update(integral.metadata())
--> 207 degree = estimate_total_polynomial_degree(integral.integrand())
208 md["estimated_polynomial_degree"] = degree
209 new_integrals.append(integral.reconstruct(metadata=md))
File ~/shared/ufl/ufl/algorithms/estimate_degrees.py:380, in estimate_total_polynomial_degree(e, default_degree, element_replace_map)
378 degrees = map_expr_dags(de, [e.integrand()])
379 else:
--> 380 degrees = map_expr_dags(de, [e])
381 degree = max(degrees) if degrees else default_degree
382 return degree
File ~/shared/ufl/ufl/corealg/map_dag.py:98, in map_expr_dags(function, expressions, compress, vcache, rcache)
96 # Cache miss: Get transformed operands, then apply transformation
97 if cutoff_types[v._ufl_typecode_]:
---> 98 r = handlers[v._ufl_typecode_](v)
99 else:
100 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File ~/shared/ufl/ufl/algorithms/estimate_degrees.py:87, in SumDegreeEstimator.coefficient(self, v)
85 e = v.ufl_element()
86 e = self.element_replace_map.get(e, e)
---> 87 d = e.embedded_superdegree # FIXME: Use component to improve accuracy for mixed elements
88 if d is None:
89 d = self.default_degree
File ~/shared/basix/python/basix/ufl.py:1540, in _QuadratureElement.embedded_superdegree(self)
1527 @property
1528 def embedded_superdegree(self) -> int:
1529 """Return the degree of the minimum degree Lagrange space that spans this element.
1530
1531 This returns the degree of the lowest degree Lagrange space such that the polynomial
(...)
1538 Lagrange space includes the degree 2 polynomial xy.
1539 """
-> 1540 raise NotImplementedError()
NotImplementedError:
Quadrature elements are currently broken with the main branches (works with 0.7.1 dolfinx/dolfinx image)
This goes into the related issue that interpolation with quadrature elements is not working: https://github.com/FEniCS/dolfinx/issues/1546 as the Quadrature elements are not real Basix-elements.