# FeenoX annotated examples

All the files needed to run theses examples are available in the examples directory of the Git repository and any of the tarballs, both source and binary.

For a presentation of these examples, see the “FeenoX Overview Presentation” (August 2021).

# 1 Hello World (and Universe)!

PRINT "Hello $1!" $ feenox hello.fee World
Hello World!
$feenox hello.fee Universe Hello Universe!$

# 2 Lorenz’ attractor (the one with the butterfly)

Solve

\begin{cases} \dot{x} &= \sigma \cdot (y - x) \\ \dot{y} &= x \cdot (r - z) - y \\ \dot{z} &= x y - b z \\ \end{cases}

for 0 < t < 40 with initial conditions

\begin{cases} x(0) = -11 \\ y(0) = -16 \\ z(0) = 22.5 \\ \end{cases}

and \sigma=10, r=28 and b=8/3, which are the classical parameters that generate the butterfly as presented by Edward Lorenz back in his seminal 1963 paper Deterministic non-periodic flow. This example’s input file ressembles the parameters, inital conditions and differential equations of the problem as naturally as possible with an ASCII file.

PHASE_SPACE x y z     # Lorenz attractor’s phase space is x-y-z
end_time = 40         # we go from t=0 to 40 non-dimensional units

sigma = 10            # the original parameters from the 1963 paper
r = 28
b = 8/3

x_0 = -11             # initial conditions
y_0 = -16
z_0 = 22.5

# the dynamical system's equations written as naturally as possible
x_dot = sigma*(y - x)
y_dot = x*(r - z) - y
z_dot = x*y - b*z

PRINT t x y z        # four-column plain-ASCII output
$feenox lorenz.fee > lorenz.dat$ gnuplot lorenz.gp
$python3 lorenz.py$ sh lorenz2x3d.sh < lorenz.dat > lorenz.html  Figure 1: The Lorenz attractor computed with FeenoX plotted with two different tools. a — Gnuplot, b — Matplotlib

# 3 NAFEMS LE10 “Thick plate pressure” benchmark

Assuming the CAD has already been created in STEP format (for instance using Gmsh with this geo file), create a tetrahedral locally-refined unstructured grid with Gmsh using the following .geo file:

// NAFEMS LE10 benchmark unstructured locally-refined tetrahedral mesh

// define physical names from the geometrical entity ids
Physical Surface("upper") = {7};
Physical Surface("DCD'C'") = {1};
Physical Surface("ABA'B'") = {3};
Physical Surface("BCB'C'") = {4, 5};
Physical Curve("midplane") = {14};
Physical Volume("bulk") = {1};

// meshing settings, read Gmsh' manual for further reference
Mesh.ElementOrder = 2;      // use second-order tetrahedra
Mesh.Algorithm = 6;         // 2D mesh algorithm:  6: Frontal Delaunayfe
Mesh.Algorithm3D = 10;      // 3D mesh algorithm: 10: HXT
Mesh.Optimize = 1;          // Optimize the mesh
Mesh.HighOrderOptimize = 1; // Optimize high-order meshes? 2: elastic+optimization

Mesh.MeshSizeMax = 80;     // main element size
Mesh.MeshSizeMin = 20;      // refined element size

// local refinement around the point D (entity 5)
Field = Distance;
Field.NodesList = {4};
Field = Threshold;
Field.IField = 1;
Field.LcMin = Mesh.MeshSizeMin;
Field.LcMax = Mesh.MeshSizeMax;
Field.DistMin = 2 * Mesh.MeshSizeMax;
Field.DistMax = 6 * Mesh.MeshSizeMax;
Background Field = {2};

and then use this pretty-straightforward input file that has a one-to-one correspondence with the original problem formulation from 1990:

# NAFEMS Benchmark LE-10: thick plate pressure
PROBLEM mechanical DIMENSIONS 3
READ_MESH nafems-le10.msh   # mesh in millimeters

BC upper    p=1      # 1 Mpa

# BOUNDARY CONDITIONS:
BC DCD'C'   v=0      # Face DCD'C' zero y-displacement
BC ABA'B'   u=0      # Face ABA'B' zero x-displacement
BC BCB'C'   u=0 v=0  # Face BCB'C' x and y displ. fixed
BC midplane w=0      #  z displacements fixed along mid-plane

# MATERIAL PROPERTIES: isotropic single-material properties
E = 210e3   # Young modulus in MPa
nu = 0.3    # Poisson's ratio

SOLVE_PROBLEM   # solve!

# print the direct stress y at D (and nothing more)
PRINT "sigma_y @ D = " sigmay(2000,0,300) "MPa"

# write post-processing data for paraview
WRITE_MESH nafems-le10.vtk sigmay VECTOR u v w
$gmsh -3 nafems-le10.geo [...]$ feenox nafems-le10.fee
sigma_y @ D =   -5.38016        MPa
$ Normal stress \sigma_y refined around point D over 5,000x-warped displacements for LE10 created with Paraview # 4 How to solve a maze without AI Say you are Homer Simpson and you want to solve a maze drawn in a restaurant’s placemat, one where both the start and end are known beforehand. In order to avoid falling into the alligator’s mouth, you can exploit the ellipticity of the Laplacian operator to solve any maze (even a hand-drawn one) without needing any fancy AI or ML algorithm. Just FeenoX and a bunch of standard open source tools to convert a bitmapped picture of the maze into an unstructured mesh. Bitmapped maze from https://www.mazegenerator.net (left) and 2D mesh (right) PROBLEM laplace 2D # pretty self-descriptive, isn't it? READ_MESH maze.msh # boundary conditions (default is homogeneous Neumann) BC start phi=0 BC end phi=1 SOLVE_PROBLEM # write the norm of gradient as a scalar field # and the gradient as a 2d vector into a .msh file WRITE_MESH maze-solved.msh \ sqrt(dphidx(x,y)^2+dphidy(x,y)^2) \ VECTOR dphidx dphidy 0  $ gmsh -2 maze.geo
[...]
$feenox maze.fee$ Solution to the maze found by FeenoX (and drawn by Gmsh)

# 5 The Fibonacci sequence

## 5.1 Using the closed-form formula as a function

When directly executing FeenoX, one gives a single argument to the executable with the path to the main input file. For example, the following input computes the first twenty numbers of the Fibonacci sequence using the closed-form formula

f(n) = \frac{\varphi^n - (1-\varphi)^n}{\sqrt{5}}

where \varphi=(1+\sqrt{5})/2 is the Golden ratio.

# the fibonacci sequence as function
phi = (1+sqrt(5))/2
f(n) = (phi^n - (1-phi)^n)/sqrt(5)
PRINT_FUNCTION f MIN 1 MAX 20 STEP 1
$feenox fibo_formula.fee | tee one 1 1 2 1 3 2 4 3 5 5 6 8 7 13 8 21 9 34 10 55 11 89 12 144 13 233 14 377 15 610 16 987 17 1597 18 2584 19 4181 20 6765$

## 5.2 Using a vector

We could also have computed these twenty numbers by using the direct definition of the sequence into a vector \vec{f} of size 20.

# the fibonacci sequence as a vector
VECTOR f SIZE 20

f[i]<1:2> = 1
f[i]<3:vecsize(f)> = f[i-2] + f[i-1]

PRINT_VECTOR i f
$feenox fibo_vector.fee > two$

## 5.3 Solving an iterative problem

Finally, we print the sequence as an iterative problem and check that the three outputs are the same.

static_steps = 20
#static_iterations = 1476  # limit of doubles

IF step_static=1|step_static=2
f_n = 1
f_nminus1 = 1
f_nminus2 = 1
ELSE
f_n = f_nminus1 + f_nminus2
f_nminus2 = f_nminus1
f_nminus1 = f_n
ENDIF

PRINT step_static f_n
$feenox fibo_iterative.fee > three$ diff one two
$diff two three$

# 6 Parametric study on a cantilevered beam

If an external loop successively calls FeenoX with extra command-line arguments, a parametric run is obtained. This file cantilever.fee fixes the face called “left” and sets a load in the negative z direction of a mesh called cantilever-$1-$2.msh, where $1 is the first argument after the inpt file and $2 the second one. The output is a single line containing the number of nodes of the mesh and the displacement in the vertical direction w(500,0,0) at the center of the cantilever’s free face.

The following Bash script first calls Gmsh to create the meshes. To do so, it first starts with a base cantilever.geo file that creates the CAD:

// https://autofem.com/examples/determining_natural_frequencie.html

L = 0.5;
b = 0.05;
h = 0.02;

Box(1) = {0,-b/2,-h/2, L, b, h};

Physical Surface("left") = {1};
Physical Surface("right") = {2};
Physical Surface("top") = {4};
Physical Volume("bulk") = {1};

Transfinite Curve {1, 3, 5, 7} = 1/(Mesh.MeshSizeFactor*Mesh.ElementOrder) + 1;
Transfinite Curve {2, 4, 6, 8} = 2/(Mesh.MeshSizeFactor*Mesh.ElementOrder) + 1;
Transfinite Curve {9, 10, 11, 12} = 16/(Mesh.MeshSizeFactor*Mesh.ElementOrder) + 1;

Transfinite Surface "*";
Transfinite Volume "*";

Then another .geo file is merged to build cantilever-${element}-${c}.msh where  Figure 2: Cantilevered beam meshed with structured tetrahedra and hexahedra. a — Tetrahedra, b — Hexahedra

It then calls FeenoX with the input cantilever.fee and passes ${element} and ${c} as extra arguments, which then are expanded as $1 and $2 respectively.

#!/bin/bash

rm -f *.dat
for element in tet4 tet10 hex8 hex20 hex27; do
for c in $(seq 1 10); do # create mesh if not alreay cached mesh=cantilever-${element}-${c} if [ ! -e${mesh}.msh ]; then
scale=$(echo "PRINT 1/${c}" | feenox -)
gmsh -3 -v 0 cantilever-${element}.geo -clscale${scale} -o ${mesh}.msh fi # call FeenoX feenox cantilever.fee${element} ${c} | tee -a cantilever-${element}.dat

done
done

After the execution of the Bash script, thanks to the design decision that output is 100% defined by the user (in this case with the PRINT instruction), one has several files cantilever-${element}.dat files. When plotted, these show the shear locking effect of fully-integrated first-order elements. The theoretical Euler-Bernoulli result is just a reference as, among other things, it does not take into account the effect of the material’s Poisson’s ratio. Note that the abscissa shows the number of nodes, which are proportional to the number of degrees of freedom (i.e. the size of the problem matrix) and not the number of elements, which is irrelevant here and in most problems. PROBLEM elastic 3D READ_MESH cantilever-$1-$2.msh # in meters E = 2.1e11 # Young modulus in Pascals nu = 0.3 # Poisson's ratio BC left fixed BC right tz=-1e5 # traction in Pascals, negative z SOLVE_PROBLEM # z-displacement (components are u,v,w) at the tip vs. number of nodes PRINT nodes %e w(500,0,0) "\#$1 $2" $ ./cantilever.sh
102     -7.641572e-05   # tet4 1
495     -2.047389e-04   # tet4 2
1372    -3.149658e-04   # tet4 3
[...]
19737   -5.916234e-04   # hex27 8
24795   -5.916724e-04   # hex27 9
37191   -5.917163e-04   # hex27 10
$pyxplot cantilever.ppl Displacement at the free tip of a cantilevered beam vs. number of nodes for different element types # 7 Optimizing the length of a tuning fork To illustrate how to use FeenoX in an optimization loop, let us consider the problem of finding the length \ell_1 of a tuning fork (fig. 3) such that the fundamental frequency on a free-free oscillation is equal to the base A frequency at 440 Hz. Figure 3: What length \ell_1 is needed so the fork vibrates at 440 Hz? The FeenoX input is extremely simple input file, since it has to solve the free-free mechanical modal problem (i.e. without any Dirichlet boundary condition). All it has to do is to print the fundamental frequency. To find the length \ell_1, FeenoX is sucessively called from a Python driving script called fork.py. This script uses Gmsh’s Python API to create the CAD and the mesh of the tuning fork given the geometrical arguments r, w, \ell_1 and \ell_2. The parameter n controls the number of elements through the fork’s thickness. Here is the driving script without the CAD & mesh details (the full implementation of the function is available in the examples directory of the FeenoX distribution): import math import gmsh import subprocess # to call FeenoX and read back def create_mesh(r, w, l1, l2, n): gmsh.initialize() ... gmsh.write("fork.msh") gmsh.finalize() return len(nodes) def main(): target = 440 # target frequency eps = 1e-2 # tolerance r = 4.2e-3 # geometric parameters w = 3e-3 l1 = 30e-3 l2 = 60e-3 for n in range(1,7): # mesh refinement level l1 = 60e-3 # restart l1 & error error = 60 while abs(error) > eps: # loop l1 = l1 - 1e-4*error # mesh with Gmsh Python API nodes = create_mesh(r, w, l1, l2, n) # call FeenoX and read scalar back # TODO: FeenoX Python API (like Gmsh) result = subprocess.run(['feenox', 'fork.fee'], stdout=subprocess.PIPE) freq = float(result.stdout.decode('utf-8')) error = target - freq print(nodes, l1, freq) Note that in this particular case, the FeenoX input files does not expand any command-line argument. The trick is that the mesh file fork.msh is overwritten in each call of the optimization loop. The detailed steps between gmsh.initialize() and gmsh.finalize() are not shown here, Since the computed frequency depends both on the length \ell_1 and on the mesh refinement level n, there are actually two nested loops: one parametric over n=1,2\dots,7 and the optimization loop itself that tries to find \ell_1 so as to obtain a frequency equal to 440 Hz within 0.01% of error. PROBLEM modal 3D MODES 1 # only one mode needed READ_MESH fork.msh # in [m] E = 2.07e11 # in [Pa] nu = 0.33 rho = 7829 # in [kg/m^2] # no BCs! It is a free-free vibration problem SOLVE_PROBLEM # write back the fundamental frequency to stdout PRINT f(1) $ python fork.py > fork.dat
$pyxplot fork.ppl$ Figure 4: Estimated length \ell_1 needed to get 440 Hz for different mesh refinement levels n