Mixed formulation of Poisson problem in dolfinx

dolfinx
1

Hi,
I am trying to rewrite the Mixed Poisson problem using dolfinx (fenicsx). Attached is my unsuccessful attempt. Any help is highly appreciated. Two issues:

  1. The BDM element does not seem to work, therefore using “CG” element in this code.
  2. Unable to properly apply the Dirichlet BC.
from mpi4py import MPI
from dolfinx import fem, mesh, plot
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem


domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                 nx=32, ny=32,
                                 cell_type=mesh.CellType.triangle)
P1 = ufl.FiniteElement('CG', domain.ufl_cell(), 1)
P2 = ufl.FiniteElement('DG', domain.ufl_cell(), 0)
element = ufl.MixedElement([P1, P2])
V = fem.FunctionSpace(domain, element)

# test and trial functions
tau, v = ufl.TestFunctions(V)
sigma, u = fem.Function(V).split()
u.name = "u"


# function determining if a node is on the tray top
def on_top_boundary(x):
    return np.logical_or(np.isclose(x[0], 0),
                         np.isclose(x[0], 1))


V0, submap = V.sub(0).collapse()

# determine boundary DOFs
boundary_dofs = fem.locate_dofs_geometrical((V.sub(0), V0), on_top_boundary)

# apply dirichlet BC to boundary DOFs
C_0 = 0.0
bc = fem.dirichletbc(ScalarType(C_0), boundary_dofs[0], V.sub(0))


x = ufl.SpatialCoordinate(domain=domain)
f = 10.0 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = ufl.sin(5 * x[0])
a = (ufl.dot(sigma, tau) + ufl.div(tau) * u + ufl.div(sigma) * v) * ufl.dx
L = -ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
2

What do you mean by this? Do you get an error?

Your DirichletBC is scalar valued. The BDM space is a vector valued function, and thus the boundary condition should be a vector.

Here is how I would approach the start of your cpde with BDM:

from mpi4py import MPI
from dolfinx import fem, mesh, plot
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem


domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                 nx=32, ny=32,
                                 cell_type=mesh.CellType.triangle)
P1 = ufl.FiniteElement('BDM', domain.ufl_cell(), 1)
P2 = ufl.FiniteElement('DG', domain.ufl_cell(), 0)
element = ufl.MixedElement([P1, P2])
V = fem.FunctionSpace(domain, element)

# test and trial functions
tau, v = ufl.TestFunctions(V)
sigma, u = ufl.TrialFunctions(V)


# function determining if a node is on the tray top
def on_top_boundary(x):
    return np.logical_or(np.isclose(x[0], 0),
                         np.isclose(x[0], 1))


V0, submap = V.sub(0).collapse()

# determine boundary DOFs
boundary_facets = mesh.locate_entities_boundary(
    domain, domain.topology.dim-1, on_top_boundary)
boundary_dofs = fem.locate_dofs_topological(
    (V.sub(0), V0), domain.topology.dim-1, boundary_facets)

Note that I am using trial-functions to define sigma and u, and locate_dofs_topological, as the BDM dofs are not defined as point evaluations: https://defelement.com/elements/examples/triangle-N2div-1.html

3

Dear @dokken,
I appreciate your quick help.
I tried applying a vector valued BC on the dofs (hope that is the correct way).
The error trace throws a type error now.

from mpi4py import MPI
from dolfinx import fem, mesh, plot
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem
import pyvista as pv


domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                 nx=32, ny=32,
                                 cell_type=mesh.CellType.triangle)
P1 = ufl.FiniteElement('BDM', domain.ufl_cell(), 1)
P2 = ufl.FiniteElement('DG', domain.ufl_cell(), 0)
element = ufl.MixedElement([P1, P2])
V = fem.FunctionSpace(domain, element)

# test and trial functions
tau, v = ufl.TestFunctions(V)
sigma, u = ufl.TrialFunctions(V)


# function determining if a node is on the tray top
def on_top_boundary(x):
    return np.logical_or(np.isclose(x[0], 0),
                         np.isclose(x[0], 1))


V0, submap = V.sub(0).collapse()

# determine boundary DOFs
boundary_facets = mesh.locate_entities_boundary(domain, domain.topology.dim-1, on_top_boundary)
boundary_dofs = fem.locate_dofs_topological((V.sub(0), V0), domain.topology.dim-1, boundary_facets)

# apply Dirichlet (vector type) on identified dofs
u_zero = np.array((0,) * domain.geometry.dim, dtype=ScalarType)  # vector type
bc = fem.dirichletbc(value=u_zero, dofs=boundary_dofs, V=V.sub(0))  # V.sub(0): Corresponds to BDM element

# apply Neumann BC
x = ufl.SpatialCoordinate(domain=domain)
f = 10.0 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = ufl.sin(5 * x[0])
a = (ufl.dot(sigma, tau) + ufl.div(tau) * u + ufl.div(sigma) * v) * ufl.dx
L = -ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

# Solve
problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

Error occurs at line starting with “bc = fem.dirichletbc(…”

Traceback (most recent call last):
  File "/Users/anil/Apps/anaconda3/envs/fenics-env/lib/python3.8/site-packages/dolfinx/fem/bcs.py", line 125, in __init__
    super().__init__(_value, dofs, V)  # type: ignore
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)
Invoked with: array([0., 0.]), [array([   0,    1,   17,   18,   40,   41,   71,   72,  110,  111,  157,
        158,  212,  213,  275,  276,  346,  347,  425,  426,  512,  513,
        607,  608,  710,  711,  821,  822,  940,  941, 1067, 1068, 1202,
       1203, 1345, 1346, 1496, 1497, 1655, 1656, 1822, 1823, 1997, 1998,
       2180, 2181, 2371, 2372, 2570, 2571, 2777, 2778, 2992, 2993, 3215,
       3216, 3446, 3447, 3685, 3686, 3932, 3933, 4187, 4188, 4194, 4195,
       4447, 4448, 4692, 4693, 4929, 4930, 5158, 5159, 5379, 5380, 5592,
       5593, 5797, 5798, 5994, 5995, 6183, 6184, 6364, 6365, 6537, 6538,
       6702, 6703, 6859, 6860, 7008, 7009, 7149, 7150, 7282, 7283, 7407,
       7408, 7524, 7525, 7633, 7634, 7734, 7735, 7827, 7828, 7912, 7913,
       7989, 7990, 8058, 8059, 8119, 8120, 8172, 8173, 8217, 8218, 8254,
       8255, 8283, 8284, 8304, 8305, 8317, 8318], dtype=int32), array([   0,    1,   14,   15,   32,   33,   56,   57,   86,   87,  122,
        123,  164,  165,  212,  213,  266,  267,  326,  327,  392,  393,
        464,  465,  542,  543,  626,  627,  716,  717,  812,  813,  914,
        915, 1022, 1023, 1136, 1137, 1256, 1257, 1382, 1383, 1514, 1515,
       1652, 1653, 1796, 1797, 1946, 1947, 2102, 2103, 2264, 2265, 2432,
       2433, 2606, 2607, 2786, 2787, 2972, 2973, 3164, 3165, 3170, 3171,
       3360, 3361, 3544, 3545, 3722, 3723, 3894, 3895, 4060, 4061, 4220,
       4221, 4374, 4375, 4522, 4523, 4664, 4665, 4800, 4801, 4930, 4931,
       5054, 5055, 5172, 5173, 5284, 5285, 5390, 5391, 5490, 5491, 5584,
       5585, 5672, 5673, 5754, 5755, 5830, 5831, 5900, 5901, 5964, 5965,
       6022, 6023, 6074, 6075, 6120, 6121, 6160, 6161, 6194, 6195, 6222,
       6223, 6244, 6245, 6260, 6261, 6270, 6271], dtype=int32)], FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Brezzi-Douglas-Marini', triangle, 1))
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
  File "/Users/anil/.local/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3378, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-ee23d4688714>", line 1, in <module>
    runfile('/Users/anil/Py/FEM/mixed_poisson.py', wdir='/Users/anil/Py/FEM')
  File "/Applications/PyCharm CE.app/Contents/plugins/python-ce/helpers/pydev/_pydev_bundle/pydev_umd.py", line 198, in runfile
    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script
  File "/Applications/PyCharm CE.app/Contents/plugins/python-ce/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "/Users/anil/Py/FEM/mixed_poisson.py", line 99, in <module>
    bc = fem.dirichletbc(value=u_zero, dofs=boundary_dofs, V=V.sub(0))
  File "/Users/anil/Apps/anaconda3/envs/fenics-env/lib/python3.8/site-packages/dolfinx/fem/bcs.py", line 182, in dirichletbc
    return formcls(value, dofs, V)
  File "/Users/anil/Apps/anaconda3/envs/fenics-env/lib/python3.8/site-packages/dolfinx/fem/bcs.py", line 127, in __init__
    super().__init__(_value, dofs, V._cpp_object)  # type: ignore
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)
Invoked with: array([0., 0.]), [array([   0,    1,   17,   18,   40,   41,   71,   72,  110,  111,  157,
        158,  212,  213,  275,  276,  346,  347,  425,  426,  512,  513,
        607,  608,  710,  711,  821,  822,  940,  941, 1067, 1068, 1202,
       1203, 1345, 1346, 1496, 1497, 1655, 1656, 1822, 1823, 1997, 1998,
       2180, 2181, 2371, 2372, 2570, 2571, 2777, 2778, 2992, 2993, 3215,
       3216, 3446, 3447, 3685, 3686, 3932, 3933, 4187, 4188, 4194, 4195,
       4447, 4448, 4692, 4693, 4929, 4930, 5158, 5159, 5379, 5380, 5592,
       5593, 5797, 5798, 5994, 5995, 6183, 6184, 6364, 6365, 6537, 6538,
       6702, 6703, 6859, 6860, 7008, 7009, 7149, 7150, 7282, 7283, 7407,
       7408, 7524, 7525, 7633, 7634, 7734, 7735, 7827, 7828, 7912, 7913,
       7989, 7990, 8058, 8059, 8119, 8120, 8172, 8173, 8217, 8218, 8254,
       8255, 8283, 8284, 8304, 8305, 8317, 8318], dtype=int32), array([   0,    1,   14,   15,   32,   33,   56,   57,   86,   87,  122,
        123,  164,  165,  212,  213,  266,  267,  326,  327,  392,  393,
        464,  465,  542,  543,  626,  627,  716,  717,  812,  813,  914,
        915, 1022, 1023, 1136, 1137, 1256, 1257, 1382, 1383, 1514, 1515,
       1652, 1653, 1796, 1797, 1946, 1947, 2102, 2103, 2264, 2265, 2432,
       2433, 2606, 2607, 2786, 2787, 2972, 2973, 3164, 3165, 3170, 3171,
       3360, 3361, 3544, 3545, 3722, 3723, 3894, 3895, 4060, 4061, 4220,
       4221, 4374, 4375, 4522, 4523, 4664, 4665, 4800, 4801, 4930, 4931,
       5054, 5055, 5172, 5173, 5284, 5285, 5390, 5391, 5490, 5491, 5584,
       5585, 5672, 5673, 5754, 5755, 5830, 5831, 5900, 5901, 5964, 5965,
       6022, 6023, 6074, 6075, 6120, 6121, 6160, 6161, 6194, 6195, 6222,
       6223, 6244, 6245, 6260, 6261, 6270, 6271], dtype=int32)], <dolfinx.cpp.fem.FunctionSpace object at 0x7ff3462c9870>
4

Hi @dokken,
I could resolve the above mentioned error by following your suggestion like this.
I am trying to reproduce the solution to Poisson equation using the mixed formulation as done using legacy FEniCS.
However, the solution is erroneous as it contains only inf values.
Attached is the new code.
Any help is highly appreciated.

from mpi4py import MPI
from dolfinx import fem, mesh, plot
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem
import pyvista as pv


domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                 nx=32, ny=32,
                                 cell_type=mesh.CellType.triangle)
P1 = ufl.FiniteElement(family='BDM', cell=domain.ufl_cell(), degree=1)
P2 = ufl.FiniteElement(family='DG', cell=domain.ufl_cell(), degree=0)
W = fem.FunctionSpace(mesh=domain, element=ufl.MixedElement([P1, P2]))

# test and trial functions
tau, v = ufl.TestFunctions(function_space=W)
sigma, u = ufl.TrialFunctions(function_space=W)


# function determining if a node is on the tray top
def on_left_and_right(x):
    return np.logical_or(np.isclose(x[0], 0),
                         np.isclose(x[0], 1))


W0, _ = W.sub(0).collapse()
u_zero = fem.Function(V=W0)
u_zero.x.array[:] = 0.0  # ----------------------------------------------------   Changed here
boundary_facets = mesh.locate_entities_boundary(mesh=domain,
                                                dim=domain.topology.dim-1,
                                                marker=on_left_and_right)
boundary_dofs = fem.locate_dofs_topological(V=(W.sub(0), W0),
                                            entity_dim=domain.topology.dim-1,
                                            entities=boundary_facets)

# apply Dirichlet BC
bc = fem.dirichletbc(value=u_zero,
                     dofs=boundary_dofs,
                     V=W.sub(0))  # V.sub(0): Corresponds to BDM element


# apply Neumann BC
x = ufl.SpatialCoordinate(domain=domain)
f = 10.0 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = ufl.sin(5 * x[0])

# Variational form
a = (ufl.dot(sigma, tau) + ufl.div(tau) * u + ufl.div(sigma) * v) * ufl.dx
L = -ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

# Solve
problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
5

I am trying to write the basic code of mixed formulation for poison using dolphinx. Could you please share with known solution and computed rate of convergence?

6

Have you found anything? I am also trying to do write the mixed formulation for poisson.

7

The problem with the code above is the usage of the default petsc “lu” solver.
Changing it to mumps or umfpack works:

from mpi4py import MPI
from dolfinx import fem, mesh
import basix.ufl
import numpy as np
import ufl
from dolfinx.fem.petsc import LinearProblem


domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                 nx=32, ny=32,
                                 cell_type=mesh.CellType.triangle)
el0 = basix.ufl.element("BDM", domain.topology.cell_name(),1)
el1 = basix.ufl.element("DG", domain.topology.cell_name(), 0)
W = fem.functionspace(mesh=domain, element=basix.ufl.mixed_element([el0, el1]))

# test and trial functions
tau, v = ufl.TestFunctions(function_space=W)
sigma, u = ufl.TrialFunctions(function_space=W)


# function determining if a node is on the tray top
def on_left_and_right(x):
    return np.isclose(x[0], 0) | np.isclose(x[0], 1)


W0, _ = W.sub(0).collapse()
u_zero = fem.Function(W0)
u_zero.x.array[:] = 0.0  # ----------------------------------------------------   Changed here
boundary_facets = mesh.locate_entities_boundary(mesh=domain,
                                                dim=domain.topology.dim-1,
                                                marker=on_left_and_right)
boundary_dofs = fem.locate_dofs_topological(V=(W.sub(0), W0),
                                            entity_dim=domain.topology.dim-1,
                                            entities=boundary_facets)

# apply Dirichlet BC
bc = fem.dirichletbc(value=u_zero,
                     dofs=boundary_dofs,
                     V=W.sub(0))


# apply Neumann BC
x = ufl.SpatialCoordinate(domain=domain)
f = 10.0 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = ufl.sin(5 * x[0])

# Variational form
a = (ufl.dot(sigma, tau) + ufl.div(tau) * u + ufl.div(sigma) * v) * ufl.dx
L = -ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

# Solve
problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                    "pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()
print(uh.x.array, problem.solver.getConvergedReason())

This code is written from dolfinx 0.7.0. The only difference is the usage of

el0 = basix.ufl.element("BDM", domain.topology.cell_name(),1)
el1 = basix.ufl.element("DG", domain.topology.cell_name(), 0)
W = fem.functionspace(mesh=domain, element=basix.ufl.mixed_element([el0, el1]))

instead of

I’ll aim to add a singular Poisson demo to the FEniCSx tutorial over the weekend.

8

Hi dokken thank you so much for your quick reply. I could run this code with minor adjustment for DOLFINx 0.6.0. This looks similar to add_demo, dolfinx_070_demo_mixed_poisson but I don’t know what is the best way to get the function \sigma(x,y) and u(x,y). I tried

sigma_, u_ = uh.split()

but it gave me the same function as uh. Suppose the problem is solved and the solutions are stored in uh. Next I want to have sigma and u functions respectively so that I can do further manipulation on them say integrate over a subdomain, and also visualize the functions respectively. I found a very powerful package based on FEniCSx for visualization called ‘Viskex’ but it looks only support DOLFINx 0.7.0… The PyVista method to create_vtk_mesh seems to only support Lagrange element. Is there any efficient way to extract both functions and do something on them?

9

Replace this with sigma_ = uh.sub(0).collapse()
and u_ = uh.sub(1).collapse().

The reason for this is that BDM element has degrees of freedom that are defined through integral moments, rather than point-evaluations, making them quite tricky to visualize in a sensible way with VTK. Thus you should interpolate the BDM function into an appropriate vector DG space (DG-1), which BDM is embedded in. You can use create_vtk_mesh with both continuous and discontinuous Lagrange functions.

10

Thank you for your guidance. While I don’t really know how to do that and didn’t find relevant materials…

I naively tried

W0_DG = fem.FunctionSpace(domain, ("DG", 1))
sigma_DG = fem.Function(W0_DG)
sigma_DG.interpolate(sigma_)

but it retuned me RuntimeError: Interpolation: elements have different value dimensions. Maybe this is because len(sigma_DG.x.array) is different from len(sigma_.x.array).

After interpolation, I should have a DG function sigma_DG the same value as sigma_ in BDM space. But I don’t know how to plot it as well. This page said it was better to export the data points and use external tools say ParaView to visualize. I wonder can I do it through something similar to create_vtk_mesh and just plot in python. So far it seems this kind of things only happens for Lagrange function. Should I further interpolate from DG to Lagrange? This page said I could use plot(Function(DG, indicators)) but Im not sure if it works for current version of FEniCSx .

I also dont understand why you say

vector DG space (DG-1).

I thought in the mixed poisson, the function value is always scalar. Where does this “vector” come from and what is (DG-1).

11

Sigma is a vector. See the original mathematics in:

Thus you use BDM, which I linked to the basis functions of here: DefElement

Has to be changed to

W0_DG = fem.VectorFunctionSpace(domain, ("DG", 1))
sigma_DG = fem.Function(W0_DG)
sigma_DG.interpolate(sigma_)

Look at the basis functions that I referred you to. You will see that if you set two elements next to eachother, they allow for discontinuities in the tangential direction.this is for instance stated in the continuity section of DefElement: Brezzi–Douglas–Marini

A vector DG space is a discontinuous Lagrange space, i.e. each cell has its own degrees of freedom, and the coefficients are not shared across cells, leading to the possibility of representing discontinuous functions

12

Hi Dokken, thank you very much for you help. Now I know how to split the mixed function appropriately . For people who may be interested this is really the mixed poisson demo for DOLFINx 0.7.0.