fem2d_poisson.cpp

c++ for poisson

	  {
			for ( int nnd = 0; nnd<NNODES; nnd++ )
			{
				  int cn=element_node[nnd+nele*NNODES];
				  if (cn % 2 ==0)
				  {
					  m_elementCount[cn+nnd*2][0] += 1;
					  m_elementCount[cn+1+nnd*2][0] += 1;
				  }
				  else
				  {
					  m_elementCount[cn+nnd*2][1] += 1;
					  m_elementCount[cn+1+nnd*2][1] += 1;
				  }
			}
	  }

	  for (int i=0; i<NODE_NUM; ++i)
	  {
		  m_data[3*i] = node_xy[2*i] ;
		  m_data[3*i+1] = node_xy[2*i+1];
		  m_data[3*i+2] = 0;

          m_countdata[4*i] = m_elementCount[2*i][0] ;
		  m_countdata[4*i+1] = m_elementCount[2*i][1] ;
		  m_countdata[4*i+2] = m_elementCount[2*i+1][0] ;
		  m_countdata[4*i+3] = m_elementCount[2*i+1][1] ;

	  }
	  
	  
// 	  for ( nele=0; nele<ELEMENT_NUM; nele++ )
// 	  {
// 			for ( int nnd = 0; nnd<NNODES; nnd++ )
// 				{
// 					int cn=element_node[nnd+nele*NNODES];
// 					m_elementCount[2*cn+nnd*2] += 1;
// 					m_elementCount[2*cn+1+nnd*2] += 1;
// 			}
// 	  }
}

void createAllTextureParameters(void) 
{
	m_stTexPara.name	= "TEXRECT - float_ARB - RGB - 32";
	m_stTexPara.texTarget = GL_TEXTURE_RECTANGLE_NV;
	m_stTexPara.texInternalFormat= GL_FLOAT_RGB32_NV;
	m_stTexPara.texFormat= GL_RGB;
	m_stTexPara.texType = GL_FLOAT;

	m_patchTexPara.name	= "TEXRECT - RECTANGLE";
	m_patchTexPara.texTarget = GL_TEXTURE_RECTANGLE_NV;
	m_patchTexPara.texInternalFormat= GL_FLOAT_RGBA16_NV;
	m_patchTexPara.texFormat= GL_RGBA;  //四个色彩分量
	m_patchTexPara.texType = GL_FLOAT;
}

// Sets up a floating point texture with NEAREST filtering.
void SetupTexturePara(const GLuint texID, int iWidth, int iHeight, struct_TextureTarameters texPara)
{
	CheckGLErrors ("SetupTexturePara()");
	// make active and bind
	glBindTexture(texPara.texTarget, texID);
	// turn off filtering and wrap modes
	glTexParameteri(texPara.texTarget, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
	glTexParameteri(texPara.texTarget, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameteri(texPara.texTarget, GL_TEXTURE_WRAP_S, GL_CLAMP);
	glTexParameteri(texPara.texTarget, GL_TEXTURE_WRAP_T, GL_CLAMP);
	// define texture with floating point format
	
	glTexImage2D(texPara.texTarget, 0, texPara.texInternalFormat,
		iWidth, iHeight, 0, texPara.texFormat, texPara.texType, 0);
	
	//	CheckGLErrors ("SetupTexturePara2()");
	
	// check error
    if (glGetError() != GL_NO_ERROR) 
	{
		exit (1);
	} 	
}

/**
* send data to texture.
*/
void SendDataToTexture (GLuint texID, float* data, int iWidth, int iHeight, struct_TextureTarameters texPara)
{
	//Transfer data to texture
	glBindTexture(texPara.texTarget, texID);
	glTexSubImage2D(texPara.texTarget, 0, 0, 0, iWidth, iHeight, 
		texPara.texFormat,GL_FLOAT, data);
	CheckGLErrors ("SendDataToTexture()");
}


/**
* send data to texture.
*/
void SendByteDataToTexture (GLuint texID, byte* data, int iWidth, int iHeight)
{
	//Transfer data to texture
	glBindTexture(m_patchTexPara.texTarget, texID);
	glTexSubImage2D(m_patchTexPara.texTarget, 0, 0, 0, iWidth, iHeight, 
		m_patchTexPara.texFormat,GL_UNSIGNED_BYTE,data);
	CheckGLErrors ("SendDataToTexture()");
}


/**
* Get data from currently texture
*/
void GetDataFromTexture(double * data)
{
    glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
    glReadPixels(0, 0, m_iWidth, m_iHeight, m_stTexPara.texFormat,
		GL_FLOAT, data);
	CheckGLErrors("GetDataFromTexture()");
	
}

void CreateTextures()
{
	createAllTextureParameters();

	CheckGLErrors ("CreateTextures()");	
	// enable texturing
    glEnable(m_stTexPara.texTarget);
	
	// create render target texture
	glGenTextures(1, &m_iTargetTextureID);
	
    // set up target textures
	int iTargetWidth = m_iWidth;
	int iTargetHeight = m_iHeight;
    SetupTexturePara(m_iTargetTextureID, iTargetWidth, iTargetHeight, m_stTexPara);
	

	// 生成输入纹理, 纹理大小与target纹理是完全一致的，而且需要每次绘制完毕之后与target纹理切换
	glGenTextures(1, &m_iOriginDataTexID);
	
    SetupTexturePara(m_iOriginDataTexID, iTargetWidth, iTargetHeight, m_stTexPara);
	SendDataToTexture(m_iOriginDataTexID, m_data, iTargetWidth, iTargetHeight, m_stTexPara);


    // create patch textures ID for element_node
    glGenTextures(1, &m_iPatchTextureID);	
	int iWidth = m_iWidth;
	int iHeight = m_iHeight;
    SetupTexturePara(m_iPatchTextureID, iWidth, iHeight, m_patchTexPara);
//	SendByteDataToTexture(m_iPatchTextureID, m_countdata, iWidth, iHeight);
	SendDataToTexture(m_iPatchTextureID, m_countdata, iWidth, iHeight, m_patchTexPara);
	

    // check if something went completely wrong
    CheckGLErrors("CreateTextures()");
}


void DrawBuffer(int num)
{
	
	// enable fragment profile
	cgGLEnableProfile(m_cgFragmentProfile);
	// bind frament program
	cgGLBindProgram(m_cgFragmentProgram);
	
	CGparameter patchParam;    //
	// and get parameter handles by name
	patchParam = cgGetNamedParameter(m_cgFragmentProgram, "count");
	
	// Set and enable the index texture and patch texture
	cgGLSetTextureParameter(patchParam, m_iPatchTextureID);
	cgGLEnableTextureParameter(patchParam);
	

	for (int t=0; t<num; t++)
	{
		if (t % 2 ==0)
		{
			// creating the framebuffer object
			glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, m_iFBO);
			
			// bind renderbuffers to framebuffer object
			glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_RECTANGLE_ARB, m_iTargetTextureID, 0);
			
			
			CGparameter originTexParam;
			// and get parameter handles by name
			originTexParam = cgGetNamedParameter(m_cgFragmentProgram, "baseinput");
			
			// Set and enable the index texture and patch texture
			cgGLSetTextureParameter(originTexParam, m_iOriginDataTexID );
			cgGLEnableTextureParameter(originTexParam);
			
			
			glPolygonMode(GL_FRONT, GL_FILL);
			// and render the quad
			{
				// render with unnormalized texcoords
				glBegin(GL_QUADS);
				glTexCoord2f(0.0, 0.0); 
				glVertex2f(0.0, 0.0);
				glTexCoord2f(m_iWidth, 0.0); 
				glVertex2f(m_iWidth, 0.0);
				glTexCoord2f(m_iWidth, m_iHeight); 
				glVertex2f(m_iWidth, m_iHeight);
				glTexCoord2f(0.0, m_iHeight); 
				glVertex2f(0.0, m_iHeight);
				glEnd();
			}
			// end time 
			glFlush();
			
			
			//	cgGLDisableTextureParameter(patchParam);
			cgGLDisableTextureParameter(originTexParam);
			
//			cgGLDisableProfile(m_cgFragmentProfile);
			

// 			glReadPixels(0, 0, m_iWidth, m_iHeight, GL_RGB, GL_FLOAT, m_data);
// 			CheckGLErrors("glReadPixels");	
// 			
// 			 
// 			 	  for (int i=0; i<NODE_NUM; ++i)
// 			 	  {
// 			 		  node_xy[2*i] = m_data[3*i] ;
// 			 		  node_xy[2*i+1] = m_data[3*i+1] ;
// 			 	  }
			
			///////////////////////////////////
			//	m_ResultData = (float * )glMapBuffer(GL_PIXEL_PACK_BUFFER_ARB, GL_READ_ONLY); 
			//	CheckGLErrors("glMapBuffer");
			//	for( int i = 0; i < m_iTargetPixelsNbr * 4; i++)
			//	{
			//		TRACE("The %d  points position result is :%f\n", i, m_ResultData[i]);
			//	}
			//   glUnmapBuffer(GL_PIXEL_PACK_BUFFER_ARB) ;	
			////////////////////////////////////////
			
			// deactivate the framebuffer object (reactivate screen as framebuffer)
			glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
			
			//	glBindBufferARB(GL_ARRAY_BUFFER_ARB, *m_iPBO) ;
			
			//    glReadPixels(0, 0, texWidth, texHeight,GL_RGB,GL_FLOAT,patch);
			// done, just do some checks if everything went smoothly.
			//    CheckFramebufferStatus();
			//   CheckGLErrors("Draw()");


		}
		else
		{
			// creating the framebuffer object
			glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, m_iFBO);
			
			// bind renderbuffers to framebuffer object
			glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_RECTANGLE_ARB, m_iOriginDataTexID, 0);
			
			
			
			// enable fragment profile
			cgGLEnableProfile(m_cgFragmentProfile);
			// bind frament program
			cgGLBindProgram(m_cgFragmentProgram);
			
			// 	CGparameter patchParam;    //
			// 	// and get parameter handles by name
			// 	patchParam = cgGetNamedParameter(m_cgFragmentProgram, "count");
			// 	
			// 	// Set and enable the index texture and patch texture
			// 	cgGLSetTextureParameter(patchParam, m_iPatchTextureID);
			// 	cgGLEnableTextureParameter(patchParam);
			
			
			CGparameter originTexParam;
			// and get parameter handles by name
			originTexParam = cgGetNamedParameter(m_cgFragmentProgram, "baseinput");
			
			
			// Set and enable the index texture and patch texture
			cgGLSetTextureParameter(originTexParam, m_iTargetTextureID );
			cgGLEnableTextureParameter(originTexParam);
			
			
			// make quad filled to hit every pixel/texel 
			glPolygonMode(GL_FRONT, GL_FILL);
			// and render the quad
			{
				// render with unnormalized texcoords
				glBegin(GL_QUADS);
				glTexCoord2f(0.0, 0.0); 
				glVertex2f(0.0, 0.0);
				glTexCoord2f(m_iWidth, 0.0); 
				glVertex2f(m_iWidth, 0.0);
				glTexCoord2f(m_iWidth, m_iHeight); 
				glVertex2f(m_iWidth, m_iHeight);
				glTexCoord2f(0.0, m_iHeight); 
				glVertex2f(0.0, m_iHeight);
				glEnd();
			}
			// end time 
			glFlush();
			
			
			//	cgGLDisableTextureParameter(patchParam);
			cgGLDisableTextureParameter(originTexParam);
			
//			cgGLDisableProfile(m_cgFragmentProfile);
			

			if (t == num-1)
			{
				glReadPixels(0, 0, m_iWidth, m_iHeight, GL_RGB, GL_FLOAT, m_data);
				CheckGLErrors("glReadPixels");	
				
				for (int i=0; i<NODE_NUM; ++i)
				{
					node_xy[2*i] = m_data[3*i] ;
					node_xy[2*i+1] = m_data[3*i+1] ;
				}
			}
			
			glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
		}
		cout<<"\n";
		cout<< t;	
	}

	cgGLDisableTextureParameter(patchParam);
	cgGLDisableProfile(m_cgFragmentProfile);

}


void GPU_computing()
{
	InitCG();
    InitGPU();
	prepareData();
	CreateTextures();

	DrawBuffer(20000);


//	GetDataFromTexture();

}

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

void element_write ( int nnodes, int element_num, int element_node[],
 char *output_filename ) 

//****************************************************************************80*
//
//  Purpose:
//
//    ELEMENT_WRITE writes the elements to a file.
//
//  Modified:
//
//    22 March 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NNODES, the number of nodes used to form one element.
//
//    Input, int ELEMENT_NUM, the number of elements.
//
//    Input, int ELEMENT_NODE[NNODES*ELEMENT_NUM]; ELEMENT_NODE(I,J) is the global 
//    index of local node I in element J.
//
//    Input, char *OUTPUT_FILENAME, the name of the file
//    in which the data should be stored.
//
{
  int element;
  int i;
  ofstream output;

  output.open ( output_filename );

  if ( !output )
  {
    cout << "\n";
    cout << "ELEMENT_WRITE - Warning!\n";
    cout << "  Could not write the node file.\n";
    return;
  }

  for ( element = 0; element < element_num; element++ )
  {
    for ( i = 0; i < nnodes; i++ )
    {
      output << setw(8)  << element_node[i+element*nnodes] << "  ";
    }
    output << "\n";
  }

  output.close ( );

  return;
}
//****************************************************************************80*



void exact ( double x, double y, double *u, double *dudx, double *dudy )

//****************************************************************************80
//
//  Purpose:
//
//    EXACT calculates the exact solution and its first derivatives.
//
//  Discussion:
//
//    The function specified here depends on the problem being
//    solved.  This is one of the routines that a user will
//    normally want to change.
//
//  Modified:
//
//    07 April 2004
//
//  Parameters:
//
//    Input, double X, Y, the coordinates of a point
//    in the region, at which the right hand side of the 
//    differential equation is to be evaluated.
//
//    Output, double *U, *DUDX, *DUDY, the value of 
//    the exact solution U and its derivatives dUdX
//    and dUdY at the point (X,Y).
//
{
# define PI 3.14159265358979323846264338327950288419716939937510

  *u    =      sin ( PI * x ) * sin ( PI * y ) + x;
  *dudx = PI * cos ( PI * x ) * sin ( PI * y ) + 1.0E+00;
  *dudy = PI * sin ( PI * x ) * cos ( PI * y );
       
  return;
# undef PI
}
//****************************************************************************80

void grid_t6 ( int nx, int ny, int nnodes, int element_num, int element_node[] )

//****************************************************************************80
//
//  Purpose:
//
//    GRID_T6 produces a grid of pairs of 6 node triangles.
//
//  Example:
//
//    Input:
//
//      NX = 4, NY = 3
//
//    Output:
//
//      ELEMENT_NODE =
//         1,  3, 15,  2,  9,  8;
//        17, 15,  3, 16,  9, 10;
//         3,  5, 17,  4, 11, 10;
//        19, 17,  5, 18, 11, 12;
//         5,  7, 19,  6, 13, 12;
//        21, 19,  7, 20, 13, 14;
//        15, 17, 29, 16, 23, 22;
//        31, 29, 17, 30, 23, 24;
//        17, 19, 31, 18, 25, 24;
//        33, 31, 19, 32, 25, 26;
//        19, 21, 33, 20, 27, 26;
//        35, 33, 21, 34, 27, 28.
//
//  Diagram:
//
//   29-30-31-32-33-34-35
//    |\ 8  |\10  |\12  |
//    | \   | \   | \   |
//   22 23 24 25 26 27 28
//    |   \ |   \ |   \ |
//    |  7 \|  9 \| 11 \|
//   15-16-17-18-19-20-21
//    |\ 2  |\ 4  |\ 6  |
//    | \   | \   | \   |
//    8  9 10 11 12 13 14
//    |   \ |   \ |   \ |
//    |  1 \|  3 \|  5 \|
//    1--2--3--4--5--6--7
//
//  Modified:
//
//    07 April 2004
//
//  Author:
//
//    John Burkardt
//