Click here to Skip to main content
12,065,045 members (46,522 online)
Click here to Skip to main content

Tagged as

Stats

14K views
1.2K downloads
12 bookmarked
Posted

A simple graphics library in C with BMP and WMF export

, 15 Jul 2011 CPOL
A graphics library to export graphical output to BMP or WMF.
BMPWMF_C_demo.zip
potential flow
potential flow.vcproj.uhwang-PC.uhwang.user
uniform-dipole.bmp
uniform-rnkovl.bmp
uniform-source.bmp
uniform-vortex.bmp
simple motion
simple motion.vcproj.uhwang-PC.uhwang.user
plate
plate.vcproj.uhwang-PC.uhwang.user
u-plate.bmp
libtg.zip
/*
	SIMPLE FLOW SIMULATION AROUND A PLATE
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#include "libtg.h"

#define PI				3.141592653589793239
#define U				1.0
#define X2(x)			((x)*(x))
#define	TUPI			(2.0*PI)
#define FOPI			(2*TUPI)
#define R2PI			(1/TUPI)
#define R4PI			(1/FOPI)
#define	ALPHA			0.0
#define toradian(x)		((x)*PI/180)

/*----------------- Lumpted Vortex Method -------------------------*/
/* Gaussian elimination */

typedef float** Array2d;
typedef float* Array1d;

void solveMatrix(Array2d A, Array1d x, Array1d b, int nx)
{
	int i, j, count, n = nx;
	float ratio=0, temp;

	for(i=0;i<(n-1);i++){
		for(j=(i+1);j<n;j++){
			ratio = A[j][i] / A[i][i];
				for(count=i;count<n;count++) {
					A[j][count] -= (ratio * A[i][count]);
                }
			b[j] -= (ratio * b[i]);
		}
	}

	/* Back substitution */
	x[n-1] = b[n-1] / A[n-1][n-1];

	for(i=(n-2);i>=0;i--){
		temp = b[i];
		for(j=(i+1);j<n;j++){
			temp -= (A[i][j] * x[j]);
		}
		x[i] = temp / A[i][i];
	}
}

Array2d createArray2D(int row, int col)
{
	int error = 0;
	int i,ii;
	Array2d p=(Array2d)malloc(row*sizeof(Array1d));

	if(!p) return p;

	for(i=0; i<row; i++) {
		p[i] = (Array1d)malloc(col*sizeof(float));
		if(!p[i]) { error = 1; break; }
	}

	if(error) {
		for(ii = 0; ii < i; ii++) free(p[ii]);
		free(p);
		p=NULL;
	}
	return p;
}

void freeArray2D(Array2d p, int row)
{
	int i;
	for(i = 0; i < row; i++) free(p[i]);
	free(p); p = NULL;
}

Array1d createArray1d(int nelem)
{
	Array1d p = (Array1d)malloc(nelem*sizeof(float));
	return p;
}

void freeArray1d(Array1d p)
{
	if(p) free(p); p = NULL;
}

void getVelocity_UPlate(Array1d gam, Array1d xv, float uu, float vv, float xs, float ys, int npan, float* velu, float* velv)
{
	int j=0;
	float r2=0;

	*velu = *velv = 0;

	for(j = 0; j < npan; j++)
	{
		r2 = (xs-xv[j])*(xs-xv[j])+ys*ys;
		*velu += -gam[j]*ys/r2;
		*velv += -gam[j]*(xs-xv[j])/r2;
	}
	*velu = (float)(uu + *velu*R2PI);
	*velv = (float)(vv - *velv*R2PI);
}

void Euler_UPlate(Array1d gam, Array1d xv, float uu, float vv, float xs, float ys, int npan, float dell, float* xm1p, float* ym1p)
{
	float um,vm,qm,ump,vmp,qmp;

	um = vm = 0;
	getVelocity_UPlate(gam, xv, uu, vv, xs, ys, npan, &um, &vm);
	qm = (float)sqrt(um*um + vm*vm);

	/* 1st order forward difference*/
	*xm1p = xs + um/qm * dell;
	*ym1p = ys + vm/qm * dell;

	ump = vmp = 0;
	getVelocity_UPlate(gam, xv, uu, vv, *xm1p, *ym1p, npan, &ump, &vmp);
	qmp = (float)sqrt(ump*ump + vmp*vmp);

	/* central difference method : improved method*/
	*xm1p = (float)(xs + 0.5 * ( um / qm + ump / qmp ) * dell);
	*ym1p = (float)(ys + 0.5 * ( vm / qm + vmp / qmp ) * dell);
}

/*
	dy/dx = f(x,y) ... x = x0, y = y0
	
	k1 = h f(xi, yi)
	k2 = h f(xi + h/2, yi + k1/2)
	k3 = h f(xi + h/2, yi + k2/2)
	k4 = h f(xi + h  , yi + k3)
	delyi = 1/6 ( k1 + 2k2 + 2k3 + k4)
	
	xi+1 = xi + h
	yi+1 = yi + delyi
*/

void RungeKutta4_UPlate(Array1d gam, Array1d xv, float uu, float vv, float xs, float ys, int npan, float dell, float* xm1p, float* ym1p)
{
	float um=0,vm=0,k1=0, k2=0, k3=0, k4=0, delyi;
	float h=dell, xi=xs,yi=ys;

	getVelocity_UPlate(gam, xv, uu, vv, xi, yi, npan, &um, &vm);
	k1 = h*vm/um;

	um=vm=0;
	getVelocity_UPlate(gam, xv, uu, vv, xi+h/2, yi+k1/2, npan, &um, &vm);
	k2 = h*vm/um;

	um=vm=0;
	getVelocity_UPlate(gam, xv, uu, vv, xi+h/2, yi+k2/2, npan, &um, &vm);
	k3 = h*vm/um;

	um=vm=0;
	getVelocity_UPlate(gam, xv, uu, vv, xi+h, yi+k3, npan, &um, &vm);
	k4 = h*vm/um;

	delyi = (float)(1./6.*(k1 + 2*k2 + 2*k3 + k4));
	*xm1p = xi+h;
	*ym1p = yi+delyi;
}

void U_Plate(TG_PlotZone_Ptr pz, float alpha, float minx, float miny,
	float maxx, float maxy, int nsgx, int nsgy)
{
	int npan = 20;
	float chord = 1, rhs=0;
	float alrad = (float)toradian(alpha);
	int i=0,j=0,k=0;
	float xs=0, ys=0, velu=0, velv=0, r2=0, slop=0, uu=0, vv=0;
	float delx = (maxx - minx)/(float)(nsgx-1);
	float dely = (maxy - miny)/(float)(nsgy-1);
	float dell = delx, xm1p=0, ym1p=0;
	TG_Plot_Ptr p = &pz->plot;
	
	Array1d xv, xc;
	Array2d A;
	Array1d b, gam, gridy;

	xv = createArray1d(npan);
	xc = createArray1d(npan);
	A = createArray2D(npan, npan);
	b = createArray1d(npan);
	gam = createArray1d(npan);
	gridy = createArray1d(nsgy);

	rhs = (float)(-2*PI*U*sin(alrad));
	for(i = 0; i < npan; i++) 
	{
		xv[i] = (float)(1+4*i)*chord/(npan*4);
		xc[i] = (float)(3+4*i)*chord/(npan*4);
		b[i]  = (float)rhs;
	}
	
	for(i = 0; i < npan; i++)
	{
		for(j = 0; j < npan; j++)
		{
			A[i][j] = 1/(xc[i]-xv[j]);
		}
	}
	
	solveMatrix(A, gam, b, npan);
	
	for(i = 0; i < nsgy; i++) gridy[i] = miny+dely*i;
	
	uu = (float)(U*cos(alrad));
	vv = (float)(U*sin(alrad));

	p->gsave();
	p->makepen(pz, TG_GetColor(TG_MAGENTA), 0.01);
	p->moveto(pz, xs,ys);
	p->lineto(pz, (float)chord,ys);
	p->deletepen();
	p->grestore();

	for(j = 0; j < nsgy; j++) 
	{
		xs = minx;
		ys = gridy[j];
		p->moveto(pz,xs,ys);

		while(1) 
		{
			Euler_UPlate(gam, xv, uu, vv, xs, ys, npan, dell, &xm1p, &ym1p);
			/*RungeKutta4_UPlate(gam,xv,uu,vv,xs,ys,npan, dell, &xm1p, &ym1p);*/
			if(xm1p < minx || xm1p > maxx || ym1p < miny || ym1p > maxy) break;
			p->lineto(pz,xm1p, ym1p);
			xs = xm1p, ys = ym1p;
		}
	}

	freeArray1d(xv);
	freeArray1d(xc);
	freeArray2D(A, npan);
	freeArray1d(b);
	freeArray1d(gam);
	freeArray1d(gridy);
}

#define BMP

int main()
{
	int xwid=500, xhgt=500;
	TG_PlotZone_Ptr zone=NULL;
	TG_Workspace_Ptr workspace = TG_InitTGLibrary("Letter");
	float minx=-1.0f, maxx=1.0f, miny=-1.0f, maxy=1.0f;
	
	if(!workspace) return 1;
	
	zone = &workspace->zone;
	TG_SetFrame(&zone->frame, 0, 0, 3, 3, minx, maxx, miny, maxy);
	
#ifdef BMP
	TG_BMP_OpenExport(zone, "u-plate.bmp", xwid, xhgt, 1);
#else
	TG_WMF_OpenExport(zone, "u-plate.wmf");
#endif

	U_Plate(zone, 10, minx,miny,maxx,maxy, 50, 40);

#ifdef BMP
	TG_BMP_CloseExport();
#else
	TG_WMF_CloseExport();	
#endif
	
	return 0;
}

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

hus
United States United States
No Biography provided

You may also be interested in...

| Advertise | Privacy | Terms of Use | Mobile
Web03 | 2.8.160204.4 | Last Updated 15 Jul 2011
Article Copyright 2011 by hus
Everything else Copyright © CodeProject, 1999-2016
Layout: fixed | fluid