Learn by examples 1: Membrane


This example is taken from the "learning by examples" section from the FreeFem++ manual. If you are new to FreeFem++ we recommend to look at the FreeFem++ manual first. Then you can come back here to study how this simple example needs to be redesigned to fit in the StabFem drivers.

Interfacing with StabFem ; principles

The example in the freefem manual is a single program which first builds the mesh and secondly solves the problem. In StabFem these two steps are generally split in two programs. The reasons are that we generally use a same mesh for a number of calculations ; moreover we are often using mesh adaptations during the course of a study, so it is better to design "mesh generation/adaptation" and "problem resolution" as independent bricks.

Hence we have two programs : - Membrane_Mesh.edp is used to build the mesh - Membrane.edp is solving the problem. you may look at explanations directly in these programs to see what has to be redisigned to use through the interface.


close all;
SF_Start('verbosity',4); % verbosity=4 will display output from the FreeFem program

NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead
WARNING - Disabling database management. Some advanced functionalities may not work properly

In this very first example the DataBase management for simplicity. see next example "SCRIPT_HeatExchanger.m" to see the difference.

Generation of the mesh

We use here the program 'Membrane_Mesh.edp' which is designed to build a mesh for an ellipical domain.

ffmesh = SF_Mesh('Membrane_Mesh.edp')
ffmesh = 

  struct with fields:

             points: [3x1360 double]
             bounds: [3x150 double]
                tri: [4x2568 double]
                 np: 1360
                nbe: 150
                 nt: 2568
             labels: [1 2]
        problemtype: 'unspecified'
           filename: './mesh.msh'
           meshtype: '2D'
           datatype: 'mesh'
    DataDescription: 'Mesh (generated by generic macro)'
     generationmode: 'initial'
    datastoragemode: 'unknown'
       TriangleSize: [1360x1 double]
                np2: 5287
          Vh_P2P2P1: [38520x1 double]
        Vh_P2P2P2P1: [53928x1 double]
              Vh_P2: [15408x1 double]
     meshgeneration: 0
           symmetry: 'N'

This is basically equivalent to launching directly "Freefem++ Membrane_Mesh.edp" in a terminal.

Note that the program Membrane_Mesh.edp produces a couple of files 'mesh.msh' and 'mesh.ff2m'. Here DataBase management is disabled hence these files are in the woking directory './'. The driver SF_Mesh imports these files as an object 'ffmesh' which can be used for subsequent calculations and plots.

Plot the mesh

figure(1);SF_Plot(ffmesh,'mesh','title','Mesh for an elliptical membrane');

Solve the Membrane Problem

VertD = SF_Launch('Membrane.edp', 'Mesh', ffmesh)
VertD = 

  struct with fields:

               mesh: [1x1 struct]
           filename: 'Data.txt'
    DataDescription: 'Vertical displacement in an elliptical domain'
           datatype: 'Membrane'
    datastoragemode: 'ReP2'
    datadescriptors: 'phi'
       meshfilename: './mesh.msh'
                phi: [5287x1 double]

This is basically equivalent to launching directly "Freefem++ Membrane.edp" in a terminal.

Note that the program Membrane.edp produces a couple of files 'Data.txt' and 'Data.ff2m'. Here DataBase management is disabled hence these files are in the woking directory './'.

The driver SF_Launch imports these files as an object 'VertD'. Note that this object contains a field 'phi' which can be plotted.

Plotting the result

figure(2);SF_Plot(VertD,'phi','title', 'Vertical displacement phi of spherical membrane','contour', 'on');
% add one line on the figure
hold on; plot([-1 1],[.5,.5],'r--')
hold off;

Post-processing : interpolation along a line

Here is how to extract data along a 1D line and plot this data

Ycut = 0.5;Xcut = -1:.01:1;
Cut = SF_ExtractData(VertD,Xcut,Ycut);
% Here the object 'Cut' has a unique field 'phi'

title('Vertical Displacement along a line : Phi(x,0.5)');xlabel('x');ylabel('Phi');
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Exercices :

  (You may either modify directly the value of b in 'Membrane_Mesh.edp',
   or learn how to pass a parameter "b" from the driver ;
   see next example <https://stabfem.gitlab.io/StabFem/TUTORIAL_CASES/Heat_Exchanger/SCRIPT_HeatExchanger.html Heat exchanger> to see how to do this !