load "msh3"

// Parameters
int nn = 20; // Mesh quality

// Mesh
int[int] labs = [1, 2, 2, 1, 1, 2]; // Label numbering
mesh3 Th = cube(nn, nn, nn, label=labs);
// Remove the ]0.5,1[^3 domain of the cube
Th = trunc(Th, (x < 0.5) | (y < 0.5) | (z < 0.5), label=1);

// Fespace
fespace Vh(Th, P1);
Vh u, v;

// Macro
macro Grad(u) [dx(u), dy(u), dz(u)] //

// Define the weak form and solve
solve Poisson(u, v, solver=CG)
    = int3d(Th)(
          Grad(u)' * Grad(v)
          1 * v
    + on(1, u=0)

// Plot
plot(u, nbiso=15);

A high level multiphysics finite element software

FreeFEM offers a fast interpolation algorithm and a language for the manipulation of data on multiple meshes.

Easy to use PDE solver

FreeFEM is a popular 2D and 3D partial differential equations (PDE) solver used by thousands of researchers across the world.

It allows you to easily implement your own physics modules using the provided FreeFEM language. FreeFEM offers a large list of finite elements, like the Lagrange, Taylor-Hood, etc., usable in the continuous and discontinuous Galerkin method framework.

Pre-built physics

Strong mesh and parallel capabilities

FreeFEM has it own internal mesher, called BAMG, and is compatible with the best open-source mesh and visualization software like Tetgen, Gmsh, Mmg and ParaView.

Written in C++ to optimize for speed, FreeFEM is interfaced with the popular mumps, PETSc and HPDDM solvers.

HPC in the cloud integration

With Qarnot's HPC platform, 7 lines of python code is all you need to run a FreeFEM simulation in the cloud. Learn how to run FreeFEM with Qarnot's sustainable HPC platform on Qarnot's blog.

FreeFEM is also available on Rescale's ScaleX® Pro. Rescale offers academic users up to 500 core hours on their HPC cloud.

Video tutorials

Thanks to Mojtaba Barzegari

Latest Articles

April 02, 2024

Topology optimization of three-dimensional structures subject to self-weight loading

Jorge Morvan Marotte Luz Filho, Antonio Andre Novotny

PurposeTopology optimization of structures under self-weight loading is a challenging problem which has received increasing attention in the past years. The use of standard formulations based on compliance minimization under volume constraint suffers from numerous difficulties for self-weight dominant scenarios, such as non-monotonic behaviour of the compliance, possible unconstrained character of the optimum and parasitic effects for low densities in density-based approaches. This paper aims to propose an alternative approach for dealing with topology design optimization of structures into three spatial dimensions subject to self-weight loading.Design/methodology/approachIn order to overcome the above first two issues, a regularized formulation of the classical compliance minimization problem under volume constraint is adopted, which enjoys two important features: (a) it allows for imposing any feasible volume constraint and (b) the standard (original) formulation is recovered once the regularizing parameter vanishes. The resulting topology optimization problem is solved with the help of the topological derivative method, which naturally overcomes the above last issue since no intermediate densities (grey-scale) approach is necessary.FindingsA novel and simple approach for dealing with topology design optimization of structures into three spatial dimensions subject to self-weight loading is proposed. A set of benchmark examples is presented, showing not only the effectiveness of the proposed approach but also highlighting the role of the self-weight loading in the final design, which are: (1) a bridge structure is subject to pure self-weight loading; (2) a truss-like structure is submitted to an external horizontal force (free of self-weight loading) and also to the combination of self-weight and the external horizontal loading; and (3) a tower structure is under dominant self-weight loading.Originality/valueAn alternative regularized formulation of the compliance minimization problem that naturally overcomes the difficulties of dealing with self-weight dominant scenarios; a rigorous derivation of the associated topological derivative; computational aspects of a simple FreeFEM implementation; and three-dimensional numerical benchmarks of bridge, truss-like and tower structures.

April 01, 2024

Parallel finite-element codes for the Bogoliubov-de Gennes stability analysis of Bose-Einstein condensates

Georges Sadaka, Pierre Jolivet, E. Charalampidis, I. Danaila

We present and distribute a parallel finite-element toolbox written in the free software FreeFem for computing the Bogoliubov-de Gennes (BdG) spectrum of stationary solutions to one- and two-component Gross-Pitaevskii (GP) equations, in two or three spatial dimensions. The parallelization of the toolbox relies exclusively upon the recent interfacing of FreeFem with the PETSc library. The latter contains itself a wide palette of state-of-the-art linear algebra libraries, graph partitioners, mesh generation and domain decomposition tools, as well as a suite of eigenvalue solvers that are embodied in the SLEPc library. Within the present toolbox, stationary states of the GP equations are computed by a Newton method. Branches of solutions are constructed using an adaptive step-size continuation algorithm. The combination of mesh adaptivity tools from FreeFem with the parallelization features from PETSc makes the toolbox efficient and reliable for the computation of stationary states. Their BdG spectrum is computed using the SLEPc eigenvalue solver. We perform extensive tests and validate our programs by comparing the toolbox's results with known theoretical and numerical findings that have been reported in the literature.

March 29, 2024

Stationary Flow of Incompressible Viscous Fluid Through a Two-Dimensional Branched Channel

Elena V. Shiryaeva, Alina S. Shokareva, Valentina P. Sibil

The results of a numerical experiment for the problem of an incompressible viscous fluid stationary flow through a branched planar (two-dimensional) channel are presented. The region in which the flow occurs simulates either blood vessels or a river delta. The finite element method and a modification of the penalty method, as well as the splitting method, are used for calculations. The implementation of the calculation algorithm is using with the help the package FreeFem++. The main goal, in addition, of course, to study the properties and structure of the stationary flow, is to demonstrate the effectiveness of the proposed modification of the penalty method. It is assumed that the region has one input boundary section through which the liquid flows into the region, and several (five) boundary sections through which the liquid flows out of the region. The remaining sections of the region boundary are considered impermeable to liquid. The boundary condition corresponding to the Poiseuille flow in a plane rectilinear channel is set at the input section. Three types of boundary conditions are considered on the output boundary secctions. 1. The boundary conditions correspond to the conditions of conservation of motion of luid particles i. e. the material derivative of the velocity is zero. 2. The boundary conditions correspond to the setting of the flow velocity. 3. The boundary conditions correspond to the setting of the same pressure. The stationary solution is constructed by the relaxation method. In fact, a non-stationary problem is solved over a sufficiently large time interval. As an initial condition for a non-stationary problem, a flow is chosen that is a Poiseuille flow in some region near the input boundary. The dependence of the convergence rate of the relaxation method on the initial data is investigated. It is found that the Poiseuille flow, given at the input section of the boundary of the region, induces similar Poiseuille flows at the output boundaries sections of the region in all these cases of boundary conditions. For a region with five sections of the output boundaries for some configuration of the region, the presence of stationary vortex flows (‘maelstorm’) is found.

December 01, 2023

Density-Based Topological Optimization of 3D-Printed Casts for Fracture Treatment with Freefem Software

K. Kokars, A. Krauze, K. Muiznieks, J. Virbulis, P. Verners, A. Gutcaits, J. Olins

Abstract 3D printed plastic casts can be used for healing bone fractures. The main requirements for these cases are: they should be light, require little printing time, have good mechanical properties, and ensure proper skin ventilation. We present a density-based topology optimization algorithm for obtaining optimal cast shapes that fulfil these requirements. The algorithm uses a linear stress model and simplified boundary conditions to model the contact problems. The cast shapes were optimized against the influence of several sharp corners. The parametric studies showed that the mass of optimized casts was reduced by 20 %–25 % in comparison with original industrial casts, and the printing time is reduced by 1.4–1.7 h for the largest cast. A major model drawback is the use of 3D numerical volume to model the density distribution. The density distribution should be homogenized across the cast layer. The overhang problem should also be addressed. We also suggest that the cast producers collect more experimental data on the cast breakages for a better calibration of the numerical model.

November 24, 2023

On the instability of travelling wave solutions for the transport-Stokes equation

M. Bonnivard, Amina Mecherbet

In this paper, we investigate the instability of the spherical travelling wave solutions for the Transport-Stokes system in $\mathbb{R}^3$. First, a classical scaling argument ensures instability among all probability measures for the Wasserstein metric and the $L^1$ norm. Secondly, we investigate the instability among patch solutions with a perturbed surface. To this end, we study the linearized system of a contour dynamics equation derived in [18] in the case where the support of the patch is axisymmetric and is described by spherical parametrization and show well-posedness of the linearized system. The actual version does not include the investigation of the eigenvalue and will be completed in the future. Eventually we investigate numerically the instability of the travelling wave by solving the Transport-Stokes equation using a finite element method on FreeFem.

November 23, 2023

Optimal Control Strategies for Mitigating Urban Heat Island Intensity in Porous Urban Environments

Nacer Sellila, M. Louaked, Waleed Mouhali, Houari Mechkour

This work is intended as an attempt to explore the use of optimal control techniques for designing green spaces and for dealing with the environmental problems related to urban heat islands appearing in cities. A three-dimensional model is established for numerical studies of the effects of urban anthropogenic heat and wind velocity in urban and rural regions. The transport mechanism of fluid in the cities is governed by the Navier–Stokes–Forschheimer porous media system. We introduce the penalty approximation method to overcome the difficulty induced by the incompressibility constraint. The partial differential equation optimal control problem is solved by using a Spectral Projected Gradient algorithm. To validate the method, we implement a numerical scheme, based on a finite element method, employing the free software FreeFem++ 14.3. We show the results for the optimized and non-optimized situations to compare the two cases.

November 19, 2023

Structure-preserving semi-convex-splitting numerical scheme for a Cahn-Hilliard cross-diffusion system in lymphangiogenesis

A. Jüngel, Boyi Wang

A fully discrete semi-convex-splitting finite-element scheme with stabilization for a degenerate Cahn-Hilliard cross-diffusion system is analyzed. The system consists of parabolic fourth-order equations for the volume fraction of the fiber phase and the solute concentration, modeling pre-patterning of lymphatic vessel morphology. The existence of discrete solutions is proved, and it is shown that the numerical scheme is energy stable up to stabilization, conserves the solute mass, and preserves the lower and upper bounds of the fiber phase fraction. Numerical experiments in two space dimensions using FreeFEM illustrate the phase segregation and pattern formation.

November 17, 2023


S. Chabbar, A. Habbal, R. Aboulaich, Nabil Ismaili, Sanaa El Majjaoui

Prostate cancer is a hormone-dependent cancer characterized by two types of cancer cells, androgen-dependent cancer cells and androgen-resistant ones. The objective of this paper is to present a novel mathematical model for the treatment of prostate cancer under combined hormone therapy and brachytherapy. Using a system of partial differential equations, we quantify and study the evolution of the different cell densities involved in prostate cancer and their responses to the two treatments. Numerical simulations of tumor growth under different therapeutic strategies are explored and presented. The numerical simulations are carried out on FreeFem++ using a 2D finite element method.

November 09, 2023

Numerical Modeling for Shoulder Injury Detection Using Microwave Imaging

Sahar Borzooei, Pierre-Henri Tournier, V. Dolean, Christian Pichot, N. Joachimowicz, Hélène Roussel, C. Migliaccio

A portable imaging system for the on-site detection of shoulder injury is necessary to identify its extent and avoid its development to severe condition. Here, firstly a microwave tomography system is introduced using state-of-the-art numerical modeling and parallel computing for imaging different tissues in the shoulder. The results show that the proposed method is capable of accurately detecting and localizing rotator cuff tears of different size. In the next step, an efficient design in terms of computing time and complexity is proposed to detect the variations in the injured model with respect to the healthy model. The method is based on finite element discretization and uses parallel preconditioners from the domain decomposition method to accelerate computations. It is implemented using the open source FreeFEM software.

September 25, 2023

Numerical Study of a Two-Dimensional Symmetric Flow of an Incompressible Fluid with a Coordinate-Dependent Viscosity Between Partially Irregular Planes

Natalia M. Polyakova, Victoria I. Tsvetkova

The structure of a symmetric two-dimensional unsteady flow between two planes having a periodic step-like irregularity is investigated by numerical methods using the finite element package Freefem++. The kinematic viscosity of the fluid is chosen as typical for turbulent flows, depending on the coordinates as 𝜇 = 𝜇0(𝜂2(𝑥) − 𝑧2), 𝜂(𝑥) = ℎ(𝑥) + 𝛿, where 𝜇0 is the characteristic viscosity, 𝛿 is the roughness, ℎ(𝑥) is the function defining the boundary 𝑧 = 𝜂(𝑥) (profiles of planes mirror-symmetric with respect to 𝑧 = 0), 𝑥, 𝑧 are longitudinal and transverse coordinates along the main flow. In contrast to the usual parabolic profile of the Poiseuille flow between the planes, with the specified choice of viscosity, the flow profile is logarithmic and has a singularity at the boundary of the region, to prevent which the roughness of the boundary is introduced. A distinctive feature of the problem under study is the setting of boundary conditions on an irregular part of the boundary. On plane sections of the boundary, we maintain the usual non-slip conditions for a viscous liquid, but on the irregular part of the boundary we set only the impermeability conditions of the boundary for the liquid (no sticking!). Computational experiment - numerical solution of the Navier-Stokes equations for a viscous incompressible fluid using a modified penalty method, showed that for a wide set of parameters, the structure of the flow described by the initial non-stationary complete problem and the quasi-stationary flow described by the asymptotic model constructed in the presented paper, for which there is an exact solution, are in good agreement - in the non-stationary flow, a stable regular system of vortices is established concentrated near those surface irregularities for which the boundary has a negative curvature (`deepening’of boundary).


on first wednesday of the month

Open Visio discussion

with Zoom form 10h to 12h paris time see for detail?

7-8 DECEMBER 2023

FreeFEM Days

Save te date and joint us for the 15th FreeFEM Days 2023 edition !
Paris, France

You are in good company

Sorbonne université INRIA ANR Genci CNRS