Implémentation

We start by setting the Feel++ environment

    using namespace Feel;
    po::options_description laplacianoptions( "Elasticity options" );
    laplacianoptions.add_options()
        ( "E", po::value<double>()->default_value( 1.0e6 ), "Young modulus" )
        ( "nu", po::value<double>()->default_value( 0.3 ), "Poisson ratio" )
        ( "no-solve", po::value<bool>()->default_value( false ), "No solve" )
        ( "weakdir", po::value<bool>()->default_value( false ), "use weak dirichlet" )
        ( "gamma", po::value<double>()->default_value( 100 ), "penalisation term" )
        ( "nullspace", po::value<bool>()->default_value( false ), "add null space" )
        ;

    Environment env( _argc=argc, _argv=argv,
                   _desc=laplacianoptions,
                   _about=about(_name="qs_elasticity",
                                _author="Feel++ Consortium",
                                _email="feelpp-devel@feelpp.org"));

Next, we define the Mesh \(\cal{T}_h\)

    tic();
    auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
    toc("loadMesh");

Next, we define the function space \(X^1_h\)

    tic();
    auto Vh = Pchv<1>( mesh );
    toc("Vh");

    auto u = Vh->element("u");
    auto v = Vh->element("v");
    auto nu = doption(_name="nu");
    auto E = doption(_name="E");
    auto lambda = E*nu/( (1+nu)*(1-2*nu) );
    auto mu = E/(2*(1+nu));
    auto deft = sym(gradt(u));
    auto def = sym(grad(u));
    auto Id = eye<FEELPP_DIM,FEELPP_DIM>();
    auto sigmat = lambda*trace(deft)*Id + 2*mu*deft;
    auto sigma = lambda*trace(def)*Id + 2*mu*def;
    auto f = expr<FEELPP_DIM,1>( soption(_name="functions.f"), "f" );
    auto g = expr<FEELPP_DIM,1>( soption(_name="functions.g"), "g" );

and the linear \(\ell\) and bilinear form \(a\)

    tic();
    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=inner(f,id(v)));
    toc("l");

    tic();
    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=inner( sigmat, grad(v) ) );

    if ( boption(_name="weakdir") )
    {
        double penaldir = doption(_name="gamma");
        a += integrate(_range=markedfaces(mesh,"Dirichlet"),
                       _expr=-inner(sigmat*N(),id(u)) + inner(-sigma*N()+std::max(2*mu,lambda)*penaldir*id(u)/hFace(),idt(u)) );

    }
    else
    {
        a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
    }
    toc("a");

We solve the algebraic system \(A \disp{d} = b\)

    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
    {
        tic();
        boost::shared_ptr<NullSpace<double> > myNullSpace( new NullSpace<double>(backend(),qsNullSpace(Vh,mpl::int_<FEELPP_DIM>())) );
        backend()->attachNearNullSpace( myNullSpace );
        if ( boption(_name="nullspace") )
            backend()->attachNearNullSpace( myNullSpace );

        a.solve(_rhs=l,_solution=u);
        toc("a.solve");
    }

We attach to the algebraic backend the near null space that may be useful for Multigrid methods and optionally the null space. To do that we build a basis spanning the (near) null space in 2D and 3D.

template <typename SpaceType>
NullSpace<double> qsNullSpace( SpaceType const& space, mpl::int_<2> /**/ ) (1)
{
    auto mode1 = space->element( oneX() ); (2)
    auto mode2 = space->element( oneY() ); (2)
    auto mode3 = space->element( vec(Py(),-Px()) ); (3)
    NullSpace<double> userNullSpace( { mode1,mode2,mode3 } ); (4)
    return userNullSpace;
}
template <typename SpaceType>
NullSpace<double> qsNullSpace( SpaceType const& space, mpl::int_<3> /**/ ) (5)
{
    auto mode1 = space->element( oneX() ); (6)
    auto mode2 = space->element( oneY() ); (6)
    auto mode3 = space->element( oneZ() ); (6)
    auto mode4 = space->element( vec(Py(),-Px(),cst(0.)) ); (7)
    auto mode5 = space->element( vec(-Pz(),cst(0.),Px()) ); (7)
    auto mode6 = space->element( vec(cst(0.),Pz(),-Py()) ); (7)
    NullSpace<double> userNullSpace( { mode1,mode2,mode3,mode4,mode5,mode6 } ); (9)
    return userNullSpace;
}

1	basis for translation \(\{(1,0,0),(0,1,0),(0,0,1)\}\)
2	basis for rotation \(\{(y,-x,0),(-z,0,x),(0,z,-y)\}\)
3	build the null space spanned by the translation and rotation basis

Finally, we export the solution for visualization in Paraview

    tic();
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->save();
    toc("Exporter");