Tips and tricks

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 04:30, 30 January 2008
Rohtav (Talk | contribs)
(Paraview Snapshots)
← Previous diff
Revision as of 04:32, 30 January 2008
Rohtav (Talk | contribs)
(Tecplot Snapshots)
Next diff →
Line 219: Line 219:
-[[Image:tecplot1.jpg|thumb|400px]]+[[Image:tecplot1.jpg|thumb|300px]]
-[[Image:tecplot2.jpg|thumb|400px]]+[[Image:tecplot2.jpg|thumb|300px]]
- +
=== Paraview Snapshots === === Paraview Snapshots ===

Revision as of 04:32, 30 January 2008

Contents

Emacs mode for Gerris files

Well, not really but something approaching. Add the following to your .emacs

(setq auto-mode-alist (cons '("\\.gfs\\'" . shell-script-mode) auto-mode-alist))

Generating several movies on-the-fly

While it is fairly simple to use the scripting mode of gfsview and unix pipes to generate a movie on the fly from a running simulation, how does one generate several movies simultaneously?

Using named unix fifos and the tee utility it is fairly easy too. For example if one has three gfsview files called wide.gfv, closeup.gfv and overview.gfv and want to generate the three corresponding movies wide.mpg, closeup.mpg and overview.mpg in one go, one could use the following script:

#!/bin/sh

movies="wide closeup overview"
rm -f $movies
mkfifo $movies

gerris3D mysimulation.gfs | tee $movies > /dev/null &
for movie in $movies; do
    gfsview-batch3D $movie.gfv < $movie | ppm2mpeg > $movie.mpg &
done
sleep 10
rm -f $movies

of course the simulation file mysimulation.gfs should contain lines looking like:

OutputSimulation { istep = 10 } stdout
EventScript { istep = 10 } { echo "Save stdout { width = 1024 height = 768 }" }

Compressing simulation files

When it is useful to save simulation results at regular intervals, the size of the files can be reduced by using on-the-fly compression. This can be done like this:

OutputSimulation { istep = 100 } sim-%ld.gfs
EventScript { istep = 100 } { gzip -f -q sim-*.gfs }

GfsView can read compressed GFS files directly.


Writing generic or customized gerris output


Since most of visualization package on the market (or as open-source) do not support quadtree/octree data format, to use benefits of third party visualization packages rather than GFSView we should convert Gerris results to a general unstructured data.

The function "gfs_write_generic_output" is written for this purpose.

This function could write output data in Tecplot or VTK format. We mean generic output because Tecplot format is almost a general unstructured data type and you could easily modify its (code or resulted file) for your own need. Alternatively you could contact me and I could do its for you (just describe your format). All data are written in ASCII format.

VTK format could be visualized by powerful free package: Paraview.

Tecplot format format could be visualized by commersial package Tecplot

Tecplot format is briefly as follows:

1) Header (include list of field variables, number of vertexes, number of elements and elements type)

2) Spatial coordinates plus filed variables related to listed variables in header

3) Element connectivity

For more details regarding to VTK format refer to VTK file format: PDF or HTML


Brief algorithm

1) Construct mesh connectivity 4/8 vertex per cell.

2) Optimize connectivity by removing repeated vertexes, this considerably reduce data size and improve considerably speed of rendering (very important for 3d or large scale simulation).

3) Interpolation of cell-wise data to vertexes.

4) Writing results on a text file.


Function Prototype

Prototype of this function is as follows (6 input arguments):

void gfs_write_generic_output (gchar *type, 
gint n_filed,
gchar **filed_name,
gint plot_depth,
GtsBBox * box,
GfsDomain * domain);


1. type: type of desired output, "VTK" or "TECPLOT" at the moment are supported.

2. n_filed: number of filed data to be visualized, 0 mean just visualization of mesh.

3. filed_name: string array, each item is correspond to each filed, you should include your desired field name, i.e., "U", "V", "W", "P" and desired tracer, e.g., "T", ..., also you could include "Vorticity" to compute and write vorticity field. If for any field, function can not find data pointer from name of variable (except for vorticity), this field is skiped.

4. plot_depth: depth of cell which are desired to visualize, -1 mean visiting all levels.

5. box: if box=NULL, all domain are considered for visualization, but user could define a box to visualize only portion of domain located in the box, to define box you should specify its two corner (by six coordinate, for 2d simulation define two z-value arbitrary).

6. domain: simulation domain.


At the moment this function could be used to write data on-the-fly but it is easy to write a short code to convert "*.sim" data (result of Gerris) to described format. Maybe in future somebody write a tool to convert "*.sim" files to desired output, you could write its yourself, for this purpose follow instruction same as "gfs2oogl.c", to load simulation object from "*.sim" file).


How to Use


To call this function on-the-fly it is sufficient to include this file in "simulation.c" file, i.e. add statement

#include "generic_output.c"

after include statements in the "simulation.c" file of your gerris and add its desired function call in time step loop of "simulation_run" function e.g. at a place like this (determined by "######")

static void simulation_run (GfsSimulation * sim)
{
...
...
...
 
while (sim->time.t < sim->time.end &&
sim->time.i < sim->time.iend) {
 
...
...
...
 
sim->time.t = sim->tnext;
sim->time.i++;
 
gts_range_add_value (&domain->timestep,
gfs_clock_elapsed(domain->timer) - tstart);
gts_range_update (&domain->timestep);
gts_range_add_value (&domain->size,
gfs_domain_size (domain, FTT_TRAVERSE_LEAFS, -1));
gts_range_update (&domain->size);
 
"######"
}

assume we like to run "heated.gfs" example

use the following statements at the specified place (writing output file after each 100 time steps):

{
gint n_field = 5, plot_depth = -1;
gchar *filed_name[5] = {"U", "V", "P", "T", "Vorticity"};
GtsBBox *box= NULL;
 
if(sim->time.i % 100 == 0)
gfs_write_generic_output ("TECPLOT", n_field,
filed_name,plot_depth,
box, domain);
 
// or
 
if(sim->time.i % 100 == 0)
gfs_write_generic_output ("VTK", n_field,
filed_name, plot_depth,
box, domain);
 
}

or if you like to limit visulization within a box use this:

{
gint n_field = 5, plot_depth = -1;
gchar *filed_name[5] = {"U", "V", "P", "T", "Vorticity"};
 
GtsBBox *box= gts_bbox_new (gts_bbox_class (), NULL,
-.5, -.5, -.5,
.5, .5, .5);
 
if(sim->time.i % 100 == 0)
gfs_write_generic_output ("TECPLOT", n_field,
filed_name, plot_depth,
box, domain);
 
// or
 
if(sim->time.i % 100 == 0)
gfs_write_generic_output ("VTK", n_field,
filed_name, plot_depth,
box, domain);
 
}

note that resulted files are names as:

   tecplot%i.plt  or paraview%i.vtk

where %i determine file number, it is automatically incremented after each output

Note that you could write various type or a type with various proberties (e.g. writing both hole data and limited data by a box) simultaneously. For this for this purpose you just include suitable call of function in simulation loop.

After modification of code, you should recompile your gerris code to activate this added feature. Of course it is better to include such data in input "*.gfs", i.e. treat its like an event, but prefarably this job should be performed by Gerris mantainer, if he detect that this feature is useful for a portion of Gerris users.

(in the case of interest, you could announce your interest to this feature in GSF user mailing list).


Tecplot Snapshots


Enlarge


Enlarge

Paraview Snapshots


Enlarge


Enlarge


Enlarge

Source Code Listing

#include <stdlib.h>
#include <math.h>
#include <gmodule.h>
#include "config.h"
#include "gfsconfig.h"
#include "simulation.h"
#include "output.h"
#include "refine.h"
#include "solid.h"
#include "version.h"
 
 
 
#define SQR(a) ((a)*(a))
 
#define my_alloc(pt, n, type, message) {\
if(((pt) = (type *) calloc(n,sizeof(type)))==NULL) {\
printf("\n\n Assertion: There is not sufficient memory to allocate: %s \n\n",message);\
exit(1);\
}\
}
 
 
typedef struct {
gint node[(FTT_DIMENSION - 1)*4];
} Conectivity;
 
 
static gint tecplot_counter = 0;
static gint vtk_counter = 0;
 
void gfs_write_generic_output (gchar *type, gint n_filed, gchar **filed_name, gint plot_depth, GtsBBox * box, GfsDomain * domain);
 
/********************************************************************
 
count number of leaf cells and compute some spatial data
 
********************************************************************/

 
void count_life_cells (FttCell * cell, gpointer * data)
{
gint *ncells;
gdouble *scene_data, h, x_corner, y_corner, z_corner;
FttVector pos;
 
g_return_if_fail (cell != NULL);
g_return_if_fail (data != NULL);
 
ncells = data[0];
scene_data = data[1];
 
h = ftt_cell_size (cell);
 
scene_data[0] = MIN(h, scene_data[0]);
scene_data[1] = MAX(h, scene_data[1]);
 
ftt_cell_pos (cell, &pos);
 
h *= 0.5;
 
x_corner = pos.x - h;
scene_data[3] = MIN(x_corner, scene_data[3]);
scene_data[6] = MAX(x_corner, scene_data[6]);
 
y_corner = pos.y - h;
scene_data[4] = MIN(y_corner, scene_data[4]);
scene_data[7] = MAX(y_corner, scene_data[7]);
 
#if FTT_3D
 
z_corner = pos.z - h;
scene_data[5] = MIN(z_corner, scene_data[5]);
scene_data[8] = MAX(z_corner, scene_data[8]);
 
#endif
 
++ *ncells;
}
 
/********************************************************************
 
fill mesh coordinates
 
********************************************************************/

 
void fill_mesh_coordinate (FttCell * cell, gpointer * data)
{
gchar *type;
gint i, id, *cell_id;
FttVector pos, *coordinate;
gdouble h, *scene_data, *distance;
gdouble scale_y, scale_z, offset_x, offset_y, offset_z;
 
g_return_if_fail (cell != NULL);
g_return_if_fail (data != NULL);
 
cell_id = data[0];
type = data[1];
scene_data = data[2];
coordinate = data[3];
distance = data[4];
 
/* scaling and offseting is performed to have only
same value of distance for exactly same points */

offset_x = scene_data[3];
offset_y = scene_data[4];
offset_z = scene_data[5];
 
scale_y = scene_data[9];
scale_z = scene_data[9] * scene_data[10];
 
id = *cell_id * (FTT_DIMENSION - 1)*4;
 
h = 0.5 * ftt_cell_size (cell);
 
ftt_cell_pos (cell, &pos);
 
#if FTT_2D
 
if(strcmp(type,"TECPLOT") == 0) {
coordinate[id+0].x = pos.x + h; coordinate[id+0].y = pos.y + h;
coordinate[id+1].x = pos.x - h; coordinate[id+1].y = pos.y + h;
coordinate[id+2].x = pos.x - h; coordinate[id+2].y = pos.y - h;
coordinate[id+3].x = pos.x + h; coordinate[id+3].y = pos.y - h;
}
 
if(strcmp(type,"VTK") == 0) {
coordinate[id+0].x = pos.x - h; coordinate[id+0].y = pos.y - h;
coordinate[id+1].x = pos.x + h; coordinate[id+1].y = pos.y - h;
coordinate[id+2].x = pos.x - h; coordinate[id+2].y = pos.y + h;
coordinate[id+3].x = pos.x + h; coordinate[id+3].y = pos.y + h;
}
 
for (i = 0; i < (FTT_DIMENSION - 1)*4; i++)
distance[id+i] = (coordinate[id+i].x - offset_x) +
(coordinate[id+i].y - offset_y) * ( scale_y + 1 ) ;
 
#else /* 3D */
 
if(strcmp(type,"TECPLOT") == 0) {
coordinate[id+0].x = pos.x + h; coordinate[id+0].y = pos.y + h; coordinate[id+0].z = pos.z + h;
coordinate[id+1].x = pos.x - h; coordinate[id+1].y = pos.y + h; coordinate[id+1].z = pos.z + h;
coordinate[id+2].x = pos.x - h; coordinate[id+2].y = pos.y - h; coordinate[id+2].z = pos.z + h;
coordinate[id+3].x = pos.x + h; coordinate[id+3].y = pos.y - h; coordinate[id+3].z = pos.z + h;
coordinate[id+4].x = pos.x + h; coordinate[id+4].y = pos.y + h; coordinate[id+4].z = pos.z - h;
coordinate[id+5].x = pos.x - h; coordinate[id+5].y = pos.y + h; coordinate[id+5].z = pos.z - h;
coordinate[id+6].x = pos.x - h; coordinate[id+6].y = pos.y - h; coordinate[id+6].z = pos.z - h;
coordinate[id+7].x = pos.x + h; coordinate[id+7].y = pos.y - h; coordinate[id+7].z = pos.z - h;
}
 
if(strcmp(type,"VTK") == 0) {
coordinate[id+0].x = pos.x - h; coordinate[id+0].y = pos.y - h; coordinate[id+0].z = pos.z - h;
coordinate[id+1].x = pos.x + h; coordinate[id+1].y = pos.y - h; coordinate[id+1].z = pos.z - h;
coordinate[id+2].x = pos.x - h; coordinate[id+2].y = pos.y + h; coordinate[id+2].z = pos.z - h;
coordinate[id+3].x = pos.x + h; coordinate[id+3].y = pos.y + h; coordinate[id+3].z = pos.z - h;
coordinate[id+4].x = pos.x - h; coordinate[id+4].y = pos.y - h; coordinate[id+4].z = pos.z + h;
coordinate[id+5].x = pos.x + h; coordinate[id+5].y = pos.y - h; coordinate[id+5].z = pos.z + h;
coordinate[id+6].x = pos.x - h; coordinate[id+6].y = pos.y + h; coordinate[id+6].z = pos.z + h;
coordinate[id+7].x = pos.x + h; coordinate[id+7].y = pos.y + h; coordinate[id+7].z = pos.z + h;
}
 
for (i = 0; i < (FTT_DIMENSION - 1)*4; i++)
distance[id+i] = (coordinate[id+i].x - offset_x) +
(coordinate[id+i].y - offset_y) * ( scale_y + 1 ) +
(coordinate[id+i].z - offset_z) * ( scale_z + 1 ) ;
 
#endif /* 3D */
 
++ *cell_id;
}
 
/********************************************************************
 
compute cell vorticity
 
********************************************************************/

 
void compute_cell_vorticity (FttCell * cell, gpointer * data) {
 
GfsVariable **velocity, *vorticity;
 
g_return_if_fail (cell != NULL);
g_return_if_fail (data != NULL);
 
velocity = data[0];
vorticity = data[1];
 
GFS_VARIABLE (cell, vorticity->i) = gfs_vorticity (cell, velocity);
 
}
 
/********************************************************************
 
interpolate cell data to vertexes and compute interpolation weight
 
********************************************************************/

 
void cell_to_vertex (FttCell * cell, gpointer * data) {
 
static FttDirection d[(FTT_DIMENSION - 1)*4][FTT_DIMENSION] = {
#if FTT_2D
{FTT_RIGHT, FTT_TOP}, {FTT_LEFT, FTT_TOP}, {FTT_LEFT, FTT_BOTTOM}, {FTT_RIGHT, FTT_BOTTOM}
#else /* 3D */
{FTT_RIGHT, FTT_TOP, FTT_FRONT}, {FTT_LEFT, FTT_TOP, FTT_FRONT},
{FTT_LEFT, FTT_BOTTOM, FTT_FRONT}, {FTT_RIGHT, FTT_BOTTOM, FTT_FRONT},
{FTT_RIGHT, FTT_TOP, FTT_BACK}, {FTT_LEFT, FTT_TOP, FTT_BACK},
{FTT_LEFT, FTT_BOTTOM, FTT_BACK}, {FTT_RIGHT, FTT_BOTTOM, FTT_BACK},
#endif /* 3D */
};
 
gint i, j, k, n, *cell_id, *n_real_filed;
gdouble **field_data;
gint id_vertex_cell[(FTT_DIMENSION - 1)*4];
 
Conectivity *conectivity;
 
GfsVariable *variable;
 
g_return_if_fail (cell != NULL);
g_return_if_fail (data != NULL);
 
cell_id = data[0];
n_real_filed = data[1];
conectivity = data[2];
field_data = data[4];
k = 5;
 
n = *n_real_filed;
 
for (i = 0; i < (FTT_DIMENSION - 1)*4; i++)
id_vertex_cell[i] = conectivity[*cell_id].node[i];
 
for (j = 0; j < n; j++) {
variable = data[j+k];
for (i = 0; i < (FTT_DIMENSION - 1)*4; i++)
field_data[j][id_vertex_cell[i]] = gfs_cell_corner_value (cell, d[i], variable, -1);
}
 
++ *cell_id;
}
 
/****************************************************
 
quick sort of double array "a" with length "n"
and return integer permutation array "t"
 
*****************************************************/

 
void my_fqsort (gint n, gdouble *a, gint *t) {
 
static gint i, j, k, l, r, s, stackl[100], stackr[100], ww, flg1, flg2;
static gdouble w, x;
 
--a; --t;
s = 1; stackl[1] = 1; stackr[1] = n;
 
L10: l = stackl[s]; r = stackr[s--];
 
L20: i = l; j = r; k = ( l + r ) / 2; x = a[k];
 
/* repeat until i > j */
flg1=1;
while (flg1) {
flg2=1;
while(flg2)
if (a[i]<x) ++i; /* search from lower end */
else flg2=0;
flg2=1;
while (flg2)
if (x<a[j]) --j; /* search from upper end */
else flg2=0;
 
if (i<=j) { // swap positions i, j
w = a[i]; ww = t[i]; a[i] = a[j]; t[i++] = t[j]; a[j] = w; t[j--] = ww;
if (i>j) flg1=0;
} else flg1=0;
};
 
if ( j - l >= r - i ) {
if (l < j ) {
stackl[++s] = l; stackr[s] = j;
}
l = i;
} else {
if ( i < r ) {
stackl[++s]= i; stackr[s] = r;
}
r = j;
}
if (l < r) goto L20;
if ( s != 0 ) goto L10;
return;
}
 
/****************************************************
 
quick sort of integer array "a" with length "n"
and return integer permutation array "t"
 
****************************************************/

 
void my_iqsort (gint n, gint *a, gint *t) {
 
static gint i, j, k, l, r, s, stackl[100], stackr[100], ww, flg1, flg2;
static gint w, x;
 
--a; --t;
s = 1; stackl[1] = 1; stackr[1] = n;
 
L10: l = stackl[s]; r = stackr[s--];
 
L20: i = l; j = r; k = ( l + r ) / 2; x = a[k];
 
/* repeat until i > j */
flg1=1;
while (flg1) {
flg2=1;
while(flg2)
if (a[i]<x) ++i; /* search from lower end */
else flg2=0;
flg2=1;
while (flg2)
if (x<a[j]) --j; /* search from upper end */
else flg2=0;
 
if (i<=j) { // swap positions i, j
w = a[i]; ww = t[i]; a[i] = a[j]; t[i++] = t[j]; a[j] = w; t[j--] = ww;
if (i>j) flg1=0;
} else flg1=0;
};
 
if ( j - l >= r - i ) {
if (l < j ) {
stackl[++s] = l; stackr[s] = j;
}
l = i;
} else {
if ( i < r ) {
stackl[++s]= i; stackr[s] = r;
}
r = j;
}
if (l < r) goto L20;
if ( s != 0 ) goto L10;
return;
}
 
/********************************************************************
 
remove repeated vertex from conectivity and re-new conectivity
 
********************************************************************/

 
void remove_repeated_vertex ( gint n_cells,
gint *n_vertex,
gdouble *scene_data,
gdouble *distance,
FttVector *coordinate,
Conectivity *conectivity
) {
 
gint i, j, n, n_vertexes_new, nv, *perm, *iperm, *label, *id_vertex;
gdouble small = 0.25 * scene_data[0];
gchar *flg;
FttVector *new_coordinate;
 
nv = 4*(FTT_DIMENSION - 1);
n = nv*n_cells;
 
/* allocate array for permutation array and label of repeated vertexes */
 
my_alloc(perm, n, gint, "temporary permutation array");
my_alloc(iperm, n, gint, "temporary ipermutation array");
my_alloc(label, n, gint, "temporary label array");
my_alloc(flg, n, gchar, "temporary flag array");
 
for (i=0; i<n; i++) {
perm [i] = label [i]= i;
flg [i] = 1;
}
 
/*
quick sort algorithm, note that qsort of ANCI C is not suitable
for our purpose, since we need permutation array rather than sorted list
*/

 
my_fqsort (n, distance, perm);
 
n_vertexes_new = 1;
 
for (i=0; i<(n-1); i++) {
 
if( fabs( distance[i] - distance[i+1] ) < small ) {
 
gint iv1 = perm[i];
gint iv2 = perm[i+1];
 
if (
fabs( coordinate[iv1].x -
coordinate[iv2].x ) < small
&&
fabs( coordinate[iv1].y -
coordinate[iv2].y ) < small
#if FTT_3D
&&
fabs( coordinate[iv1].z -
coordinate[iv2].z ) < small
#endif
) {
label[iv2] = label[iv1];
flg[iv2] = 0;
} else ++n_vertexes_new;
} else ++n_vertexes_new;
}
 
my_alloc(id_vertex, n_vertexes_new, gint, "temporary id vertex array");
 
n_vertexes_new = 0;
 
for (i=0; i<n; i++) if( flg [i] > 0 ) {
perm [n_vertexes_new] = i;
id_vertex [n_vertexes_new++] = label[i];
}
 
my_iqsort (n_vertexes_new, id_vertex, perm);
 
/* re-new conectivity */
 
for (i=0; i<n_vertexes_new; i++) iperm [perm[i]] = i;
 
for (i=0; i<n_cells; i++)
for (j=0; j<nv; j++)
conectivity[i].node[j] = iperm [label[i*nv + j]];
 
g_free(label);
g_free(perm);
g_free(iperm);
g_free(id_vertex);
 
my_alloc(new_coordinate, n_vertexes_new, FttVector, "temporary array for new_coordinate");
 
j=0;
for(i=0; i<n; i++) if(flg[i]>0) {
new_coordinate[j].x = coordinate[i].x;
new_coordinate[j].y = coordinate[i].y;
#if FTT_3D
new_coordinate[j].z = coordinate[i].z;
#endif
++j;
}
 
*n_vertex = n_vertexes_new;
 
for(i=0; i<n_vertexes_new; i++) {
coordinate[i].x = new_coordinate[i].x;
coordinate[i].y = new_coordinate[i].y;
#if FTT_3D
coordinate[i].z = new_coordinate[i].z;
#endif
}
 
g_free(flg);
g_free(new_coordinate);
}
 
/********************************************************************
 
write a tecplot file
 
*********************************************************************/

 
void gfs_write_generic_output (gchar *type, gint n_filed, gchar **filed_name, gint plot_depth, GtsBBox * box, GfsDomain * domain)
{
gchar file_name[100], temp_string[60];
FILE * fp;
gint n, n_cells, n_vertex, cell_id, i, j, k, n_real_filed, id_real_filed[20], vort, have_vel_vector[4];
gdouble *distance, scene_data[20], *field_data[20];
 
GfsVariable *vorticity;
 
Conectivity *conectivity;
FttVector *coordinate;
 
gpointer data[20];
 
/* count number of leaf cells and some useful spatial data */
 
n_cells = 0;
 
scene_data[0] = G_MAXDOUBLE; /* h_min */
scene_data[1] = 0.0; /* h_max */
 
scene_data[3] = G_MAXDOUBLE; /* x_min */
scene_data[4] = G_MAXDOUBLE; /* y_min */
scene_data[5] = G_MAXDOUBLE; /* z_min */
 
scene_data[6] = G_MINDOUBLE; /* x_max */
scene_data[7] = G_MINDOUBLE; /* y_max */
scene_data[8] = G_MINDOUBLE; /* z_max */
 
data[0] = &n_cells;
data[1] = scene_data;
 
if(box != NULL)
gfs_domain_cell_traverse_box (domain, box, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, plot_depth,
(FttCellTraverseFunc) count_life_cells, data);
else
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, plot_depth,
(FttCellTraverseFunc) count_life_cells, data);
 
scene_data[9] = scene_data[6] - scene_data[3]; /* lenght_x */
scene_data[10] = scene_data[7] - scene_data[4]; /* lenght_y */
scene_data[11] = scene_data[8] - scene_data[5]; /* lenght_z */
 
 
/* allocate temporary memory to save mesh coordinates and its conectivity */
 
n = n_cells * 4*(FTT_DIMENSION - 1);
 
my_alloc(conectivity, n_cells, Conectivity, "temporary mesh conectivity");
my_alloc(coordinate, n, FttVector, "temporary mesh coordinates");
my_alloc(distance, n, gdouble, "temporary distance of mesh nodes");
 
/* fill coordinates of leaf cells and their distance from a reference point */
/* note that in this stage connectivity is priori known */
 
cell_id = 0;
 
data[0] = &cell_id;
data[1] = type;
data[2] = scene_data;
data[3] = coordinate;
data[4] = distance;
 
if(box != NULL)
gfs_domain_cell_traverse_box (domain, box, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, plot_depth,
(FttCellTraverseFunc) fill_mesh_coordinate, data);
else
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, plot_depth,
(FttCellTraverseFunc) fill_mesh_coordinate, data);
 
/*
removing repeated nodes from global connectivity
for this: we sort nodes based on the distance from
a refrence ponit and then detect repeated ones, remove
and renew connectivity. Note that to be compatible with formal
post-processors, like Tecplot, VTK, Ensight, we have to convert
cell-based data to vertex-based data and use weighted interpolation
of cell-wise data to compute data at vertexes. Also precence of
repeated vertexes decrease rendering speed of post-processor
*/

 
remove_repeated_vertex ( n_cells, &n_vertex, scene_data, distance, coordinate, conectivity);
 
/* allocate memory to store vertex interpolation of cell-wise data */
k = 5;
vort = -1;
have_vel_vector[0] = 0;
n_real_filed = 0;
for(i=0; i<n_filed; i++) {
data[k+n_real_filed] = gfs_variable_from_name (domain->variables, filed_name[i]);
if(data[k+n_real_filed] == NULL) {
if(strcmp( filed_name[i], "Vorticity") == 0 ) vort = i;
}
else {
if(strcmp( filed_name[i], "U") == 0 ) {++have_vel_vector[0]; have_vel_vector[1] = i;}
if(strcmp( filed_name[i], "V") == 0 ) {++have_vel_vector[0]; have_vel_vector[2] = i;}
if(strcmp( filed_name[i], "W") == 0 ) {++have_vel_vector[0]; have_vel_vector[3] = i;}
 
my_alloc(field_data[n_real_filed], n_vertex, gdouble, "temporary array to store vertex interpolated data");
id_real_filed[n_real_filed++] = i;
}
}
 
if(vort>0) {
 
vorticity = gfs_temporary_variable (domain);
 
data[0] = gfs_domain_velocity (domain);
data[1] = vorticity;
 
if(box != NULL)
gfs_domain_cell_traverse_box (domain, box, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, (plot_depth == -1 ? (-1):(plot_depth+1)),
(FttCellTraverseFunc) compute_cell_vorticity, data);
else
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, (plot_depth == -1 ? (-1):(plot_depth+1)),
(FttCellTraverseFunc) compute_cell_vorticity, data);
 
data[k+n_real_filed] = vorticity;
 
my_alloc(field_data[n_real_filed], n_vertex, gdouble, "temporary array to store vertex interpolated data");
 
id_real_filed[n_real_filed++] = vort;
}
 
if(n_real_filed>0) {
for(j=0; j<n_real_filed; j++)
for(i=0; i<n_vertex; i++)
field_data[j][i] = 0.;
cell_id = 0;
 
data[0] = &cell_id;
data[1] = &n_real_filed;
data[2] = conectivity;
data[4] = field_data;
 
if(box != NULL)
gfs_domain_cell_traverse_box (domain, box, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, plot_depth,
(FttCellTraverseFunc) cell_to_vertex, data);
else
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, plot_depth,
(FttCellTraverseFunc) cell_to_vertex, data);
 
}
 
 
/* write data on a file */
 
/* ---------------------------------------------------
Tecplot output, it could be used as a general output
you could easily modify sort of data based on your
specific data format
------------------------------------------------- */

 
if(strcmp(type,"TECPLOT") == 0) {
 
int nn = (FTT_DIMENSION - 1)*4;
 
/* name of output file */
 
++ tecplot_counter;
 
strcpy(file_name,"tecplot");
sprintf(temp_string,"%i", tecplot_counter);
strcat(temp_string, ".plt");
strcat(file_name,temp_string);
 
fp = fopen (file_name, "w");
 
/* write header of vtk file */
 
fprintf (fp, " TITLE = 'Gerris Flow Solver Output' \n");
 
#if FTT_2D
fprintf (fp, " VARIABLES = 'X', 'Y'");
#else /* 3D */
fprintf (fp, " VARIABLES = 'X', 'Y', 'Z'");
#endif
 
for(j=0; j<n_real_filed; j++)
fprintf (fp, ", '%s' ",filed_name[id_real_filed[j]]);
fprintf (fp, " \n");
 
#if FTT_2D
fprintf (fp, " ZONE N=%i, E=%i, F=FEPOINT, ET=QUADRILATERAL \n", n_vertex, n_cells);
#else /* 3D */
fprintf (fp, " ZONE N=%i, E=%i, F=FEPOINT, ET=BRICK \n", n_vertex, n_cells);
#endif
 
 
/* write vertex coordinates and related field data */
 
for(i=0; i<n_vertex; i++) {
 
#if FTT_2D
fprintf (fp, " %f %f ", coordinate[i].x, coordinate[i].y);
#else /* 3D */
fprintf (fp, " %f %f %f ", coordinate[i].x, coordinate[i].y, coordinate[i].z);
#endif
for(j=0; j<n_real_filed; j++)
fprintf (fp, " %f ",field_data[j][i]);
fprintf (fp, " \n");
}
 
/* write conectivity data */
 
for(i=0; i<n_cells; i++) {
for (j = 0; j < nn; j++)
fprintf (fp, " %i ", conectivity[i].node[j]+1);
fprintf (fp, " \n");
}
fclose(fp);
fprintf (stderr, "\n %s file %i is dumped n_cells = %i \n", type, tecplot_counter, n_cells);
}
 
 
/* ---------------------------------------------------
VTK output, it could be imported and visualized by
free visualization package ParaView, due to
consistancy with various platform ASCII format is used
--------------------------------------------------- */

 
if(strcmp(type,"VTK") == 0) {
 
int nn = (FTT_DIMENSION - 1)*4;
 
/* name of output file */
 
++ vtk_counter;
 
strcpy(file_name,"paraview");
sprintf(temp_string,"%i", vtk_counter);
strcat(temp_string, ".vtk");
strcat(file_name,temp_string);
 
fp = fopen (file_name, "w");
 
 
/* write header of vtk file */
 
fprintf (fp, "# vtk DataFile Version 2.0\n");
fprintf (fp, "Unstructured Grid Example\n");
fprintf (fp, "ASCII\n");
fprintf (fp, "DATASET UNSTRUCTURED_GRID\n");
fprintf (fp, "\n");
 
/* write spatial coordinates */
 
fprintf (fp, "POINTS %i float\n", n_vertex);
 
#if FTT_2D
for(i=0; i<n_vertex; i++)
fprintf (fp, "%f %f 0 \n", coordinate[i].x, coordinate[i].y);
#else /* 3D */
for(i=0; i<n_vertex; i++)
fprintf (fp, "%f %f %f\n", coordinate[i].x, coordinate[i].y, coordinate[i].z);
#endif
fprintf (fp, "\n");
 
/* write conectivity data */
 
fprintf (fp, "CELLS %i %i \n", n_cells, n_cells*(nn+1));
 
for(i=0; i<n_cells; i++) {
fprintf (fp, " %i ",nn);
for (j = 0; j < nn; j++)
fprintf (fp, " %i ", conectivity[i].node[j]);
fprintf (fp, " \n");
}
fprintf (fp, "\n");
 
fprintf (fp, "CELL_TYPES %i \n",n_cells);
 
for(i=0; i<n_cells; i++)
#if FTT_2D
fprintf (fp, "8 \n");
#else
fprintf (fp, "11 \n");
#endif
fprintf (fp, "\n");
 
/* write scaler fields */
 
if(n_real_filed>0) {
 
fprintf (fp, "POINT_DATA %i\n",n_vertex);
fprintf (fp, "\n");
 
for(j=0; j<n_real_filed; j++) {
 
fprintf (fp, "SCALARS %s float\n",filed_name[id_real_filed[j]]);
fprintf (fp, "LOOKUP_TABLE default \n");
 
for(i=0; i<n_vertex; i++)
fprintf(fp, "%f\n",field_data[j][i]);
fprintf (fp, "\n");
}
}
 
/* write velosity fields as a vector field,
if user include all velosity components */

 
if(have_vel_vector[0] == FTT_DIMENSION) {
 
fprintf (fp, "VECTORS Velocity_vector float\n");
 
#if FTT_2D
for(i=0; i<n_vertex; i++)
fprintf (fp, "%f %f 0\n",field_data[have_vel_vector[1]][i], field_data[have_vel_vector[2]][i]);
#else
for(i=0; i<n_vertex; i++)
fprintf (fp, "%f %f %f\n",field_data[have_vel_vector[1]][i], field_data[have_vel_vector[2]][i], field_data[have_vel_vector[3]][i]);
#endif
}
fclose(fp);
fprintf (stderr, "\n %s file %i is dumped n_cells = %i \n", type, vtk_counter, n_cells);
}
 
g_free(conectivity);
g_free(coordinate);
g_free(distance);
for(j=0; j<n_real_filed; j++) g_free(field_data[j]);
}
Personal tools
communication