Exportación e importación de datos

Resumen: En esta sección se muestran algunos de los útiles de FreeFEM para exportar datos a ficheros.

Exportación de mallados

La función savemesh permite salvar a un fichero la información relativa a un mallado:

• savemesh(Th, "fichero.msh") salva toda la información relativa al mallado 2D Th a un fichero de texto de nombre fichero.msh
• savemesh(Th, "fichero.mesh") salva la información relativa al mallado 3D Th a un fichero de texto de nombre fichero.mesh.
• También es posible salvar los datos a ficheros con estructuras adecuadas para ser usados por otros programas. Una descripción más exhaustiva, se puede consultar aquí.

Importación de mallados .msh y .mesh

Para recuperar un mallado de un fichero msh o mesh se puede usar la función readmesh:

• readmesh("fichero.msh") para un mallado 2D.
• readmesh("fichero.mesh") para un mallado 3D.

Observación: los nombres dados a las curvas que definen el contorno no son almacenados en el fichero. Por lo tanto todas las referencias tendrán que hacerse mediante las etiquetas.

Ejemplo

//  construccion del mallado
int[int] lab = [3,5,8,12];
mesh Th = square(20, 20, flags=3, label = lab); 
plot(Th, cmm = "The original mesh", wait=1);
//  salvaguarda
savemesh(Th, "prueba.msh");
// recuperacion
mesh Mh;
Mh = readmesh("prueba.msh");
plot(Mh, cmm = "The read mesh", wait=true);

Escritura y lectura de datos en ficheros

En la Introducción ya se describieron las instrucciones básicas de lectura y escritura de FreeFEM (ver aquí): cin y cout. Estas instrucciones utilizan los dispositivos de entrada y salida por defecto definidos en el sistema (normalmente, son el teclado y la pantalla

Para escribir/leer en un fichero, lo primero que hay que hacer es abrirlo y darle un nombre (interno al programa):

• ofstream file("fichero.sufijo") abre el archivo de nombre fichero.sufijo para escritura y le da el nombre interno file. Si el archivo fichero.sufijo ya existe, lo abre, borra su contenido y se queda listo para escribir en él. Si no existe, se crea.
• ifstream file("fichero.sufijo") abre el archivo fichero.sufijo para lectura y le da el nombre interno file. Si el archivo no existe, se emite un mensaje de error.

Después habrá que dirigir los datos a escribir/leer a ese nombre.

Ejemplo Este código crea/abre un archivo de nombre hola.txt y escribe en él tres líneas:

{
  ofstream mifich("hola.txt");
  mifich << "Hola" << endl << pi << endl << "fin" << endl;
}

Observación: En FreeFEM, el código que se encierra entre llaves forma un bloque . El fichero se cierra de forma automática cuando se cierra el bloque o se termina el programa.

Ejemplo Este código abre el fichero creado en el ejemplo anterior y lee los datos almacenados en él, asignándolos, respectivamente, a las variables var1, num y var2.

{
  ifstream file("hola.txt");
  string var1, var2;
  real num;
  file >> var1 >> num >> var2;
  cout << var1 << endl << var2 << endl << num << endl;
}

Observación: más información y opciones para estas instrucciones aquí.

Ejemplo En este ejemplo se salva a un fichero la solución de un problema variacional calculada previamente:

//
mesh Th = square(20, 20, flags=3); 
fespace Vh(Th, P2);
Vh u, v;
fespace Vh0(Th, P0);
Vh0 f = 60*(abs(x-0.8) < 0.1)*(abs(y-0.8) < 0.1);
solve aDir(u, v) = 
              int2d(Th)(u*v+0.5*(dx(u)*dx(v)+dy(u)*dy(v)))
            - int2d(Th)(f*v)
            + on(1,2, u=1.) + on(3, u=x) + on(4, u=(y-1)^2);
savemesh(Th, "Malla.msh");
{
  ofstream file("Solucion.txt");
  file << u[] << endl;
}

Ejemplo En este ejemplo se recuperan los datos (mallado y solución) guardados en ficheros en el ejemplo anterior y se dibuja la solución.

//
mesh Mh;
Mh = readmesh("Malla.msh");
fespace Vh(Mh, P2);
Vh u;
{
  ifstream file("Solucion.txt");
  file >> u[];
}
plot(Mh, u, cmm = "Solucion", fill = true, wait = true);

Exportación para MATLAB

Los macros ffSaveVh y ffSaveData permiten salvar a un fichero, respectivamente, la información de un espacio de elementos finitos y los datos de una función de dicho espacio. Dichos macros están definidos, junto con otros, en el fichero ffmatlib.idp. Para poder utilizarlos es preciso "incluir" dicho fichero en nuestro código.

Ejemplo

include "ffmatlib.idp"                   // incluir la ruta si es preciso
ffSaveVh(Th, Vh, "export_vh.txt");       // salva el espacio de EF
ffSaveData(u,    "export_data.txt");     // salva la solucion

Para posteriormente utilizar estos datos con MATLAB, es preciso usar una serie de M-funciones de la librería ffmatlib , que se encuentra en la carpeta del mismo nombre. En el programa MATLAB hay que añadir la ruta de dicha carpeta al PATH de MATLAB:

addpath('ffmatlib')                      % incluir la ruta si es preciso

Observación: Tanto el fichero ffmatlib.idp con los macros para FreeFEM como la carpeta ffmatlib con la librería para MATLAB se pueden descargar aquí.

Ejemplo

Añadir el código siguiente al final del programa Elliptic3D.edp:

savemesh(Th3d, "mallado3D.mesh");
//  library for exporting data for MATLAB visualization
include "ffmatlib.idp"
ffSaveVh(Th3d, Vh, "export_vh.txt");
ffSaveData(u,      "export_data.txt");

El siguiente código MATLAB permite visualizar la solución de dicho problema:

xxxxxxxxxx
%%  Visualizacion elemental
addpath('ffmatlib')
[p, b, t] = ffreadmesh('mallado3D.mesh');
[vh] = ffreaddata('export_vh.txt');
[u]  = ffreaddata('export_data.txt');
ffpdeplot3D(p,b,t, 'VhSeq', vh, 'XYZData', u)
xlabel('X')
ylabel('Y')
zlabel('Z')
view(125, 25)
shg

El código MATLAB siguiente permite ver la solución en un corte del objetpo 3D paralelo al plano OXZ. Los puntos S1, S2 y S3 determinan el plano de corte:

xxxxxxxxxx
%%  Corte paralelo al plano OXZ
addpath('ffmatlib')
[p, b, t] = ffreadmesh('mallado3D.mesh');
[vh] = ffreaddata('export_vh.txt');
[u]  = ffreaddata('export_data.txt');
%  Corte paralelo al plano OXZ
y = 0.7;
S1 = [ 0, y, 0];
S2 = [ 0, y, 0.5];
S3 = [ 1, y, 0];
ffpdeplot3D(p,b,t, 'VhSeq', vh, 'XYZData', u, ...
                   'XYZStyle', 'noface', ...
                   'Slice', S1, S2, S3)
xlabel('X')
ylabel('Y')
zlabel('Z')
view(125, 25)
colorbar
shg
%  Otras opciones posibles
%   'Mesh', 'off'
%   'BoundingBox', 'on'
%   Probar con otras coordenadas de S1, S2, S3

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