I am using a Windows machine and have installed FEniCS 2019.1.0 through docker.
I am solving 1D axially loaded bar with BCs: u(0) = 0 and EAdu/dx(x=1)=F.
I have created a mesh with 2 elements and solve the problem with usual approach, but I am not getting the correct solution.
I think I am not applying the Neuman BC correctly. I wrote g = 5*x[0]…(5=F/EA in my problem).
I think this is a wrong approach because it will also impose constrain on the left boundary, but my left boundary has a DIrichlet BC.
I am not understanding how to give BC in this problem.
Please post your minimal example as actual code and not as a screenshot. Please also encapsulate it with 3x` to make it render properly.
You say that you are solving a 1D problem, however, your plot is 2 dimensional (and your definititon of the mesh is not included). This means that you are applying a Neumann boundary condition \partial u /\partial n = 5 x on the top and bottom boundary as well, (not just at x=1).
See for instance: Poisson equation with multiple subdomains — DOLFIN documentation
on how to define integration measures over a subset of the boundary.
Sir, I have created 2d mesh because I didn’t find any way of creating 1d mesh. If there is any way, please tell me. I am following this approach because, in one tutorial from the FEniCS book, the 1d problem is solved as 2d by considering the zero gradients at the top and bottom boundary (du/dn = 0).
As you told that I am specifying the Neuman BC on top and bottom boundary as well, I agree that I am doing a mistake in the code.
The way I want to solve the problem is by any of the following way.
- Either creating a 1d mesh of only 2 elements, if its possible then my problem got solved.
- Or creating a 2d mesh and specifying u=0 on left boundary, du/dn=5 on right boundary, and du/dn=0 on top and bottom boundary.
Please help in any of the cases.
I am attaching the complete code here
from fenics import *
import matplotlib.pyplot as plt
Problem parameters
L = 1
F = 10
EA = 2
p = F/EA
Defining mesh
mesh = UnitSquareMesh(2, 1)
V = FunctionSpace(mesh, ‘P’, 1)
plot(mesh)
Defining Dirichlet boundary condition
tol = 1E-14 # tolerance for coordinate comparisons
def Dirichlet_boundary0(x, on_boundary):
return on_boundary and abs(x[0]) < tol
bc = DirichletBC(V, Constant(0), Dirichlet_boundary0)
Defining the problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = -inner(grad(u), grad(v))dx
g = Expression('5x[0]’, degree=1)
L = fvdx - gvds
#problem = VariationalProblem(a, L, bc)
#u = problem.solve()
u = Function(V)
solve(a == L, u, bc)
I am also attaching the problem statement. Please find it.
You can create an IntervalMesh as shown in: Bitbucket
In the previous post I’ve linked you to a code that shows how to apply different Neumann conditions on different parts of the boundary.
I created an interval mesh now only and got the correct and.
Thank you very much sir