1. Spaces and elements

You’ve learned how to discretize the space you want to compute on. You now have to learn how to define and use function spaces and elements of functions spaces. For advanced informations on this subject, you can look in the Function Space documentation.

1.1. Constructing a function space

  • Loading a Mesh in 2D

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

  • For basic function spaces, we have predetermined constructors:

auto Xh = Pch<2>( mesh );

  • Defining an element

auto u = Xh->element( "u" );
auto w = Xh->element( "w" );

One can also use :

  • Pdh<ORDER>(mesh) : Polynomial Discontinuous

  • Pvh<ORDER>(mesh) : Polynomial Continuous Vectorial

  • Pdhv<ORDER>(mesh) : Polynomial Discontinuous Vectorial

  • Pchm<ORDER>(mesh) : Polynomial Continuous Matrix

  • Dh<ORDER>(mesh) : RaviartThomas function space

  • Ned1h<ORDER>(mesh) : Nedelec function spaces

1.2. Implementation

The implementation reads are follows

#include <feel/feelcore/environment.hpp>
#include <feel/feeldiscr/pch.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/vf.hpp>

using namespace Feel;

int main( int argc, char** argv )
{

    //Initialize Feel++ Environment
    Environment env( _argc=argc, _argv=argv,
                     _about=about( _name="myfunctionspace",
                                   _author="Feel++ Consortium",
                                   _email="feelpp-devel@feelpp.org" )  );

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

    // function space \( X_h \) using order 2 Lagrange basis functions
    auto Xh = Pch<2>( mesh );

    auto g = expr<4>( soption(_name="functions.g"));
    auto gradg = grad<3>(g);

    // elements of \( u,w \in X_h \)
    auto u = Xh->element( "u" );
    auto w = Xh->element( "w" );

		// build the interpolant of u
    u.on( _range=elements( mesh ), _expr=g );

		// build the interpolant of the interpolation error
    w.on( _range=elements( mesh ), _expr=idv( u )-g );

    // compute L2 norms \(||\cdot||_{L^2}\)
    double L2g = normL2( elements( mesh ), g );
    double H1g = normL2( elements( mesh ), _expr=g,_grad_expr=gradg );
    double L2uerror = normL2( elements( mesh ), ( idv( u )-g ) );
    double H1uerror = normH1( elements( mesh ), _expr=( idv( u )-g ),
                              _grad_expr=( gradv( u )-gradg ) );
    
    Feel::cout << "||u-g||_0 = " << L2uerror/L2g << std::endl;
    Feel::cout << "||u-g||_1 = " << H1uerror/H1g << std::endl;

    // export for post-processing
    auto e = exporter( _mesh=mesh );

	 	// save interpolant
    e->add( "g", u );

	 	// save interpolant of interpolation error
    e->add( "u-g", w );

    e->save();

}

and the associated config file

[gmsh]
# default is hypercube
shape=hypercube
# you can try ellipsoid
#shape=ellipsoid

[functions]
g=x*y*(y*2.2+4.3)*(2.1*x-1.3):x:y