Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Pages
Harmonic extension of a boundary displacement

Table of Contents

Author
Christophe Prud'homme
Date
2013-03-13



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 \(N\) then the geometric transformation associated to the mesh should be of degree \(N\) 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(),
_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();