This is a follow up to an older topic. Consider the following 2d Poisson Equation problem -\nabla ^2\mathbf{u}=\left( -1,-2 \right) with BC \mathbf{u}\left( \mathbf{x} \right) =\left( 1,2 \right) ,\mathbf{x}\in \partial \Omega where I simplify the problem somewhat since the BC won’t be 1 or 2 but results from another calculation on the mesh boundary points. The code is taken from the older topic Fenicsx: use an array of data for boundary conditions with some changes
from IPython import embed
import dolfinx
from mpi4py import MPI
import numpy as np
import random
import ufl
from dolfinx import fem
from petsc4py.PETSc import ScalarType
Nx = 5
Ny = 6
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, Nx, Ny)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
u_bc_x = dolfinx.fem.Function(V)
u_bc_y = dolfinx.fem.Function(V)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)
dof_layout = V.dofmap.dof_layout
coords = V.tabulate_dof_coordinates()
num_dofs = 0
bc_dofs = []
for facet in boundary_facets:
# Find local index for each facet
cells = f_to_c.links(facet)
assert len(cells) == 1
facets = c_to_f.links(cells[0])
local_index = np.flatnonzero(facets == facet)
# Find local index for dofs in cell
closure_dofs = dof_layout.entity_closure_dofs(
mesh.topology.dim-1, local_index)
cell_dofs = V.dofmap.cell_dofs(cells[0])
for dof in closure_dofs:
local_dof = cell_dofs[dof]
dof_coordinate = coords[local_dof]
print(local_dof, dof_coordinate)
for b in range(V.dofmap.bs):
num_dofs += 1
u_bc_x.x.array[local_dof * V.dofmap.bs + b] = 1
u_bc_y.x.array[local_dof * V.dofmap.bs + b] = 2
bc = dolfinx.fem.dirichletbc(u_bc_x, boundary_dofs)
u_x = ufl.TrialFunction(V)
v_x = ufl.TestFunction(V)
f_x = fem.Constant(mesh, ScalarType(-1))
a_x = ufl.dot(ufl.grad(u_x), ufl.grad(v_x)) * ufl.dx
L_x = f_x * v_x * ufl.dx
problem = fem.petsc.LinearProblem(a_x, L_x, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh_x = problem.solve()
ux=uh_x.x.array
bc = dolfinx.fem.dirichletbc(u_bc_y, boundary_dofs)
u_y = ufl.TrialFunction(V)
v_y = ufl.TestFunction(V)
f_y = fem.Constant(mesh, ScalarType(-2))
a_y = ufl.dot(ufl.grad(u_x), ufl.grad(v_y)) * ufl.dx
L_y = f_y * v_y * ufl.dx
problem = fem.petsc.LinearProblem(a_y, L_y, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh_y = problem.solve()
uy=uh_y.x.array
What we want to do know if solve the problem as a vector function problem basically to maybe add some constraints later and going in the relevant 3d problem.
What is the easier way to include the boundary conditions in the above sense in a multidimensional problem? Any input is greatly appreciated!
Thanks, everything works, sorry for asking something somewhat trivial. I will post the relevant code here, including calculation of the solution in the mesh points in case somebody finds this discussion useful.
Thanks again!
from IPython import embed
import dolfinx
from mpi4py import MPI
import numpy as np
import random
import ufl
from dolfinx import fem
from petsc4py.PETSc import ScalarType
Nx = 5
Ny = 6
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, Nx, Ny)
V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 1))
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
u_bc = dolfinx.fem.Function(V)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)
values = [1,2]
dof_layout = V.dofmap.dof_layout
coords = V.tabulate_dof_coordinates()
num_dofs = 0
bc_dofs = []
for facet in boundary_facets:
# Find local index for each facet
cells = f_to_c.links(facet)
assert len(cells) == 1
facets = c_to_f.links(cells[0])
local_index = np.flatnonzero(facets == facet)
# Find local index for dofs in cell
closure_dofs = dof_layout.entity_closure_dofs(
mesh.topology.dim-1, local_index)
cell_dofs = V.dofmap.cell_dofs(cells[0])
for dof in closure_dofs:
local_dof = cell_dofs[dof]
dof_coordinate = coords[local_dof]
for b in range(V.dofmap.bs):
num_dofs += 1
print(b,local_dof, dof_coordinate)
u_bc.x.array[local_dof * V.dofmap.bs + b] = values[b]
#print(u_bc.x.array)
bc = dolfinx.fem.dirichletbc(u_bc, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(mesh, ScalarType((-1,-2)))
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh_x = problem.solve()
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
points = mesh.geometry.x
u_values = []
from dolfinx import geometry
bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points)
for i, point in enumerate(points):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
# 3d array x_points has mesh points (including the trivial z axis) and u_values=(u_x,u_y) the relevant 2d solution of the Poisson equation in each point
x_points=points_on_proc
u_values = uh_x.eval(points_on_proc, cells)