1. Evaluating function

Once you have created an element, you may want to give it a value, that can depends on a lot of parameters ( mainly spaces, but others may apply ).

To do so, Feel++ relies on expressions. We may use various kind of expressions :

1.1. Built-in

Let’s begin with the evaluation of the expression \(sin(\pi x)\) on a unit circle.

First at all, we define the unit circle and its function space :

  // circle - geometrical order: 2
  auto mesh = unitCircle<2>();

  // circle - geometrical order: 1
  auto meshp1 = unitCircle<1>();

  auto Xh = Pch<2>( mesh );

Then the expression we would like to evaluate :

  auto myExpr = sin(pi*Px());

Px() refers to the variable x of our space.

With this,we can project it on our function space :

  auto v = project( _space=Xh, _range=elements(mesh), _expr=myExpr);

The expression will be evaluated on each point of our mesh.

In order to visualize the result, we create an exporter, named exhi, and add to it the projection.

  auto v = project( _space=Xh, _range=elements(mesh), _expr=myExpr);

1.1.1. Code

#include <feel/feelcore/environment.hpp>
#include <feel/feeldiscr/pch.hpp>
#include <feel/feeldiscr/operatorlagrangep1.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelfilters/unitcircle.hpp>

int main(int argc, char**argv )
{
  using namespace Feel;
  Environment env( _argc=argc, _argv=argv,
      _about=about(_name="myexporter",
        _author="Christophe Prud'homme",
        _email="christophe.prudhomme@feelpp.org"));

  // circle - geometrical order: 2
  auto mesh = unitCircle<2>();

  // circle - geometrical order: 1
  auto meshp1 = unitCircle<1>();

  auto Xh = Pch<2>( mesh );

  auto myExpr = sin(pi*Px());

  auto v = project( _space=Xh, _range=elements(mesh), _expr=myExpr);

  auto exhi = exporter( _mesh=mesh, _name="exhi" );
  auto exlo = exporter( _mesh=meshp1, _name="exlo" );
//  auto exhilo = exporter( _mesh=lagrangeP1(_space=Xh)->mesh(),_name="exhilo");

  int max = 10; double dt = 0.1;
  double time = 0;
  for (int i = 0; i<max; i++)
  {
    //exhilo->step( time )->add( "vhilo", v );
    exlo->step( time )->add( "vlo", idv(v) );
    exhi->step( time )->add( "vhi", v );
    time += dt;
    exhi->save();
    exlo->save();
    //exhilo->save();
  }
}

The list of the Feel++ Keyword is here.

1.2. Hard Coded

In this second method, we will use Functor :

struct Functor
{
    static const size_type context = vm::JACOBIAN|vm::POINT;
    typedef double value_type;
    typedef Feel::uint16_type uint16_type;
    static const uint16_type rank = 0;
    static const uint16_type imorder = 1;
    static const bool imIsPoly = true;
    myFunctor( int i ): coord(i){}
    double operator()( uint16_type a, uint16_type b, ublas::vector<double> const& x, ublas::vector<double> const& n ) const
        {
            return xtag::coord[];
        }
    int coord = 0;
};

We create a unit square meshed by triangles and we define the associated function space :

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Xh= Pch<2>(mesh);

From this space, we can define two elements, here one equals to the variable \(x\) and the other to the variable \(y\), obtain from Functor class.

    auto u = Xh.element(idf(Functor(0))); // Will contain x
    auto v = Xh.element(idf(Functor(1))); // Will contain y

The data exportation is the final step to visualize our expression \(x\) and \(y\) on the defined mesh.

    auto ex1 = exporter(_mesh=mesh,_name="ex1");
    ex1->add("u",u);
    ex1->add("v",v);
    ex1->save();

1.2.1. Code

The complete code reads as follows

#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feeldiscr/pch.hpp>
#include <feel/feeldiscr/pchv.hpp>
#include <feel/feelvf/projectors.hpp>
#include <feel/feeldiscr/operatorlagrangep1.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/function.hpp>
using namespace Feel;

namespace Feel {
struct Functor
{
    static const size_type context = vm::JACOBIAN|vm::POINT;
    typedef double value_type;
    typedef Feel::uint16_type uint16_type;
    static const uint16_type rank = 0;
    static const uint16_type imorder = 1;
    static const bool imIsPoly = true;
    myFunctor( int i ): coord(i){}
    double operator()( uint16_type a, uint16_type b, ublas::vector<double> const& x, ublas::vector<double> const& n ) const
        {
            return xtag::coord[];
        }
    int coord = 0;
};
}

int main(int argc, char**argv )
{
    Environment env( _argc=argc, _argv=argv,
                     _about=about(_name="myfunctor",
                                  _author="Feel++ Consortium",
                                  _email="feelpp-devel@feelpp.org"));


    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Xh= Pch<2>(mesh);


    auto u = Xh.element(idf(Functor(0))); // Will contain x
    auto v = Xh.element(idf(Functor(1))); // Will contain y

    auto ex1 = exporter(_mesh=mesh,_name="ex1");
    ex1->add("u",u);
    ex1->add("v",v);
    ex1->save();
}

Consult the list of Feel++ Mathematics keywords for more details.