I am trying to understand how FEniCSx processes the discontinuous functions used in defining a region dependent conductivity in
https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html
So in the following example, I define a region dependent variable H, which is 1 in region 0 and 0 in region 1. Then I solve the variational formulation of the equation u=H for u, and I print out the results and make a plot. Here is the code
import numpy as np
import dolfinx
from mpi4py import MPI
from dolfinx import fem, mesh, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import ufl
from ufl import dx
from dolfinx import default_scalar_type
import matplotlib.pyplot as plt
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.mesh import locate_entities
nx = 50
domain = mesh.create_interval(MPI.COMM_WORLD, nx,[0., 1.])
Q = functionspace(domain, ("DG", 0))
def Omega_0(x):
return x[0] <= 0.6
def Omega_1(x):
return x[0] >= 0.6
H = Function(Q) # a function that is 0 for y>=0.6 and 1 for y<=0.6
cells_0 = locate_entities(domain, domain.topology.dim, Omega_0)
cells_1 = locate_entities(domain, domain.topology.dim, Omega_1)
H.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
H.x.array[cells_1] = np.full_like(cells_1, 0, dtype=default_scalar_type)
V=fem.functionspace(domain, ("Lagrange", 1))
u=fem.Function(V)
v=ufl.TestFunction(V)
F=((u-H)*v)*dx
problem = NonlinearProblem(F, u, bcs=[])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True
log.set_log_level(dolfinx.log.LogLevel.OFF)
n, converged = solver.solve(u)
assert (converged)
print(f"Number of interations: {n:d}")
print(u.x.array)
x=u.function_space.tabulate_dof_coordinates()[:,0]
plt.figure(1)
plt.plot(x,u.x.array)
Here is the plot
I suspect that the oscillations are coming from the fact that, in my variational formulation, I have equated H (of type (“DG”,0) ) to u (of type (“Lagrange”,1) ). Can someone explain this a little better to me?
Returning to the original tutorial, why don’t I see these oscillations near the boundary of the two regions when I plot uh in the tutorial example? Is it because somehow grad(uh) has the same function type as kappa? How can I print out and make a plot of grad(uh) using matplotlib? I know how to extract values of uh using uh.x.array, but how can I do this with ufl.grad(uh) or its dot product with some other vector to make it a scalar? Below is a slightly modified version of the tuturial, in which I try to make these plots in matplotlib and fail. Can someone tell me how to correct this? My failed attempt appears in the second figure. The first figure shows uh, when I have made the conductivities in the two regions very different. Note that there are no oscillations appearing in the first plot.
# %%
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import create_unit_square, locate_entities
from dolfinx.plot import vtk_mesh
import dolfinx.fem as df
from ufl import dot, dx, grad
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
dx, grad, inner)
from mpi4py import MPI
import meshio
import gmsh
import numpy as np
import pyvista
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
Q = functionspace(mesh, ("DG", 0))
# %%
def Omega_0(x):
return x[1] <= 0.5
def Omega_1(x):
return x[1] >= 0.5
# %%
kappa = Function(Q)
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)
# %%
kappa.x.array[cells_0] = np.full_like(cells_0, 1000, dtype=default_scalar_type)
kappa.x.array[cells_1] = np.full_like(cells_1, 0.01, dtype=default_scalar_type)
# %%
V = functionspace(mesh, ("Lagrange", 1))
u, v = TrialFunction(V), TestFunction(V)
a = inner(kappa * grad(u), grad(v)) * dx
x = SpatialCoordinate(mesh)
L = Constant(mesh, default_scalar_type(1)) * v * dx
dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
bcs = [dirichletbc(default_scalar_type(1), dofs, V)]
# %%
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# %%
import matplotlib.pyplot as plt
upnts =uh.function_space.tabulate_dof_coordinates()
xx, yy = upnts[:,0], upnts[:,1]
x_line = xx[np.isclose(yy,0.4)] # if this does not work then use x_line = x[np.isclose(y, 0)]
y_line = yy[np.isclose(xx,0.4)]
u_line = uh.x.array[np.isclose(xx,0.4)]
plt.ion()
plt.figure(1)
plt.plot(y_line,u_line)
WW = functionspace(mesh, ("Lagrange", 1))
graduh = Function(WW)
graduh_expr = df.Expression(uh.dx(0), WW.element.interpolation_points())
graduh.interpolate(graduh_expr)
gradd=graduh.x.array[np.isclose(xx,0.4)]
plt.figure(2)
plt.plot(y_line,gradd)