INTRODUCCIÓN A FreeFEM

Introducción al lenguaje de FreeFEM

FreeFem es un software y un lenguaje de programación orientados a la resolución de ecuaciones en derivadas parciales usando el método de los elementos finitos. FreeFem está escrito en C++ y ha sido desarrollado y mantenido por el Laboratoire Jacques-Louis Lions de la Universidad Pierre et Marie Curie (Paris VI; Sorbonne Université). Corre sobre Linux, Solaris, macOS y Microsoft Windows systems. FreeFem es software libre (LGPL, GNU Lesser General Public License)

Resumen: En estas notas se hace una introducción (no exhaustiva) a los elementos básicos del lenguaje de FreeFEM.

Instrucciones I/O elementales

Con objeto de poder observar los resultados de los ejemplos que se vayan viendo, presentamos en primer lugar instrucciones básicas de lectura y escritura.

	Instrucciones básicas de lectura/escritura
`cout`	Escritura en el dispositivo de salida por defecto (normalmente es la consola)
`cin`	Lectura desde el dispositivo de entrada por defecto (normalmente es el teclado)
`endl`	Fin de línea

Ejemplo 1


// Se escribe en la pantalla el valor de pi, seguido de un salto de linea
cout << pi << endl;

Tipos de datos

Tipos de datos básicos

	Cómo se escriben	Cómo se almacenan	Declaración
Enteros	Números sin punto decimal	Se almacenan en una palabra	`int`
Reales	Con punto decimal, en coma fija o en coma flotante (formato científico)	Norma IEEE 74, doble precisión	`real`
Complejos	Con el símbolo `i` (unidad imaginaria)	Dos reales doble precisión	`complex`
Booleanos	Datos lógicos: `true` y `false`	En un bit	`bool`
String	(Cadenas de caracteres) Entre comillas dobles.	Código ASCII	`string`

Las declaraciones se pueden hacer en cualquier lugar del programa y en ellas se puede inicializar el valor de la variable.

Ejemplos 2


// EJEMPLOS DE IMPRESION EN PANTALLA
// Se escribe un numero en la pantalla
cout << pi << endl;
// Se escribe una concatenacion  texto + numero
cout << "El valor de pi es: " + pi << endl;
// Se deja una linea en blanco
cout << " " << endl;


// EJEMPLOS DE VARIABLES ENTERAS
// Se declara entera la variable var y se lee su valor del teclado
int mivar;
// Se escribe un mensaje indicando que se leera un numero
cout << "Escribe un numero entero" << endl;
cin >> mivar;
// Se imprime el valor de mivar con una explicacion y se deja una linea en blanco
cout << "mivar = " + mivar << endl << endl;


// EJEMPLOS DE VARIABLES DE CADENAS DE CARACTERES
// Se declara string la variable nom
string nom;
// Se escribe un mensaje pidiendo un valor (string)
cout << "Escribe tu nombre: " << endl;
// Se lee un valor (string) y se almacena en nom
cin >> nom;
// Se escribe un texto concatenado con nom y luego un texto en otra linea
cout << "Hola, " + nom << endl << "Todo bien? "<< endl << endl;


// EJEMPLOS DE VARIABLES REALES
// Se declaran real las variables a, b y c
// Se puede asignar un valor en la misma declaracion
real  a = 1.2345, b = -0.321e-1, c;
// Para concatenar una expresion con un texto, el valor de la expresion se pone entre parentesis. Observa la diferencia entre las siguientes instrucciones. Observa, en particular que la instruccion entre parentesis (c=a+b) se ha ejecutado (se ha asignado a c el valor de la suma a+b)
cout << a+b << endl;
cout << " a + b = " + a + b << endl;
cout << " a + b = " +  (a+b) << endl;
cout << " a + b = " +  (c = a+b) << endl;
cout << c << endl << endl;


// EJEMPLO DE VARIABLE COMPLEJA
// Se declara complex la variable cc, se le asigna un valor y se imprime
complex cc;
cc = 2+3i;
cout << cc << endl << endl;


// EJEMPLO DE VARIABLES LOGICAS
// Se declaran bool las variables log1 y log2
// true y false son mostrados como 1 y 0 resp.
// false y 0 son equivalentes
// true y cualquier numero distinto de cero son equivalentes
bool log1 = true, log2 = 3;
cout << log1 << endl << log2 << endl;
cout << (log1 = false) << endl;
cout << (7<6) << endl;

Otros tipos de datos

	Para qué sirven	Declaración
Matrices	Almacenar matrices huecas sparse	`matrix`
Contornos	Definir la frontera del dominio a mallar	`border`
Mallados 2D	Almacenar el mallado de un dominio plano	`mesh`
Mallados 3D	Almacenar el mallado de un dominio tridimensional	`mesh3`
Espacio de elementos finitos	Definir el tipo de elemento finito que se va a utilizar	`fespace`
Funciones	Definir funciones. Para definirlas se pueden usar las variables reservadas `x`, `y`, `z` que representan las coordenadas del punto en que se evaluará la función (solo `x` e `y` en el caso 2D)	`func`
Problema	Definir el problema variacional que se quiere resolver	`problem`

Algunas variables reservadas

Algunos nombres de variables están reservados, por ser variables globales del lenguaje. Se citan aquí los que pueden provocar más confusiones:

Variable reservada (global)	Lo que representa
`N`	Vector normal exterior. `N.x` , `N.y` , `N.z` son sus tres componentes
`P`	Punto actual (cuando se hacen cálculos elementales). `P.x` , `P.y` , `P.z` son sus tres coordenadas
`pi`	$\pi$
`x`, `y`, `z`	Representan las coordenadas de un punto en los cálculos elementales ( solo `x` e `y` en el caso 2D)
`verbosity`	Determina el nivel de verbosidad (RAE: Abundancia de palabras en la elocución) en la información que FreeFEM proporciona durante la ejecución de un programa. Un valor entero entre 0 (ninguna información) y 10 (mucha información).
`CG`, `Cholesky` , `Crout`, `GMRES`, `LU`, `sparsesolver`	Son nombres de solvers de sistemas lineales (algoritmos/rutinas para resover sistemas lineales de ecuaciones)

Formatos de escritura de números reales

En construcción ...

Operaciones entre números

Tipo de las operaciones

Las operaciones aritméticas entre datos de tipos distintos se hacen siempre en el tipo más "alto" de los que intervienen. La jerarquía es la siguiente: booleanos < enteros < reales < complejos.

Ejemplos 3


cout << " 1/2 (entero/entero) = " + (1/2)  << endl; 
cout << " 1./2 (real/entero)    " + (1./2) << endl << endl;

Operadores aritméticos compuestos

Los operadores compuestos += -= *= /= deben entenderse en el sentido siguiente: h *= a es lo mismo que h = h * a. Análogamente para los demás.

Ejemplo


xxxxxxxxxx
int h;
h = 5;
cout << " h = " + h << endl;
cout << " h += 2  (h = h+2) --> " + (h += 2)  + ", h =" + h << endl;
cout << " h -= 3  (h = h-3) --> " + (h += 3)  + ", h =" + h << endl;
cout << " h *= 2  (h = h*2) --> " + (h *= 2)  + ", h =" + h << endl;
cout << " h /= 4  (h = h/4) --> " + (h /= 4)  + ", h =" + h << endl << endl;

Operadores de comparación

Operador	Operación
==	igual a
!=	no igual a
<	menor estricto que
<=	menor o igual que
>	mayor estricto que
>=	mayor o igual que

Resto de la división entera (módulo)

Operador	Operación
`n % m`	$\dfrac{n}{m}$
`a % b`	Si `a``b` $\dfrac{\text{Parte entera}(a)}{ \text{Parte entera}(b)}$

Operadores incremento y decremento

Operador	Operación
`++i`, `--i`	Incrementa (decrementa) en una unidad el valor de `i` (es lo mismo que `i=i+1` o `i=i-1`)
`j = ++i`, `j = --i`	Incrementa (decrementa) en una unidad el valor de `i` y lo almacena en `j`
`j = i++`, `j = i--`	Incrementa (decrementa) en una unidad el valor de `i` y almacena en `j` el valor de `i` antes del incremento (decremento )

Vectores y matrices

En FreeFEM existen dos tipos de vectores:

con índice entero (atención: los índices comienzan en 0)
con índice string (no los tratamos aquí)

Declaración de vectores

Es preciso declarar el tipo y el tamaño de los vectores. Se pueden declarar vectores dobles (con dos subíndices), que no son lo mismo que una matriz (ver Declaración de matrices).

Ejemplos


xxxxxxxxxx
// Se declara un vector real con indice entero, de longitud 5
// A_0, A_1, A_2, A_3, A_4
real[int] A(5);
// Se declara un vector entero con indice entero, de longitud 200
// TRI_0, TRI_1, ........, TRI_199
int[int] TRI(200);
// Se declara un vector doble de números reales con subindices enteros, con 3 filas y 2 columnas
real[int, int]  D(3, 2);

Asignación de valores a un vector unidimensional

Se pueden asignar valores a los vectores (de una dimension) de varias formas:

Escribiendo sus componentes, separadas por comas, entre corchetes

Ejemplo


xxxxxxxxxx
real[int] A(5);
A = [1, 0, 2, 0, 3];
cout << A << endl;

Asignando el mismo valor a todas las componentes

Ejemplo


xxxxxxxxxx
real[int] A(5);
A = 0;
cout << A << endl;

Asignando valor a una componente aislada

Ejemplo


xxxxxxxxxx
// Observa que las componentes aisladas se designan mediante el nombre del vector,
// seguido del indice ENTRE CORCHETES
real[int] A(5);
A = 0;
A[2] = 22;
cout << A << endl;

Asignando valor a un sub-vector (tiene que ser conexo)

Ejemplo


xxxxxxxxxx
real[int] A(5);
A = 0;
A(2:3) = [22, 33];
cout << A << endl;
// Observacion
// para designar índices no consecutivos se usa la sintaxis a : b : paso
//  Ejemplos
//  A(0:5:2) --> [A(0), A(2), A(4)]
A(0:4:2) = 77;
cout << A << endl;

El valor final de A en este ejemplo es:

\begin{matrix} A = (\begin{matrix} 77 & 0 & 77 & 33 & 0 \end{matrix}) \end{matrix}

Asignando valores regularmente espaciados

Ejemplos


xxxxxxxxxx
real[int] A(5);
A = 6:10;
cout << A << endl;


xxxxxxxxxx
real[int] A(5);
A = 10:-1:6;
cout << A << endl;


xxxxxxxxxx
int n = 20;
real[int] A(n);
A = 0:n-1;
cout << A << endl;

Asignando valores en la declaración

Ejemplos


xxxxxxxxxx
real[int] A([5, 10, 15]);
cout << A << endl;


xxxxxxxxxx
real[int] A(1:2:9);
cout << A << endl;

Asignación de valores a un vector doble

Escribiendo sus componentes, por filas, formando un "vector de filas"

Ejemplo


xxxxxxxxxx
real[int, int] D(3, 2);
D = [ [1, 2], [3,4], [5,6] ];
cout << D << endl;

Asignando el mismo valor a todas las componentes

Ejemplo


xxxxxxxxxx
real[int, int] D(3, 2);
D = 1;
cout << D << endl;

Asignando valor a una componente aislada

Ejemplo


xxxxxxxxxx
real[int, int] D(3, 2);
D = 1;
D(0, 1) = 22;
cout << D << endl;
// OBSERVACION
// En el caso de los vectores dobles, las componentes aisladas no se designan con corchetes [ ], sino con parentesis ( )

Asignando valores de otro vector

Ejemplos


xxxxxxxxxx
 real[int, int] D(3, 2);
 D = 1:6;
 cout << D << endl;
// Observa que aqui la asignacion es por columnas


xxxxxxxxxx
real[int, int] D(3, 2);
real[int] A(6);
A = 1:6;
D = A;
cout << D << endl;


xxxxxxxxxx
real[int, int] D(3, 2);
 real[int] A(20);
 A = 1:20;
 D = A(11:16);
 cout << D << endl;

Declaración de matrices

Como se ha indicado antes, hay que distinguir entre vectores dobles y matrices, ya que son objetos con distintas propiedades y métodos. Las matrices son matrices sparse, esto es, solo se almacenan los elementos que no son ceros y sus subíndices (fila y columna). Por ejemplo, la matriz

\begin{matrix} (1) & A = (\begin{matrix} 1 & 2 & 0 \\ 0 & - 1 & 0 \\ - 2 & 0 & 5 \end{matrix}) \end{matrix}

se almacena utilizando dos vectores enteros y uno real: los dos primeros contienen los índices de las filas y columnas de las componentes no nulas, y el tercero contiene los valores de dichas componentes:

\begin{matrix} (2) & I F = (\begin{matrix} 0 \\ 0 \\ 1 \\ 2 \\ 2 \end{matrix}), I C = (\begin{matrix} 0 \\ 1 \\ 1 \\ 0 \\ 2 \end{matrix}), A = (\begin{matrix} 1 \\ 2 \\ - 1 \\ - 2 \\ 5 \end{matrix}) \end{matrix}

$A$ pondremos A(i,j)

Las matrices se declaran con la instrucción matrix y son siempre matrices bidimensionales y, por defecto, reales.

Ejemplos


xxxxxxxxxx
// Se declara una matriz real con tres filas y dos columnas
matrix A(3,2);
// Se declara una matrix compleja con 4 filas y 4 columnas
matrix<complex> C(4,4);

Tamaño dinámico de vectores y matrices

En realidad, en la declaración de vectores y matrices solo hay que indicar el tipo de datos y el número de dimensiones (en el caso de los vectores). El tamaño de ambos es dinámico: se "adapta" a lo que se le asigne.

Ejemplos


xxxxxxxxxx
// Se declara un vector real con índice entero, de longitud 5
real[int] v(5);
// Se asigna 0 a todas sus componentes
v = 0;
cout << v << endl << endl;
// Se asigna a v un vector con 10 componentes
v = 1:10;
cout << v << endl << endl;
// Se asigna a v un vector con 3 componentes
v = [-1, -2, -3];
cout << v << endl << endl;


xxxxxxxxxx
// Se declara un vector doble con 2 filas y dos columnas
real[int,int] D(2,2);
// Se asigna 0 a todas las componentes de D
D = 0;
cout << D << endl;
// Se asigna a D una "matriz" 2x5
D = [ [1,2,3,4,5], [0,1,0,1,0] ];
cout << D << endl;
// Se asigna el valor 7 a todas las componentes de D
// (con su tamaño actual)
D = 7;
cout << D << endl;
// Se asigna a D una "matriz" 1x1
D = [[7]];
cout << D << endl;

Asignación de valores a una matriz

Se pueden asignar valores a los vectores (de una dimension) de varias formas:

Escribiendo sus filas, separadas por comas, entre corchetes

Ejemplo


xxxxxxxxxx
matrix A(3, 2);
A = [ [1, 0], [5, 2], [-2, 0]];
cout << A << endl;

Cuando se imprime una matriz (matrix) se obtiene un resultado como el de la figura siguiente:

sparse_matrix

Asignando un vector como diagonal

Ejemplo


xxxxxxxxxx
real[int] V(1:6);
matrix A(6,6);
A = V;
cout << A << endl;

El resultado de este ejemplo es:

\begin{matrix} (3) & A = (\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 & 0 & 0 \\ 0 & 0 & 0 & 0 & 5 & 0 \\ 0 & 0 & 0 & 0 & 0 & 6 \end{matrix}) \end{matrix}

Asignando un vector doble Ejemplo


xxxxxxxxxx
real[int,int] D(3,2);
matrix A;
D = [ [1, 2], [3,4], [5,6] ];
A = D;
cout << A << endl;

Definiendo la matriz por bloques

Ejemplo


xxxxxxxxxx
matrix M1(2,2), M2(2,2), M3(2,2), M4(2,2);
matrix A;
M1 = [[1,1],[1,1]];
M2 = [[2,2],[2,2]];
M3 = [[3,3],[3,3]];
M4 = [[4,4],[4,4]];
A = [ [M1,M2], [M3,M4] ];
cout << A << endl;

El resultado de este ejemplo es:

\begin{matrix} (4) & A = (\begin{matrix} 1 & 1 & 2 & 2 \\ 1 & 1 & 2 & 2 \\ 3 & 3 & 4 & 4 \\ 3 & 3 & 4 & 4 \end{matrix}) \end{matrix}

Definiéndola como la matriz asociada a una forma bilineal de una formulación variacional (se verá más adelante).

Longitud, máximos y mínimos de vectores y matrices

	Vector `real[int] v(n);`	Vector doble `real[int] v(n,m);`	Matriz `matrix v(n,m);`
`v.n`	Longitud de `v`	Número de filas de `v`	Número de filas de `v`
`v.m`		Número de columnas de `v`	Número de columnas de `v`
`v.nnz`			Número de elementos no nulos de `v`
`v.min`	Mínimo de `v`	Mínimo de `v`
`v.max`	Máximo de `v`	Máximo de `v`
`v'`		Traspuesta de `v`	Traspuesta de `v`
v.sort	Ordena `v` (ver ejemplo)

Ejemplos


xxxxxxxxxx
real[int, int] D(3,2);
D = [[1,2],[3,4],[5,6]];
cout << D' << endl;
// Observacion: la instruccion anterior muestra la traspuesta de D, pero NO CAMBIA D


xxxxxxxxxx
// Estas instrucciones producen un error
matrix D(3,2);
D = [[1,2],[3,4],[5,6]];
cout << D' << endl;

Explicación Si D es una matriz (no un vector doble), no se puede "imprimir" su traspuesta. Se puede asignar a otra matriz e imprimir ésta, como en el ejemplo siguiente


xxxxxxxxxx
matrix D(3,2);
D = [[1,2],[3,4],[5,6]];
D = D';
cout << D << endl;


xxxxxxxxxx
real[int] v(10:-2:2);   // -->v = [10,8,6,4,2]
v.sort;
cout << v << endl;


xxxxxxxxxx
// ATENCION: las siguientes instrucciones no solo muestran el vector v ordenado, 
// sino que también lo dejan ordenado
real[int] v(10:-2:2);
cout << v.sort << endl;
cout << v << endl;

Operaciones con vectores y matrices

En construcción

Pasos en la resolución de un problema con FreeFEM

La resolución de un Problema de EDP's utilizando el software FreeFEM pasa por la realización de las siguientes etapas:

Disponer de la formulación variacional del problema.
Definir el dominio en el que se plantea el problema, a partir de su frontera.
Contruir un mallado del dominio y su frontera (en realidad de una aproximación de ambos).
Definir el espacio de elementos finitos en el que se plantea el problema aproximado.
Definir el problema aproximado.
Resolver el problema aproximado.
Representar gráficamente la solución aproximada.

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla