libmesh5.c

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

/*----------------------------------------------------------*/
/*                                                          */
/*                       LIBMESH V 5.0                      */
/*                                                          */
/*----------------------------------------------------------*/
/*                                                          */
/*    Description:        handle .meshb file format I/O     */
/*    Author:                Loic MARECHAL                  */
/*    Creation date:        feb 16 2007                     */
/*    Last modification:    feb 23 2007                     */
/*                                                          */
/*----------------------------------------------------------*/


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

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <ctype.h>
#include "libmesh5.h"


/*----------------------------------------------------------*/
/* Structures                                               */
/*----------------------------------------------------------*/

typedef struct
{
    int pos, typ, SolSiz, NmbLin, NmbTyp, TypTab[ GmfMaxTyp ];
    char fmt[ GmfMaxTyp ];
}KwdSct;

typedef struct
{
    int dim, ver, iter, mod, typ, cod, NexKwdPos;
    double angle, bbox[3][2], time;
    KwdSct KwdTab[ GmfMaxKwd + 1 ];
    FILE *hdl;
    char FilNam[ GmfStrSiz ];
}GmfMshSct;


/*----------------------------------------------------------*/
/* Defines                                                    */
/*----------------------------------------------------------*/

#define Asc 1
#define Bin 2
#define MshFil 4
#define SolFil 8
#define MaxMsh 100
#define InfKwd 1
#define RegKwd 2
#define SolKwd 3
#define WrdSiz 4


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

int NmbMsh=0;
GmfMshSct *MshTab[ MaxMsh + 1 ];
char *KwdFmt[ GmfMaxKwd + 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"},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"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"},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"BoundingBox", "", "drdr"},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"Reserved", "", ""},
    {"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"},
    {"Reserved","",""}
 };


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

static void ScaWrd(GmfMshSct *, unsigned char *);
static void ScaDblWrd(GmfMshSct *, unsigned char *);
static void RecWrd(GmfMshSct *, unsigned char *);
static void RecDblWrd(GmfMshSct *, unsigned char *);
static int ScaKwdTab(GmfMshSct *);
static void ExpFmt(GmfMshSct *, int);
static void ScaKwdHdr(GmfMshSct *, int);


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

int GmfOpenMesh(const char *FilNam, int mod, ...)
{
    int KwdCod, res, *PtrVer, *PtrDim;
    char str[ GmfStrSiz ];
    va_list par;
    GmfMshSct *msh;

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

    if( !(msh = calloc(1, sizeof(GmfMshSct))) || (NmbMsh >= MaxMsh) )
        return(0);

    MshTab[ ++NmbMsh ] = msh;

    /* Copy the FilNam into the structure */

    if(strlen(FilNam) + 7 >= GmfStrSiz)
        return(0);

    strcpy(msh->FilNam, FilNam);

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

    msh->mod = mod;

    if(strstr(msh->FilNam, ".meshb"))
        msh->typ |= (Bin | MshFil);
    else if(strstr(msh->FilNam, ".mesh"))
        msh->typ |= (Asc | MshFil);
    else if(strstr(msh->FilNam, ".solb"))
        msh->typ |= (Bin | SolFil);
    else if(strstr(msh->FilNam, ".sol"))
        msh->typ |= (Asc | SolFil);
    else
        return(0);

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

    if(msh->mod == GmfRead)
    {

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

        va_start(par, mod);
        PtrVer = va_arg(par, int *);
        PtrDim = va_arg(par, int *);
        va_end(par);

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

        if(!(msh->hdl = fopen(msh->FilNam, "rb")))
            return(0);
        /* Read the endian coding tag, the mesh version and the mesh dimension (mandatory kwd) */

        if(msh->typ & Bin)
        {
            fread((unsigned char *)&msh->cod, WrdSiz, 1, msh->hdl);

            if( (msh->cod != 1) && (msh->cod != 16777216) )
                return(0);

            ScaWrd(msh, (unsigned char *)&msh->ver);
            ScaWrd(msh, (unsigned char *)&KwdCod);

            if(KwdCod != GmfDimension)
                return(0);

            ScaWrd(msh, (unsigned char *)&KwdCod);
            ScaWrd(msh, (unsigned char *)&msh->dim);
        }
        else
        {
            do
            {
                res = fscanf(msh->hdl, "%s", str);
            }while( (res != EOF) && strcmp(str, "MeshVersionFormatted") );

            if(res == EOF)
                return(0);

            fscanf(msh->hdl, "%d", &msh->ver);

            do
            {
                res = fscanf(msh->hdl, "%s", str);
            }while( (res != EOF) && strcmp(str, "Dimension") );

            if(res == EOF)
                return(0);

            fscanf(msh->hdl, "%d", &msh->dim);
        }

        if( (msh->dim != 2) && (msh->dim != 3) )
            return(0);

        (*PtrVer) = msh->ver;
        (*PtrDim) = msh->dim;

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

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

        if(!ScaKwdTab(msh))
            return(0);
    }
    else if(msh->mod == GmfWrite)
    {

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

        msh->cod = 1;

        /* Check if the user provided a valid version number and dimension */

        va_start(par, mod);
        msh->ver = va_arg(par, int);
        msh->dim = va_arg(par, int);
        va_end(par);

        if( (msh->ver != 1) && (msh->ver != 2) )
            return(0);
        if( (msh->dim != 2) && (msh->dim != 3) )
            return(0);

        /* Create the mesh file */

        if(!(msh->hdl = fopen(msh->FilNam, "wb")))
            return(0);


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

        /* Write the mesh version and dimension */

        if(msh->typ & Asc)
        {
            fprintf(msh->hdl, "%s %d\n\n", KwdFmt[ GmfVersionFormatted ][0], msh->ver);
            fprintf(msh->hdl, "%s %d\n", KwdFmt[ GmfDimension ][0], msh->dim);
        }
        else
        {
            RecWrd(msh, (unsigned char *)&msh->cod);
            RecWrd(msh, (unsigned char *)&msh->ver);
            GmfSetKwd(NmbMsh, GmfDimension, 0);
            RecWrd(msh, (unsigned char *)&msh->dim);
        }
    }
    else
        return(0);

    return(NmbMsh);
}


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

int GmfCloseMesh(int MshIdx)
{
    int ret;
    GmfMshSct *msh;

    if( (MshIdx < 1) || (MshIdx > MaxMsh) )
        return(0);

    msh = MshTab[ MshIdx ];

    /* In write down the "End" kw in write mode */

    if(msh->mod == GmfWrite)
        {
          if(msh->typ & Asc)
            fprintf(msh->hdl, "\n%s\n", KwdFmt[ GmfEnd ][0]);
          else
            GmfSetKwd(MshIdx, GmfEnd, 0);
        }

    /* Close the file and free the mesh structure */
    /* modif FH */
    ret=fclose(msh->hdl);

    free(msh);
    return ret!=0;
}


/*----------------------------------------------------------*/
/* Read the number of lines and set the position to this kwd*/
/*----------------------------------------------------------*/

int GmfStatKwd(int MshIdx, int KwdCod, ...)
{
    int i, *PtrNmbTyp, *PtrSolSiz, *TypTab;
    GmfMshSct *msh;
    KwdSct *kwd;
    va_list par;

    if( (MshIdx < 1) || (MshIdx > MaxMsh) )
        return(0);

    msh = MshTab[ MshIdx ];

    if( (KwdCod < 1) || (KwdCod > GmfMaxKwd) )
        return(0);

    kwd = &msh->KwdTab[ KwdCod ];

    if(!kwd->NmbLin)
        return(0);

    /* Read further arguments if this kw is a sol */

    if(kwd->typ == SolKwd)
    {
        va_start(par, KwdCod);

        PtrNmbTyp = va_arg(par, int *);
        *PtrNmbTyp = kwd->NmbTyp;

        PtrSolSiz = va_arg(par, int *);
        *PtrSolSiz = kwd->SolSiz;

        TypTab = va_arg(par, int *);