1. Bilinear and Linears Forms

We consider in this section bilinear and linear forms \(a: X_h \times X_h \rightarrow \mathbb{R}\) and \(\ell: X_h \rightarrow \mathbb{R}.\)

We suppose in this section that you know how to define your Mesh and your function spaces. You may need integration tools too, see Integrals.

There are Feel++ tools you need to create linear and bilinear forms in order to solve variational formulation.

from now on, u denotes an element from your trial function space (unknown function) and v an element from your test function space

1.1. Building Forms

1.1.1. Using form1

To construct a linear form \ell: X_h \rightarrow \mathbb{R}, proceed as follows

auto mesh = ...;
// build a P1/Q1 approximation space
auto Xh = Pch<1>( mesh );
auto l = form1(_test=Xh);
Table 1. Table of parameters

Name

Parameter

Description

Status

test

function space e.g. Xh

define test function space

Required

Here are some examples taken from the Feel++ tutorial.

From mylaplacian.cpp

// right hand side
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh), _expr=id(v));

From myadvection.cpp

// right hand side
auto l = form1( _test=Xh );
l+= integrate( _range=elements( mesh ), _expr=f*id( v ) );
The operators += and = are supported by linear and bilinear forms. 
auto a1 = form2(_test=Xh,_trial=Xh);
auto a2 = form2(_test=Xh,_trial=Xh);
// operations on a2 ...
// check that they have the same type and
// copy matrix associated to a2 in a1
a1 = a2;

1.1.2. Using form2

To define a bilinear form a: X_h \times X_h \rightarrow \mathbb{R}, for example a(u,v)=\int_\Omega uv

Building form2

The free-function form2 allows you to simply define such a bilinear form using the Feel++ language:

// define function space
auto Xh = ...;
// define a : Xh x Xh -> R
auto a = form2(_trial=Xh, _test=Xh );
// a(u,v) = \int_\Omega u v
a = integrate(_range=elements(mesh), _expr=idt(u)*id(v));
Table 2. Table of parameters

Name

Parameter

Description

Status

test

function space e.g. Xh

define test function space

Required

trial

function space e.g. Xh

define trial function space

Optional

Here are some examples taken from the Feel++ tutorial

From mylaplacian.cpp

// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
              _expr=gradt(u)*trans(grad(v)) );

From mystokes.cpp:

// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
              _expr=trace(gradt(u)*trans(grad(u))) );
a+= integrate(_range=elements(mesh),
              _expr=-div(u)*idt(p)-divt(u)*id(p));
see note above on operators += and =
Solving variational formulations

Once you created your linear and bilinear forms you can use the solve() member function of your bilinear form.

The following generic example solves: find u \in X_h \text{ such that } a(u,v)=l(v) \forall v \in X_h

Example
auto Xh = ...; // function space
auto u = Xh->element();
auto a = form2(_test=Xh, _trial=Xh);
auto l = form1(_test=Xh);

a.solve(_solution=u, _rhs=l, _rebuild=false, _name="");
Table 3. Table of arguments for solve()

Name

Parameter

Description

Status

_solution

element of domain function space

the solution

Required

_rhs

linear form

right hand side

Required

_rebuild

boolean(Default = false)

rebuild the solver components

Optional

_name

string(Default = "")

name of the associated Backend

Optional

Here are some examples from the Feel++ tutorial.

From laplacian.cpp
// solve the equation  a(u,v) = l(v)
a.solve(_rhs=l,_solution=u);
Using on for Dirichlet conditions

The function on() allows you to add Dirichlet conditions to your bilinear form before using the solve function.

The interface is as follows

Interface
on(_range=..., _rhs=..., _element=..., _expr=...);

Required Parameters:

  • _range domain concerned by this condition (see Integrals ).

  • _rhs right hand side. The linear form.

  • _element element concerned.

  • _expr the condition.

This function is used with += operator.

Here are some examples from the Feel++ tutorial.

From mylaplacian.cpp
// apply the boundary condition
a+=on(_range=boundaryfaces(mesh),
       _rhs=l,
       _element=u,
      _expr=expr(soption("functions.alpha")) );

There we add the condition: u = 0 \text{ on }\;\partial\Omega \;.

From mystokes.cpp
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
      _expr=expr<2,1,5>(u_exact,syms));

You can also apply boundary conditions per component:

Component-wise Dirichlet conditions
a+=on(_range=markedfaces(mesh,"top"),
      _element=u[Component::Y],
      _rhs=l,
      _expr=cst(0.))

The notation u[Component:Y] allows to access the Y component of u. Component::X and Component::Z are respectively the X and Z components.