Example of submesh and DG

General
1

Hello,

I know there are few examples of the use of submesh:

So I created a MWE below where DG is used. We solve the scalar advection domain on a 1D domain (so that the solution simply converges towards the inflow value, here 2.5). The mesh is divided into a left domain and a right domain. Lax-Friedrich flux is used (even at the interface between the two domains).

from scifem import BlockedNewtonSolver
import dolfinx

import numpy as np
import pylab as plt
#
import ufl
#
from mpi4py import MPI
from packaging.version import Version

# Mesh
msh = dolfinx.mesh.create_interval(comm=MPI.COMM_WORLD, points=(0.0, 1.0), nx=4)
gdim = msh.geometry.dim
n = ufl.dot(ufl.FacetNormal(msh), ufl.as_vector([1.0]))

# Convective velocity
velocity_norm = 1.0
velocity = ufl.as_vector([velocity_norm])

# Value at inflow
f_inflow = dolfinx.fem.Constant(msh, 2.5)

# Create submeshes (left and right)
LEFT = 1
RIGHT = 2
cells = dolfinx.mesh.locate_entities(msh, gdim, lambda x: np.full_like(x[0], 1, dtype=np.int8))
midpoints = dolfinx.mesh.compute_midpoints(msh, gdim, cells)
marker = np.full(midpoints.shape[0], -1, dtype=np.int32)
marker[(midpoints[:,0] <= 0.5)] = LEFT
marker[(midpoints[:,0] > 0.5)] = RIGHT
assert(np.all(marker!=-1))
cell_tag = dolfinx.mesh.meshtags(msh, gdim, np.arange(midpoints.shape[0], dtype=np.int32), marker)
dx = ufl.Measure("dx", domain=msh, subdomain_data=cell_tag)

left_mesh, left_mesh_cell_map, left_mesh_vertex_map, _ = dolfinx.mesh.create_submesh(msh, gdim, cell_tag.find(LEFT))
right_mesh, right_mesh_cell_map, right_mesh_vertex_map, _ = dolfinx.mesh.create_submesh(msh, gdim, cell_tag.find(RIGHT))
if (Version(dolfinx.__version__) >= Version('0.10')):
    left_mesh_cell_map = left_mesh_cell_map.sub_topology_to_topology(np.arange(np.where(cell_tag==LEFT)[0].size), False)
    right_mesh_cell_map = right_mesh_cell_map.sub_topology_to_topology(np.arange(np.where(cell_tag==RIGHT)[0].size), False)

parent_to_sub_left = np.full(midpoints.shape[0], -1, dtype=np.int32)
parent_to_sub_left[left_mesh_cell_map] = np.arange(len(left_mesh_cell_map), dtype=np.int32)
parent_to_sub_right = np.full(midpoints.shape[0], -1, dtype=np.int32)
parent_to_sub_right[right_mesh_cell_map] = np.arange(len(right_mesh_cell_map), dtype=np.int32)
entity_maps = {left_mesh: parent_to_sub_left, right_mesh: parent_to_sub_right}

# Defining the function spaces
V_left, V_right = dolfinx.fem.functionspace(left_mesh, ("DG", 1)), dolfinx.fem.functionspace(right_mesh, ("DG", 1))
W = ufl.MixedFunctionSpace(V_left, V_right)
f_left_test, f_right_test = ufl.TestFunctions(W)
f_left_state, f_right_state = dolfinx.fem.Function(V_left), dolfinx.fem.Function(V_right)

# boundary marker
LEFT_BC_MARKER = 1
RIGHT_BC_MARKER = 2
bdry_facets = dolfinx.mesh.locate_entities_boundary(msh, gdim-1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
midpoints = dolfinx.mesh.compute_midpoints(msh, gdim-1, bdry_facets)
bcMarkers = np.full(midpoints.shape[0], -1, dtype='int32')
bcMarkers[(midpoints[:,0] == 0.0)] = LEFT_BC_MARKER
bcMarkers[(midpoints[:,0] == 1.0)] = RIGHT_BC_MARKER
assert(np.where(bcMarkers==LEFT_BC_MARKER)[0].size==1)
fts = dolfinx.mesh.meshtags(msh, gdim-1, bdry_facets, bcMarkers)
ds = ufl.Measure("ds", domain=msh, subdomain_data=fts)

# Inner markers
INTERFACE = 3
inner_facets =  dolfinx.mesh.locate_entities(msh, gdim-1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
midpoints = dolfinx.mesh.compute_midpoints(msh, gdim-1, inner_facets)
inner_markers = np.zeros(midpoints.shape[0], dtype='int32')
inner_markers[(midpoints[:,0] <  0.5)] = LEFT
inner_markers[(midpoints[:,0] >  0.5)] = RIGHT
inner_markers[(midpoints[:,0] ==  0.5)] = INTERFACE
inner_facet_tags = dolfinx.mesh.meshtags(msh, gdim-1, inner_facets, inner_markers)
dS = ufl.Measure("dS", msh, subdomain_data=inner_facet_tags)

# Create from
def flux_LF(f_int, f_ext, n): # Lax Friedrich flux
    return 0.5 * (f_int + f_ext) * velocity_norm * n + 0.5 * velocity_norm * (f_int - f_ext)

F = ufl.inner(-velocity * f_left_state, ufl.grad(f_left_test))*dx(LEFT)\
  + ufl.inner(-velocity * f_right_state, ufl.grad(f_right_test))*dx(RIGHT)

F += ufl.inner(flux_LF(f_left_state('+'), f_left_state('-'), n('+')), f_left_test('+')) * dS(LEFT)
F += ufl.inner(flux_LF(f_left_state('-'), f_left_state('+'), n('-')), f_left_test('-')) * dS(LEFT)
F += ufl.inner(flux_LF(f_right_state('+'), f_right_state('-'), n('+')), f_right_test('+')) * dS(RIGHT)
F += ufl.inner(flux_LF(f_right_state('-'), f_right_state('+'), n('-')), f_right_test('-')) * dS(RIGHT)

F += ufl.inner(flux_LF(f_left_state('+'), f_right_state('-'), n('+')), f_left_test('+')) * dS(INTERFACE)
F += ufl.inner(flux_LF(f_right_test('+'), f_left_state('-'), n('-')), f_right_test('+')) * dS(INTERFACE)

F += ufl.inner(flux_LF(f_left_state, f_inflow, n), f_left_test) * ds(LEFT_BC_MARKER)
F += ufl.inner(flux_LF(f_right_state, dolfinx.fem.Constant(msh, 0.0) , n), f_right_test) * ds(RIGHT_BC_MARKER) # dummy ext value (not used)


F_block = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
J_block = dolfinx.fem.form(ufl.extract_blocks(ufl.derivative(F, [f_left_state, f_right_state], ufl.TrialFunctions(W))), entity_maps=entity_maps)
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
solver = BlockedNewtonSolver(F_block, [f_left_state, f_right_state], J=J_block, petsc_options=petsc_options)
solver.solve()

coords = V_left.tabulate_dof_coordinates()
ids = np.argsort(coords[:,0])
plt.plot(coords[ids,0], f_left_state.x.array[ids], 'b')
coords = V_right.tabulate_dof_coordinates()
ids = np.argsort(coords[:,0])
plt.plot(coords[ids,0], f_right_state.x.array[ids], 'r')
plt.show()

Unfortunately this script fails (dolfinx v0.9) with the error:

Traceback (most recent call last):
File “/stck/fpascal/Python/Code/FENICS/Generic/advection_1D_submesh.py”, line 99, in
F_block = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py”, line 337, in form
return _create_form(form)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py”, line 333, in _create_form
return list(map(lambda sub_form: _create_form(sub_form), form))
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py”, line 333, in
return list(map(lambda sub_form: _create_form(sub_form), form))
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py”, line 331, in _create_form
return _form(form)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py”, line 254, in _form
ufcx_form, module, code = jit.ffcx_jit(
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/jit.py”, line 62, in mpi_jit
return local_jit(*args, **kwargs)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/jit.py”, line 212, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ffcx/codegeneration/jit.py”, line 225, in compile_forms
raise e
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ffcx/codegeneration/jit.py”, line 205, in compile_forms
impl = _compile_objects(
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ffcx/codegeneration/jit.py”, line 330, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ffcx/compiler.py”, line 108, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options[“scalar_type”]) # type: ignore
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ffcx/analysis.py”, line 94, in analyze_ufl_objects
form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ffcx/analysis.py”, line 94, in
form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ffcx/analysis.py”, line 180, in _analyze_form
form_data = ufl.algorithms.compute_form_data(
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ufl/algorithms/compute_form_data.py”, line 427, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ufl/algorithms/check_arities.py”, line 213, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ufl/algorithms/check_arities.py”, line 194, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ufl/corealg/map_dag.py”, line 35, in map_expr_dag
(result,) = map_expr_dags(
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ufl/corealg/map_dag.py”, line 103, in map_expr_dags
r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ufl/algorithms/check_arities.py”, line 60, in sum
f"Adding expressions with non-matching form arguments {_afmt(a)} vs {_afmt(b)}."
File “/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/ufl/algorithms/check_arities.py”, line 20, in _afmt
arg, conj = atuple
ValueError: not enough values to unpack (expected 2, got 1)

If I comment the part with dS(INTERFACE) the error becomes:

Traceback (most recent call last):
  File "/stck/fpascal/Python/Code/FENICS/Generic/advection_1D_submesh.py", line 103, in <module>
    solver.solve()
  File "/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/scifem/solvers.py", line 437, in solve
    n, converged = super().solve(self._x)
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 76, Error in external library

If in addition I comment the ds part the script runs fine (but obviously a wrong solution is obtained).

So I guess I’m not far from having a working script, could you please help to me to correct what is wrong here?

Lucas

2

Here you are hitting an issue that relates to: Discontinuity at interface using mixed domains - #2 by dokken
which has several workarounds required for coupling across the submeshes.

The support for this has been greatly improved in DOLFINx main branch, see for instance:

where this coupling over interior submesh boundaries have been simplified greatly (copmared to my aforementioned implementation from 2024).

3

Thank you for you reply and for the link.

I understand that I must use the function compute_interface_data to modify the EntityMap. However in my case I use DG so I need dS to define as well the inner flux, which was not the case in your example. To this end, following your example I created dGamma that I use to define the flux at the interface (and I still have dS(LEFT) and dS(RIGHT) to define the inner flux). Unfortunately this fails with the error:

Traceback (most recent call last):
  File "/stck/fpascal/Python/Code/FENICS/Generic/advection_1D_submesh.py", line 162, in <module>
    F_block = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
  File "/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 337, in form
    return _create_form(form)
  File "/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 333, in _create_form
    return list(map(lambda sub_form: _create_form(sub_form), form))
  File "/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 333, in <lambda>
    return list(map(lambda sub_form: _create_form(sub_form), form))
  File "/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 331, in _create_form
    return _form(form)
  File "/tmp_user/juno/fpascal/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 249, in _form
    assert all([d is data[0] for d in data])
AssertionError

Should I define a single dS object? If so I do not see how I should set the subdomain_data and subdomain_id (I must admit that I do not fully understand those two parameters).

The new script is given below:

from scifem import BlockedNewtonSolver
import dolfinx
import numpy as np
import pylab as plt
#
import ufl
#
from mpi4py import MPI
from packaging.version import Version

# Mesh
msh = dolfinx.mesh.create_interval(comm=MPI.COMM_WORLD, points=(0.0, 1.0), nx=4)
gdim = msh.geometry.dim
n = ufl.dot(ufl.FacetNormal(msh), ufl.as_vector([1.0]))
num_cells_local = msh.topology.index_map(msh.topology.dim).size_local + msh.topology.index_map(msh.topology.dim).num_ghosts

# Convective velocity
velocity_norm = 1.0
velocity = ufl.as_vector([velocity_norm])

# Value at inflow
f_inflow = dolfinx.fem.Constant(msh, 2.5)

# Create submeshes (left and right)
LEFT = 1
RIGHT = 2
cells = dolfinx.mesh.locate_entities(msh, gdim, lambda x: np.full_like(x[0], 1, dtype=np.int8))
midpoints = dolfinx.mesh.compute_midpoints(msh, gdim, cells)
marker = np.full(midpoints.shape[0], -1, dtype=np.int32)
marker[(midpoints[:,0] <= 0.5)] = LEFT
marker[(midpoints[:,0] > 0.5)] = RIGHT
assert(np.all(marker!=-1))
cell_tag = dolfinx.mesh.meshtags(msh, gdim, np.arange(midpoints.shape[0], dtype=np.int32), marker)
dx = ufl.Measure("dx", domain=msh, subdomain_data=cell_tag)

left_mesh, left_mesh_cell_map, left_mesh_vertex_map, _ = dolfinx.mesh.create_submesh(msh, gdim, cell_tag.find(LEFT))
right_mesh, right_mesh_cell_map, right_mesh_vertex_map, _ = dolfinx.mesh.create_submesh(msh, gdim, cell_tag.find(RIGHT))
# if (Version(dolfinx.__version__) >= Version('0.10')):
    # left_mesh_cell_map = left_mesh_cell_map.sub_topology_to_topology(np.arange(np.where(cell_tag==LEFT)[0].size), False)
    # right_mesh_cell_map = right_mesh_cell_map.sub_topology_to_topology(np.arange(np.where(cell_tag==RIGHT)[0].size), False)


# Defining the function spaces
V_left, V_right = dolfinx.fem.functionspace(left_mesh, ("DG", 1)), dolfinx.fem.functionspace(right_mesh, ("DG", 1))
W = ufl.MixedFunctionSpace(V_left, V_right)
f_left_test, f_right_test = ufl.TestFunctions(W)
f_left_state, f_right_state = dolfinx.fem.Function(V_left), dolfinx.fem.Function(V_right)

# boundary marker
LEFT_BC_MARKER = 1
RIGHT_BC_MARKER = 2
bdry_facets = dolfinx.mesh.locate_entities_boundary(msh, gdim-1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
midpoints = dolfinx.mesh.compute_midpoints(msh, gdim-1, bdry_facets)
bcMarkers = np.full(midpoints.shape[0], -1, dtype='int32')
bcMarkers[(midpoints[:,0] == 0.0)] = LEFT_BC_MARKER
bcMarkers[(midpoints[:,0] == 1.0)] = RIGHT_BC_MARKER
assert(np.where(bcMarkers==LEFT_BC_MARKER)[0].size==1)
fts = dolfinx.mesh.meshtags(msh, gdim-1, bdry_facets, bcMarkers)
ds = ufl.Measure("ds", domain=msh, subdomain_data=fts)

# Inner markers
INTERFACE = 3
inner_facets =  dolfinx.mesh.locate_entities(msh, gdim-1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
midpoints = dolfinx.mesh.compute_midpoints(msh, gdim-1, inner_facets)
inner_markers = np.zeros(midpoints.shape[0], dtype='int32')
inner_markers[(midpoints[:,0] ==  0.0)] = -1
inner_markers[(midpoints[:,0] ==  1.0)] = -1
inner_markers[(midpoints[:,0] <  0.5)] = LEFT
inner_markers[(midpoints[:,0] >  0.5)] = RIGHT
inner_markers[(midpoints[:,0] ==  0.5)] = INTERFACE
inner_facet_tags = dolfinx.mesh.meshtags(msh, gdim-1, inner_facets, inner_markers)
dS = ufl.Measure("dS", msh, subdomain_data=inner_facet_tags)

# Create integration measure for interface
# Interior marker is considered as ("+") restriction
def compute_interface_data(
    cell_tags: dolfinx.mesh.MeshTags, facet_indices: np.typing.NDArray[np.int32]
) -> np.typing.NDArray[np.int32]:
    """
    Compute interior facet integrals that are consistently ordered according to the `cell_tags`,
    such that the data `(cell0, facet_idx0, cell1, facet_idx1)` is ordered such that
    `cell_tags[cell0]`<`cell_tags[cell1]`, i.e the cell with the lowest cell marker is considered the
    "+" restriction".

    Args:
        cell_tags: MeshTags that must contain an integer marker for all cells adjacent to the `facet_indices`
        facet_indices: List of facets (local index) that are on the interface.
    Returns:
        The integration data.
    """
    # Future compatibilty check
    integration_args: tuple[int] | tuple
    if Version("0.10.0") <= Version(dolfinx.__version__):
        integration_args = ()
    else:
        fdim = cell_tags.dim - 1
        integration_args = (fdim,)
    idata = dolfinx.cpp.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet,
        cell_tags.topology,
        facet_indices,
        *integration_args,
    )
    ordered_idata = idata.reshape(-1, 4).copy()
    switch = (
        cell_tags.values[ordered_idata[:, 0]]
        > cell_tags.values[ordered_idata[:, 2]]
    )
    if True in switch:
        ordered_idata[switch, :] = ordered_idata[switch][:, [2, 3, 0, 1]]
    return ordered_idata

ordered_integration_data = compute_interface_data(cell_tag, inner_facet_tags.find(INTERFACE))

gamma_marker = 83
dGamma = ufl.Measure(
        "dS",
        domain=msh,
        subdomain_data=[(gamma_marker, ordered_integration_data.flatten())],
        subdomain_id=gamma_marker,
    )

parent_cells_plus = ordered_integration_data[:, 0]
parent_cells_minus = ordered_integration_data[:, 2]
parent_to_sub_left = np.full(num_cells_local, -1, dtype=np.int32)
parent_to_sub_left[left_mesh_cell_map] = np.arange(len(left_mesh_cell_map), dtype=np.int32)
parent_to_sub_right = np.full(num_cells_local, -1, dtype=np.int32)
parent_to_sub_right[right_mesh_cell_map] = np.arange(len(right_mesh_cell_map), dtype=np.int32)

entity_maps = {
    left_mesh: parent_to_sub_left,
    right_mesh: parent_to_sub_right,
}
entity_maps[left_mesh][parent_cells_minus] = entity_maps[left_mesh][parent_cells_plus]
entity_maps[right_mesh][parent_cells_plus] = entity_maps[right_mesh][parent_cells_minus]

# Create from
def flux_LF(f_int, f_ext, n): # Lax Friedrich flux
    return 0.5 * (f_int + f_ext) * velocity_norm * n + 0.5 * velocity_norm * (f_int - f_ext)

F = ufl.inner(-velocity * f_left_state, ufl.grad(f_left_test))*dx(LEFT)\
  + ufl.inner(-velocity * f_right_state, ufl.grad(f_right_test))*dx(RIGHT)

F += ufl.inner(flux_LF(f_left_state('+'), f_left_state('-'), n('+')), f_left_test('+')) * dS(LEFT)
F += ufl.inner(flux_LF(f_left_state('-'), f_left_state('+'), n('-')), f_left_test('-')) * dS(LEFT)
F += ufl.inner(flux_LF(f_right_state('+'), f_right_state('-'), n('+')), f_right_test('+')) * dS(RIGHT)
F += ufl.inner(flux_LF(f_right_state('-'), f_right_state('+'), n('-')), f_right_test('-')) * dS(RIGHT)

F += ufl.inner(flux_LF(f_left_state('+'), f_right_state('-'), n('+')), f_left_test('+')) * dGamma
F += ufl.inner(flux_LF(f_right_test('+'), f_left_state('-'), n('-')), f_right_test('+')) * dGamma
# F += ufl.inner(flux_LF(f_left_state('+'), f_right_state('-'), n('+')), f_left_test('+')) * dS(INTERFACE)
# F += ufl.inner(flux_LF(f_right_test('+'), f_left_state('-'), n('-')), f_right_test('+')) * dS(INTERFACE)

F += ufl.inner(flux_LF(f_left_state, f_inflow, n), f_left_test) * ds(LEFT_BC_MARKER)
F += ufl.inner(flux_LF(f_right_state, dolfinx.fem.Constant(msh, 0.0) , n), f_right_test) * ds(RIGHT_BC_MARKER) # dummy ext value (not used)


F_block = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
J_block = dolfinx.fem.form(ufl.extract_blocks(ufl.derivative(F, [f_left_state, f_right_state], ufl.TrialFunctions(W))), entity_maps=entity_maps)
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
solver = BlockedNewtonSolver(F_block, [f_left_state, f_right_state], J=J_block, petsc_options=petsc_options)
solver.solve()

coords = V_left.tabulate_dof_coordinates()
ids = np.argsort(coords[:,0])
plt.plot(coords[ids,0], f_left_state.x.array[ids], 'b')
coords = V_right.tabulate_dof_coordinates()
ids = np.argsort(coords[:,0])
plt.plot(coords[ids,0], f_right_state.x.array[ids], 'r')
plt.show()