Equation on the boundary

On 06/20/2012 05:41 PM, ktjunior wrote:
> Question #200991 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200991
>
> Description changed to:
> Hello,
>
> I do not know how to solve the following problem:
>
> My unknowns are: u, v
>
> For u I have a laplace equation with the Dirichlet boundary condition
> du/dt = v
> and v satisfies some nonlinear PDE. I tried this
>
> bc1 = DirichletBC(W.sub(0), u0 + dt*v, boundary_parts, 0)
>
> (u0 is the result from the previous timestep) but it does not work. Can
> you help me to write such Dirichlet condition? Thanks, ktj
>

What is v? If that is a function have to add the underlying degrees of
freedoms.

u1 = Function(V)
u1.vector()[:] = u0.vector()
u1.vector().axpy(dt, v.vector())
bc1 = DirichletBC(W.sub(0), u1, boundary_parts, 0)

Johan

It does not work to me (I wrote u1=Function(U)), it says
u1.vector().axpy(dt, v.vector())
AttributeError: 'ListTensor' object has no attribute 'vector'

Thank you

Just to be clear:

U = VectorFunctionSpace(mesh, "CG", 2)
V = VectorFunctionSpace(mesh, "CG", 2)
W = MixedFunctionSpace([U, V])

w = Function(W)
(u,v) = (as_vector((w[0], w[1])),as_vector((w[2], w[3])))

u0 = Function(U)
v0 = Function(V)

as a solver I am using NonlinearVariationalSolver.

On 06/20/2012 11:06 PM, ktjunior wrote:
> Question #200991 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200991
>
> ktjunior posted a new comment:
> It does not work to me (I wrote u1=Function(U)), it says
> u1.vector().axpy(dt, v.vector())
> AttributeError: 'ListTensor' object has no attribute 'vector'
>
> Thank you
>
> Just to be clear:
>
> U = VectorFunctionSpace(mesh, "CG", 2)
> V = VectorFunctionSpace(mesh, "CG", 2)
> W = MixedFunctionSpace([U, V])
>
> w = Function(W)
> (u,v) = (as_vector((w[0], w[1])),as_vector((w[2], w[3])))
>
> u0 = Function(U)
> v0 = Function(V)
>
> as a solver I am using NonlinearVariationalSolver.

So v is not a Function. v is a UFL vector expression of dimension 2, and
u0 is a vector Function of dimension 3. I am not sure what your are
trying to do.

Johan

u0 is vector of dimension 2 (my domain is 2D), u and v are also vector Function of dimension 2, both unknowns. I just do not know how to add Dirichlet condition for the unknown u computed from the unknown v. Am I clear? Thank you.

Physical problem:

I have also pressure p as unknown, but this was not important in my problem, v is velocity and u is deformation of the domain. I am solving the problem of fluid which is being deformed and doing it with the help of ALE method (not using any Fenics ALE tool, just have my own ALE and want to implement it). For this I need the Dirichled boundary condition du/dt = v wchich means that the points on the boundary are the material points.

You then need to provide a minimal script reproducing the error. This is
always the best way of getting help.

Johan

On 06/20/2012 11:31 PM, ktjunior wrote:
> Question #200991 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200991
>
> ktjunior posted a new comment:
> u0 is vector of dimension 2 (my domain is 2D), u and v are also vector
> Function of dimension 2, both unknowns. I just do not know how to add
> Dirichlet condition for the unknown u computed from the unknown v. Am I
> clear? Thank you.
>

If you think I should prepare simplier problem, I will try. The problem is in the line
bc1 = DirichletBC(W.sub(0), u0 + dt*v, boundary_parts, 0)

If it works, the residuum should be zero (there is no force...)

Thank you.

from dolfin import *

# Create mesh and define function space
mesh = Rectangle(0, 0, 2, 1, 20, 10, "right")

# Define function spaces
U = VectorFunctionSpace(mesh, "CG", 2) #deformation
V = VectorFunctionSpace(mesh, "CG", 2) #velocity
Q = FunctionSpace(mesh, "CG", 1) #pressure
W = MixedFunctionSpace([U, V, Q])

#Define boundaries
boundary_parts = FacetFunction("uint", mesh) #MeshFunction("uint", mesh, mesh.topology().dim()-1)
boundary_parts.set_all(0)

ds = Measure("ds")[boundary_parts]

# Define variational problem
w = Function(W)

(u,v,p) = (as_vector((w[0], w[1])),as_vector((w[2], w[3])),w[4])
(u_tst, v_tst, p_tst) = TestFunctions(W)

dt = 0.01
t_end = 2.0

u0 = Function(U)
v0 = Function(V)
p0 = Function(Q)

rho = 1000
eps1 = 1.0

derivative_u = (u - u0)/dt

# Dirichlet boundary conditions
bc1 = DirichletBC(W.sub(0), u0 + dt*v, boundary_parts, 0)
bcs = [bc1]

I = Identity(U.cell().d) # Identity tensor

F = I + grad(u)
detF = det(F)
F_inv = inv(F)
gradv = grad(v)*F_inv

T = -p*I + eps1*(gradv + gradv.T)

rce_0 = (tr(gradv)*p_tst)*dx
rce_1 = (inner((rho*detF*(v-v0)/dt + rho*detF*gradv*(v-derivative_u)),v_tst))*dx + (inner(detF*T*F_inv.T,grad(v_tst)))*dx
rce_2 = (inner(grad(u),grad(u_tst)))*dx

rce = rce_0 + rce_1 + rce_2

J = derivative(rce, w)

ffc_options = {"optimize": True, "eliminate_zeros": True, "quadrature_degree": 8 }
parameters["form_compiler"]["cpp_optimize"] = True

problem=NonlinearVariationalProblem(rce,w,bcs,J,form_compiler_parameters=ffc_options)
solver=NonlinearVariationalSolver(problem)

# Create files for storing solution
ufile = File("results/displacement.pvd")
vfile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

# Time-stepping
t = dt

# Initial conditions
u0.assign(Expression(("0","0")))
v0.assign(Expression(("0","0")))
p0.assign(Expression("0"))

ufile << u0
vfile << v0
pfile << p0

while t < t_end:

    print "t =", t

    # Compute
    begin("Solving ....")
    solver.solve()
    end()

    # Extract solutions:
    (u,v,p) = w.split()

    # Save to file
    ufile << u
    vfile << v
    pfile << p

    # Move to next time step
    u0.assign(u)
    v0.assign(v)
    t += dt