dolfin team mailing list archive

[Bing Jian <bjian@xxxxxxxxxxxx>]

Hi Johan,


   Thanks for your clarification.  That pdf file and your explanation
are really helpful. Seems that in order to solve my problems, I need
to write a specific module for that, including a subclass of PDE and
a corresponding solver. Also I need to obtain the lefthandside and
righthandside by multiplying a test function and integration. I
guess one of problems would be that there might have a lot of
coeffecients in my class, each of them being an element function. :-(


-- 
Best wishes,
Bing Jian
bjian@xxxxxxxxxxxx


On Tue, 3 Aug 2004, Johan Jansson wrote:

> 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
>
>