# Second deriative on a submesh **URL:** \ **Category:** General\ **Created:** [May 22, 2024, 9:57am UTC](https://fenicsproject.discourse.group/t/second-deriative-on-a-submesh/14757 "2024-05-22T09:57:31Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![pierremrtt](https://avatars.discourse-cdn.com/v4/letter/p/e0b2c6/32.png) [@pierremrtt](https://fenicsproject.discourse.group/u/pierremrtt)\ **Post date:** [May 22, 2024, 9:57am UTC](https://fenicsproject.discourse.group/t/second-deriative-on-a-submesh/14757/1 "2024-05-22T09:57:32Z")
Dear Community, I have an issue with the second derivative of a trial function over a submesh. My project is around 10k code lines, and it is not really easy to explain you the context in a simple post. But basically, I have a differential equation on a global domain where an integral over a surface appears, in this equation p is the main function and q is the secondary function that appears under the integral over the surface, this surface belongs to the global domain, and is also considered as a subdomain. I have another equation relative to p and q that is only defined along the subdomain. So far, I’ve implemented the order 1 and 2 of my operator succesfuly, at the order 3, long story short, I need to have the second derivate of q along the subdomain, and that’s where problems appear. Here is the system for the 2nd order, I didn’t take the time to write the 3rd one with latex ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/6/6fa2da7cf4b520b0a6c5694e4bfed384c9297b96.png) I’ve build a code as small as possible to reproduce the error, this code doesn’t correspond to the problem, but just reaches the error: # Needed modules ```auto from mpi4py import MPI import gmsh from dolfinx.io import gmshio import dolfinx.mesh as msh import numpy as np from dolfinx import plot from dolfinx.fem import FunctionSpace, form from ufl import Measure, TrialFunction, TestFunction, grad, inner ``` # Geometry with gmsh ```auto gmsh.initialize() comm = MPI.COMM_WORLD model_rank = 0 model = gmsh.model() gmsh.model.add("test") side_box = 1 lc = 1e-1 # Definition of the points p1 = gmsh.model.geo.addPoint(0, 0, 0, lc) p2 = gmsh.model.geo.addPoint(side_box, 0, 0, lc) p3 = gmsh.model.geo.addPoint(side_box, side_box, 0, lc) p4 = gmsh.model.geo.addPoint(0, side_box, 0, lc) l1 = gmsh.model.geo.addLine(p1, p2) l2 = gmsh.model.geo.addLine(p2, p3) l3 = gmsh.model.geo.addLine(p3, p4) l4 = gmsh.model.geo.addLine(p4, p1) cl1 = [l1, l2, l3, l4] s1 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl1)]) gmsh.model.geo.synchronize() gmsh.model.addPhysicalGroup(1, [l1], tag=1) gmsh.model.addPhysicalGroup(2, [s1], tag=1) gmsh.model.mesh.generate(2) final_mesh, cell_tags, facet_tags = gmshio.model_to_mesh(model, comm, model_rank) gmsh.finalize() tdim = final_mesh.topology.dim fdim = tdim - 1 submesh, entity_map = msh.create_submesh(final_mesh, fdim, facet_tags.find(1))[0:2] ``` # Rest of the code ```auto deg = 3 family = "Lagrange" P = FunctionSpace(final_mesh, (family, deg)) Q = FunctionSpace(submesh, (family, deg)) p, q = TrialFunction(P), TrialFunction(Q) v, u = TestFunction(P), TestFunction(Q) dx = Measure("dx", domain=final_mesh, subdomain_data=cell_tags) ds = Measure("ds", domain=final_mesh, subdomain_data=facet_tags) ddqx = grad(grad(q)[0])[0] e0x = inner(ddqx, u)*ds(1) a = form(e0x) ``` Here is the error (Very long) # Error ```auto --------------------------------------------------------------------------- IndexError Traceback (most recent call last) Cell In[4], line 3 1 ddqx = grad(grad(q)[0])[0] 2 e0x = inner(ddqx, u)*ds(1) ----> 3 a = form(e0x) File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py:194, in form(form, dtype, form_compiler_options, jit_options, entity_maps) 191 return list(map(lambda sub_form: _create_form(sub_form), form)) 192 return form --> 194 return _create_form(form) File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py:189, in form.._create_form(form) 186 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise 187 return form argument""" 188 if isinstance(form, ufl.Form): --> 189 return _form(form) 190 elif isinstance(form, collections.abc.Iterable): 191 return list(map(lambda sub_form: _create_form(sub_form), form)) File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py:152, in form.._form(form) 150 if mesh is None: 151 raise RuntimeError("Expecting to find a Mesh in the form.") --> 152 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form, 153 form_compiler_options=form_compiler_options, 154 jit_options=jit_options) 156 # For each argument in form extract its function space 157 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()] File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/jit.py:56, in mpi_jit_decorator..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 /usr/local/dolfinx-real/lib/python/dist-packages/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.FiniteElementBase): 206 r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit) File /usr/local/lib/python3.10/dist-packages/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 /usr/local/lib/python3.10/dist-packages/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 /usr/local/lib/python3.10/dist-packages/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 /usr/local/lib/python3.10/dist-packages/ffcx/compiler.py:102, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise) 100 # Stage 2: intermediate representation 101 cpu_time = time() --> 102 ir = compute_ir(analysis, object_names, prefix, options, visualise) 103 _print_timing(2, time() - cpu_time) 105 # Stage 3: code generation File /usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py:198, in compute_ir(analysis, object_names, prefix, options, visualise) 192 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names) 193 for e in analysis.unique_elements] 195 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names) 196 for e in analysis.unique_elements] --> 198 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names, 199 options, visualise) 200 for (i, fd) in enumerate(analysis.form_data)] 201 ir_integrals = list(itertools.chain(*irs)) 203 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers, 204 finite_element_names, dofmap_names, object_names) 205 for (i, fd) in enumerate(analysis.form_data)] File /usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py:198, in (.0) 192 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names) 193 for e in analysis.unique_elements] 195 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names) 196 for e in analysis.unique_elements] --> 198 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names, 199 options, visualise) 200 for (i, fd) in enumerate(analysis.form_data)] 201 ir_integrals = list(itertools.chain(*irs)) 203 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers, 204 finite_element_names, dofmap_names, object_names) 205 for (i, fd) in enumerate(analysis.form_data)] File /usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py:493, in _compute_integral_ir(form_data, form_index, element_numbers, integral_names, finite_element_names, options, visualise) 490 integrands = {rule: integral.integrand() for rule, integral in sorted_integrals.items()} 492 # Build more specific intermediate representation --> 493 integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type, 494 ir["entitytype"], integrands, ir["tensor_shape"], 495 ir["mixed_dim"], options, visualise) 497 ir.update(integral_ir) 499 # Fetch name File /usr/local/lib/python3.10/dist-packages/ffcx/ir/integral.py:84, in compute_integral_ir(cell, integral_type, entitytype, integrands, argument_shape, is_mixed_dim, p, visualise) 74 # Build terminal_data from V here before factorization. Then we 75 # can use it to derive table properties for all modified 76 # terminals, and then use that to rebuild the scalar graph more 77 # efficiently before argument factorization. We can build 78 # terminal_data again after factorization if that's necessary. 80 initial_terminals = {i: analyse_modified_terminal(v['expression']) 81 for i, v in S.nodes.items() 82 if is_modified_terminal(v['expression'])} ---> 84 mt_table_reference = build_optimized_tables( 85 quadrature_rule, 86 cell, 87 integral_type, 88 entitytype, 89 initial_terminals.values(), 90 ir["unique_tables"], 91 is_mixed_dim, 92 rtol=p["table_rtol"], 93 atol=p["table_atol"]) 95 # Fetch unique tables for this quadrature rule 96 table_types = {v.name: v.ttype for v in mt_table_reference.values()} File /usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py:351, in build_optimized_tables(quadrature_rule, cell, integral_type, entitytype, modified_terminals, existing_tables, is_mixed_dim, rtol, atol) 349 new_table = [] 350 for ref in range(2): --> 351 new_table.append(get_ffcx_table_values( 352 permute_quadrature_interval(quadrature_rule.points, ref), cell, 353 integral_type, element, avg, entitytype, local_derivatives, flat_component, is_mixed_dim)) 355 t = new_table[0] 356 t['array'] = np.vstack([td['array'] for td in new_table]) File /usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py:144, in get_ffcx_table_values(points, cell, integral_type, element, avg, entitytype, derivative_counts, flat_component, is_mixed_dim) 141 raise RuntimeError(f"Domain cell and ufl element cell not compatible for {integral_type} integral") 143 tbl = component_element.tabulate(deriv_order, entity_points) --> 144 tbl = tbl[basix_index(derivative_counts)] 145 component_tables.append(tbl) 147 if avg in ("cell", "facet"): 148 # Compute numeric integral of the each component table IndexError: index 4 is out of bounds for axis 0 with size 3 ``` Within my project, in 3D with a 2D subdomain, the last line of the error is ```auto IndexError: index 7 is out of bounds for axis 0 with size 6 ``` I’ve been very deep in FEniCSx files in order to understand what’s going on and how to solve it. But I have no clue.. Once again, the second order derivative of p (And all derivatives of p) doesn’t provide any error, the first derivative of q doesn’t neither. Thank’s a billion for your help, Regards, Pierre. ---
**Author:** ![pierremrtt](https://avatars.discourse-cdn.com/v4/letter/p/e0b2c6/32.png) [@pierremrtt](https://fenicsproject.discourse.group/u/pierremrtt)\ **Post date:** [May 24, 2024, 12:13pm UTC](https://fenicsproject.discourse.group/t/second-deriative-on-a-submesh/14757/2 "2024-05-24T12:13:34Z")
Hi Folks, Did anyone had chance to get into this? Regards, ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 24, 2024, 12:49pm UTC](https://fenicsproject.discourse.group/t/second-deriative-on-a-submesh/14757/3 "2024-05-24T12:49:51Z")
Dear Pierre, As far as I can tell, you are looking for the support sketched out in: [Add support for mixed-domain problems of codimension 1 by jpdean · Pull Request #3224 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/pull/3224) by @jpdean . I.e. a mixed assembly of spaces defined on submeshes of different co-dimensions.