Regarding moving boundary

dolfinx
1

Hello everyone,

I’m having trouble enforcing a moving boundary condition in my simulation. One end of the boundary is supposed to follow the position vector [0,0,sin⁡(πt)] but after running the simulation, I see that the x-coordinate of the boundary is not staying at zero as expected. I can’t figure out why this is happening.

Here’s a snippet of my code:

while t < T_end:
    u_bc_left = fem.Constant(domain, PETSc.ScalarType(np.array([0.0 ,0.0, np.sin(np.pi*t)])))
    left_bc = dirichletbc(u_bc_left, fem.locate_dofs_geometrical(V, left_boundary), V)
    bcs = [left_bc]
    
    # Update r_ddot (acceleration term) based on finite difference
    solver.solve(r)

The total code is is given below if you want to reproduce

from dolfinx import mesh, fem ,default_real_type
from mpi4py import MPI
import ufl
from basix.ufl import element
import numpy as np
from petsc4py import PETSc
import dolfinx
import ufl.coefficient
from dolfinx.fem import dirichletbc
from dolfinx.nls.petsc import NewtonSolver
from dolfinx import default_real_type
from dolfinx.fem.petsc import NonlinearProblem
# Create 1D mesh for s in [0, L]
L = 10.0  # Length of the domain
num_elements = 30  # Number of elements in the 1D mesh
domain = mesh.create_interval(MPI.COMM_WORLD, num_elements, [0, L])

# Define a vector function space in 3D (Lagrange element of degree 1)
#element =ufl.element("Lagrange", domain.ufl_cell(), 1, dim=3)
typeelement = element("Lagrange", domain.basix_cell(),1,shape = (3,) ,dtype=default_real_type)
V = fem.functionspace(domain, typeelement)

# Define trial and test functions
r = fem.Function(V)  # Solution function (trial function)
w = ufl.TestFunction(V)  # Test function
r_prime = ufl.nabla_grad(r) 
norm_r_prime = ufl.sqrt(ufl.inner(r_prime, r_prime))
#print(r.x.array)
# Define constants
from dolfinx.fem import Function, assemble_scalar, form
import ufl
mu = fem.Constant(domain, 0.03)  # Mass density, example value
md = fem.Constant(domain, 10.0)
EA = fem.Constant(domain, 100000.0 ) # Stiffness, example value
ti= fem.Constant(domain, 0.0)
gr= np.transpose([0,0,-10])
g = ufl.as_vector(gr)  # Gravity vector in 3D
Fr = np.transpose([0.1,0,1.0*103])
Force = ufl.as_vector(Fr)
# Fix r at s = 0 to be zero (Dirichlet condition)
def left_boundary(x):
    return np.isclose(x[0], 0)

bcs = []
dt = 0.001  # Time step size
t = 0.0
T_end = 20 # Final time

# Time-dependent acceleration term (assuming some finite difference approximation for r_ddot)
r_n = fem.Function(V)  # Solution at previous time step
r_n1 = fem.Function(V)
r_n.interpolate(lambda x: np.vstack((0*x[0], 0 * x[0], x[0])))
#print(r_n.x.array)
r.x.array [:]= r_n.x.array
r_n1.x.array[:] = r_n.x.array  # Zero initial velocity
#r_prime = ufl.transpose(r_prime)
#print(r_prime.ufl_shape)
# Weak form definition
F = (
      ufl.inner(mu * ((r - 2 * r_n + r_n1) / dt**2),w)*ufl.dx - ufl.inner(mu*g,w)* ufl.dx
    +EA * ufl.inner(((norm_r_prime - 1) / norm_r_prime) * r_prime, ufl.nabla_grad(w)) * ufl.dx
    +ufl.inner(md* ((r - 2 * r_n + r_n1) / dt**2 ) - md*g - Force,w)*ufl.ds
)   
# Create nonlinear problem
problem = NonlinearProblem(F, r, bcs)

# Create solver
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-6
solver.rtol = 1e-6
x = r.x.array[0::3]  # Start at index 0, take every 3rd element
y = r.x.array[1::3]  # Start at index 1, take every 3rd element
z = r.x.array[2::3]  # Start at index 2, take every 3rd element
x_coord_matrix = x
y_coord_matrix = y
z_coord_matrix = z
#dolfinx.log.set_log_level(dolfinx.log.LogLevel.DEBUG)
# Time-stepping loop
while t < T_end:
    u_bc_left = fem.Constant(domain, PETSc.ScalarType(np.array([0.0 ,0.0, np.sin(np.pi*t)])))
    left_bc = dirichletbc(u_bc_left, fem.locate_dofs_geometrical(V, left_boundary), V)
    bcs = [left_bc]
    
    # Update r_ddot (acceleration term) based on finite difference
    solver.solve(r)

    # Update previous solutions
    r_n1.x.array[:], r_n.x.array[:] =r_n.x.array, r.x.array
    ti.value = t
    #Now x, y, z are arrays of length 101 containing the coordinates

    x = r.x.array[0::3]  # Start at index 0, take every 3rd element
    y = r.x.array[1::3]  # Start at index 1, take every 3rd element
    z = r.x.array[2::3]  # Start at index 2, take every 3rd element
    x_coord_matrix = np.vstack((x_coord_matrix, x))
    y_coord_matrix = np.vstack((y_coord_matrix, y))
    z_coord_matrix = np.vstack((z_coord_matrix, z))
    # Update time
    t += dt
2

You shouldn’t rely on this. There is no guarantee that DOFs are oredered in that way (they may be, or not, I can’t honestly remember and right now I don’t have time to check). Use instead r.sub(0).collapse(), and then get the array out of the collapsed function.

3

Various maps between mesh nodes, degrees of freedom and vertices have been explained to some extent in: How to apply a boundary condition to a range of DOFs? - #7 by dokken

4

Thank you for your insight.I have solved my problem by locating Dof of the boundary usinglocate_dofs_topological and assigning values directly .Here’s a code snippet of boundary condition that might help others facing a similar problem.

s0_location = mesh.locate_entities_boundary(domain,0,left_boundary)
x0 =fem.locate_dofs_topological(V.sub(0), 0,s0_location)
y0 =fem.locate_dofs_topological(V.sub(1), 0,s0_location)
z0 =fem.locate_dofs_topological(V.sub(2), 0,s0_location)
    
ti= fem.Constant(domain, 0.0)
left_bcx = dirichletbc(0.0, x0, V.sub(0))
left_bcy = dirichletbc(0.0, y0, V.sub(1))
left_bcz = dirichletbc( ti, z0, V.sub(2))
bcs = [left_bcx,left_bcy,left_bcz]

By changing the value of ti inside the loop as i wanted i have solved my problem.