Function Spaces in Feel++

1. Introduction

Function spaces are fundamental to finite element methods in Feel++. They define the mathematical framework for approximating solutions to partial differential equations using polynomial basis functions over a mesh.

2. Basic Function Space Types

Feel++ provides several pre-defined function space types:

2.1. Scalar Function Spaces

For scalar-valued functions (temperature, pressure, density):

auto Vh1 = Pch<1>(mesh);  // Continuous piecewise linear (P1)
auto Vh2 = Pch<2>(mesh);  // Continuous piecewise quadratic (P2)
auto Vh3 = Pch<3>(mesh);  // Continuous piecewise cubic (P3)

2.2. Vector Function Spaces

For vector-valued functions (velocity, displacement):

auto Vh_vec = Pchv<1>(mesh);  // Vector P1 elements
auto Vh_vec2 = Pchv<2>(mesh); // Vector P2 elements

2.3. Mixed Function Spaces

For coupled systems (Stokes: velocity + pressure):

auto Vh = Pchv<2>(mesh);     // Velocity space
auto Qh = Pch<1>(mesh);      // Pressure space
auto Wh = Vh * Qh;           // Mixed space

3. Function Space Properties

Each function space provides important information:

auto Vh = Pch<2>(mesh);

LOG(INFO) << "Degrees of freedom: " << Vh->nDof();
LOG(INFO) << "Local DoFs per element: " << Vh->nLocalDof();
LOG(INFO) << "Polynomial order: " << Vh->basis()->nOrder;
LOG(INFO) << "Space dimension: " << Vh->worldComm().globalSize();

4. Creating and Manipulating Elements

Elements are functions that belong to a function space:

4.1. Basic Element Creation

auto Vh = Pch<2>(mesh);
auto u = Vh->element("u");      // Create element named "u"
auto v = Vh->element();         // Anonymous element

4.2. Element Initialization

From Expressions:

u.on( _range=elements(mesh), _expr=sin(pi*Px())*cos(pi*Py()) );

From Functions:

u.on( _range=elements(mesh), _expr=expr("x*x+y*y:x:y") );

From Other Elements:

auto u2 = Vh->element();
u2 = u;  // Copy values

5. Projections

Projections are used to represent continuous functions in discrete function spaces:

5.1. L² Projection

auto Vh = Pch<2>(mesh);
auto u_proj = vf::project( _space=Vh,
                          _range=elements(mesh),
                          _expr=exp(-((Px()-0.5)*(Px()-0.5) + (Py()-0.5)*(Py()-0.5))/0.1) );

5.2. H¹ Projection

For functions requiring continuity of gradients:

auto u_h1 = vf::project( _space=Vh,
                        _range=elements(mesh),
                        _expr=sin(2*pi*Px())*sin(2*pi*Py()),
                        _grad_expr=vec(2*pi*cos(2*pi*Px())*sin(2*pi*Py()),
                                      2*pi*sin(2*pi*Px())*cos(2*pi*Py())) );

6. Working with Vector Functions

Vector function spaces handle multi-component fields:

auto Vh_vec = Pchv<2>(mesh);
auto u_vec = Vh_vec->element("velocity");

// Initialize with vector expression
u_vec.on( _range=elements(mesh),
          _expr=vec(cos(Px()), sin(Py())) );

// Access components
auto u_x = u_vec.comp(X);  // X-component
auto u_y = u_vec.comp(Y);  // Y-component

7. Mixed Function Spaces

Mixed spaces combine multiple function spaces for coupled problems:

7.1. Stokes System Example

auto Vh = Pchv<2>(mesh);    // Velocity (P2 vector)
auto Qh = Pch<1>(mesh);     // Pressure (P1 scalar)
auto Wh = Vh * Qh;          // Mixed space

// Create mixed element
auto U = Wh->element("(u,p)");
auto u = U.template element<0>("u");  // Velocity component
auto p = U.template element<1>("p");  // Pressure component

// Initialize
u.on( _range=elements(mesh), _expr=vec(sin(Px()), cos(Py())) );
p.on( _range=elements(mesh), _expr=Px()*Py() );

8. Advanced Function Space Features

8.1. Marked Elements

Create function spaces on subdomains:

// Function space on marked region
auto Vh_subdomain = Pch<2>( mesh, markedelements(mesh, "fluid") );

// Function space with level set
auto Vh_levelset = Pch<2>( mesh, elements(mesh, expr("x*x+y*y-0.25<0")) );

8.2. Periodic Boundary Conditions

For problems with periodic domains:

auto Vh_periodic = Pch<2>( mesh,
                          periodicity( Periodic<>(BoundaryCondition("left"), BoundaryCondition("right")) ) );

9. Norms and Error Computation

Function spaces provide various norm computations:

9.1. L² and H¹ Norms

auto u = Vh->element();
// ... initialize u ...

// L² norm: ||u||_{L²}
double l2_norm = normL2( _range=elements(mesh), _expr=idv(u) );

// H¹ norm: ||u||_{H¹} = (||u||²_{L²} + ||∇u||²_{L²})^{1/2}
double h1_norm = normH1( _range=elements(mesh),
                        _expr=idv(u),
                        _grad_expr=gradv(u) );

// H¹ semi-norm: |u|_{H¹} = ||∇u||_{L²}
double h1_semi = normH1( _range=elements(mesh),
                        _grad_expr=gradv(u) );

9.2. Error Norms

Compute errors between solutions:

auto u_exact = expr("sin(pi*x)*sin(pi*y):x:y");
auto u_numerical = Vh->element();
// ... solve for u_numerical ...

double error_l2 = normL2( _range=elements(mesh),
                         _expr=idv(u_numerical) - u_exact );

double error_h1 = normH1( _range=elements(mesh),
                         _expr=idv(u_numerical) - u_exact,
                         _grad_expr=gradv(u_numerical) - grad<2>(u_exact) );

10. Interpolation vs Projection

Understanding the difference between interpolation and projection:

10.1. Interpolation

Point-wise matching at interpolation nodes:

auto u_interp = Vh->element();
u_interp.on( _range=elements(mesh), _expr=sin(pi*Px())*cos(pi*Py()) );

10.2. Projection

Optimal approximation in the function space:

auto u_proj = vf::project( _space=Vh,
                          _range=elements(mesh),
                          _expr=sin(pi*Px())*cos(pi*Py()) );

11. Complete Example

Here’s a comprehensive example demonstrating function space usage:

#include <feel/feel.hpp>

int main( int argc, char** argv )
{
    using namespace Feel;
    using namespace Feel::vf;

    // Initialize Feel++ Environment with command line options
    po::options_description myoptions( "Function space tutorial options" );
    myoptions.add_options()
        ( "order", po::value<int>()->default_value( 2 ), "polynomial order" )
        ( "hsize", po::value<double>()->default_value( 0.1 ), "mesh size" );

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

    // Create mesh
    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    
    // Get polynomial order from command line
    int order = ioption(_name="order");
    LOG(INFO) << "Using polynomial order: " << order;

    // Create function spaces of different types
    // Scalar function space (temperature, pressure, etc.)
    auto Vh_scalar = Pch<1>( mesh );           // P1 elements
    auto Vh_scalar_ho = Pch<2>( mesh );        // P2 elements
    
    // Vector function space (velocity, displacement, etc.)
    auto Vh_vector = Pchv<1>( mesh );          // Vector P1 elements
    auto Vh_vector_ho = Pchv<2>( mesh );       // Vector P2 elements
    
    // Mixed function space (Stokes: velocity + pressure)
    auto Vh_mixed = Pchv<2>( mesh ) * Pch<1>( mesh );
    
    LOG(INFO) << "Scalar P1 space has " << Vh_scalar->nDof() << " degrees of freedom";
    LOG(INFO) << "Scalar P2 space has " << Vh_scalar_ho->nDof() << " degrees of freedom";
    LOG(INFO) << "Vector P1 space has " << Vh_vector->nDof() << " degrees of freedom";
    LOG(INFO) << "Vector P2 space has " << Vh_vector_ho->nDof() << " degrees of freedom";
    LOG(INFO) << "Mixed space has " << Vh_mixed->nDof() << " degrees of freedom";

    // Create elements (functions) in these spaces
    auto u_scalar = Vh_scalar->element("u_scalar");
    auto u_vector = Vh_vector->element("u_vector");
    
    // Initialize with expressions
    u_scalar.on( _range=elements(mesh), _expr=sin(pi*Px())*cos(pi*Py()) );
    u_vector.on( _range=elements(mesh), _expr=vec(cos(Px()), sin(Py())) );
    
    // Demonstrate projections
    auto proj_scalar = vf::project( _space=Vh_scalar_ho, 
                                   _range=elements(mesh), 
                                   _expr=exp(-((Px()-0.5)*(Px()-0.5) + (Py()-0.5)*(Py()-0.5))/0.1) );
    
    // Compute norms
    double l2_norm_scalar = normL2( _range=elements(mesh), _expr=idv(u_scalar) );
    double h1_norm_scalar = normH1( _range=elements(mesh), _expr=idv(u_scalar), _grad_expr=gradv(u_scalar) );
    
    LOG(INFO) << "L2 norm of scalar function: " << l2_norm_scalar;
    LOG(INFO) << "H1 norm of scalar function: " << h1_norm_scalar;
    
    // Export results
    auto e = exporter( _mesh=mesh, _name="functionspace" );
    e->add( "u_scalar", u_scalar );
    e->add( "u_vector", u_vector );
    e->add( "projected", proj_scalar );
    e->save();
    
    return 0;
}

12. Performance Considerations

12.1. Memory Usage

Function spaces store basis function information. For large problems: - Use appropriate polynomial orders (don’t over-refine) - Consider mixed formulations to reduce total DoFs - Use periodic boundary conditions when applicable

12.2. Parallel Distribution

In parallel computations, function spaces automatically handle: - DoF distribution across processes - Ghost element communication - Consistent global numbering

13. Best Practices

  1. Choose Appropriate Orders: Match polynomial order to solution smoothness

  2. Use Mixed Formulations: For coupled physics (Stokes, elasticity)

  3. Leverage Auto Types: Let compiler deduce function space types

  4. Initialize Properly: Always initialize elements before use

  5. Check Consistency: Ensure function spaces match your problem requirements

Function spaces are the foundation of finite element methods in Feel++. Understanding their creation, manipulation, and properties is essential for effective PDE solving.