Hi,
I don’t understand what is the meaning of “first 0, then 1”.
I know the first part of the return is to return if the cell is on the boundary, but what is the latter part, what are the 0 and 1 representing?
As clearly stated in the solution, x[0] and x[1] represents the x and y coordinates.
We are then checking first if we are on the boundary, and then if the x coordinate of your subdomain is close to 0 (in the case of BoundaryLeft) or close to 1 (in case of BoundaryRight).
Imagine that you have a unit cube (created with tetrahedral elements) spanning [0,0,0]\times[1,1,1], and you would like to mark all the facets (triangles) on the side of the cube where x is 0 (using BoundaryLeft) and similarly for the side of the cube where x is 1.
The first 0, then 1, is simply referring to the fact that you have to separate functions BoundaryLeft (x close to 0) and BoundaryRight (x close to 1).
Thank you, for a 1d problem (x in [0, L]), if I want to make the right side of the domain remain zero forever, namely
BC u(L,t)=0
From the results, this code seems like makes left equals to zero, can you help me found which part is wrong?
https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html
class Right(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],1)
boundaries = MeshFunction("size_t", mesh,0)
boundaries.set_all(0)
# Now mark your Neumann boundary
rightbound = Right()
rightbound.mark(boundaries, 1)
ds=Measure('ds',subdomain_data=boundaries)
# Time variables
dt = 0.001; t = 0; T = 1
# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)
# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)
left = CompiledSubDomain("(std::abs(x[0] - 0.0) < DOLFIN_EPS) && on_boundary")
a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx-u0*v*dx+ exp(t)*v*ds
bc = DirichletBC(V, 0, left)
A, b = assemble_system(a, L, bc)
First of all, you have not defined a complete code, as your mesh definition is missing. If L!=1, the following code
is not correct, as x[0] has to be replaced with L.
Secondly, you state that you would like to set the right side to zero, but you are defining a Zero condition on the left side:
It is also not clear to me what you are referring to in
I got this from my previous post, 1 d wave equations (x \in [0, L])
two ICs:
u(x,0)=0,
u_t(x,0)=0
two BCs are
u(L,t)=0
u_x(0,t)=exp(t)
from dolfin import *
import numpy as np
c =1000
mesh = IntervalMesh(100,0,1.0)
V=FunctionSpace(mesh, "Lagrange", 1)
class Right(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],1)
boundaries = MeshFunction("size_t", mesh,0)
boundaries.set_all(0)
# Now mark your Neumann boundary
rightbound = Right()
rightbound.mark(boundaries, 1)
ds=Measure('ds',subdomain_data=boundaries)
# Time variables
dt = 0.001; t = 0; T = 1
# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)
# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)
left = CompiledSubDomain("(std::abs(x[0] - 0.0) < DOLFIN_EPS) && on_boundary")
a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx-u0*v*dx+ exp(t)*v*ds
bc = DirichletBC(V, 0, left)
A, b = assemble_system(a, L, bc)
u=Function(V)
vtkfile=File('wave/u.pvd')
while t <= T:
A, b = assemble_system(a, L, bc)
solve(A, u.vector(), b)
u0.assign(u1)
u1.assign(u)
t += dt
vtkfile<<(u,t)
How to modify the left side to make zero, can I just delete
left = CompiledSubDomain("(std::abs(x[0] - 0.0) < DOLFIN_EPS) && on_boundary")
and
bc = DirichletBC(V, 0, right)?
Yes, you can. For future posts, please test such suggesting before asking the question, as you could easily have answered this my testing it and looking at the results
bc = DirichletBC(V, 0, Right)
A, b = assemble_system(a, L, bc)
But this is not working, I don’t know if this is due to the another BC, DirichletBC cannot take the class as argument.
You need to use an initialized version of the class. Ie.
rightbound = Right()
bc= DrichletBC(V,0, rightbound)
Thank you dokken, back to the mesh problem, I changed x[0] to L, but I still don’t understand why x[0] must be L.
Length=10.0
mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)
class Right(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(Length,1)
And a warning
Should be
class Right(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],Length)