User defined Boundary Conditions in multidimensional problems

Hi,

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!

The block_size

takes care of the blocked components in a vector or tensor space.

Thanks for the reply! That’s what I thought too but I get the error setting an array element with a sequence

The minimal code is as follows

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)

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.array[local_dof * V.dofmap.bs + b] = [1,2]
     


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()

This should be

values = [1,2]
# add all the other code here
# ….
       for b in range(V.dofmap.bs):
            num_dofs += 1
            u_bc.x.array[local_dof * V.dofmap.bs + b] = values[b]

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)