Fenicsx: use an array of data for boundary conditions

belmont · September 2, 2022, 1:39pm

Hi.

I have the following problem: Consider the following matrix in a 3d Poisson equation problem: (x_i, y_i, z_i, bc_i). The (x_i, y_i, z_i) are the bounds of a 3d mesh and bc_i is the relevant boundary conditions. So, consider now making the 3d mesh and perform a numerical calculation to find the bc_i values of each boundary point (x_i, y_i, z_i).

Using those boundary conditions, I would like to use fenicsx to solve the Poisson equation in the whole mesh.

I have solved this problem some years ago in the legacy fenics as follows, with the help of a dictionary variable. Please see the example code. First, I find the boundaries of my geometry

# Calculate the boundary data
# Define the mesh of the boundary    
boundary_mesh = BoundaryMesh(mesh, 'exterior')
# Create an array with the node points
boundary_nodes = boundary_mesh.coordinates()  
# Dimension of the array
length = np.shape(boundary_nodes)[0] 
 
# Create the points data from the mesh
x_boundary, y_boundary, z_boundary = np.zeros(length), np.zeros(length), np.zeros(length)
x_boundary, y_boundary, z_boundary = [boundary_nodes[i][0] for i in range(length)], [boundary_nodes[i][1] for i in range(length)], [boundary_nodes[i][2] for i in range(length)]

Then I ran the numerical calculation to find the values bc_i. Then I use the following functions and dictionary to “feed” the boundary conditions in legacy fenics:

def bound_dictionary(x,y,z,dictionary): 
 try:
    dictionary[(x,y,z)]
 except KeyError:
    b=0
 else:
    a=dictionary[(x,y,z)]
    return a   


class BCfunction_dictionary(UserExpression):
    def eval(self, value, x):
        value[0] = bound_dictionary(x[0],x[1],x[2],dictionary)
    def value_shape(self):
        return ()

 dictionary = {(x_boundary[i], y_boundary[i], z_boundary[i]): bc[i] for i in range(length)}
 
 u_D = BCfunction_dictionary()

 def boundary(x, on_boundary):
  return on_boundary

 bc = DirichletBC(V, u_D, boundary)

# Define variational problem
 u_x = TrialFunction(V)
 v = TestFunction(V)
 f = Constant(0)
 a = dot(grad(u_x), grad(v))*dx
 L = f*v*dx

# Compute solution
 u_x = Function(V)
 solve(a == L, u_x, bc)

 u=np.array(u_x.compute_vertex_values(mesh))

Obviously, it is not elegant, but it works.

Now to my question: what your suggestions for implementing this type of calculation in fenicsx are?

So far, I have found how to find in fenicsx the boundary points, but I am not sure of how to proceed from here. I am not sure how to write the appropriate function.

dokken · September 5, 2022, 7:39am

In DOLFINx, you deal with boundary conditions in multiple stages:

Locate all degrees of freedom associated with the topological entities of interest (here your boundary facets).
Interpolate function into appropriate function space
Use this function in dirichlet bc.

To highlight this, consider the following minimal example

import dolfinx
from mpi4py import MPI
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)


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)

dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
print(dof_coordinates)
u_bc = dolfinx.fem.Function(V)


def f(x):
    print(x)
    # Look-up x,y,z coordinates, preferrably in a vectorized fashion
    return x[0]
u_bc.interpolate(f)
bc = dolfinx.fem.dirichletbc(u_bc, boundary_dofs)

u_boundary = dolfinx.fem.Function(V)
dolfinx.fem.petsc.set_bc(u_boundary.vector, [bc])
u_boundary.x.scatter_forward()
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", u_boundary) as vtx:
    vtx.write(0.0)

which yields the following bc:

belmont · September 5, 2022, 3:46pm

Thanks for the example dokken. I have a better understanding of how boundary conditions work. In your example len(dof_coordinates) is 40 in the main code and print(len(x.T)) gives 600 i.e. it has more “points” (for lack of a better term). Shouldn’t the points be the same? If I consider my 3d problem this is the problem I am facing. I have the array dof_coordinates and the boundary values as an array in those values but in the actual user supplied function, i.e. def f(x) in your example, the vectorized argument x has points not present in the dof_coordinates matrix (in this example it has all that mesh points) and I am not sure where they are from. In my case they don’t even seem to be points from the actual mesh. The boundary function f(x) is evaluated (or interpolated) on the whole mesh? This is a trivial modification of your code that shows plots of what I try to explain

import dolfinx
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)


V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = np.flatnonzero(dolfinx.mesh.compute_boundary_facets(mesh.topology))
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)

dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
n=len(dof_coordinates)
xx,yy=np.zeros(n),np.zeros(n)
for i in range(n):
 xx[i],yy[i]=dof_coordinates[i][0],dof_coordinates[i][1]

# Plot the result
fig = plt.figure()
ax = fig.add_subplot(111)
xs,ys = xx,yy
p = ax.scatter(xs, ys)#, cmap='hot')
ax.set_title('Boundary Points 1', fontsize=15)
ax.set_xlabel('$x$', fontsize=15)
ax.set_ylabel('$y$', fontsize=15)
fig.savefig('bp1.png')

u_bc = dolfinx.fem.Function(V)


def f(x):
   n=len(x.T)
   xx,yy=np.zeros(n),np.zeros(n)
   for i in range(n):
    xx[i],yy[i]=x.T[i][0],x.T[i][1]
   fig = plt.figure()
   ax = fig.add_subplot(111)
   xs,ys = xx,yy
   p = ax.scatter(xs, ys)#, cmap='hot')
   ax.set_title('Boundary Points 2', fontsize=15)
   ax.set_xlabel('$x$', fontsize=15)
   ax.set_ylabel('$y$', fontsize=15)
   fig.savefig('bp2.png') 
   return x[0]


u_bc.interpolate(f)
bc = dolfinx.fem.dirichletbc(u_bc, boundary_dofs)

u_boundary = dolfinx.fem.Function(V)
dolfinx.fem.petsc.set_bc(u_boundary.vector, [bc])
u_boundary.x.scatter_forward()

I really apologise if I am talking about something very naive.

dokken · September 5, 2022, 5:34pm

This depends on how you do the interpolation into the appropriate space.
In the code above, i did not specify which cells to use interpolation for, i.e.

as opposed to

u_bc.interpolate(f, np.array([0,1,2], dtype=np.int32))

which would only interpolate over the cells with local index 0, 1 and 2.
In general, for the problem you are considering, where you have only data at a sub-set of the dofs in a cell, you should do something like the following:

from IPython import embed
import dolfinx
from mpi4py import MPI
import numpy as np
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)

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
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(dof, dof_coordinate)
        for b in range(V.dofmap.bs):
            num_dofs += 1
            u_bc.x.array[local_dof * V.dofmap.bs + b] = dof_coordinate[b]
print(num_dofs, 2 * 2 * Nx + 2 * 2 * Ny)
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", u_bc) as vtx:
    vtx.write(0.0)

which you can see only evaluates 44 times in the example I’ve presented, which is equal to going through each cell with a boundary facet once, and only accessing the dofs located on the boundary.

belmont · September 7, 2022, 9:04am

Thank you very much for your help till now. I can see that the actual boundaries are returnred. I tried to solve the minimum example from the first tutorial using the dof from your code. I considered a valuable close to 1 for the boundary. However, when I print the actual points that are interpolated they seem to be … different from the ones from the above example? I must be doing something wrong. Again I apologise for being somewhat naive

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 = 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] = dof_coordinate[b]
        bc_dofs.append(local_dof)    
print(num_dofs, 2 * 2 * Nx + 2 * 2 * Ny)     

def f(x):
   xT = x.T
   n=x.shape[1]
   values = np.zeros(n)
   for i in range(n):
    values[i]=random.gauss(1,0.01)
    print(xT[i][0],xT[i][1])
   return values

u_bc.interpolate(f, bc_dofs)
bc = dolfinx.fem.dirichletbc(u_bc, boundary_dofs)
u_x = ufl.TrialFunction(V)
v_x = ufl.TestFunction(V)
f_x = fem.Constant(mesh, ScalarType(-6))
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

print(ux)

dokken · September 7, 2022, 12:04pm

belmont:

def f(x):
   xT = x.T
   n=x.shape[1]
   values = np.zeros(n)
   for i in range(n):
    values[i]=random.gauss(1,0.01)
    print(xT[i][0],xT[i][1])
   return values

u_bc.interpolate(f, bc_dofs)

You should remove this code, as the whole point of

belmont:

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] = dof_coordinate[b]
        bc_dofs.append(local_dof)

was to replace the interpolate command

belmont · September 7, 2022, 1:02pm

Thanks for everything. A hopefully final question, while I hope I am not saying something obvious.

To pass the actual value I want in the boundary point, if I understood everything correctly, is to simply input something like

u_bc.x.array[local_dof * V.dofmap.bs + b] = 1

instead of

u_bc.x.array[local_dof * V.dofmap.bs + b] = dof_coordinate[b]

in a simple example when all boundary conditions are simply equal to 1. One can of course use the code above as a basis to give a different value in each point (like I am currently trying to do) or function value etc.

dokken · September 7, 2022, 1:32pm

Yes, that is correct.

belmont · September 8, 2022, 9:44pm

Thanks again for the help. Everything seems to work as indented in my code in both 2d and 3d cases. I will mark the topic as solved. Thanks again, hopefully this info helps more people.

Meena_Pargaei · May 22, 2023, 9:00am

Thank you it is really nice discussion.
I am also doing the same thing but in my case:
(1) A cuboid with a Dirichlet boundary on the bottom surface i.e. where z=0.
(2) I would like to give a .txt data of size (size(x(i)), size(y(j))) into the Dirichlet boundary.

May you please help me with the 2nd step?

Topic		Replies	Views
How to solve an equation with an array data General	6	263	July 16, 2023
Mesh in FEniCS vs mesh in gmsh	2	527	March 23, 2021
2D Poisson equation with periodic boundaries - hangs forever	1	652	March 10, 2020
What changed between using the linear solver and the nonlinear solver? General	2	344	November 7, 2022
Another question about FenicsX boundary conditions General	4	572	January 12, 2023

Fenicsx: use an array of data for boundary conditions

Related Topics