Hi thank you for your professional reply. While I am doing in layman’s view, I may miss something. Forecast: this is going to be a lengthy reply.
For the mixed Poisson problem, we have two set of boundary condition, which are both Dirichlet, in a sense that the unknown functions, instead of their derivatives, are specified on the boundary. When I was reformulating, I chose different ways to impose the boundary condition:
- Strongly imposing u=u_0 on \Gamma_D, thus v=0 on \Gamma_D, and weakly imposing \vec{\sigma}\cdot\vec{n}=g on \Gamma_N, leading to
\int_{\Omega}(\vec{\sigma}\cdot\vec{\tau}-\vec{\tau}\cdot\vec{\nabla}u-\vec{\sigma}\cdot\vec{\nabla}v)dx=-\int_{\Omega}fvdx-\int_{\Gamma_N}vgds.
For this imposition, I can have the same answer as the original one, as answered by myself.
- Strongly imposing both boundary conditions, where we can even leave the original form without doing integration by part. Or we can integrate by part for both derivatives,
\int_{\Omega}(\vec{\sigma}\cdot\vec{\tau}+u\vec{\nabla}\cdot\vec{\tau}-\vec{\sigma}\cdot\vec{\nabla}v)dx-\int_{\Gamma_D}u\vec{\tau}\cdot\vec{n}ds+\int_{\Gamma_N}v\vec{\sigma}\cdot\vec{n}=-\int_{\Omega}fvdx.
In these two method, I can have the same answer but they are a little different from the original answer.
- Weakly imposing both boundary conditions.
\int_{\Omega}(\vec{\sigma}\cdot\vec{\tau}+u\vec{\nabla}\cdot\vec{\tau}-\vec{\sigma}\cdot\vec{\nabla}v)dx-\int_{\Gamma_N}u\vec{\tau}\cdot\vec{n}ds+\int_{\Gamma_D}v\vec{\sigma}\cdot\vec{n}=-\int_{\Omega}fvdx-\int_{\Gamma_N}gvds.
In this way I have the same answer as the original one.
But you are saying one has to impose the boundary condition as described in the demo. I don’t understand this. I checked Brenner and Scott’s. In the heading of Chapter 12 they said
“one characteristic of mixed methods is that not all choices of finite element spaces will lead to convergent approximations”.
Maybe this is what you were saying or the reason I got a different answer when I impose two bc strongly? BTW their mixed element was for static stokes equation if I recognized that.
I am perfectly fine with the word strongly imposing and weakly imposing. By stating the bc specifically in the input, one determines the bc term directly. While the other is solved through the algraic equations, where numerical error may occur.
I am doing mixed element because I wanted to solve a coupled equation for research purpose. So I am studying this.
I attach the code (dolfinx 0.6.0) here.
For part 2, the bilinear form for integration by part goes with boundary terms. The NO-integration-by-part goes WITHOUT boundary terms.
import numpy as np
import pyvista
from ufl import FiniteElement, MixedElement, VectorElement
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import(Measure, SpatialCoordinate, FacetNormal, TestFunctions, TrialFunctions,
div, grad, exp, inner)
from mpi4py import MPI
from petsc4py import PETSc
domain = mesh.create_unit_square(MPI.COMM_WORLD, 12, 12, mesh.CellType.triangle)
k = 1
Q_el = VectorElement("Lagrange", domain.ufl_cell(), k)
P_el = FiniteElement("Lagrange", domain.ufl_cell(), k)
V_el = MixedElement([Q_el, P_el])
V = fem.FunctionSpace(domain, V_el)
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
dx = Measure("dx", domain)
# bilinear and linear form for integration by part for both
a = (inner(sigma, tau) + u * div(tau) - inner(sigma, grad(v))) * dx
L = -f * v * dx
# bilinear and linear form for: leave the original form, do not perform integration by part
a = (inner(sigma, tau) - inner(grad(u), tau) + div(sigma) * v) * dx
L = - f * v * dx
boundaries = [(1, lambda x: np.isclose(x[1], 0.0)),
(2, lambda x: np.isclose(x[0], 1.0)),
(3, lambda x: np.isclose(x[1], 1.0)),
(4, lambda x: np.isclose(x[0], 0.0))]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with io.XDMFFile(domain.comm, "t1mixed_poisson_alt1/facet_tags.xdmf","w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(facet_tag)
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)
n = FacetNormal(domain)
# boundary terms for integration by part for both
a += + v * inner(sigma, n) * ds(1) + v * inner(sigma, n) * ds(3) - u * inner(tau, n) * ds(2) - u * inner(tau, n) * ds(4)
facet_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facet_top)
def f1(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = np.sin(5 * x[0])
return values
f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0))
facet_bottom = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facet_bottom)
def f2(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = -np.sin(5 * x[0])
return values
f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0))
facet_left_or_right = mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)))
P, _ = V.sub(1).collapse()
dofs_left_or_right = fem.locate_dofs_topological((V.sub(1), P), fdim, facet_left_or_right)
u_left_or_right = fem.Function(P)
u_left_or_right.x.array[:] = 0.0
bc_left_or_right = fem.dirichletbc(u_left_or_right, dofs_left_or_right, V.sub(1))
P, _ = V.sub(1).collapse()
u_left_or_right = fem.Function(P)
u_left_or_right.x.array[:] = 0.0
bc = [bc_top, bc_bottom, bc_left_or_right]
problem = LinearProblem(a, L, bc, petsc_options={"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"})
w_h = problem.solve()
sigma_h, u_h = w_h.sub(0).collapse(), w_h.sub(1).collapse()
grid = pyvista.UnstructuredGrid(*plot.create_vtk_mesh(P))
grid.point_data["u"] = u_h.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.show()
For the weak imposition for both bc (part 3):
import numpy as np
import pyvista
from ufl import FiniteElement, MixedElement, VectorElement
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import(Measure, SpatialCoordinate, FacetNormal, TestFunctions, TrialFunctions,
div, grad, exp, inner, sin)
from mpi4py import MPI
from petsc4py import PETSc
domain = mesh.create_unit_square(MPI.COMM_WORLD, 12, 12, mesh.CellType.triangle)
k = 1
Q_el = VectorElement("Lagrange", domain.ufl_cell(), k)
P_el = FiniteElement("Lagrange", domain.ufl_cell(), k)
V_el = MixedElement([Q_el, P_el])
V = fem.FunctionSpace(domain, V_el)
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = sin(5 * x[0])
dx = Measure("dx", domain)
a = (inner(sigma, tau) + u * div(tau) - inner(sigma, grad(v))) * dx
L = - f * v * dx
boundaries = [(1, lambda x: np.isclose(x[1], 0.0)),
(2, lambda x: np.isclose(x[0], 1.0)),
(3, lambda x: np.isclose(x[1], 1.0)),
(4, lambda x: np.isclose(x[0], 0.0))]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
facet_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facet_top)
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)
n = FacetNormal(domain)
a += - u * inner(tau, n) * ds(1) - u * inner(tau, n) * ds(3)\
+ v * inner(sigma, n) * ds(2) + v * inner(sigma, n) * ds(4)
L += - g * v * ds(1) - g * v * ds(3)
problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"})
w_h = problem.solve()
sigma_h, u_h = w_h.sub(0).collapse(), w_h.sub(1).collapse()
P, _ = V.sub(1).collapse()
grid = pyvista.UnstructuredGrid(*plot.create_vtk_mesh(P))
grid.point_data["u"] = u_h.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.show()