Subject: [Fenics] Noslip in x-direction of the hyperelastic cbc.twist case

From: Roland Siegbert <[email-masked]>
Date: Thu, 14 Feb 2013 18:23:03 +0100

Dear Fenics-Userlist,

I want to release the twisted cube from cbc.twist such that at
x[0]==1.0 only movement in y,z-direction is possible.

Does anyone have an idea? I don't see, how I should formulate the
boundary condition using Expression() to fix only one value of the
displacement-vector. The code below produces a pretty weird result...

Best,

Roland

from cbc.twist import *

class Release(Hyperelasticity):

    def mesh(self):
        n = 8
        return UnitCube(n, n, n)

    def end_time(self):
        return 2.0

    def time_step(self):
        return 2.e-2

    def is_dynamic(self):
        return True

    def time_stepping(self):
        return "CG1"

    def reference_density(self):
        return 1.0

    def initial_conditions(self):
        """Return initial conditions for displacement field, u0, and
        velocity field, v0"""
        u0 = "twisty.txt"
        v0 = Expression(("0.0", "0.0", "0.0"))
        return u0, v0

    def dirichlet_values(self):
        clamp = Expression(("0.0", "0.0", "0.0"))
        twist = Expression(("0.0", "x[1]", "x[2]"))
        return [clamp, twist]
        #return [clamp]

    def dirichlet_boundaries(self):
        return ["x[0] == 0.0", "x[0] == 1.0"]
        #return ["x[0] == 0.0"]

    def material_model(self):
        mu    = 3.8461
        lmbda = 5.76
        material = StVenantKirchhoff([mu, lmbda])
        return material

    def __str__(self):
        return "A prestrained hyperelastic cube being let go"

# Setup and solve problem
problem = Release()
print problem
problem.parameters["solver_parameters"]["save_solution"] = True
problem.solve()
interactive()

---

From: Anders Logg <[email-masked]>
Date: Thu, 14 Feb 2013 18:27:43 +0100

I don't think this is possible with the existing interface of
CBC.Twist which assumes you want to set a Dirichlet boundary condition
for all components of the displacement. But it should be possible to
add a Dirichlet condition only for the x-displacement by modifying the
solver itself. It's fairly well documented so it should be possible to
figure out where it happens.

Use the V.sub(i) functionality for setting boundary conditions only
for a component of the full space (see demos for examples).

--
Anders

On Thu, Feb 14, 2013 at 06:23:03PM +0100, Roland Siegbert wrote:
> _______________________________________________
> Mailing list: https://launchpad.net/~fenics
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~fenics
> More help   : https://help.launchpad.net/ListHelp

---

From: Roland Siegbert <[email-masked]>
Date: Fri, 15 Feb 2013 11:10:59 +0100

> solver itself. It's fairly well documented so it should be possible to
The solver is created in cbc.twist.problem_definitions.py in the
solve() routine of the problem directly, thus I cannot access the
boundary condition list bcu defined in the CBCSolver derived class
CG1MomentumBalanceSolver because creation and solve() coincide. My
thought was to simply append my condition like
BoundaryDC(vector.sub(0), Constant(0.0), "x[0] == 1.0") to the solvers
boundary condition list. A solution I see is to remove the coupling
between problem and solver, allowing me to change the "solver" before
executing or recreate the Hyperelasticity problem using DOLFIN only,
what is quite a nasty thing to do.

> for a component of the full space (see demos for examples).
Thanks!

Roland

---

From: Anders Logg <[email-masked]>
Date: Fri, 15 Feb 2013 11:52:46 +0100

> what is quite a nasty thing to do.

My suggestion would be to edit the CBC.Twist code to make it do what
you want. If you find a good solution and a nice interface, it might
be merged back in.

--
Anders

---

From: Roland Siegbert <[email-masked]>
Date: Sun, 17 Feb 2013 11:54:12 +0100

> you want. If you find a good solution and a nice interface, it might
OK, for anybody interested in the result: I added a create_solver()
routine and modified solve() of the CBCProblem class derived CBC.Twist
class Hyperelasticity, in a way that doesn't break the other demos and
gave me the possibility to do this using release.py written by Harish
Narayanan:

<code>
problem.create_solver()

bc = DirichletBC(problem.solver.vector.sub(0), Constant(0.0), "x[0] == 1.0")
problem.solver.bcu.append(bc)

</code>

As can be seen, I also added vector to the instance attributes of the
CBCSolver derived class in cbc.Twist (CG1...).

I'll write a patch for review, after I understand how launchpad works
in that way. ;-)

Roland