diff --git a/examples/incompressible-navier-stokes/tnlNSFastBuildConfig.h b/examples/incompressible-navier-stokes/tnlNSFastBuildConfig.h new file mode 100644 index 0000000000000000000000000000000000000000..fc03ebd1ff9a5f274f89e92d99cc1b72c7c29df7 --- /dev/null +++ b/examples/incompressible-navier-stokes/tnlNSFastBuildConfig.h @@ -0,0 +1,71 @@ +/*************************************************************************** + tnlNSFastBuildConfig.h - description + ------------------- + begin : Jul 7, 2014 + copyright : (C) 2014 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#ifndef TNLNSFASTBUILDCONFIG_H_ +#define TNLNSFASTBUILDCONFIG_H_ + +#include <solvers/tnlConfigTags.h> + +class tnlNSFastBuildConfig +{ + public: + + static void print() { cerr << "tnlNSFastBuildConfig" << endl; } +}; + +/**** + * Turn off support for float and long double. + */ +template<> struct tnlConfigTagReal< tnlNSFastBuildConfig, float > { enum { enabled = false }; }; +template<> struct tnlConfigTagReal< tnlNSFastBuildConfig, long double > { enum { enabled = false }; }; + +/**** + * Turn off support for short int and long int indexing. + */ +template<> struct tnlConfigTagIndex< tnlNSFastBuildConfig, short int >{ enum { enabled = false }; }; +template<> struct tnlConfigTagIndex< tnlNSFastBuildConfig, long int >{ enum { enabled = false }; }; + +/**** + * 1, 2, and 3 dimensions are enabled by default + */ +template<> struct tnlConfigTagDimensions< tnlNSFastBuildConfig, 1 >{ enum { enabled = false }; }; +template<> struct tnlConfigTagDimensions< tnlNSFastBuildConfig, 2 >{ enum { enabled = true }; }; +template<> struct tnlConfigTagDimensions< tnlNSFastBuildConfig, 3 >{ enum { enabled = false }; }; + +/**** + * Use of tnlGrid is enabled for allowed dimensions and Real, Device and Index types. + */ +template< int Dimensions, typename Real, typename Device, typename Index > + struct tnlConfigTagMesh< tnlNSFastBuildConfig, tnlGrid< Dimensions, Real, Device, Index > > + { enum { enabled = tnlConfigTagDimensions< tnlNSFastBuildConfig, Dimensions >::enabled && + tnlConfigTagReal< tnlNSFastBuildConfig, Real >::enabled && + tnlConfigTagDevice< tnlNSFastBuildConfig, Device >::enabled && + tnlConfigTagIndex< tnlNSFastBuildConfig, Index >::enabled }; }; + +/**** + * Please, chose your preferred time discretisation here. + */ +template<> struct tnlConfigTagTimeDiscretisation< tnlNSFastBuildConfig, tnlExplicitTimeDiscretisationTag >{ enum { enabled = true }; }; +template<> struct tnlConfigTagTimeDiscretisation< tnlNSFastBuildConfig, tnlSemiImplicitTimeDiscretisationTag >{ enum { enabled = true }; }; +template<> struct tnlConfigTagTimeDiscretisation< tnlNSFastBuildConfig, tnlImplicitTimeDiscretisationTag >{ enum { enabled = false }; }; + +/**** + * Only the Runge-Kutta-Merson solver is enabled by default. + */ +template<> struct tnlConfigTagExplicitSolver< tnlNSFastBuildConfig, tnlExplicitEulerSolverTag >{ enum { enabled = false }; }; + +#endif /* TNLNSFASTBUILDCONFIG_H_ */ diff --git a/examples/incompressible-navier-stokes/tnlPoissonProblem.h b/examples/incompressible-navier-stokes/tnlPoissonProblem.h new file mode 100644 index 0000000000000000000000000000000000000000..74d2ff515eb46f0f07772c9509cfd4134835ee5d --- /dev/null +++ b/examples/incompressible-navier-stokes/tnlPoissonProblem.h @@ -0,0 +1,110 @@ +#ifndef TNLPOISSONPROBLEM_H +#define TNLPOISSONPROBLEM_H + +#include <core/vectors/tnlVector.h> +#include <mesh/tnlGrid.h> + +template< typename Mesh, + typename Real = typename Mesh::RealType, + typename Index = typename Mesh::IndexType > +class tnlPoissonProblem +{ + +}; + + +template< typename MeshReal, + typename Device, + typename MeshIndex, + typename Real, + typename Index, + typename MatrixType> +class tnlPoissonProblem< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index > +{ + public: + + typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType; + typedef typename MeshType::CoordinatesType CoordinatesType; + typedef Real RealType; + typedef Device DeviceType; + typedef Index IndexType; + + static tnlString getType() + { + return tnlString( "tnlPoisssonMatrix< " ) + + MeshType::getType() + ", " + + ::getType< Real >() + ", " + + ::getType< Index >() + " >"; + } + + template< typename Vector > +#ifdef HAVE_CUDA + __device__ __host__ +#endif + Real getValue( const MeshType& mesh, + const IndexType cellIndex, + const CoordinatesType& coordinates, + const Vector& u, + const Real& time ) const; + +#ifdef HAVE_CUDA + __device__ __host__ +#endif + Index getLinearSystemRowLength( const MeshType& mesh, + const IndexType& index, + const CoordinatesType& coordinates ) const + { + return 5; + } + +#ifdef HAVE_CUDA + __device__ __host__ +#endif + void initLinearSystem( + const MeshType& mesh, + const IndexType& index, + ) + { + _matrix.setElement(index,index,4); + _matrix.setElement(index,mesh.getCellYPredecessor( index ),-1); + _matrix.setElement(index,mesh.getCellXPredecessor( index ),-1); + _matrix.setElement(index,mesh.getCellXSuccessor( index ),-1); + _matrix.setElement(index,mesh.getCellYSuccessor( index ),-1); + + /*matrixRow.setElement( 0, mesh.getCellYPredecessor( index ), -1 ); + matrixRow.setElement( 1, mesh.getCellXPredecessor( index ), -1 ); + matrixRow.setElement( 2, index, 4.0 ); + matrixRow.setElement( 3, mesh.getCellXSuccessor( index ), -1 ); + matrixRow.setElement( 4, mesh.getCellYSuccessor( index ), -1 );*/ + } + + void init(const MeshType& mesh) + { + typename MatrixType::RowLengthsVector rowLenghts; + IndexType N = mesh.getNumberOfCells(); + rowLenghts.setSize(N); + rowLenghts.setValue(5); + _matrix.setDimensions( N, N ); + _matrix.setRowLengths(rowLenghts); + + for (IndexType i = 0; i < N; i++) + initLinearSystem(mesh, i); + + } + template< typename Vector, typename Solver > + void solve(Solver & solver, const Vector & rhs, Vector & result) + { + solver.setMatrix(_matrix); + solver.solve(rhs,result); + } + + MatrixType _matrix; +}; + + + + +#include "tnlPoissonProblem_impl.h" + + +#endif /* TNLPOISSONPROBLEM_H */ diff --git a/examples/incompressible-navier-stokes/tnlPoissonProblem_impl.h b/examples/incompressible-navier-stokes/tnlPoissonProblem_impl.h new file mode 100644 index 0000000000000000000000000000000000000000..43f8c29d951edc1a6aaf05c39191f965ff473838 --- /dev/null +++ b/examples/incompressible-navier-stokes/tnlPoissonProblem_impl.h @@ -0,0 +1,37 @@ + +#ifndef TNLPOISSONPROBLEM_IMP_H +#define TNLPOISSONPROBLEM_IMP_H + +#include "tnlPoissonProblem.h" +#include <mesh/tnlGrid.h> + + + +template< typename MeshReal, + typename Device, + typename MeshIndex, + typename Real, + typename Index > +template< typename Vector > +#ifdef HAVE_CUDA +__device__ __host__ +#endif +Real +tnlPoissonProblem< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: +getValue( const MeshType& mesh, + const IndexType cellIndex, + const CoordinatesType& coordinates, + const Vector& u, + const Real& time ) const +{ + return ( u[ mesh.getCellXPredecessor( cellIndex ) ] + - 2.0 * u[ cellIndex ] + + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() + + ( u[ mesh.getCellYPredecessor( cellIndex ) ] + - 2.0 * u[ cellIndex ] + + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse(); +} + + + +#endif /* TNLPOISSONPROBLEM_IMP_H */ diff --git a/examples/incompressible-navier-stokes/visit_writer.cpp b/examples/incompressible-navier-stokes/visit_writer.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c27478eacd5849a2cd2c8425007c5c27cff3977a --- /dev/null +++ b/examples/incompressible-navier-stokes/visit_writer.cpp @@ -0,0 +1,1066 @@ +/* ************************************************************************* // +// visit_writer.c // +// ************************************************************************* */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +#include "visit_writer.h" //mozna uvozovky + + +/* + * Globals. + */ + +static FILE *fp = NULL; +static int useBinary = 0; +static int numInColumn = 0; + + +/* **************************************************************************** + * Function: end_line + * + * Purpose: + * If floats or ints have been written using the write_float or write_int + * functions, this will issue a newline (if necessary) so that a new + * heading can be placed. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void end_line(void) +{ + if (!useBinary) + { + char str2[8] = "\n"; + fprintf(fp, str2); + numInColumn = 0; + } +} + + +/* **************************************************************************** + * Function: open_file + * + * Purpose: + * Opens a file for writing and assigns the handle to the global variable + * "fp". + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void open_file(const char *filename) +{ + char full_filename[1024]; + if (strstr(filename, ".vtk") != NULL) + { + strcpy(full_filename, filename); + } + else + { + sprintf(full_filename, "%s.vtk", filename); + } + + fp = fopen(full_filename, "w+"); +} + + +/* **************************************************************************** + * Function: close_file + * + * Purpose: + * Closes the file with handle "fp" (a global variable). + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void close_file(void) +{ + end_line(); + fclose(fp); + fp = NULL; +} + + +/* **************************************************************************** + * Function: force_big_endian + * + * Purpose: + * Determines if the machine is little-endian. If so, then, for binary + * data, it will force the data to be big-endian. + * + * Note: This assumes that all inputs are 4 bytes long. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void force_big_endian(unsigned char *bytes) +{ + static int doneTest = 0; + static int shouldSwap = 0; + if (!doneTest) + { + int tmp1 = 1; + unsigned char *tmp2 = (unsigned char *) &tmp1; + if (*tmp2 != 0) + shouldSwap = 1; + doneTest = 1; + } + + if (shouldSwap & useBinary) + { + unsigned char tmp = bytes[0]; + bytes[0] = bytes[3]; + bytes[3] = tmp; + tmp = bytes[1]; + bytes[1] = bytes[2]; + bytes[2] = tmp; + } +} + +/* **************************************************************************** + * Function: force_double_big_endian + * + * Purpose: + * Determines if the machine is little-endian. If so, then, for binary + * data, it will force the data to be big-endian. + * + * Note: This assumes that all inputs are 8 bytes long. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void force_double_big_endian(unsigned char *bytes) +{ + static int doneTest = 0; + static int shouldSwap = 0; + if (!doneTest) + { + int tmp1 = 1; + unsigned char *tmp2 = (unsigned char *) &tmp1; + if (*tmp2 != 0) + shouldSwap = 1; + doneTest = 1; + } + + if (shouldSwap & useBinary) + { + unsigned char tmp = bytes[0]; + bytes[0] = bytes[7]; + bytes[7] = tmp; + tmp = bytes[1]; + bytes[1] = bytes[6]; + bytes[6] = tmp; + tmp = bytes[2]; + bytes[2] = bytes[5]; + bytes[5] = tmp; + tmp = bytes[3]; + bytes[3] = bytes[4]; + bytes[4] = tmp; + } +} + + +/* **************************************************************************** + * Function: write_string + * + * Purpose: + * Writes a character to the open file. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void write_string(const char *str) +{ + fprintf(fp, str); +} + + +/* **************************************************************************** + * Function: new_section + * + * Purpose: + * Adds a new line, provided we didn't already just do so and we are + * writing an ASCII file. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void new_section(void) +{ + if (numInColumn != 0) + end_line(); + numInColumn = 0; +} + + +/* **************************************************************************** + * Function: write_int + * + * Purpose: + * Writes an integer to the currently open file. This routine takes + * care of ASCII vs binary issues. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void write_int(int val) +{ + if (useBinary) + { + force_big_endian((unsigned char *) &val); + fwrite(&val, sizeof(int), 1, fp); + } + else + { + char str[128]; + sprintf(str, "%d ", val); + fprintf(fp, str); + if (((numInColumn++) % 9) == 8) + { + char str2[8] = "\n"; + fprintf(fp, str2); + numInColumn = 0; + } + } +} + + +/* **************************************************************************** + * Function: write_float + * + * Purpose: + * Writes an float to the currently open file. This routine takes + * care of ASCII vs binary issues. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * Modifications: + * + * Hank Childs, Fri Apr 22 09:14:44 PDT 2005 + * Make precision changes suggested by Jeff McAninch + * + * ************************************************************************* */ + +static void write_float(float val) +{ + if (useBinary) + { + force_big_endian((unsigned char *) &val); + fwrite(&val, sizeof(float), 1, fp); + } + else + { + char str[128]; + sprintf(str, "%20.12e ", val); + fprintf(fp, str); + if (((numInColumn++) % 9) == 8) + { + end_line(); + } + } +} +/* **************************************************************************** + * Function: write_double + * + * Purpose: + * Writes a double to the currently open file. This routine takes + * care of ASCII vs binary issues. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * Modifications: + * + * Hank Childs, Fri Apr 22 09:14:44 PDT 2005 + * Make precision changes suggested by Jeff McAninch + * + * ************************************************************************* */ + +static void write_double(double val) +{ + if (useBinary) + { + force_double_big_endian((unsigned char *) &val); + fwrite(&val, sizeof(double), 1, fp); + } + else + { + char str[128]; + sprintf(str, "%20.12e ", val); + fprintf(fp, str); + if (((numInColumn++) % 9) == 8) + { + end_line(); + } + } +} + + +/* **************************************************************************** + * Function: write_header + * + * Purpose: + * Writes the standard VTK header to the file. This should be the first + * thing written to the file. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static void write_header(void) +{ + fprintf(fp, "# vtk DataFile Version 2.0\n"); + fprintf(fp, "Written using VisIt writer\n"); + if (useBinary) + fprintf(fp, "BINARY\n"); + else + fprintf(fp, "ASCII\n"); +} + + +/* **************************************************************************** + * Function: write_variables + * + * Purpose: + * Writes the variables to the file. This can be a bit tricky. The + * cell data must be written first, followed by the point data. When + * writing the [point|cell] data, one variable must be declared the + * primary scalar and another the primary vector (provided scalars + * or vectors exist). The rest of the arrays are added through the + * "field data" mechanism. Field data should support groups of arrays + * with different numbers of components (ie a scalar and a vector), but + * there is a failure with the VTK reader. So the scalars are all written + * one group of field data and then the vectors as another. If you don't + * write it this way, the vectors do not show up. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +void write_variables(int nvars, int *vardim, int *centering, + const char * const * varname, double **vars, + int npts, int ncells) +{ + char str[1024]; + int i, j, first_scalar, first_vector; + int num_scalars, num_vectors; + int num_field = 0; + + new_section(); + sprintf(str, "CELL_DATA %d\n", ncells); + write_string(str); + + first_scalar = 0; + first_vector = 0; + num_scalars = 0; + num_vectors = 0; + /* The field data is where the non-primary scalars and vectors are + * stored. They must all be grouped together at the end of the point + * data. So write out the primary scalars and vectors first. + */ + for (i = 0 ; i < nvars ; i++) + { + if (centering[i] == 0) + { + int num_to_write = 0; + int should_write = 0; + + if (vardim[i] == 1) + { + if (first_scalar == 0) + { + should_write = 1; + sprintf(str, "SCALARS %s double\n", varname[i]); + write_string(str); + write_string("LOOKUP_TABLE default\n"); + first_scalar = 1; + } + else + num_scalars++; + } + else if (vardim[i] == 3) + { + if (first_vector == 0) + { + should_write = 1; + sprintf(str, "VECTORS %s double\n", varname[i]); + write_string(str); + first_vector = 1; + } + else + num_vectors++; + } + else + { + printf("Only supported variable dimensions are 1 and 3.\n"); + printf("Ignoring variable %s.\n", varname[i]); + continue; + } + + if (should_write) + { + num_to_write = ncells*vardim[i]; + for (j = 0 ; j < num_to_write ; j++) + { + write_double(vars[i][j]); + } + end_line(); + } + } + } + + first_scalar = 0; + if (num_scalars > 0) + { + sprintf(str, "FIELD FieldData %d\n", num_scalars); + write_string(str); + for (i = 0 ; i < nvars ; i++) + { + int should_write = 0; + if (centering[i] == 0) + { + if (vardim[i] == 1) + { + if (first_scalar == 0) + { + first_scalar = 1; + } + else + { + should_write = 1; + sprintf(str, "%s 1 %d double\n", varname[i], ncells); + write_string(str); + } + } + } + + if (should_write) + { + int num_to_write = ncells*vardim[i]; + for (j = 0 ; j < num_to_write ; j++) + { + write_double(vars[i][j]); + } + end_line(); + } + } + } + + first_vector = 0; + if (num_vectors > 0) + { + sprintf(str, "FIELD FieldData %d\n", num_vectors); + write_string(str); + for (i = 0 ; i < nvars ; i++) + { + int should_write = 0; + if (centering[i] == 0) + { + int num_to_write = 0; + + if (vardim[i] == 3) + { + if (first_vector == 0) + { + first_vector = 1; + } + else + { + should_write = 1; + sprintf(str, "%s 3 %d double\n", varname[i], ncells); + write_string(str); + } + } + } + + if (should_write) + { + int num_to_write = ncells*vardim[i]; + for (j = 0 ; j < num_to_write ; j++) + { + write_double(vars[i][j]); + } + end_line(); + } + } + } + + new_section(); + sprintf(str, "POINT_DATA %d\n", npts); + write_string(str); + + first_scalar = 0; + first_vector = 0; + num_scalars = 0; + num_vectors = 0; + /* The field data is where the non-primary scalars and vectors are + * stored. They must all be grouped together at the end of the point + * data. So write out the primary scalars and vectors first. + */ + for (i = 0 ; i < nvars ; i++) + { + if (centering[i] != 0) + { + int num_to_write = 0; + int should_write = 0; + + if (vardim[i] == 1) + { + if (first_scalar == 0) + { + should_write = 1; + sprintf(str, "SCALARS %s double\n", varname[i]); + write_string(str); + write_string("LOOKUP_TABLE default\n"); + first_scalar = 1; + } + else + num_scalars++; + } + else if (vardim[i] == 3) + { + if (first_vector == 0) + { + should_write = 1; + sprintf(str, "VECTORS %s double\n", varname[i]); + write_string(str); + first_vector = 1; + } + else + num_vectors++; + } + else + { + printf("Only supported variable dimensions are 1 and 3.\n"); + printf("Ignoring variable %s.\n", varname[i]); + continue; + } + + if (should_write) + { + num_to_write = npts*vardim[i]; + for (j = 0 ; j < num_to_write ; j++) + { + write_double(vars[i][j]); + } + end_line(); + } + } + } + + first_scalar = 0; + if (num_scalars > 0) + { + sprintf(str, "FIELD FieldData %d\n", num_scalars); + write_string(str); + for (i = 0 ; i < nvars ; i++) + { + int should_write = 0; + if (centering[i] != 0) + { + if (vardim[i] == 1) + { + if (first_scalar == 0) + { + first_scalar = 1; + } + else + { + should_write = 1; + sprintf(str, "%s 1 %d double\n", varname[i], npts); + write_string(str); + } + } + } + + if (should_write) + { + int num_to_write = npts*vardim[i]; + for (j = 0 ; j < num_to_write ; j++) + { + write_double(vars[i][j]); + } + end_line(); + } + } + } + + first_vector = 0; + if (num_vectors > 0) + { + sprintf(str, "FIELD FieldData %d\n", num_vectors); + write_string(str); + for (i = 0 ; i < nvars ; i++) + { + int should_write = 0; + if (centering[i] != 0) + { + int num_to_write = 0; + + if (vardim[i] == 3) + { + if (first_vector == 0) + { + first_vector = 1; + } + else + { + should_write = 1; + sprintf(str, "%s 3 %d double\n", varname[i], npts); + write_string(str); + } + } + } + + if (should_write) + { + int num_to_write = npts*vardim[i]; + for (j = 0 ; j < num_to_write ; j++) + { + write_double(vars[i][j]); + } + end_line(); + } + } + } +} + + +/* **************************************************************************** +// Function: write_point_mesh +// +// Purpose: +// Writes out a point mesh. +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// npts The number of points in the mesh. +// pts The spatial locations of the points. This array should +// be size 3*npts. The points should be encoded as: +// <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn> +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +void write_point_mesh(const char *filename, int ub, int npts, float *pts, + int nvars, int *vardim, const char * const *varnames, + double **vars) +{ + int i; + char str[128]; + int *centering = NULL; + + useBinary = ub; + open_file(filename); + write_header(); + + write_string("DATASET UNSTRUCTURED_GRID\n"); + sprintf(str, "POINTS %d float\n", npts); + write_string(str); + for (i = 0 ; i < 3*npts ; i++) + { + write_float(pts[i]); + } + + new_section(); + sprintf(str, "CELLS %d %d\n", npts, 2*npts); + write_string(str); + for (i = 0 ; i < npts ; i++) + { + write_int(1); + write_int(i); + end_line(); + } + + new_section(); + sprintf(str, "CELL_TYPES %d\n", npts); + write_string(str); + for (i = 0 ; i < npts ; i++) + { + write_int(VISIT_VERTEX); + end_line(); + } + + centering = (int *) malloc(nvars*sizeof(int)); + for (i = 0 ; i < nvars ; i++) + centering[i] = 1; + write_variables(nvars, vardim, centering, varnames, vars, npts, npts); + free(centering); + + close_file(); +} + + +/* **************************************************************************** + * Function: num_points_for_cell + * + * Purpose: + * Determines the number of points for the type of cell. + * + * Programmer: Hank Childs + * Creation: September 3, 2004 + * + * ************************************************************************* */ + +static int num_points_for_cell(int celltype) +{ + int npts = 0; + switch (celltype) + { + case VISIT_VERTEX: + npts = 1; + break; + case VISIT_LINE: + npts = 2; + break; + case VISIT_TRIANGLE: + npts = 3; + break; + case VISIT_QUAD: + npts = 4; + break; + case VISIT_TETRA: + npts = 4; + break; + case VISIT_HEXAHEDRON: + npts = 8; + break; + case VISIT_WEDGE: + npts = 6; + break; + case VISIT_PYRAMID: + npts = 5; + break; + } + return npts; +} + + +/* **************************************************************************** +// Function: write_unstructured_mesh +// +// Purpose: +// Writes out a unstructured mesh. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// npts The number of points in the mesh. +// pts The spatial locations of the points. This array should +// be size 3*npts. The points should be encoded as: +// <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn> +// ncells The number of cells. +// celltypes The type of each cell. +// conn The connectivity array. +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +void write_unstructured_mesh(const char *filename, int ub, int npts, + float *pts, int ncells, int *celltypes, int *conn, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars) +{ + int i, j; + char str[128]; + int conn_size = 0; + int *curr_conn = conn; + + useBinary = ub; + open_file(filename); + write_header(); + + write_string("DATASET UNSTRUCTURED_GRID\n"); + sprintf(str, "POINTS %d float\n", npts); + write_string(str); + for (i = 0 ; i < 3*npts ; i++) + { + write_float(pts[i]); + } + + new_section(); + for (i = 0 ; i < ncells ; i++) + { + int npts = num_points_for_cell(celltypes[i]); + + conn_size += npts+1; + } + sprintf(str, "CELLS %d %d\n", ncells, conn_size); + write_string(str); + for (i = 0 ; i < ncells ; i++) + { + int npts = num_points_for_cell(celltypes[i]); + write_int(npts); + for (j = 0 ; j < npts ; j++) + write_int(*curr_conn++); + end_line(); + } + + new_section(); + sprintf(str, "CELL_TYPES %d\n", ncells); + write_string(str); + for (i = 0 ; i < ncells ; i++) + { + write_int(celltypes[i]); + end_line(); + } + + write_variables(nvars, vardim, centering, varnames, vars, npts, ncells); + + close_file(); +} + + +/* **************************************************************************** +// Function: write_rectilinear_mesh +// +// Purpose: +// Writes out a rectilinear mesh. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// dims An array of size 3 = { nX, nY, nZ }, where nX is the +// number of points in the X-dimension, etc. +// x An array of size dims[0] that contains the x-coordinates. +// y An array of size dims[1] that contains the x-coordinates. +// z An array of size dims[2] that contains the x-coordinates. +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// Modifications: +// +// Hank Childs, Wed Apr 6 16:22:57 PDT 2005 +// Fix problem with 2D structured meshes and assessing cell count. +// +// ***************************************************************************/ + +void write_rectilinear_mesh(const char *filename, int ub, int *dims, + float *x, float *y, float *z, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars) +{ + int i, j; + char str[128]; + int npts = dims[0]*dims[1]*dims[2]; + int ncX = (dims[0] - 1 < 1 ? 1 : dims[0] - 1); + int ncY = (dims[1] - 1 < 1 ? 1 : dims[1] - 1); + int ncZ = (dims[2] - 1 < 1 ? 1 : dims[2] - 1); + int ncells = ncX*ncY*ncZ; + + useBinary = ub; + open_file(filename); + write_header(); + + write_string("DATASET RECTILINEAR_GRID\n"); + sprintf(str, "DIMENSIONS %d %d %d\n", dims[0], dims[1], dims[2]); + write_string(str); + sprintf(str, "X_COORDINATES %d float\n", dims[0]); + write_string(str); + for (i = 0 ; i < dims[0] ; i++) + write_float(x[i]); + new_section(); + sprintf(str, "Y_COORDINATES %d float\n", dims[1]); + write_string(str); + for (i = 0 ; i < dims[1] ; i++) + write_float(y[i]); + new_section(); + sprintf(str, "Z_COORDINATES %d float\n", dims[2]); + write_string(str); + for (i = 0 ; i < dims[2] ; i++) + write_float(z[i]); + + write_variables(nvars, vardim, centering, varnames, vars, npts, ncells); + + close_file(); +} + + +/* **************************************************************************** +// Function: write_regular_mesh +// +// Purpose: +// Writes out a regular mesh. A regular mesh is one where the data lies +// along regular intervals. "Brick of bytes/doubles", +// "Block of bytes/double", and MRI data all are examples of data that +// lie on regular meshes. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// dims An array of size 3 = { nX, nY, nZ }, where nX is the +// number of points in the X-dimension, etc. +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +void write_regular_mesh(const char *filename, int ub, int *dims, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars) +{ + int i; + + float *x = (float *) malloc(sizeof(float)*dims[0]); + float *y = (float *) malloc(sizeof(float)*dims[1]); + float *z = (float *) malloc(sizeof(float)*dims[2]); + + for (i = 0 ; i < dims[0] ; i++) + x[i] = (float) i; + for (i = 0 ; i < dims[1] ; i++) + y[i] = (float) i; + for (i = 0 ; i < dims[2] ; i++) + z[i] = (float) i; + + write_rectilinear_mesh(filename, ub, dims, x, y, z, nvars, vardim, + centering, varnames, vars); + + free(x); + free(y); + free(z); +} + + +/* **************************************************************************** +// Function: write_curvilinear_mesh +// +// Purpose: +// Writes out a curvilinear mesh. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// dims An array of size 3 = { nI, nJ, nK }, where nI is the +// number of points in the logical I dimension, etc. +// pts An array of size nI*nJ*nK*3. The array should be layed +// out as (pt(i=0,j=0,k=0), pt(i=1,j=0,k=0), ... +// pt(i=nI-1,j=0,k=0), pt(i=0,j=1,k=0), ...). +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// Modifications: +// +// Hank Childs, Wed Apr 6 16:22:57 PDT 2005 +// Fix problem with 2D structured meshes and assessing cell count. +// +// ***************************************************************************/ + +void write_curvilinear_mesh(const char *filename, int ub, int *dims,float *pts, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars) +{ + int i, j; + char str[128]; + int npts = dims[0]*dims[1]*dims[2]; + int ncX = (dims[0] - 1 < 1 ? 1 : dims[0] - 1); + int ncY = (dims[1] - 1 < 1 ? 1 : dims[1] - 1); + int ncZ = (dims[2] - 1 < 1 ? 1 : dims[2] - 1); + int ncells = ncX*ncY*ncZ; + + useBinary = ub; + open_file(filename); + write_header(); + + write_string("DATASET STRUCTURED_GRID\n"); + sprintf(str, "DIMENSIONS %d %d %d\n", dims[0], dims[1], dims[2]); + write_string(str); + sprintf(str, "POINTS %d float\n", npts); + write_string(str); + for (i = 0 ; i < 3*npts ; i++) + { + write_float(pts[i]); + } + + write_variables(nvars, vardim, centering, varnames, vars, npts, ncells); + + close_file(); +} + + diff --git a/examples/incompressible-navier-stokes/visit_writer.h b/examples/incompressible-navier-stokes/visit_writer.h new file mode 100644 index 0000000000000000000000000000000000000000..ae67e3ccc8de624ccf7f027cddf223cc32cd476c --- /dev/null +++ b/examples/incompressible-navier-stokes/visit_writer.h @@ -0,0 +1,271 @@ +/* ************************************************************************* // +// visit_writer.h // +// ************************************************************************* */ + +/* +// This file contains function prototypes for writing out point meshes, +// unstructured meshes, rectilinear meshes, regular meshes, and +// structured/curvilinear meshes into files that can later be read by VisIt. +// +// Each routine assumes that the data being written is three-dimensional. +// If the data is two-dimensional, you must still write out the data +// as three-dimensional (ie pad arrays so that they are the correct size, etc). +// However: the VisIt reader will determine that the data is truly two- +// dimensional and visualize it as a two-dimensional dataset. +// +// All writers have an ASCII vs Binary decision. The tradeoffs are the +// standard ones: ASCII is human readable, but slow. The +// binary is much faster, but not human readable. Note: the binary format +// is portable, since it converts all data to be big-endian (this was a +// design decision for the format the visit_writer writes to -- the VTK +// format). +// +// If you have multiple grids, you can write out one file for each grid. +// There are potential pitfalls in doing this, where extra geometry and +// interpolation problems appear along grid boundaries. For additional +// help with this issue, e-mail visit-help@llnl.gov +*/ + + +/* **************************************************************************** +// Function: write_point_mesh +// +// Purpose: +// Writes out a point mesh. +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// npts The number of points in the mesh. +// pts The spatial locations of the points. This array should +// be size 3*npts. The points should be encoded as: +// <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn> +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +void write_point_mesh(const char *filename, int useBinary, int npts, + float *pts, int nvars, int *vardim, + const char * const *varnames, double **vars); + + + +/* **************************************************************************** +// Function: write_unstructured_mesh +// +// Purpose: +// Writes out a unstructured mesh. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// npts The number of points in the mesh. +// pts The spatial locations of the points. This array should +// be size 3*npts. The points should be encoded as: +// <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn> +// ncells The number of cells. +// celltypes The type of each cell. +// conn The connectivity array. +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// Example: +// You have two triangles. The first has points (0,0,0), (0,1,0), and +// (1,1,0). The second has points (0,0,0), (1,1,0), and (1,0,0). +// +// There are four unique points. +// +// float pts[12] = { 0,0,0, 0,1,0, 1,1,0, 1,0,0 }; +// +// It is important the points list contain only unique points, +// because VisIt is not able to correctly determine the connectivity of a +// dataset when points are duplicated. +// +// There are two triangles. +// int ncells = 2; +// +// The cells are both triangles. +// int celltypes[2] = { VISIT_TRIANGLE, VISIT_TRIANGLE }; +// +// The connectivity contains indices into the points list. The indexing +// assumes that each point has size 3 (x,y,z). +// +// int conn[6] = { 0, 1, 2, 0, 2, 3 }; +// +// Hint: +// When writing an unstructured mesh, it is easy to get the orientation +// of a cell backwards. VisIt typically does okay with this, but it +// can cause problems. To test if this is happening, bring up VisIt on +// your newly outputted dataset and make a Pseudocolor plot of +// "mesh_quality/volume" for 3D datasets or "mesh_quality/area" for 2D +// datasets. If the cells are inside-out, the volumes or areas will be +// negative. +// +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +#define VISIT_VERTEX 1 +#define VISIT_LINE 3 +#define VISIT_TRIANGLE 5 +#define VISIT_QUAD 9 +#define VISIT_TETRA 10 +#define VISIT_HEXAHEDRON 12 +#define VISIT_WEDGE 13 +#define VISIT_PYRAMID 14 + +void write_unstructured_mesh(const char *filename, int useBinary, int npts, + float *pts, int ncells, int *celltypes, int *conn, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars); + + + +/* **************************************************************************** +// Function: write_regular_mesh +// +// Purpose: +// Writes out a regular mesh. A regular mesh is one where the data lies +// along regular intervals. "Brick of bytes/floats", +// "Block of bytes/floats", and MRI data all are examples of data that +// lie on regular meshes. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// dims An array of size 3 = { nX, nY, nZ }, where nX is the +// number of points in the X-dimension, etc. +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +void write_regular_mesh(const char *filename, int useBinary, int *dims, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars); + + + + +/* **************************************************************************** +// Function: write_rectilinear_mesh +// +// Purpose: +// Writes out a rectilinear mesh. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// dims An array of size 3 = { nX, nY, nZ }, where nX is the +// number of points in the X-dimension, etc. +// x An array of size dims[0] that contains the x-coordinates. +// y An array of size dims[1] that contains the x-coordinates. +// z An array of size dims[2] that contains the x-coordinates. +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// +// Example: +// You have a rectilinear mesh with x = { 0, 1, 2}, y = { 1, 1.5, 2, 3 }, +// and z = { 2.5, 3.5 }. +// +// Then dims = { 3, 4, 2 }. +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +void write_rectilinear_mesh(const char *filename, int useBinary, + int *dims, float *x, float *y, float *z, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars); + + + + +/* **************************************************************************** +// Function: write_curvilinear_mesh +// +// Purpose: +// Writes out a curvilinear mesh. +// +// +// Arguments: +// filename The name of the file to write. If the extension ".vtk" is +// not present, it will be added. +// useBinary '0' to write ASCII, !0 to write binary +// dims An array of size 3 = { nI, nJ, nK }, where nI is the +// number of points in the logical I dimension, etc. +// pts An array of size nI*nJ*nK*3. The array should be layed +// out as (pt(i=0,j=0,k=0), pt(i=1,j=0,k=0), ... +// pt(i=nI-1,j=0,k=0), pt(i=0,j=1,k=0), ...). +// nvars The number of variables. +// vardim The dimension of each variable. The size of vardim should +// be nvars. If var i is a scalar, then vardim[i] = 1. +// If var i is a vector, then vardim[i] = 3. +// centering The centering of each variable. The size of centering +// should be nvars. If centering[i] == 0, then the variable +// is cell-based. If centering[i] != 0, then the variable +// is point-based. +// vars An array of variables. The size of vars should be nvars. +// The size of vars[i] should be npts*vardim[i]. +// +// +// Programmer: Hank Childs +// Creation: September 2, 2004 +// +// ***************************************************************************/ + +void write_curvilinear_mesh(const char *filename, int useBinary, + int *dims, float *pts, + int nvars, int *vardim, int *centering, + const char * const *varnames, double **vars); + + +