libmesh3.c

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

/*----------------------------------------------------------*/
/*															*/
/*						LIBMESH V 3.0						*/
/*															*/
/*----------------------------------------------------------*/
/*															*/
/*	Description:		handle .meshb file format I/O		*/
/*	Author:				Loic MARECHAL						*/
/*	Creation date:		aug  2 2003							*/
/*	Last modification:	jan 25 2006							*/
/*															*/
/*----------------------------------------------------------*/


/*----------------------------------------------------------*/
/* Includes													*/
/*----------------------------------------------------------*/

#include <string.h>
#include <float.h>
#include <math.h>
#include <ctype.h>
#include "libmesh3.h"


/*----------------------------------------------------------*/
/* Global variables											*/
/*----------------------------------------------------------*/

char *LM_kw_table[ LM_NBKW + 1 ][3] = 
{	{"Reserved", "", ""},
	{"MeshVersionFormatted", "", "i"},
	{"Reserved", "", ""},
	{"Dimension", "", "i"},
	{"Vertices", "i", "dri"},
	{"Edges", "i", "iii"},
	{"Triangles", "i", "iiii"},
	{"Quadrilaterals", "i", "iiiii"},
	{"Tetrahedra", "i", "iiiii"},
	{"Pentahedra", "i", "iiiiiii"},
	{"Hexahedra", "i", "iiiiiiiii"},
	{"SubDomainFromGeom", "i", "iiii"},
	{"SubDomainFromMesh", "i", "iiii"},
	{"Corners", "i", "i"},
	{"Ridges", "i", "i"},
	{"RequiredVertices", "i", "i"},
	{"RequiredEdges", "i", "i"},
	{"RequiredTriangles", "i", "i"},
	{"RequiredQuadrilaterals", "i", "i"},
	{"TangentAtEdgeVertices", "i", "iii"},
	{"NormalAtVertices", "i", "ii"},
	{"NormalAtTriangleVertices", "i", "iii"},
	{"NormalAtQuadrilateralVertices", "i", "iiii"},
	{"AngleOfCornerBound", "", "r"},
	{"Geometry", "", "c"},
	{"VertexOnGeometricVertex", "i", "ii"},
	{"VertexOnGeometricEdge", "i", "iir"},
	{"VertexOnGeometricTriangle", "i", "iirr"},
	{"VertexOnGeometricQuadrilateral", "i", "iirr"},
	{"EdgeOnGeometricEdge", "i", "ii"},
	{"TriangleOnGeometricTriangle", "i", "ii"},
	{"TriangleOnGeometricQuadrilateral", "i", "ii"},
	{"QuadrilateralOnGeometricTriangle", "i", "ii"},
	{"QuadrilateralOnGeometricQuadrilateral", "i", "ii"},
	{"MeshSupportOfVertices", "", "c"},
	{"VertexOnSupportVertex", "i", "ii"},
	{"VertexOnSupportEdge", "i", "iir"},
	{"VertexOnSupportTriangle", "i", "iirr"},
	{"VertexOnSupportQuadrilateral", "i", "iirr"},
	{"VertexOnSupportTetrahedron", "i", "iirrr"},
	{"VertexOnSupportPentahedron", "i", "iirrr"},
	{"VertexOnSupportHexahedron", "i", "iirrr"},
	{"CrackedEdges", "i", "ii"},
	{"CrackedTriangles", "i", "ii"},
	{"CrackedQuadrilaterals", "i", "ii"},
	{"EquivalentEdges", "i", "ii"},
	{"EquivalentTriangles", "i", "ii"},
	{"EquivalentQuadrilaterals", "i", "ii"},
	{"PhysicsReference", "i", "ic"},
	{"IncludeFile", "", "c"},
	{"BoundingBox", "", "drdr"},
	{"Identifier", "", "c"},
	{"IdentityOfGeometry", "", "c"},
	{"IdentityOfMeshSupport", "", "c"},
	{"End", "", ""},
	{"Reserved", "", ""},
	{"Reserved", "", ""},
	{"Reserved", "", ""},
	{"Reserved", "", ""},
	{"Tangents", "i", "dr"},
	{"Normals", "i", "dr"},
	{"TangentAtVertices", "i", "ii"},
	{"SolAtVertices", "i", "sr"},
	{"SolAtEdges", "i", "sr"},
	{"SolAtTriangles", "i", "sr"},
	{"SolAtQuadrilaterals", "i", "sr"},
	{"SolAtTetrahedra", "i", "sr"},
	{"SolAtPentahedra", "i", "sr"},
	{"SolAtHexahedra", "i", "sr"},
	{"DSolAtVertices", "i", "sr"},
	{"ISolAtVertices", "i", "i"},
	{"ISolAtEdges", "i", "ii"},
	{"ISolAtTriangles", "i", "iii"},
	{"ISolAtQuadrilaterals", "i", "iiii"},
	{"ISolAtTetrahedra", "i", "iiii"},
	{"ISolAtPentahedra", "i", "iiiiii"},
	{"ISolAtHexahedra", "i", "iiiiiiii"},
	{"Iterations","","i"},
	{"Time","","r"},
	{"VertexHack","","drdr"}
 };


/*----------------------------------------------------------*/
/* Prototypes of local procedures							*/
/*----------------------------------------------------------*/

static void write_kw(LM_mesh_struct *, int);
static int read_int(LM_mesh_struct *);
static void write_int(LM_mesh_struct *, int);
static void file2kw_tab(LM_mesh_struct *);
static void kw_tab2file(LM_mesh_struct *);
static int expand_format(LM_mesh_struct *, int, char *);
static void swap_bytes(void *, void *, int);
static void read_sol_headers(LM_mesh_struct *);


/*----------------------------------------------------------*/
/* Open a mesh file in read or write mode					*/
/*----------------------------------------------------------*/

int LM_open_mesh(const char *filename, int mode, LM_mesh_struct *mesh, ...)
{
	int i;
	va_list pa;
	

	/*---------------------*/
	/* MESH STRUCTURE INIT */
	/*---------------------*/

	/* Init the kw table */

	for(i=0;i<=LM_NBKW;i++)
	{
		mesh->kw_counters[i] = 0;
		mesh->kw_pos[i][0] = 0;
		mesh->kw_pos[i][1] = 0;
		mesh->kw_pos[i][2] = 0;
		mesh->sol_headers[i] = NULL;
	}

	/* Allocate a string large enough to contain the full filename and path plus the ".meshb" extension. */

	mesh->filename = (char *)calloc((strlen(filename) + 7), sizeof(char));
	strcpy(mesh->filename, filename);

	/* Store the opening mode (read or write) and guess the filetype (binary or ascii) depending on the extension */

	mesh->mode = mode;
	mesh->current_kw = mesh->type = 0;
	mesh->endian = 1;

	if(strstr(mesh->filename, ".meshb"))
		mesh->type |= (LM_BINARY | LM_MESH);
	else if(strstr(mesh->filename, ".mesh"))
		mesh->type |= (LM_ASCII | LM_MESH);
	else if(strstr(mesh->filename, ".solb"))
		mesh->type |= (LM_BINARY | LM_SOL);
	else if(strstr(mesh->filename, ".sol"))
		mesh->type |= (LM_ASCII | LM_SOL);
	else
		return(0);

	/* Open the file in the required mode and initialyse the mesh structure */

	if(mesh->mode == LM_READ)
	{

		/*-----------------------*/
		/* OPEN FILE FOR READING */
		/*-----------------------*/

		/* Create the name string and open the file */

		if(!(mesh->handle = fopen(mesh->filename, "rb")))
			return(0);

		/* Read the endian tag and the mesh version in binary */

		if(mesh->type & LM_BINARY)
		{
			mesh->endian = read_int(mesh);

			if( (mesh->endian != 1) && (mesh->endian != 16777216) )
				return(0);

			mesh->version = read_int(mesh);
		}


		/*------------*/
		/* KW READING */
		/*------------*/

		/* Read the list of kw present in the file */

		file2kw_tab(mesh);

		/* Check the mesh dimension */

		if(!mesh->kw_counters[ LM_Dimension ])
			return(0);

		LM_read_field(mesh, LM_Dimension, 1, &mesh->dimension);

		if( (mesh->dimension != 2) && (mesh->dimension != 3) )
			return(0);

		/* Read the meshversion in ascii case */

		if(mesh->type & LM_ASCII)
			LM_read_field(mesh, LM_MeshVersionFormatted, 1, &mesh->version);

		/* And read the extended sol headers */

		read_sol_headers(mesh);
	}
	else if(mesh->mode == LM_WRITE)
	{

		/*-----------------------*/
		/* OPEN FILE FOR WRITING */
		/*-----------------------*/

		/* Check if the user provided a valid dimension */

		va_start(pa, mesh);
		mesh->dimension = va_arg(pa, int);
		va_end(pa);

		if( (mesh->dimension != 2) && (mesh->dimension != 3) )
			return(0);

		/* If no extension has been provided, create a binary file */

		if(!(mesh->handle = fopen(mesh->filename, "wb")))
			return(0);


		/*------------*/
		/* KW WRITING */
		/*------------*/

		/* Initialyse the required fields. The kw will be stored afterward. */

		mesh->version = LM_MESH_VERSION;
		mesh->endian = 1;

		/* Write the mesh version */

		if(mesh->type & LM_ASCII)
			LM_write_field(mesh, LM_MeshVersionFormatted, 1, &mesh->version);
		else
		{
			write_int(mesh, mesh->endian);
			write_int(mesh, mesh->version);
		}

		/* Write the mesh dimension */

		LM_write_field(mesh, LM_Dimension, 1, &mesh->dimension);
	}
	else
		return(0);

	return(1);
}


/*----------------------------------------------------------*/
/* Close a meshfile in the right way						*/
/*----------------------------------------------------------*/

int LM_close_mesh(LM_mesh_struct *mesh)
{
	if(mesh->mode == LM_WRITE)
	{
		/* Test if the user wrote the "End" kw */

		if(!mesh->kw_counters[ LM_End ])
			LM_write_field(mesh, LM_End, 0, NULL);

		/* Write down the number lines written in each field to the file */

		kw_tab2file(mesh);
	}

	if(fclose(mesh->handle))
		return(0);
	else
		return(1);
}


/*----------------------------------------------------------*/
/* Bufferized read of a whole orpart of a field				*/
/*----------------------------------------------------------*/

int LM_read_field(LM_mesh_struct *mesh, int kw_code, int nbl, void *buffer)
{
	int i, j, swaped, size, *int_buffer = (int *)buffer, string_size;
	float *flt_buffer = (float *)buffer;
	char format[256], letter, str_buf[256];

	/* Check if the kw code is valid */

	if( (kw_code < 1) || (kw_code > LM_NBKW) )
		return(0);

	/* Check if this kw has a format */

	if(!strlen(LM_kw_table[ kw_code ][2]))
		return(0);

	/* If this kw is only a header, the number of lines to be read is set to one */

	if(!strlen(LM_kw_table[ kw_code ][1]))
		nbl = 1;

	/* Check if the user is not asking more lines than the remaining lines in the file */

	if(nbl > mesh->kw_counters[ kw_code ] - mesh->kw_pos[ kw_code ][2])
		nbl = mesh->kw_counters[ kw_code ] - mesh->kw_pos[ kw_code ][2];

	if(!nbl)
		return(0);

	/* Set the curent position in file to the begining of this kw's data */

	fseek(mesh->handle, mesh->kw_pos[ kw_code ][1], SEEK_SET);

	/* Transform the internal format into a "c" format string for the scanf and compute the size of 
		field's line in order to compute the right adresses in the buffer */

	size = expand_format(mesh, kw_code, format);

	if(mesh->type & LM_ASCII)
	{
		for(i=0;i<nbl;i++)
			for(j=0;j<size;j++)
				if(format[j] == 'i')
					fscanf(mesh->handle, "%d", &int_buffer[ i * size + j ]);
				else if(format[j] == 'r')
					fscanf(mesh->handle, "%g", &flt_buffer[ i * size + j ]);
				else if(format[j] == 'c')
				{
					string_size = 0;

					do
					{
						fscanf(mesh->handle, "%c", &letter);
					}while(letter != '"');

					do
					{
						fscanf(mesh->handle, "%c", &letter);
						str_buf[ string_size++ ] = letter;
					}while( (letter != '"') && (string_size <= 256) );

					str_buf[ string_size-1 ] = 0;

					memset(&flt_buffer[ i * size + j ], 0, 256);
					strcpy((char *)&flt_buffer[ i * size + j ], str_buf);
				}
	}
	else
	{
		fread(buffer, nbl * size * 4, 1, mesh->handle);

		/* Swap the bytes in the whole buffer in case of different endian */

		if(mesh->endian != 1)
			for(i=0;i<nbl*size;i++)
			{
				swap_bytes((void *)&int_buffer[i], (void *)&swaped, 4);
				int_buffer[i] = swaped;
			}
	}

	/* Then store the curent position and the total number of read lines in case we didn't read the whole data */

	mesh->kw_pos[ kw_code ][1] = ftell(mesh->handle);
	mesh->kw_pos[ kw_code ][2] += nbl;

	return(nbl);
}


/*----------------------------------------------------------*/
/* Bufferized write of a whole field or part of it			*/
/*----------------------------------------------------------*/

int LM_write_field(LM_mesh_struct *mesh, int kw_code, int nbl, void *buffer, ...)
{
	int i, j, size, *int_buffer = (int *)buffer, nbsol;
	float *flt_buffer = (float *)buffer;
	char format[256];
	va_list pa;

	/* Check if the kw code is valid */

	if( (kw_code < 1) || (kw_code > LM_NBKW) )
		return(0);

	/* Read further arguments if this kw is solution field and the extra header was not provided by the user */

	if(!mesh->sol_headers[ kw_code ] && !strcmp(LM_kw_table[ kw_code ][2], "sr"))
	{
		va_start(pa, buffer);

		nbsol = va_arg(pa, int);

		if(!(mesh->sol_headers[ kw_code ] = malloc((nbsol+2) * sizeof(int))))
			return(0);

		mesh->sol_headers[ kw_code ][0] = nbsol;
		mesh->sol_headers[ kw_code ][1] = 0;

		for(i=1;i<=nbsol;i++)
		{
			mesh->sol_headers[ kw_code ][i+1] = va_arg(pa, int);

			switch(mesh->sol_headers[ kw_code ][i+1])
			{
				case 1 : mesh->sol_headers[ kw_code ][1] += 1; break;
				case 2 : mesh->sol_headers[ kw_code ][1] += mesh->dimension; break;
				case 3 : mesh->sol_headers[ kw_code ][1] += (mesh->dimension * (mesh->dimension+1)) / 2; break;
				case 4 : mesh->sol_headers[ kw_code ][1] += mesh->dimension * mesh->dimension; break;
			}
		}

		va_end(pa);
	}

	/* If this kw is only a header, the number of lines to be read is set to one */

	if(!strlen(LM_kw_table[ kw_code ][1]))
		nbl = 1;

	if(!mesh->kw_counters[ kw_code ])
		write_kw(mesh, kw_code);

	mesh->kw_counters[ kw_code ] += nbl;

	/* Check if this kw has a format */

	if(!strlen(LM_kw_table[ kw_code ][2]))
		return(0);

	size = expand_format(mesh, kw_code, format);

	if(mesh->type & LM_ASCII)
	{
		for(i=0;i<nbl;i++)
		{
			for(j=0;j<size;j++)