1. Defining and using expressions

The next step is to construct a function space over the mesh.

1.1. Step by step explanations

  • We start by loading a Mesh in 2D.

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

  • then we define some expression through the command line of config file: g is a scalar field and f is a vector field, here is an example how to enter them :

./feelpp_tut_myexpression --a=3 --functions.g="a*x*y:x:y:a" --functions.f="{sin(pi*x),cos(pi*y)}:x:y"

You can print back the expression to the screen to check that everything is ok. You want to use as expression a*x+b*y, you have to define a and b as option (either in your code, either in the library).

  • then we compute the gradient of g and f.

    auto grad_g=grad<2>(g);
    auto grad_f=grad(f);
    std::cout << "grad(g)=" << grad_g << std::endl;
    std::cout << "grad(f)=" << grad_f << std::endl;
template argument are given to grad to specify the shape of the gradient: in the case of \(\nabla g\), it is \(1\times2\) and \(2\times 2\) for \(\nabla f\) since we are in 2D.

  • then we compute the laplacian of g and f.

    auto laplacian_g=laplacian(g);
    std::cout << "laplacian(g)=" << laplacian_g << std::endl;

    auto laplacian_f=laplacian(f);
    std::cout << "laplacian(f)=" << laplacian_f << std::endl;

  • then we compute the divergence of f.

    auto div_f=div(f);
    std::cout << "div(f)=" << div_f << std::endl;

  • and the curl of f

    auto curl_f=curl(f);
    std::cout << "curl(f)=" << curl_f << std::endl;

  • Finally we evaluate these expressions at one point given by the option x and y.

1.2. Implementation

An implementation reads as follows:

#include <feel/feelcore/environment.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelvf/ginac.hpp>
using namespace Feel;

inline
po::options_description
makeOptions()
{
    po::options_description EXPRoptions( "DAR options" );
    EXPRoptions.add_options()
    ( "a", po::value<double>()->default_value( 1 ), "a parameter" )
    ( "b", po::value<double>()->default_value( 2 ), "a parameter" )
    ;
    return EXPRoptions;
}


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

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

    auto g = expr(soption(_name="functions.g"));
    std::cout << "g=" << g << std::endl;

    auto f = expr<2,1>(soption(_name="functions.f"));
    std::cout << "f=" << f << std::endl;

    double aVal = doption("a")+doption("b");
    std::map<std::string,double> myMap; myMap["aVal"]=aVal;
    auto i = expr(soption("functions.i"),myMap);
    std::cout << "i=" << i << std::endl;

    auto grad_g=grad<2>(g);
    auto grad_f=grad(f);
    std::cout << "grad(g)=" << grad_g << std::endl;
    std::cout << "grad(f)=" << grad_f << std::endl;

    auto laplacian_g=laplacian(g);
    std::cout << "laplacian(g)=" << laplacian_g << std::endl;

    auto laplacian_f=laplacian(f);
    std::cout << "laplacian(f)=" << laplacian_f << std::endl;


    auto div_f=div(f);
    std::cout << "div(f)=" << div_f << std::endl;

    auto curl_f=curl(f);
    std::cout << "curl(f)=" << curl_f << std::endl;

    std::cout << "Evaluation  at  (" << doption("x") << "," << doption("y") << "):" << std::endl;
    std::cout << "           g(x,y)=" << g.evaluate() << std::endl;
    std::cout << "           f(x,y)=" << f.evaluate() << std::endl;
    std::cout << "           i(x,y)=" << i.evaluate() << std::endl;
    std::cout << "Gradient:\n";
    std::cout << "     grad(g)(x,y)=" << grad_g.evaluate() << std::endl;
    std::cout << "     grad(f)(x,y)=" << grad_f.evaluate() << std::endl;
    std::cout << "Divergence:\n";
    std::cout << "      div(f)(x,y)=" << div_f.evaluate() << std::endl;
    std::cout << "Curl:\n";
    std::cout << "     curl(f)(x,y)=" << curl_f.evaluate() << std::endl;
    std::cout << "Laplacian:\n";
    std::cout << "laplacian(g)(x,y)=" << laplacian_g.evaluate() << std::endl;
    std::cout << "laplacian(f)(x,y)=" << laplacian_f.evaluate() << std::endl;
}

and the associated config file

a=12
b=-1
[functions]
g=(a-x)*x+(a/b)*y^3:x:y:a:b
f={1,1}:x:y
i=(x-aVal)*y:x:y:aVal

1.3. Execution

$ ./feelpp_tut_myexpression

or

$ ./feelpp_tut_myexpression --a=3 --functions.g="<your_function>" --functions.f="<your_function>"

We start with the following function g=1 and f=(1,1).

 $./feelpp_tut_myexpression --functions.g=1:x:y --functions.f="{1,1}:x:y"

and get something like this

g=1
f={1,1}
i=(x-aVal)*y
grad(g)=[[0,0]]
grad(f)=[[0,0],[0,0]]
laplacian(g)=[[0]]
laplacian(f)=[[0],[0]]
div(f)=[[0]]
curl(f)=[[0]]
Evaluation  at  (0,0):
           g(x,y)=1
           f(x,y)=1
1
           i(x,y)=-0
Gradient:
     grad(g)(x,y)=0 0
     grad(f)(x,y)=0 0
0 0
Divergence:
      div(f)(x,y)=0
Curl:
     curl(f)(x,y)=0
Laplacian:
laplacian(g)(x,y)=0
laplacian(f)(x,y)=0
0

The symbolic calculus system worked as expected.