dolfin team mailing list archive

[Johan Jansson <johanjan@xxxxxxxxxxxxxxxx>]

On Mon, Aug 02, 2004 at 03:28:00PM -0400, Bing Jian wrote:
> Hi,
> 
>    I have a couple of questions on elasticity stationary module.
> 
>    First of all, what's elasticity stationary equation? I googled
> "Elasticity Stationary Equation" but did not find a clear answer.

I have put up a definition (from a book by Peter Hansbo) here (page 12):

http://www.math.chalmers.se/~johanjan/Hansbo.pdf

Here is the ASCII version:

sigma = lambda * div(u) * I + 2 * mu * epsilon(u) (1)

-div(sigma) = f (2)

Where lambda and mu are positive coefficients (material coefficients)
and epsilon is defined as:

epsilon_ij(u) = 0.5 * (du_i / dx_j + du_j / dx_i)

and div(sigma) is defined as:

(sum_j=1^3 dsigma_ij / dx_j)_i=1^3

Note that equation 1 (the definition of sigma) has 9 components (the
terms are matrices) while equation 2 (the actual equation we use in
the end) has 3 components.

We use some intermediate expressions as "learning wheels" in the
DOLFIN module which I think are helpful for clarity:

sigma00 = lambda * (u(0).ddx() + u(1).ddy() + u(2).ddz()) + 2 * mu * epsilon00;
sigma01 = 2 * mu * epsilon01;
sigma02 = 2 * mu * epsilon02;
sigma10 = 2 * mu * epsilon10;
sigma11 = lambda * (u(0).ddx() + u(1).ddy() + u(2).ddz()) + 2 * mu * epsilon11;
sigma12 = 2 * mu * epsilon12;
sigma20 = 2 * mu * epsilon20;
sigma21 = 2 * mu * epsilon21;
sigma22 = lambda * (u(0).ddx() + u(1).ddy() + u(2).ddz()) + 2 * mu * epsilon22;

Those are the components of sigma. Then just apply the divergence as
defined and you have the left hand side of the equation (multiply by
-1 also).

>
>    From jj's reply to my previous questions, it seems to me that
> elasticity-stationary is also a system of PDEs. However, according
> to the code in src/demo/solvers/elasticity-stationary/main.cpp,
> 
> int main(int argc, char **argv)
> {
>   Mesh mesh("tetmesh-4.xml.gz");
>   Problem elasticitystationary("elasticity-stationary", mesh);
> 
>   elasticitystationary.set("source", f);
>   elasticitystationary.set("diffusivity", a);
>   elasticitystationary.set("convection", b);
>   elasticitystationary.set("boundary condition", mybc);
>   elasticitystationary.set("final time", 2.0);
>   elasticitystationary.set("time step", 0.1);
> 
>   elasticitystationary.solve();
> 
>   return 0;
> }
> 
>    It looks like this equation depends on the parameters f(),a(),b(),
> right?  What's the exact formula of it in forms of something like
>         -div(c*grad(u))+a*u=f

See above. a and b are not used, they have been now been removed from
the file.

> 
>    By reading the source files, I assume the statement
>        elasticitystationary.solve();
>    will call the function void ElasticityStationarySolver::solve()
> defined in
> modules/elasticity-stationary/ElasticityStationarySolver.cpp
> 
> void ElasticityStationarySolver::solve()
> {
>   Matrix A;
>   Vector x10, x11, x20, x21, b, xcomp;
> 
>   Function::Vector u0(mesh, x10, 3);
>   Function::Vector u1(mesh, x11, 3);
>   Function::Vector w0(mesh, x20, 3);
>   Function::Vector w1(mesh, x21, 3);
> 
>   Function::Vector f("source", 3);
> 
>   ElasticityStationary   elasticity(f, u0, w0);
>   KrylovSolver solver;
>   File         file("elasticitystationary.m");
> 
>   // Time independent
> 
>   // Assemble matrix and vector
>   FEM::assemble(elasticity, mesh, A, b);
> 
>   // Solve the linear system
>   //solver.setMethod(KrylovSolver::CG);
>   //solver.setPreconditioner(KrylovSolver::ILU0);
>   solver.solve(A, x11, b);
>   //x1.show();
> 
>   //A.show();
>   //b.show();
> 
>   //cout << "dot(b, b): " << (b * b) << endl;
>   //cout << "dot(x, x): " << (x11 * x11) << endl;
>   //printf("dot(b, b): %30.30lf\n", (b * b));
>   //printf("dot(x, x): %30.30lf\n", (x11 * x11));
> 
>   file << u1;
> }
> 
>    But I am really confused by above code. I can not see
> (1) how the parameters f,a,b set in main.cpp are used here,
>    seems the only thing passed in is the mesh.

The key function is this:

FEM::assemble(elasticity, mesh, A, b);

which does this (for the left hand side, corresponding to A):

(for every cell)

// Iterate over test and trial functions
for (FiniteElement::Vector::TestFunctionIterator v(element); !v.end(); ++v)
  for (FiniteElement::Vector::TrialFunctionIterator u(element); !u.end(); ++u)
    A(v.dof(cell), u.dof(cell)) += pde.lhs(*u, *v);

For the stationary elasticity case, pde.lhs() is defined in
ElasticityStationary.h as:

    real lhs(ShapeFunction::Vector& u, ShapeFunction::Vector& v)
    {

      ...

      return 
        (sigma00 * epsilonv00 + sigma01 * epsilonv01 + sigma02 * epsilonv02 +
         sigma10 * epsilonv10 + sigma11 * epsilonv11 + sigma12 * epsilonv12 +
         sigma20 * epsilonv20 + sigma21 * epsilonv21 + sigma22 * epsilonv22)
        * dx;
    }

and is exactly the left hand side of the weak form of the elasticity
equation, which is derived on pages 12-13 in the Hansbo document.

The right hand side is defined as:

    real rhs(ShapeFunction::Vector& v)
    {
      return (f(0) * v(0) + f(1) * v(1) + f(2) * v(2)) * dx;
    }

and there you have the f().

> (2) what are the meanings of those local variables A,u0,u1,w0,w1,
> x10,x11,x20,x21,f,b? seems they are initialized by mesh only
> (3) how the execution of solver.solve(A,x11,b) changes the
> value of u1? looks like the result is stored in the u1.

In the stationary case, only A, b, f, u1 and x11 (we could call them
just u and x instead) are used. u1 is the solution function and x11 is
the vector containing the nodal values. They are coupled through:

Function::Vector u1(mesh, x11, 3);

If you have the mesh and the nodal values x11, then u1 (the solution)
is known, since it is defined as a linear interpolation of the nodal
values on the mesh (in the implementation, u1 becomes a DofFunction IIRC).

> 
>    I would be rather thankful if anybody can clarify my above
> questions.
> 

Hope this helps,
  Johan