dirichletBC doesn't seems to work

dolfinx
1

Hi, I’m trying to do a simple experiment but with kind of weird boundary condition. It’s a unit square with x[0]=0 on the left face, x[0]=0.1 on the right face and x[1]=0 on the [0,0.5] node.

The problem is that the Paraview plotted result doesn’t seems to follow this constraints, the x[1] displacement seems to have it minimum at the bottom and not in the middle of the square.
Sorry for my english and thanks in advance for the help!

image
image1193×704 10.6 KB

Code:

from mpi4py import MPI
from dolfinx import mesh, fem, plot
import numpy as np
from petsc4py.PETSc import ScalarType
import ufl
import dolfinx

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8,8, mesh.CellType.triangle)
V = fem.VectorFunctionSpace(domain, ("CG", 1)) 

def left(x):
    return np.isclose(x[0], 0)
def right(x):
    return np.isclose(x[0], 1)

boundary_facets = mesh.locate_entities_boundary(domain, domain.topology.dim-1, left)
boundary_dofs_x = fem.locate_dofs_topological(V.sub(0), domain.topology.dim-1, boundary_facets)
boundary_dofs_x

bcx = fem.dirichletbc(ScalarType(0), boundary_dofs_x, V.sub(0))
boundary_facets2 = mesh.locate_entities_boundary(domain, domain.topology.dim-1, right)
boundary_dofs_x2 = fem.locate_dofs_topological(V.sub(0), domain.topology.dim-1, boundary_facets2)
bcx2 = fem.dirichletbc(ScalarType(0.1), boundary_dofs_x2, V.sub(0))
def ambos(x):
    return np.logical_and(np.isclose(x[0],0),np.isclose(x[1],0.5))
boundary_facets3 = mesh.locate_entities_boundary(domain, domain.topology.dim-2, ambos)
boundary_dofs_x3 = fem.locate_dofs_topological(V.sub(1), domain.topology.dim-2, boundary_facets3)

bcx3 = fem.dirichletbc(ScalarType(0), boundary_dofs_x3,V.sub(1))
boundary_dofs_x3

bcs = [bcx,bcx2,bcx3]

E = 210000 #MPa
lambda_ = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))
T = fem.Constant(domain, ScalarType((0, 0)))
mt = dolfinx.mesh.meshtags(domain, 1, boundary_facets, 1)
ds = ufl.Measure("ds", subdomain_data=mt)

def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, 0)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

from dolfinx import io
with io.VTKFile(domain.comm, "output.pvd", "w") as vtk:
    vtk.write([uh._cpp_object])
with io.XDMFFile(domain.comm, "output.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(uh)
2

Two comments:

  • After adding nu=0.3, your code runs fine on my machine and I can’t reproduce your error (the dirichlet BC is correctly applied to the central node).
  • Applying a dirichlet BC at a single coordinate is not guaranteed to always work, it depends on your mesh having a node/dof at the specified location or not. For more general applications, I recommend instead to perform a bounding box tree search from your specified coordinate, look for the closest node/dof and apply the BC to it.
    In the specific code you are sharing, it should work though as a 8x8 discretization of the unit square happens to have a node at the coordinate (0, 0.5). But if you chose instead, e.g. a 7x7 discretization, then you can no longer expect your BC to be applied.
2 Likes
3

Thanks for your answer Calva.
Yes, I didn’t copy the nu = 0.3, my mistake.
I knew about the mesh sensible bc, I was just testing it like that but thanks for the suggestion, i didn’t know about the bounding box tree search, I’ll implement it in the future for sure.
About the error, that’s so weird, I keep getting the result shown in the uploaded image, any idea why I’m getting it wrong and you right?
Could you upload your output file please? So I can see if I’m interpreting the results correctly
Thank you again!

4

Here is the output file I get from your code. About why you’re get the result wrong, I’m not sure. If it can help, I’m using version 0.4.1 of dolfinx.

1 Like
5

As a side note to @Calva’s great response, you can:

  1. Check the value of the given degree of freedom with:
x = V.tabulate_dof_coordinates()
print(uh.x.array[boundary_dofs_x3], x[boundary_dofs_x3//V.dofmap.bs])
  1. Use glyphs in Paraview to properly visualize the data, using magnitudal plots with a non-zero min/max is probably what confuses you.
1 Like