Feel++ 0.91.0
|
In this example, we would like to solve for the following problem in 2D find such that
with
and is the exact solution
The following boundary conditions apply
We propose here two possible variational formulations. The first one, handles the Dirichlet boundary conditions strongly, that is to say the condition is {incorporated} into the function space definitions. The second one handles the Dirichlet condition {weakly} and hence we have a uniform treatment for all types of boundary conditions.
The variational formulation reads as follows, we introduce the spaces
and
We multiply (eq:1}) by then integrate over
and obtain
We integrate by parts and reformulate the problem as follows: we look for such that for all
In the present space setting~(eq:12}) and boundary conditions~(eq:4}), we have the boundary term from the integration by parts which is dropped being equal to 0
recalling that
where is the outward normal to
by convention.We now discretize the problem, we create a mesh out of
, we have
where can be segments, triangles or tetrahedra depending on
and we have
of them. We introduce the finite dimensional spaces of continuous piecewise polynomial of degree
functions
and
which are out trial and test function spaces respectively. We now have the problem we seek to solve which reads in our continuous Galerkin framework
we look for such that for all
There is an alternative formulation which allows to treat weakly Dirichlet(Essential) boundary conditions similarly to Neumann(Natural) and Robin conditions. Following a similar development as in the previous section, the problem reads
We look for such that for all
where
In (eq:16}), is defined by (eq:3}).
serves as a penalisation parameter which should be
, e.g. between 2 and 10, and
is the size of the face. The inconvenient of this formulation is the introduction of the parameter
, but the advantage is the {weak} treatment of the Dirichlet condition.
{sec:feel-formulation-1}
First we define the and
. To do that we use the new C++
auto
keyword and associate to f
and g
their corresponding expressions
//# marker1 # value_type pi = M_PI; auto g = sin(pi*Px())*cos(pi*Py())*cos(pi*Pz()); gproj = vf::project( Xh, elements(mesh), g ); auto f = pi*pi*Dim*g; //# endmarker1 #
where M_PI
is defined in the header cmath
. Using auto
allows to defined f
and g
--- which are moderately complex object --- without having to know the actual type. auto
determines automatically the type of the expression using the __typeof__
keyword internally.
Then we form the right hand side by defining a linear form whose algebraic representation will be stored in a vector_ptrtype
which is provided by the chosen linear algebra backend. The linear form is equated with an integral expression defining our right hand side.
//# marker2 # vector_ptrtype F( M_backend->newVector( Xh ) ); form1( _test=Xh, _vector=F, _init=true ) = integrate( elements(mesh), f*id(v) )+ integrate( markedfaces( mesh, mesh->markerName("Neumann") ), nu*gradv(gproj)*vf::N()*id(v) ); //# endmarker2 #
form1
generates an instance of the object representing linear forms, that is to say it mimics the mathematical object such that
which is represented algebraically in the code by the vector F
using the argument _vector
. The last argument _init
, if set to true
, will zero-out the entries of the vector F
.
_init
is set to false
by default.We now turn to the left hand side and define the bilinear form using the form2
helper function which is passed (i) the trial function space using the _trial
option, (ii) the test function space using the _test
option, (iii) the algebraic representation using _matrix
, i.e. a sparse matrix whose type is derived from one of the linear algebra backends and (iv) whether the associated matrix should initialized using _init
.
//# marker3 # sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) ); form2( Xh, Xh, D, _init=true ) = integrate( elements(mesh), nu*gradt(u)*trans(grad(v)) ); //# endmarker3 #
Finally, we deal with the boundary condition, we implement both formulation described in appendix~sec:vari-form-1}. For a {strong} treatment of the Dirichlet condition, we use the on()
keyword of FEEL++ as follows
//# marker5 # D->close(); form2( Xh, Xh, D ) += on( markedfaces(mesh, mesh->markerName("Dirichlet")), u, F, g ); //# endmarker5 #
Notice that we add, using +=
, the Dirichlet contribution for the bilinear form. The first argument is the set of boundary faces to apply the condition: in gmsh the points satisfying are marked using the flags
and
(
and
respectively.)
To implement the weak Dirichlet boundary condition, we add the following contributions to the left and right hand side:
//# marker41 # form1( _test=Xh, _vector=F ) += integrate( markedfaces(mesh,mesh->markerName("Dirichlet")), g*(-grad(v)*vf::N()+penaldir*id(v)/hFace()) ); //# endmarker41 #
//# marker10 # form2( Xh, Xh, D ) += integrate( markedfaces(mesh,mesh->markerName("Dirichlet")), -(gradt(u)*vf::N())*id(v) -(grad(v)*vf::N())*idt(u) +penaldir*id(v)*idt(u)/hFace()); D->close(); //# endmarker10 #
Note that we use the command line option --weakdir
set to 1 by default to decide between weak/strong Dirichlet handling. Apart the uniform treatment of boundary conditions, the weak Dirichlet formulation has the advantage to work also in a parallel environment.
Next we solve the linear system
where the solve
function is implemented as follows
Finally we check for the error in our approximation by computing
where is the exact solution and is equal to
and
is the numerical solution of the problem and the components of
in the
Lagrange basis are given by solving the linear system above.
The code reads
//# marker7 # double L2error2 =integrate(elements(mesh), (idv(u)-g)*(idv(u)-g) ).evaluate()(0,0); double L2error = math::sqrt( L2error2 ); Log() << "||error||_L2=" << L2error << "\n"; //# endmarker7 #
You can now verify that the error norm behaves like
where
is the mesh size and
the polynomial order. The
error norm would be checked similarly in
. The figure~fig:2} displays the results using Paraview.
![]()
Solution of the Laplacian problem | ![]()
Warped solution of the Laplacian problem |
We are now interested in solving the Stokes equations, we would like to solve for the following problem in 2D
find such that {