Problem

This program feelpp_doc_harmonic shows how to solve for the harmonic extension of a displacement given at a boundary of a domain. Given $\eta$ on $\partial \Omega$ , find $\eta^H$ such that

$\begin{split} -\Delta \eta^H &= 0 \mbox{ in } \Omega\\ \eta^H &=\eta \mbox{ on } \partial \Omega. \end{split}$

$\eta^H$ is the harmonic extension of $\eta$ in $\Omega$ , then we move the mesh vertices using $\eta^H$ .

Remarks: Note that the degrees of freedom of $\eta^H$ must coincide with the mesh vertices which implies that the mesh order (the order of the geometric transformation) should be the same as $\eta^H$ . If $\eta^H$ is piecewise polynomial of degree then the geometric transformation associated to the mesh should be of degree too, i.e. Mesh<Simplex<2,N>>.

Results

We consider the following displacement

$\eta = ( 0, 0.08 (x+0.5) (x-1) (x^2-1) )^T$

on the bottom boundary and $\eta$ being 0 on the remaining borders and we look for $\eta^H$ as a $P_2$ piecewise polynomial function in $\Omega$ and the mesh associated is of the same order i.e. Mesh<Simplex<2,2>>.

Implementation

The implementation is as follows

    using namespace Feel;
    Environment env( _argc=argc, _argv=argv,
                     _desc=feel_options(),
                     _directory=".",
                     _about=about(_name="harmonic",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));
    auto mesh = createGMSHMesh( _mesh=new Mesh<Simplex<2,2> >,
                                _desc=domain( _name="harmonic",
                                              _usenames=false,
                                              _shape="hypercube",
                                              _h=option(_name="gmsh.hsize").as<double>(),
                                              _xmin=-0.5, _xmax=1,
                                              _ymin=-0.5, _ymax=1.5 ) );
    auto Vh = Pchv<2>( mesh );
    auto u = Vh->element();
    auto v = Vh->element();
    auto l = form1( _test=Vh );
    auto a = form2( _trial=Vh, _test=Vh );
    a = integrate(_range=elements(mesh),
                  _expr=trace(gradt(u)*trans(grad(v))) );
    // boundary conditions
    a+=on(_range=markedfaces(mesh,1), _rhs=l, _element=u,
          _expr=zero<2,1>() );
    a+=on(_range=markedfaces(mesh,3), _rhs=l, _element=u,
          _expr=zero<2,1>() );
    a+=on(_range=markedfaces(mesh,4), _rhs=l, _element=u,
          _expr=zero<2,1>() );
    a+=on(_range=markedfaces(mesh,2), _rhs=l, _element=u,
          _expr=vec(cst(0.),0.08*(Px()+0.5)*(Px()-1)*(Px()*Px()-1)));
    a.solve(_rhs=l,_solution=u);
    auto m1 = lagrangeP1(_space=Vh)->mesh();
    auto XhVisu = Pchv<1>(m1);
    auto opIVisu = opInterpolation(_domainSpace=Vh,
                                   _imageSpace=XhVisu,
                                   _type=InterpolationNonConforme(false,true,false) );
    auto uVisu = opIVisu->operator()(u);
    // exporter mesh and harmonic extension
    auto e = exporter( _mesh=m1, _name="initial" );
    e->step(0)->setMesh( m1 );
    e->step(0)->add( "uinit", uVisu );
    e->save();
    // move the mesh vertices
    meshMove( m1, uVisu );
    // export mesh after moving the vertices
    auto e1 = exporter( _mesh=m1, _name="moved" );
    e1->step(0)->setMesh( m1  );
    e1->step(0)->add( "umoved", uVisu );
    e1->save();

Table of Contents

Problem

Results

Implementation