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

Table of Contents

Author
Christophe Prud'homme
Date
2013-05-01

Theory

We are interested in solving $ -\nu \Delta \bvec{u} + \nabla p = \bvec{f},\quad \nabla \cdot \bvec{u} = 0 $ in curl-curl formulation in 2D.

First we introduce the following notations:

\[ u \times n = u_1 n_2 - u_2 n_!, \quad u = (u_1, u_2),\ n = (n_1, n_2) \]

The curl of a vectorial field $\phi$

\[ \mathrm{curl} \phi = (-\frac{\partial \phi_2}{\partial x_1}, -\frac{\partial \phi_1}{\partial x_2} ) \]

The curl of a scalar field $ $\psi$

\[ \nabla \times \psi = \frac{\partial \psi}{\partial x_1} -\frac{\partial \psi}{\partial x_2} \]

Implementation

int main(int argc, char**argv )
{
Environment env( _argc=argc, _argv=argv,
_desc=feel_options(),
_about=about(_name="stokes_curl",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));
// create the mesh
auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );
// function space
auto Vh = THch<2>( mesh );
// element U=(u,p) in Vh
auto U = Vh->element();
auto u = U.element<0>();
auto p = U.element<1>();
// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
_expr=curlxt(u)*curlx(u) + divt(u)*div(u) );
a+= integrate(_range=elements(mesh),
_expr=-div(u)*idt(p)-divt(u)*id(p));
a+= integrate(_range=markedfaces(mesh,"inlet"),
_expr=( curlxt(u)*cross(id(u),N()) +
curlx(u)*cross(idt(u),N()) +
30*cross(id(u),N())*cross(idt(u),N())/hFace() ) );
a+= integrate(_range=markedfaces(mesh,"outlet"),
_expr=( curlxt(u)*cross(id(u),N()) +
curlx(u)*cross(idt(u),N()) +
30*cross(id(u),N())*cross(idt(u),N())/hFace() ) );
// right hand side
auto l = form1( _test=Vh );
l += integrate(_range=markedfaces(mesh,"inlet"),
_expr=2*trans(id(u))*N() );
l += integrate(_range=markedfaces(mesh,"outlet"),
_expr=1*trans(id(u))*N() );
// Dirichlet condition
a+=on(_range=markedfaces(mesh,"wall"), _rhs=l, _element=u,
_expr=vec(cst(0.),cst(0.)));
// solve a(u,v)=l(v)
a.solve(_rhs=l,_solution=U);
// save results
auto e = exporter( _mesh=mesh );
e->add( "u", u );
e->add( "p", p );
e->save();
}

Generated on Fri Oct 25 2013 14:24:27 for Feel++ by doxygen 1.8.4