Learn by examples 1: Membrane
Contents
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.
Initialization
close all; addpath('../../SOURCES_MATLAB'); SF_Start('verbosity',4); % verbosity=4 will display output from the FreeFem program SF_DataBase('disable');
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: [0x1 double] status: 'loaded' 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'); pause(0.1);
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] status: 'loaded'
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' figure(3);plot(Xcut,real(Cut.phi),'r-'); title('Vertical Displacement along a line : Phi(x,0.5)');xlabel('x');ylabel('Phi'); pause(0.1);
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Exercices :
- Modify the example so that the membrane is circular (a=1,b=1).
(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 !
% [[PUBLISH]] % fclose(fopen('SCRIPT_membrane.success','w+'));