fem2d_poisson.cpp

c++ for poisson

# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <fstream>
# include <cmath>
# include <ctime>
# include <windows.h>
# include <stdlib.h>


#include "GL/glew.h"
#include "GL/glut.h"
#include "Cg/Cg.h"
#include "Cg/cgGL.h"


using namespace std;

int main (int argc, char **argv);
void area_set ( int node_num, double node_xy[], int nnodes,
int element_num, int element_node[], double element_area[] );
int element_num, int element_node[], int nq, 
double wq[], double xq[], double yq[], double element_area[], int indx[], 

int bandwidth ( int nnodes, int element_num, int element_node[], 
int node_num, int indx[] );
void boundary ( int nx, int ny, int node_num, double node_xy[], int indx[], 
int ib, int nunk, double a[], double f[] );
void compare ( int node_num, double node_xy[], int indx[], int nunk, double f[] );


void element_write ( int nnodes, int element_num, int element_node[],
char *triangulation_txt_file_name );
void errors ( double element_area[], int element_node[], int indx[], 
double node_xy[], double f[], int element_num, int nnodes, 
int nunk, int node_num, double *el2, double *eh1 );
void exact ( double x, double y, double *u, double *dudx, double *dudy );
void grid_t6 ( int nx, int ny, int nnodes, int element_num, int element_node[] );
int i4_max ( int i1, int i2 );
int i4_min ( int i1, int i2 );

void indx_set ( int nx, int ny, int node_num, int indx[], int *nunk );
void nodes_plot ( char *file_name, int node_num, double node_xy[], 
bool node_label );
void nodes_write ( int node_num, double node_xy[], char *output_filename );
void qbf ( double x, double y, int element, int inode, double node_xy[], 
int element_node[], int element_num, int nnodes, 
int node_num, double *bb, double *bx, double *by );
void quad_a ( double node_xy[], int element_node[], 
int element_num, int node_num, int nnodes, double wq[], double xq[], 
double yq[] );
void quad_e ( double node_xy[], int element_node[], 
int element, int element_num, int nnodes, int node_num, int nqe, 
double wqe[], double xqe[], double yqe[] );
double r8_huge ( void );
double r8_max ( double x, double y );
double r8_min ( double x, double y );
int r8_nint ( double x );
void r8vec_print_some ( int n, double a[], int max_print, char *title );
double rhs ( double x, double y );
int s_len_trim ( char *s );
void solution_write ( double f[], int indx[], int node_num, int nunk, 
char *output_filename, double node_xy[] );
void timestamp ( void );
char *timestring ( void );
void triangulation_order6_plot ( char *file_name, int node_num, double node_xy[], 
int tri_num, int triangle_node[], int node_show, int triangle_show );
void xy_set ( int nx, int ny, int node_num, double xl, double xr, double yb, 
double yt, double node_xy[] );



# define NNODES 4             // 单元的节点数目
# define NQ 3                 // 积分点的数目
# define NX 500                // X方向上单元的数目(可以进行设置，确定计算规模的大小)
# define NY 500                // Y方向上单元的数目(可以进行设置，确定计算规模的大小)
# define ELEMENT_NUM ( NX - 1 ) * ( NY - 1 ) 
# define NODE_NUM ( 2 * NX - 1 ) * ( 2 * NY - 1 )


int element_node[NNODES*ELEMENT_NUM];
int ib;
int indx[NODE_NUM];


char *node_eps_file_name = "fem2d_poisson_nodes.eps";
char *node_txt_file_name = "fem2d_poisson_nodes.txt";
bool node_label;
int node_show;
double node_xy[2*NODE_NUM];
int nunk;

char *solution_txt_file_name = "fem2d_poisson_solution.txt";
int triangle_show;
char *triangulation_eps_file_name = "fem2d_poisson_elements.eps";
char *triangulation_txt_file_name = "fem2d_poisson_elements.txt";



// 以下为GPU计算所使用的结构和数据
void GPU_computing();

struct struct_TextureTarameters 
{
	char* name;
	GLenum texTarget;
	GLenum texInternalFormat;
	GLenum texFormat;
	GLenum texType;
	char* shader_source;
};

struct_TextureTarameters m_stTexPara, m_patchTexPara;	


// 	 fragment program
CGprogram m_cgFragmentProgram;//计算顶点位置的片元程序
CGprofile m_cgFragmentProfile;
GLuint m_iFBO; // The ID of frame buffer object
int m_iWidth;  // The width of render target
int m_iHeight; // The height of render target
GLuint m_iTargetTextureID; // the ID of the render target texture
GLuint m_iPatchTextureID; // the ID of the Patch texture
GLuint m_iOriginDataTexID ;  // 保存原始node_xy数据的纹理

	
GLuint  frameBufferObject;

byte m_elementCount[2*NODE_NUM][2];  // 用于保存计数

float  m_data[3*NODE_NUM];       // 用于保存RGB色彩形式的node_xy数据 
float  m_countdata[4*NODE_NUM];  // 用于保存RGBA色彩形式的计数

HWND   g_hWnd          = NULL;
HDC	   g_hDC           = NULL;
HGLRC  g_hRC           = NULL;



//****************************************************************************80

int main (int argc, char **argv)

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for FEM2D_POISSON.
//
//  Discussion:
//
//    FEM2D_POISSON solves 
//
//      -Laplacian U(X,Y) = F(X,Y)
//
//    in a rectangular region in the plane.  Along the boundary,
//    Dirichlet boundary conditions are imposed.  
//
//      U(X,Y) = G(X,Y)
//
//    The code uses continuous piecewise quadratic basis functions on 
//    triangles determined by a uniform grid of NX by NY points.
//
//  Local parameters:
//
//    Local, double A[(3*IB+1)*NUNK], the coefficient matrix.
//
//    Local, double ELEMENT_AREA[ELEMENT_NUM], the area of each element.
//
//    Local, double C[NUNK], the finite element coefficients, solution of A * C = F.
//
//    Local, double EH1, the H1 seminorm error.
//
//    Local, double EL2, the L2 error.
//
//    Local, int ELEMENT_NODE[ELEMENT_NUM*NNODES]; ELEMENT_NODE(I,J) is the global node index of
//    the local node J in element I.
//
//    Local, int ELEMENT_NUM, the number of elements.
//
//    Local, double F[NUNK], the right hand side.
//
//    Local, int IB, the half-bandwidth of the matrix.
//
//    Local, int INDX[NODE_NUM], gives the index of the unknown quantity
//    associated with the given node.
//
//    Local, int NNODES, the number of nodes used to form one element.
//
//    Local, double NODE_XY[2*NODE_NUM], the X and Y coordinates of nodes.
//
//    Local, int NQ, the number of quadrature points used for assembly.
//
//    Local, int NUNK, the number of unknowns.
//
//    Local, int NX, the number of points in the X direction.
//
//    Local, int NY, the number of points in the Y direction.
//
//    Local, double WQ[NQ], quadrature weights.
//
//    Local, double XL, XR, YB, YT, the X coordinates of
//    the left and right sides of the rectangle, and the Y coordinates
//    of the bottom and top of the rectangle.
//
//    Local, double XQ[NQ*ELEMENT_NUM], YQ[NQ*ELEMENT_NUM], the X and Y 
//    coordinates of the quadrature points in each element.
//
{


	/* GLUT initialization */
	glutInit(&argc, argv);
	/* set the displaying-window size */
	glutInitWindowSize(640, 480);
	/* set the display mode							*/
	/* --------------------							*/
	/*	GLUT_DOUBLE:	double buffer				*/
	/*	GLUT_DEPTH:		enable depth buffer			*/
	/*	GLUT_RGB:		the RGB pixel format		*/
	/*	GLUT_ALPHA:		enable the alpha channel	*/
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
	/* create the displaying window with the desired caption */
	glutCreateWindow("GPU computing for FEM");

  double xl = 0.0E+00;
  double xr = 1.0E+00;
  double yb = 0.0E+00;
  double yt = 1.0E+00;

  timestamp ( );

  cout << "\n";
  cout << "FEM2D_POISSON (C++ version):\n";
  cout << "\n";
  cout << "  Compiled on " << __DATE__ << " at " << __TIME__ << ".\n";
  cout << "\n";
  cout << "  Solution of the Poisson equation on a unit box\n";
  cout << "  in 2 dimensions.\n";
  cout << "\n";
  cout << "  - Uxx - Uyy = F(x,y) in the box\n";
  cout << "       U(x,y) = G(x,y) on the boundary.\n";
  cout << "\n";
  cout << "  The finite element method is used, with piecewise\n";
  cout << "  quadratic basis functions on 6 node triangular\n";
  cout << "  elements.\n";
  cout << "\n";
  cout << "  The corner nodes of the triangles are generated by an\n";
  cout << "  underlying grid whose dimensions are\n";
  cout << "\n";
  cout << "  NX =                 " << NX << "\n";
  cout << "  NY =                 " << NY << "\n";
  cout << "\n";
  cout << "  Number of nodes    = " << NODE_NUM << "\n";
  cout << "  Number of elements = " << ELEMENT_NUM << "\n";
//
//  Set the coordinates of the nodes.
//
  xy_set ( NX, NY, NODE_NUM, xl, xr, yb, yt, node_xy );
//
//  Organize the nodes into a grid of 6-node triangles.
//
  grid_t6 ( NX, NY, NNODES, ELEMENT_NUM, element_node );

//
//  Determine which nodes are boundary nodes and which have a
//  finite element unknown.  Then set the boundary values.
//
  indx_set ( NX, NY, NODE_NUM, indx, &nunk );

  cout << "  Number of unknowns =       " << nunk << "\n";

//
//  Make an EPS picture of the nodes.
//
  if ( NX <= 10 && NY <= 10 )
  {
    node_label = true;
    nodes_plot ( node_eps_file_name, NODE_NUM, node_xy, node_label );

    cout << "\n";
    cout << "FEM2D_POISSON:\n";
    cout << "  Wrote an EPS file\n"; 
    cout << "    \"" << node_eps_file_name << "\".\n";
    cout << "  containing a picture of the nodes.\n";
  }
//
//  Write the nodes to an ASCII file that can be read into MATLAB.
//

  nodes_write ( NODE_NUM, node_xy, node_txt_file_name );

  cout << "\n";
  cout << "FEM2D_POISSON:\n";
  cout << "  Wrote an ASCII node file\n"; 
  cout << "    " << node_txt_file_name << "\n";
  cout << "  of the form\n";
  cout << "    X(I), Y(I)\n";
  cout << "  which can be used for plotting.\n";
//
//  Make a picture of the elements.
//
  if ( NX <= 10 && NY <= 10 )
  {
    node_show = 1;
    triangle_show = 2;

    triangulation_order6_plot ( triangulation_eps_file_name, NODE_NUM, 
      node_xy, ELEMENT_NUM, element_node, node_show, triangle_show );

    cout << "\n";
    cout << "FEM2D_POISSON:\n";
    cout << "  Wrote an EPS file\n"; 
    cout << "    \"" << triangulation_eps_file_name << "\".\n";
    cout << "  containing a picture of the elements.\n";
  }
//
//  Write the elements to a file that can be read into MATLAB.
//
  element_write ( NNODES, ELEMENT_NUM, element_node, 
    triangulation_txt_file_name );

//********************************************************************
// start to test mode for explicit algorithm 
// This procedure can be regarded as explicit algorithm
// The current and previous steps are relevant
// 以上的代码实质上是网格的划分，相当于前处理程序，可以不予理会，以下代码
// 则作为显式算法的一种模式，t代表时间，20000代表进行迭代的数目。首先对每
// 个单元进行循环，更新单元中所有节点的相关信息，包括速度、加速度、应力和
// 应变。为了减少程序的复杂性，所有的力学计算过程我们都忽略不计，而采用挂
// 起的方式来控制其时间。因此，只要对以下循环进行并行化就可以了，非常简单
// 但是，程序的顺序是不能改变的，就是说程序在每一个迭代步是可以并行处理，
// 但在总的时间历程内，是顺序的，不能任意调整。

//   for (int t=0; t<20000; t++)
//   {
// 	  for ( int nele=0; nele<ELEMENT_NUM; nele++ )
// 	  {
// 			for ( int nnd = 0; nnd<NNODES; nnd++ )
// 			{
// 				//********************************************************************
// 				// caculate the acc, vel and strain, stress internal,external force
// 				// This procedure can be ignored due to test mode
// 				// Thus, we set time suspend for a while here.
// 				  int cn=element_node[nnd+nele*NNODES];
// 				  node_xy[2*cn+nnd*2]=node_xy[2*cn+nnd*2]+0.00012;
// 				  node_xy[2*cn+1+nnd*2]=node_xy[2*cn+1+nnd*2]+0.00015;
// //				  Sleep(5);
// 				//********************************************************************
// 			}
// 	  }
// 	  cout<<"\n";
// 	  cout<< t;
//   }


  GPU_computing();


// end
//********************************************************************
  nodes_write ( NODE_NUM, node_xy, node_txt_file_name );


  cout << "\n";
  cout << "FEM2D_POISSON:\n";
  cout << "  Normal end of execution.\n";

  cout << "\n";
  timestamp ( );

  system("PAUSE");

  return 0;

}
//****************************************************************************80


void CheckGLErrors (const char *label) 
{
    GLenum errCode;
    const GLubyte *errStr;
    
    if ((errCode = glGetError()) != GL_NO_ERROR) 
	{
		errStr = gluErrorString(errCode);
        return;
    }
}

void CheckCGError()
{
	CGerror error = cgGetError(); 

	if (error != CG_NO_ERROR)
	{
		const char* errorString = cgGetErrorString(error);
		return;
	}

}


void InitGPU()
{

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
	int iWidth = 2*NX - 1 ;
	int iHeight = 2*NY - 1 ;
	
	m_iWidth = iWidth;
	m_iHeight = iHeight ;
	
    gluOrtho2D(0.0, iWidth, 0.0, iHeight);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glViewport(0, 0, iWidth, iHeight);
	
	
	// create FBO and bind it as render target;
	//(that is, use offscreen render target)
	
	// create FBO (off-screen framebuffer)
    glGenFramebuffersEXT(1, &m_iFBO); 
    // bind offscreen framebuffer (that is, skip the window-specific render target)
    glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, m_iFBO);

	glEnable(GL_TEXTURE_RECTANGLE_NV);
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
	
}


void InitCG() 
{
	
	//Sets up GLEW to initialise OpenGL extensions
    int err = glewInit();
	// error check
    if (GLEW_OK != err) 
	{
		return;
    } 


	// Cg vars
	CGcontext cgContext;
	
	// set up Cg
//	cgSetErrorCallback(cgErrorCallback);
	cgContext = cgCreateContext();
	
	//
	// Search for a valid pixel shader profile in this order:
	//
	// CG_PROFILE_ARBFP1 - GL_ARB_fragment_program
	// CG_PROFILE_FP30   - GL_NV_fragment_program
	// CG_PROFILE_FP20   - NV_texture_shader & NV_register_combiners
	//
	

    if( cgGLIsProfileSupported(CG_PROFILE_FP30) )
	{
        m_cgFragmentProfile = CG_PROFILE_FP30;
	}
	else if( cgGLIsProfileSupported(CG_PROFILE_FP20) )
	{
        m_cgFragmentProfile = CG_PROFILE_FP20;
	}
	else if( cgGLIsProfileSupported(CG_PROFILE_ARBFP1) )
	{
        m_cgFragmentProfile = CG_PROFILE_ARBFP1;
	}
    else
    {
        MessageBox( NULL,"Failed to initialize pixel shader! Hardware doesn't "
			"support any of the required pixel shading extensions!",
			"ERROR",MB_OK|MB_ICONEXCLAMATION );
		return;
    }

	m_cgFragmentProfile = cgGLGetLatestProfile(CG_GL_FRAGMENT);
//	m_cgFragmentProfile = CG_PROFILE_FP30;
	cgGLSetOptimalOptions(m_cgFragmentProfile);
//	CheckCGError();
	
	
	// create fragment program
	m_cgFragmentProgram = cgCreateProgramFromFile( cgContext,
		CG_SOURCE,
		".\\FEM.cg",

		m_cgFragmentProfile,
		"main", 
		NULL );
	// load program
	CheckCGError();
	cgGLLoadProgram(m_cgFragmentProgram);
	CheckCGError();
}

void prepareData()
{
	  memset(m_elementCount, 0, 4*NODE_NUM);
	  for ( int nele=0; nele<ELEMENT_NUM; nele++ )