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} \]
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>>
.
using namespace Feel;
Environment env( _argc=argc, _argv=argv,
_desc=feel_options(),
_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))) );
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);
auto e = exporter( _mesh=m1, _name="initial" );
e->step(0)->setMesh( m1 );
e->step(0)->add( "uinit", uVisu );
e->save();
meshMove( m1, uVisu );
auto e1 = exporter( _mesh=m1, _name="moved" );
e1->step(0)->setMesh( m1 );
e1->step(0)->add( "umoved", uVisu );
e1->save();