libmesh5.c

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

					if(kwd->fmt[i] == 'r')
					{
						DblPtr = va_arg(par, double *);
						fscanf(msh->hdl, "%lf", DblPtr);
					}
					else
					{
						IntPtr = va_arg(par, int *);
						fscanf(msh->hdl, "%d", IntPtr);
					}
				}
			}
			else
			{
				for(i=0;i<kwd->SolSiz;i++)
				{
					if(kwd->fmt[i] == 'r')
					{
						DblPtr = va_arg(par, double *);
						ScaDblWrd(msh, (unsigned char *)DblPtr);
					}
					else
					{
						IntPtr = va_arg(par, int *);
						ScaWrd(msh, (unsigned char *)IntPtr);
					}
				}
			}
		}
	}
	else
	{
		if(msh->ver == 1)
		{
			FltSolTab = va_arg(par, float *);

			if(msh->typ & Asc)
				for(j=0;j<kwd->SolSiz;j++)
					fscanf(msh->hdl, "%f", &FltSolTab[j]);
			else
				for(j=0;j<kwd->SolSiz;j++)
					ScaWrd(msh, (unsigned char *)&FltSolTab[j]);
		}
		else if(msh->ver == 2)
		{
			DblSolTab = va_arg(par, double *);

			if(msh->typ & Asc)
				for(j=0;j<kwd->SolSiz;j++)
					fscanf(msh->hdl, "%lf", &DblSolTab[j]);
			else
				for(j=0;j<kwd->SolSiz;j++)
					ScaDblWrd(msh, (unsigned char *)&DblSolTab[j]);
		}
	}

	va_end(par);
}


/*----------------------------------------------------------*/
/* Write a full line from the current kwd					*/
/*----------------------------------------------------------*/

void GmfSetLin(int MshIdx, int KwdCod, ...)
{
	double d, *DblSolTab;
	float f, *FltSolTab;
	int i, j;
	va_list par;
	GmfMshSct *msh = MshTab[ MshIdx ];
	KwdSct *kwd = &msh->KwdTab[ KwdCod ];

	/* Start decoding the arguments */

	va_start(par, KwdCod);

	if(kwd->typ != SolKwd)
	{
		if(msh->ver == 1)
		{
			if(msh->typ & Asc)
			{
				for(i=0;i<kwd->SolSiz;i++)
				{
					if(kwd->fmt[i] == 'r')
					{
						d = va_arg(par, double);
						fprintf(msh->hdl, "%g ", (float)d);
					}
					else
					{
						j = va_arg(par, int);
						fprintf(msh->hdl, "%d ", j);
					}
				}
			}
			else
			{
				for(i=0;i<kwd->SolSiz;i++)
				{
					if(kwd->fmt[i] == 'r')
					{
						d = va_arg(par, double);
						f = d;
						RecWrd(msh, (unsigned char *)&f);
					}
					else
					{
						j = va_arg(par, int);
						RecWrd(msh, (unsigned char *)&j);
					}
				}
			}
		}
		else
		{
			if(msh->typ & Asc)
			{
				for(i=0;i<kwd->SolSiz;i++)
				{
					if(kwd->fmt[i] == 'r')
					{
						d = va_arg(par, double);
						fprintf(msh->hdl, "%.15lg ", d);
					}
					else
					{
						j = va_arg(par, int);
						fprintf(msh->hdl, "%d ", j);
					}
				}
			}
			else
			{
				for(i=0;i<kwd->SolSiz;i++)
				{
					if(kwd->fmt[i] == 'r')
					{
						d = va_arg(par, double);
						RecDblWrd(msh, (unsigned char *)&d);
					}
					else
					{
						j = va_arg(par, int);
						RecWrd(msh, (unsigned char *)&j);
					}
				}
			}
		}
	}
	else
	{
		if(msh->ver == 1)
		{
			FltSolTab = va_arg(par, float *);

			if(msh->typ & Asc)
				for(j=0;j<kwd->SolSiz;j++)
					fprintf(msh->hdl, "%g ", FltSolTab[j]);
			else
				for(j=0;j<kwd->SolSiz;j++)
					RecWrd(msh, (unsigned char *)&FltSolTab[j]);
		}
		else if(msh->ver == 2)
		{
			DblSolTab = va_arg(par, double *);

			if(msh->typ & Asc)
				for(j=0;j<kwd->SolSiz;j++)
					fprintf(msh->hdl, "%.15lg ", DblSolTab[j]);
			else
				for(j=0;j<kwd->SolSiz;j++)
					RecDblWrd(msh, (unsigned char *)&DblSolTab[j]);
		}
	}

	va_end(par);

	if(msh->typ & Asc)
		fprintf(msh->hdl, "\n");
}


/*----------------------------------------------------------*/
/* Private procedure for transmesh : copy a whole line		*/
/*----------------------------------------------------------*/

void GmfCpyLin(int InpIdx, int OutIdx, int KwdCod)
{
	double d;
	float f;
	int i, a;
	GmfMshSct *InpMsh = MshTab[ InpIdx ], *OutMsh = MshTab[ OutIdx ];
	KwdSct *kwd = &InpMsh->KwdTab[ KwdCod ];

	for(i=0;i<kwd->SolSiz;i++)
	{
		if(kwd->fmt[i] == 'r')
		{
			if(InpMsh->ver == 1)
			{
				if(InpMsh->typ & Asc)
					fscanf(InpMsh->hdl, "%f", &f);
				else
					ScaWrd(InpMsh, (unsigned char *)&f);

				d = f;
			}
			else
			{
				if(InpMsh->typ & Asc)
					fscanf(InpMsh->hdl, "%lf", &d);
				else
					ScaDblWrd(InpMsh, (unsigned char *)&d);

				f = (float)d;
			}

			if(OutMsh->ver == 1)
				if(OutMsh->typ & Asc)
					fprintf(OutMsh->hdl, "%g ", f);
				else
					RecWrd(OutMsh, (unsigned char *)&f);
			else
				if(OutMsh->typ & Asc)
					fprintf(OutMsh->hdl, "%.15g ", d);
				else
					RecDblWrd(OutMsh, (unsigned char *)&d);
		}
		else
		{
			if(InpMsh->typ & Asc)
				fscanf(InpMsh->hdl, "%d", &a);
			else
				ScaWrd(InpMsh, (unsigned char *)&a);

			if(OutMsh->typ & Asc)
				fprintf(OutMsh->hdl, "%d ", a);
			else
				RecWrd(OutMsh, (unsigned char *)&a);
		}
	}

	if(OutMsh->typ & Asc)
		fprintf(OutMsh->hdl, "\n");
}


/*----------------------------------------------------------*/
/* Find every kw present in a meshfile						*/
/*----------------------------------------------------------*/

static int ScaKwdTab(GmfMshSct *msh)
{
	int KwdCod, NexPos, CurPos, EndPos;
	char str[ GmfStrSiz ];

	if(msh->typ & Asc)
	{
		/* Scan each string in the file until the end */

		while(fscanf(msh->hdl, "%s", str) != EOF)
		{
			/* Fast test in order to reject quickly the numeric values */

			if(isalpha(str[0]))
			{
				/* Search which kwd code this string is associated with, 
					then get its header and save the curent position in file (just before the data) */

				for(KwdCod=1; KwdCod<= GmfMaxKwd; KwdCod++)
					if(!strcmp(str, KwdFmt[ KwdCod ][0]))
					{
						ScaKwdHdr(msh, KwdCod);
						break;
					}
			}
			else if(str[0] == '#')
				while(fgetc(msh->hdl) != '\n');
		}
	}
	else
	{
		/* Get file size */

		CurPos = ftell(msh->hdl);
		fseek(msh->hdl, 0, SEEK_END);
		EndPos = ftell(msh->hdl);
		fseek(msh->hdl, CurPos, SEEK_SET);

		/* Jump through kwd positions in the file */

		do
		{
			/* Get the kwd code and the next kwd position */

			ScaWrd(msh, (unsigned char *)&KwdCod);
			ScaWrd(msh, (unsigned char *)&NexPos);

			if(NexPos > EndPos)
				return(0);

			/* Check if this kwd belongs to this mesh version */

			if( (KwdCod >= 1) && (KwdCod <= GmfMaxKwd) )
				ScaKwdHdr(msh, KwdCod);

			/* Go to the next kwd */

			if(NexPos)
				fseek(msh->hdl, NexPos, SEEK_SET);
		}while(NexPos && (KwdCod != GmfEnd));
	}

	return(1);
}


/*----------------------------------------------------------*/
/* Read and setup the keyword's header						*/
/*----------------------------------------------------------*/

static void ScaKwdHdr(GmfMshSct *msh, int KwdCod)
{
	int i;
	KwdSct *kwd = &msh->KwdTab[ KwdCod ];

	if(!strcmp("i", KwdFmt[ KwdCod ][1]))
	{
		if(msh->typ & Asc)
			fscanf(msh->hdl, "%d", &kwd->NmbLin);
		else
			ScaWrd(msh, (unsigned char *)&kwd->NmbLin);
	}
	else
		kwd->NmbLin = 1;

	if(!strcmp("sr", KwdFmt[ KwdCod ][2]))
	{
		if(msh->typ & Asc)
		{
			fscanf(msh->hdl, "%d", &kwd->NmbTyp);

			for(i=0;i<kwd->NmbTyp;i++)
				fscanf(msh->hdl, "%d", &kwd->TypTab[i]);
		}
		else
		{
			ScaWrd(msh, (unsigned char *)&kwd->NmbTyp);

			for(i=0;i<kwd->NmbTyp;i++)
				ScaWrd(msh, (unsigned char *)&kwd->TypTab[i]);
		}
	}

	ExpFmt(msh, KwdCod);
	kwd->pos = ftell(msh->hdl);
}


/*----------------------------------------------------------*/
/* Expand the compacted format and compute the line size	*/
/*----------------------------------------------------------*/

static void ExpFmt(GmfMshSct *msh, int KwdCod)
{
	int i, j, TmpSiz=0;
	char chr, *InpFmt = KwdFmt[ KwdCod ][2];
	KwdSct *kwd = &msh->KwdTab[ KwdCod ];

	/* Set the kwd's type */

	if(!strlen(KwdFmt[ KwdCod ][1]))
		kwd->typ = InfKwd;
	else if(!strcmp(InpFmt, "sr"))
		kwd->typ = SolKwd;
	else
		kwd->typ = RegKwd;

	/* Get the solution-field's size */

	if(kwd->typ == SolKwd)
		for(i=0;i<kwd->NmbTyp;i++)
			switch(kwd->TypTab[i])
			{
				case GmfSca    : TmpSiz += 1; break;
				case GmfVec    : TmpSiz += msh->dim; break;
				case GmfSymMat : TmpSiz += (msh->dim * (msh->dim+1)) / 2; break;
				case GmfMat    : TmpSiz += msh->dim * msh->dim; break;
			}

	/* Scan each character from the format string */

	i = 0;

	while(i < (int) strlen(InpFmt))
	{
		chr = InpFmt[ i++ ];

		if(chr == 'd')
		{
			chr = InpFmt[i++];

			for(j=0;j<msh->dim;j++)
				kwd->fmt[ kwd->SolSiz++ ] = chr;
		}
		else if(chr == 's')
		{
			chr = InpFmt[i++];

			for(j=0;j<TmpSiz;j++)
				kwd->fmt[ kwd->SolSiz++ ] = chr;
		}
		else
			kwd->fmt[ kwd->SolSiz++ ] = chr;
	}
}


/*----------------------------------------------------------*/
/* Read a four bytes word in a mesh file					*/
/*----------------------------------------------------------*/

static void ScaWrd(GmfMshSct *msh, unsigned char *wrd)
{
	unsigned char swp;

	fread(wrd, WrdSiz, 1, msh->hdl);

	if(msh->cod == 1)
		return;

	swp = wrd[3];
	wrd[3] = wrd[0];
	wrd[0] = swp;

	swp = wrd[2];
	wrd[2] = wrd[1];
	wrd[1] = swp;
}


/*----------------------------------------------------------*/
/* Read an eight bytes word in a mesh file					*/
/*----------------------------------------------------------*/

static void ScaDblWrd(GmfMshSct *msh, unsigned char *wrd)
{
	int i;
	unsigned char swp;

	fread(wrd, WrdSiz, 2, msh->hdl);

	if(msh->cod == 1)
		return;

	for(i=0;i<4;i++)
	{
		swp = wrd[7-i];
		wrd[7-i] = wrd[i];
		wrd[i] = swp;
	}
}


/*----------------------------------------------------------*/
/* Write a four bytes word in a mesh file					*/
/*----------------------------------------------------------*/

static void RecWrd(GmfMshSct *msh, unsigned char *wrd)
{
	fwrite(wrd, WrdSiz, 1, msh->hdl);
}


/*----------------------------------------------------------*/
/* Write an eight bytes word in a mesh file					*/
/*----------------------------------------------------------*/

static void RecDblWrd(GmfMshSct *msh, unsigned char *wrd)
{
	fwrite(wrd, WrdSiz, 2, msh->hdl);
}


/*----------------------------------------------------------*/
/* Fortran 77 API											*/
/*----------------------------------------------------------*/

int call(gmfopenmeshf77)(char *FilNam, int *mod, int *ver, int *dim, int StrSiz)
{
	int i;
	char TmpNam[ GmfStrSiz ];

	for(i=0;i<StrSiz;i++)
		TmpNam[i] = FilNam[i];

	TmpNam[ StrSiz ] = 0;

	if(*mod == GmfRead)
		return(GmfOpenMesh(TmpNam, *mod, ver, dim));
	else
		return(GmfOpenMesh(TmpNam, *mod, *ver, *dim));
}

int call(gmfclosemeshf77)(int *idx)
{
	return(GmfCloseMesh(*idx));
}

int call(gmfstatkwdf77)(int *MshIdx, int *KwdIdx, int *NmbTyp, int *SolSiz, int *TypTab)
{
	if(!strcmp(KwdFmt[ *KwdIdx ][2], "sr"))
		return(GmfStatKwd(*MshIdx, *KwdIdx, NmbTyp, SolSiz, TypTab));
	else
		return(GmfStatKwd(*MshIdx, *KwdIdx));
}

int call(gmfgotokwdf77)(int *MshIdx, int *KwdIdx)
{
	return(GmfGotoKwd(*MshIdx, *KwdIdx));
}

int call(gmfsetkwdf77)(int *MshIdx, int *KwdIdx, int *NmbLin, int *NmbTyp, int *TypTab)
{
	if(!strcmp(KwdFmt[ *KwdIdx ][2], "sr"))
		return(GmfSetKwd(*MshIdx, *KwdIdx, *NmbLin, *NmbTyp, TypTab));
	else if(strlen(KwdFmt[ *KwdIdx ][1]))
		return(GmfSetKwd(*MshIdx, *KwdIdx, *NmbLin));
	else
		return(GmfSetKwd(*MshIdx, *KwdIdx));
}

int call(gmfgetvertex2df77)(int *MshIdx, float *x, float *y, int *ref)
{
	GmfGetLin(*MshIdx, GmfVertices, x, y, ref);
	return(1);
}

int call(gmfgetvertex3df77)(int *MshIdx, float *x, float *y, float *z, int *ref)
{
	GmfGetLin(*MshIdx, GmfVertices, x, y, z, ref);
	return(1);
}

int call(gmfsetvertex2df77)(int *MshIdx, float *x, float *y, int *ref)
{
	GmfSetLin(*MshIdx, GmfVertices, *x, *y, *ref);
	return(1);
}

int call(gmfsetvertex3df77)(int *MshIdx, float *x, float *y, float *z, int *ref)
{
	GmfSetLin(*MshIdx, GmfVertices, *x, *y, *z, *ref);
	return(1);
}

int call(gmfgettrianglef77)(int *MshIdx, int *p1, int *p2, int *p3, int *ref)
{
	GmfGetLin(*MshIdx, GmfTriangles, p1, p2, p3, ref);
	return(1);
}

int call(gmfsettrianglef77)(int *MshIdx, int *p1, int *p2, int *p3, int *ref)
{
	GmfSetLin(*MshIdx, GmfTriangles, *p1, *p2, *p3, *ref);
	return(1);
}

int call(gmfgettetrahedronf77)(int *MshIdx, int *p1, int *p2, int *p3, int *p4, int *ref)
{
	GmfGetLin(*MshIdx, GmfTetrahedra, p1, p2, p3, p4, ref);
	return(1);
}

int call(gmfsettetrahedronf77)(int *MshIdx, int *p1, int *p2, int *p3, int *p4, int *ref)
{
	GmfSetLin(*MshIdx, GmfTetrahedra, *p1, *p2, *p3, *p4, *ref);
	return(1);
}

int call(gmfgetedgef77)(int *MshIdx, int *p1, int *p2, int *ref)
{
	GmfGetLin(*MshIdx, GmfEdges, p1, p2, ref);
	return(1);
}

int call(gmfsetedgef77)(int *MshIdx, int *p1, int *p2, int *ref)
{
	GmfSetLin(*MshIdx, GmfEdges, *p1, *p2, *ref);
	return(1);
}

int call(gmfgetquadrilateralf77)(int *MshIdx, int *p1, int *p2, int *p3, int *p4, int *ref)
{
	GmfGetLin(*MshIdx, GmfQuadrilaterals, p1, p2, p3, p4, ref);
	return(1);
}

int call(gmfsetquadrilateralf77)(int *MshIdx, int *p1, int *p2, int *p3, int *p4, int *ref)
{
	GmfSetLin(*MshIdx, GmfQuadrilaterals, *p1, *p2, *p3, *p4, *ref);
	return(1);
}

int call(gmfgethexahedronf77)(int *MshIdx, int *p1, int *p2, int *p3, int *p4, int *p5, int *p6, int *p7, int *p8, int *ref)
{
	GmfGetLin(*MshIdx, GmfHexahedra, p1, p2, p3, p4, p5, p6, p7, p8, ref);
	return(1);
}

int call(gmfsethexahedronf77)(int *MshIdx, int *p1, int *p2, int *p3, int *p4, int *p5, int *p6, int *p7, int *p8, int *ref)
{
	GmfSetLin(*MshIdx, GmfHexahedra, *p1, *p2, *p3, *p4, *p5, *p6, *p7, *p8, *ref);
	return(1);
}

int call(gmfsetsolf77)(int *MshIdx, int *kwd, int *SolTab)
{
	GmfSetLin(*MshIdx, *kwd, SolTab);
	return(1);
}

int call(gmfgetsolf77)(int *MshIdx, int *kwd, int *SolTab)
{
	GmfGetLin(*MshIdx, *kwd, SolTab);
	return(1);
}