Tutorial example demonstrating the basic functionalities of the StabFem Software

This tutorial script is designed as a support for chapter 2 of the StabFem manual available here

We show how to compute, plot and manipulate the solution of three simple problems formulated in a L-shaped domain, namely steady and unsteady heat conduction, and Navier-Stokes flow with bottom-driven boundary conditions.

Source code for this program : SCRIPT_Lshape.m

NB in this simple tutorial the DataBase management is disabled. After studying the present example you may have a look at the related example SCRIPT_Lshape_DataBase.m

Contents

Initialization

close all;
addpath([fileparts(fileparts(pwd)) '/SOURCES_MATLAB']);

SF_Start('verbosity',3);
% higest levels of verbosity results in more messages.
% recommended values are 2 (minimal notice) and 4 (notice + freefem output). 6 and more is debug mode.

SF_DataBase('create','./WORK/');
% this is to initialize the database management. All results will be placed
% in a folder "WORK" and subfolders.

mkdir('FIGURES');
set(0,'defaulttextfontsize',16); % set font size for figures
set(0,'defaultaxesfontsize',14); % idem
NOTICE  - Initialization already done
NOTICE  - Verbosity set to 3 : Notice messages
NOTICE  - Database directory :./WORK/
NOTICE  - Working dir does not exist : create it
NOTICE  - Database directory :./WORK/
NOTICE  - This directory already exists
NOTICE  - Please type "SF_Status" to see what is available in this dir
NOTICE  -     or type "SF_core_arborescence('cleanall')" to erase any previous data
NOTICE  - CLEANALL in SF_core_arborescence ; working dir ./WORK/ is now empty
NOTICE  - Cleaning ALL arborescence ./WORK/
Warning: Directory already exists. 

Generation of the mesh

The program Lshape_Mesh.edp is a simple mesh generator producing a mesh.msh file. We launch this program and import the mesh as a structure using the SF_Mesh.m driver :

Ndensity =40;
ffmesh = SF_Mesh('Lshape_Mesh.edp','Options',{'meshdensity',Ndensity})
NOTICE  -       ### INITIAL MESH CREATED WITH np = 1429 vertices

ffmesh = 

  struct with fields:

             points: [3x1429 double]
             bounds: [3x160 double]
                tri: [4x2696 double]
                 np: 1429
                nbe: 160
                 nt: 2696
             labels: [1 2]
        problemtype: 'unspecified'
           filename: './WORK/MESHES/FFMESH_000001.msh'
           meshtype: '2D'
           datatype: 'mesh'
    DataDescription: 'Mesh (generated by generic macro)'
     generationmode: 'initial'
    datastoragemode: 'unknown'
       TriangleSize: [1429x1 double]
                np2: 5553
          Vh_P2P2P1: [40440x1 double]
        Vh_P2P2P2P1: [56616x1 double]
              Vh_P2: [16176x1 double]
     meshgeneration: 0
           symmetry: 'N'

NB this sequence is equivalent to executing the program with FreeFem++ in a shell terminal as follows :

> FreeFem++ Lshape_Mesh.edp -meshdensity 40

(but in addition the function returns a "handle' to the mesh which can be used for subsequent computations, plotted, etc...)

How to plot the mesh

figure();SF_Plot(ffmesh,'mesh','title','Mesh for L-shape body');
saveas(gca,'FIGURES/Lshape_Mesh','png'); % to save the figure
pause(0.1);

NB the SF_Plot.m function, which is built upon the ffpdeplot.m function developped by Markus Meisinger, allows a variety of plots, including mesh plots, data plots using color levels, contour plots, or quiver plots. A large number of options are available, see the function or type "help SF_Plot".

First problem : steady conduction

Here is how to launch a simple FreeFem++ program Lshape_Steady.edp and import the data as a structure using the SF_Launch.m driver. Have a look at these programs to understand how it works !

heatS=SF_Launch('Lshape_Steady.edp','Mesh',ffmesh,'Store','HeatSteady')
WARNING -  Folder ./WORK//HeatSteady mentioned for database storing does not exist ; creating it
NOTICE  - Storing result in folder HeatSteady

heatS = 

  struct with fields:

               mesh: [1x1 struct]
           filename: './WORK/HeatSteady/FFDATA_000001.txt'
    DataDescription: 'Temperature in a L-Shape domain'
           datatype: 'HeatSteady'
    datastoragemode: 'ReP2'
    datadescriptors: 'T'
       meshfilename: './WORK/MESHES/FFMESH_000001.msh'
                  T: [5553x1 double]

NB this sequence is equivalent to executing the program with FreeFem++ in a shell terminal as follows :

> FreeFem++ Lshape_Steady.edp

(but in addition the data is imported as an object which can be used for subsequent computations, plotted, etc...)

Note that the data is imported as a structure witch includes plottable fields (here "T", "dxT", "dyT"), numerical fields (here "HeatFlux") along as a number of string-valued descriptor fields. All this data is read from the file Data.ff2m generated by the FreeFem program ; the .ff2m format is designed for easy customization.

Note the option 'store' which tells the driver to store the results in a folder 'HeatSteady' (see last section about database management)

Plot the temperature field

figure(1);SF_Plot(heatS,'T','xlim',[0 1],'ylim',[0 1],'title', 'Steady conduction on a L-shaped domain: Temperature field ');
% add one line on the figure
hold on; plot([0 1],[.25,.25],'r--')
saveas(gca,'FIGURES/Lshape_T0','png');
hold off;

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

Ycut = 0.25;Xcut = 0:.01:1;
Tcut = SF_ExtractData(heatS,'T',Xcut,Ycut);

figure(2);plot(Xcut,real(Tcut),'r-',Xcut,imag(Tcut),'b--');
title('Temperature along a line : T(x,0.25)');xlabel('x');ylabel('T');
saveas(gca,'FIGURES/Lshape_T0_Cut','png')
pause(0.1);
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Second problem : unsteady conduction

Assuming harmonic dependance as $T(x,y,t) = Re(\hat{T}(x,y) e^{-i \omega t})$, this leads to:

$- i \omega \hat{T} = \kappa \Delta \hat{T}$ on domain, $\hat{T} = 1$ along boundary

Have a look at the program Lshape_Unsteady.edp to see how this problem is solved by FreeFem++ !

The program computes the field $\hat{T}(x,y)$ as well as the complex heat flux such as $\phi(t) = Re(\hat{\phi} e^{-i \omega t})$ ; namely : $\hat{\phi} = \int_{\partial \Omega} \nabla \hat{T} \cdot {\bf n}$.

First we solve the problem for a single value of $\omega$:

omega = 100;
heatU=SF_Launch('Lshape_Unsteady.edp','Options',{'omega',omega,'kappa',1},'Mesh',ffmesh,'DataFile','Heat_Unsteady.ff2m','Store','HeatSteady')
NOTICE  - Storing result in folder HeatSteady

heatU = 

  struct with fields:

               mesh: [1x1 struct]
           filename: './WORK/HeatSteady/FFDATA_000002.txt'
    DataDescription: 'Temperature field in a L-shaped region ; Unsteady conduction with harmonic forcing at boundaries'
           datatype: 'HeatUnsteady'
    datastoragemode: 'CxP2'
    datadescriptors: 'T'
       meshfilename: './WORK/MESHES/FFMESH_000001.msh'
              omega: 100
           INDEXING: [1x1 struct]
           HeatFlux: 17.1650 -21.3108i
          normGradT: [1429x1 double]
                  T: [5553x1 double]

NB this sequence is equivalent to executing the program with FreeFem++ in a shell terminal as follows :

> FreeFem++ Lshape_Unsteady.edp -omega 40 -kappa 1

the parameters passed as "-omega 40 -kappa 1" are detected within the program using the "getARGV" Freefem function.

We plot the results

figure(); SF_Plot(heatU,'T.re','colormap','redblue','title',{'Unsteady heat condition problem ' , 'real(colors) and imaginary(levels) parts of T_1(x)'});
hold on;  SF_Plot(heatU,'T.im','contour','on','xystyle','off');
saveas(gca,'FIGURES/Lshape_Tc','png');

Note that resulting value 'heatU' has a field 'HeatFlux' corresponding to $\hat{\phi}$

heatU.HeatFlux
ans =

  17.1650 -21.3108i

Adapts the mesh to both steady and unsteady solutions

We now use the SF_Adapt.m function to adapt the mesh to the results. This function (which is a wrapper for the FreeFem++ program AdaptMesh.edp) allows to perform mesh adaptation using up to 8 fields with various data type (scalar P1 or P2 data ; composite data for instance [P2xP2xP1] ,...) Have a look at these programs to understand how it works !

ffmesh_adapted = SF_Adapt(ffmesh,heatS,heatU)
NOTICE  - SF_Adapt : created new mesh ; np = 3331 vertices 

ffmesh_adapted = 

  struct with fields:

             points: [3x3331 double]
             bounds: [3x236 double]
                tri: [4x6424 double]
                 np: 3331
                nbe: 236
                 nt: 6424
             labels: [1 2]
        problemtype: 'unspecified'
           filename: './WORK/MESHES/FFMESH_000002.msh'
           meshtype: '2D'
           datatype: 'mesh'
    DataDescription: 'Mesh (generated by generic macro)'
     generationmode: 'adapted'
    datastoragemode: 'unknown'
       TriangleSize: [3331x1 double]
                np2: 13085
          Vh_P2P2P1: [96360x1 double]
        Vh_P2P2P2P1: [134904x1 double]
              Vh_P2: [38544x1 double]
     meshgeneration: 0
           symmetry: 'N'

NB here the first argument is the mesh handle and the return object is a handle to the new mesh. We could also have used the alternative syntax: [heatS,heatU] = SF_Adapt(heatS,heatU) which, in addition, will project back the two fields on the new mesh and return a handle to them. It is not mandatory to specify the mesh as first argument since the mesh is also a sub-object of the fields.

figure;
subplot(1,2,1); SF_Plot(ffmesh,'mesh','title','Initial mesh');
subplot(1,2,2); SF_Plot(ffmesh_adapted,'mesh','title','Adapted mesh');

recomputes the solution on new mesh and plot the sum of both temperature fields

heatS=SF_Launch('Lshape_Steady.edp','Mesh',ffmesh_adapted,'store','HeatSteady');
heatU=SF_Launch('Lshape_Unsteady.edp','Options',{'omega',omega},'Mesh',ffmesh_adapted,'DataFile','Heat_Unsteady.ff2m','store','HeatUnsteady');
heatSum = SF_Add({heatS,heatU},'coefs',[1 1])
NOTICE  - Storing result in folder HeatSteady
WARNING -  Folder ./WORK//HeatUnsteady mentioned for database storing does not exist ; creating it
NOTICE  - Storing result in folder HeatUnsteady

heatSum = 

  struct with fields:

               mesh: [1x1 struct]
           filename: './WORK/MISC/FFDATA_000001.txt'
    DataDescription: '(no description provided when creating file)'
           datatype: 'Addition'
    datastoragemode: 'ReP2'
    datadescriptors: 'T'
       meshfilename: './WORK/MESHES/FFMESH_000002.msh'
                  T: [13085x1 double]

plot the results

figure();SF_Plot(heatSum,'T','colormap','redblue','title','Temperature (sum of two problems)');
saveas(gca,'FIGURES/Lshape_Tc','png');

Third problem : Loop over a range of omega

The next program Lshape_Impedance.edp solves the unsteady heat conduction problem over a range of $\omega$ and returns the corresponding values of $\hat{\phi}$.

res = SF_Launch('Lshape_Impedance.edp','Mesh',ffmesh_adapted,'Options',{'kappa',1},'DataFile','Heat_Unsteady_loop.ff2m','Store','ThermalImpedances')
WARNING -  Folder ./WORK//ThermalImpedances mentioned for database storing does not exist ; creating it
NOTICE  - Storing result in folder ThermalImpedances

res = 

  struct with fields:

               mesh: [1x1 struct]
           filename: './WORK/ThermalImpedances/FFDATA_000001.txt'
    DataDescription: 'Solution of the unsteady thermal problem over a range of omega'
           datatype: 'EXAMPLE'
    datastoragemode: 'columns'
    datadescriptors: 'omega,Flux_r,Flux_i'
       meshfilename: './WORK/MESHES/FFMESH_000002.msh'
              kappa: 1
           INDEXING: [1x1 struct]
              omega: [50x1 double]
               Flux: [50x1 double]

We plot the 'thermal impedance' defined as $Z_{\theta} = 1/\hat{\phi}$. We compare with the asymptotic solutions in the low-frequency and high-frequency limits.

figure;loglog(res.omega,real(1./res.Flux),'r','LineWidth',2);
hold on;loglog(res.omega,imag(1./res.Flux),'b','LineWidth',2);
rangeomegaLOW = logspace(-1,0,10);
loglog(rangeomegaLOW,2./rangeomegaLOW,'k--','LineWidth',2);
rangeomegaHIGH = logspace(2,3,10);
loglog(rangeomegaHIGH,.25.*rangeomegaHIGH.^(-.5),'k:','LineWidth',2);
xlabel('\omega');ylabel('Re(Z_\theta), Im(Z_\theta)');
title('Thermal impedance Z_\theta(\omega)')
legend('$Re(Z_\theta$)', '$Im(Z_\theta)$','Low-$\omega$ asymptotics','High-$\omega$ asymptotics','Interpreter','latex');
saveas(gca,'FIGURES/Z_omega','png');

Fourth problem : Navier-Stokes flow

Eventually we call a third program Lshape_NavierStokes.edp which computes the flow in the L-shape domain (with bottom-driven boundary conditions). by solving the Navier-Stokes Flow using Newton iteration.

Ndensity =30;
ffmesh = SF_Mesh('Lshape_Mesh.edp','Options',{'meshdensity',Ndensity});
viscosity = 1;Umax = 1;
Flow = SF_Launch('Lshape_NavierStokes.edp','Mesh',ffmesh,'Options',{'Re',10},'Store','NavierStokes')
NOTICE  -       ### INITIAL MESH CREATED WITH np = 826 vertices
WARNING -  Folder ./WORK//NavierStokes mentioned for database storing does not exist ; creating it
NOTICE  - Storing result in folder NavierStokes

Flow = 

  struct with fields:

               mesh: [1x1 struct]
           filename: './WORK/NavierStokes/FFDATA_000001.txt'
    DataDescription: 'Navier-Stokes flow in a L-shaped bottom-driven cavity'
           datatype: 'StokesFlow'
    datastoragemode: 'ReP2P2P1.1'
    datadescriptors: 'ux,uy,p,Re'
       meshfilename: './WORK/MESHES/FFMESH_000003.msh'
                 Re: 10
           INDEXING: [1x1 struct]
                 Fx: -1.7808
          vorticity: [826x1 double]
                 ux: [9180x1 double]
                 uy: [9180x1 double]
                  p: [4590x1 double]

Remarks : - Here one option is specified. The FreeFem program is launched as follows: > FreeFem++ Lshape_NavierStokes.edpc -Re 10

- In this case the data is a vectorial dataset (ux,uy,p) defined using (P2,P2,P1) mixed elements by FreeFem. Note that the SF_Launch driver is able to recognize this and import the fields (ux,uy,p) along with the vorticity.

We adapt the mesh and recompute, this time for Re=20

ffmesh = SF_Adapt(ffmesh,Flow);
Flow = SF_Launch('Lshape_NavierStokes.edp','Mesh',ffmesh,'Options',{'Re',20},'Store','NavierStokes')
NOTICE  - SF_Adapt : created new mesh ; np = 1933 vertices 
NOTICE  - Storing result in folder NavierStokes

Flow = 

  struct with fields:

               mesh: [1x1 struct]
           filename: './WORK/NavierStokes/FFDATA_000002.txt'
    DataDescription: 'Navier-Stokes flow in a L-shaped bottom-driven cavity'
           datatype: 'StokesFlow'
    datastoragemode: 'ReP2P2P1.1'
    datadescriptors: 'ux,uy,p,Re'
       meshfilename: './WORK/MESHES/FFMESH_000004.msh'
                 Re: 20
           INDEXING: [1x1 struct]
                 Fx: -1.3219
          vorticity: [1933x1 double]
                 ux: [21984x1 double]
                 uy: [21984x1 double]
                  p: [10992x1 double]

We plot the vorticity fields with isocontours; We then add the velocity field using quiver plot

figure;SF_Plot(Flow,'vorticity','colorrange',[-20 20]);
hold on;SF_Plot(Flow,{'ux','uy'},'ylim',[0 1],'xlim',[0 1],'title','Navier-Stokes flow (vorticity and velocity fields)');
saveas(gca,'FIGURES/StokesFlow','png');

Appendix : Exploring the database

In addition to providing a "handle" to the meshes and datasets, the drivers SF_Mesh and SF_Launch store all the output files produced by the FreeFem programs in specific subfolders of the "DataBase" folder declared at the beginning. Son everything computed so far has been stored and sorted. We can reaccess the results even after closing Matlab/Octave and reinitializing StabFem :

clear all;
SF_Start('verbosity',3);
SF_DataBase('open','./WORK/');
 Warning  : SF_core_start was not previously lanched ; launching right now
WARNING - librairy ARPACK NOT available
NOTICE  - Verbosity set to 3 : Notice messages
NOTICE  - Database directory :./WORK/
NOTICE  - This directory already exists
NOTICE  - Please type "SF_Status" to see what is available in this dir
NOTICE  -     or type "SF_core_arborescence('cleanall')" to erase any previous data
NOTICE  - Database directory :./WORK/
NOTICE  - This directory already exists
NOTICE  - Please type "SF_Status" to see what is available in this dir
NOTICE  -     or type "SF_core_arborescence('cleanall')" to erase any previous data

To see what is available in the database, just use "SF_Status" to get a summary :

sfs = SF_Status
################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER :    ./WORK/
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/MESHES :
     (list of meshes previously created/adapted ; couples of .msh/.ff2m files )
 
Index | Name              | generation mode | Date                 | Nv      
1     | FFMESH_000001.msh | initial         | 11-juil.-2021 19:59:42 | 1429    
2     | FFMESH_000002.msh | adapted         | 11-juil.-2021 19:59:54 | 3331    
3     | FFMESH_000003.msh | initial         | 11-juil.-2021 20:00:28 | 826     
4     | FFMESH_000004.msh | adapted         | 11-juil.-2021 20:00:34 | 1933    
 
#################################################################################
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/HeatSteady
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file        
 1     | FFDATA_000001.txt | HeatSteady   | 11-juil.-2021 19:59:45 | FFMESH_000001.msh
 2     | FFDATA_000002.txt | HeatUnsteady | 11-juil.-2021 19:59:50 | FFMESH_000001.msh
 3     | FFDATA_000003.txt | HeatSteady   | 11-juil.-2021 19:59:58 | FFMESH_000002.msh
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/HeatUnsteady
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | omega      | HeatFlux        
 1     | FFDATA_000001.txt | HeatUnsteady | 11-juil.-2021 19:59:59 | FFMESH_000002.msh | 100        | 17.2-21.5117i   
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/MISC
     (couples of .txt/.ff2m files )
   NB this directory may contain various secondary files produced by StabFem,
      such as flowfields projected after adaptation but not recomputed, adaptation masks, etc... 
 
 Index | Name              | Type         | Date                 | Mesh file        
 1     | FFDATA_000001.txt | Addition     | 11-juil.-2021 19:59:59 | FFMESH_000002.msh
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/NavierStokes
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Fx        
 1     | FFDATA_000001.txt | StokesFlow   | 11-juil.-2021 20:00:33 | FFMESH_000003.msh | 10         | -1.7808   
 2     | FFDATA_000002.txt | StokesFlow   | 11-juil.-2021 20:00:44 | FFMESH_000004.msh | 20         | -1.3219   
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/ThermalImpedances
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | kappa     
 1     | FFDATA_000001.txt | EXAMPLE      | 11-juil.-2021 20:00:26 | FFMESH_000002.msh | 1         
 
#################################################################################

sfs = 

  struct with fields:

               MESHES: [1x4 struct]
           HeatSteady: [1x3 struct]
         HeatUnsteady: [1x1 struct]
                 MISC: [1x1 struct]
         NavierStokes: [1x2 struct]
    ThermalImpedances: [1x1 struct]

Everything can be reaccessed using SF_Load. For instance Here is how to load and replot the second field stored in folder 'Steady' (as computed at line 178)

heatS = SF_Load('HeatSteady',2);
figure(1);SF_Plot(heatS,'T','xlim',[0 1],'ylim',[0 1],'title', 'Temperature field in a L-shaped domain');
NOTICE  - Loading field number 2 from folder HeatSteady

Notes :

- All Meshes (initial one and adapted ones) are stored in the "MESH" subfolder" - solutions of 1st and 2nd problems are stored in 'ThermalSteady' and 'ThermalUnsteady' folders as requested by the option 'store' when calling SF_Launch - Intermediate calculations (here, the addition of two fields) are stored in a 'MISC' folder - The return object 'sfs' contains summaries of all the folders, returned as an array containing all detected metadata. For instance the objet 'sfs.HeatUnsteady' is an array of dimension 2, and 'omega' and 'HeatFlux' are detected as metadata.

sfs.HeatUnsteady
ans = 

  struct with fields:

        filename: './WORK/HeatUnsteady/FFDATA_000001.txt'
           omega: 100
        HeatFlux: 17.2000 -21.5117i
    meshfilename: './WORK/MESHES/FFMESH_000002.msh'

Final remark :

To understand what the FreeFem++ programs actually do behind the interface, you may have a look to the file ".stabfem_log.bash" generated at the end of this script, which regroups all the calls to FreeFem programs and file manipulations !

% [[PUBLISH]] -> this tag is here to tell the git runner to automatically generate a web page from this script


%
fclose(fopen('SCRIPT_Lshape.success','w+'));