Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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(),
_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();

Generated on Sun Dec 22 2013 13:11:16 for Feel++ by doxygen 1.8.5