wxWidget写的FEM的网格划分程序

文件1：fem.h

// Download by http://www.codefans.net
#ifndef FEM__H
#define FEM__H
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

int Gauss(double a[],double b[],int n);    //全选主元高斯消去法

/*double GaussIntegral(int n,int js[],    //计算多重积分的高斯方法        
             double(*fn)(int n,double x[]),
             void(*fnGaussLimit)(int j, int n,double x[],double y[]));  */

struct ETNode{                            //单元结点结构体
    double x,y;                            //单元结点坐标
    int number;                            //单元结点在总体区域划分中的编号
};

struct ElementTriangle{                    //三角形单元结构体    
    ETNode nd[3];                        //存储相对应的总体结点号
    double    a[3],b[3],c[3];                //基函数的系数
    double    A;                            //单元面积
    double    Aij[3][3];                    //单元有限元特征式系数矩阵
    double    fi[3];                        //单元上的积分值
    
        ElementTriangle *rightNeighbor;  //右是邻居
        
            ElementTriangle *leftChild;         //左是内层三角形
};


//    int    i,j,k,                            
//        id;                                
//     int nx=21,ny=21;                
//     double Lx=1,Ly=1;                
//    double    D;                        
//     int iNode=(nx+1)*(nx+1);        
//    double* pMatrix;                    
//    double* pMf;                        
//    ElementTriangle* pE;                
//    double ai,bi,ci;                    
//-----------------------------------------------------

//double fn(int n,double x[2])            
//{
//    return(ai+bi*x[0]+ci*x[1]);
//}


/*
void run(void)
{
    
//    pMatrix=(double*)malloc(iNode*iNode*sizeof(double));
//    pE=(ElementTriangle*)malloc(nx*ny*2*sizeof(ElementTriangle));
//    pMf=(double*)malloc(iNode*sizeof(double));

    try{
    //---------------------------------------------------------

    for(id=0;id<nx*ny*2;id++)
    {
        for(i=0;i<3;i++)
        {
            if(i==0)        j=1,k=2;
            else if(i==1)    j=2,k=0;
            else if(i==2)    j=0,k=1;
            
            pE[id].A=( (pE[id].nd[j].x-pE[id].nd[i].x)*(pE[id].nd[k].y-pE[id].nd[i].y)-
                (pE[id].nd[j].y-pE[id].nd[i].y)*(pE[id].nd[k].x-pE[id].nd[i].x) )/2.0;
            D=2.0*pE[id].A;
            pE[id].a[i]=( pE[id].nd[j].x*pE[id].nd[k].y- pE[id].nd[k].x*pE[id].nd[j].y )/D;
            pE[id].b[i]=( pE[id].nd[j].y-pE[id].nd[k].y )/D;
            pE[id].c[i]=( pE[id].nd[k].x-pE[id].nd[j].x )/D;
        }
    }printf("OK!\n");
    
    int l,m;
    for(i=0;i<nx*ny*2;i++)
    {    for(l=0;l<3;l++)                                   
            for(m=0;m<3;m++)                                     // Respaired
            {    
                pE[i].Aij[l][m]=( pE[i].b[l]*pE[i].b[m] +
                    pE[i].c[l]*pE[i].c[m] ) * pE[i].A;
            }
    }

    printf("OK!\n");



 
    int idx=0;                                    
    for(i=0;i<2*nx*ny;i++)
    {   for(idx=0;idx<3;idx++)
       {
         pE[i].fi[idx]=(1.0/3.0)*(0.5)*(Lx/nx)*(Ly/ny);         // Respaired
    }
    }

        printf("OK!\n");

    for(idx=0;idx<nx*ny*2;idx++)        
        for(i=0;i<3;i++)
        {
            for(j=0;j<3;j++)
                pMatrix[pE[idx].nd[i].number*iNode+pE[idx ].nd[j].number] +=pE[idx].Aij[i][j];

            pMf[ pE[idx].nd[i].number ]+=pE[idx].fi[i];
        }
    printf("OK!\n");

    double dBig=pow(10,20);            
    double Ur=0.0;                    
    for(i=0;i<nx+1;i++)
    {    j=nx+1;
        pMatrix[(j*nx+i)*iNode+(j*nx+i)]*=dBig;
        pMf[(j*nx+i)]=pMatrix[(j*nx+i)*iNode+(j*nx+i)]*Ur;
    }
    for(j=0;j<nx+1;j++)
    {    i=(nx+1)*(j+1)-1;
        pMatrix[i*iNode+i]*=dBig;
        pMf[i]=pMatrix[i*iNode+i]*Ur;
    }

    Gauss(pMatrix,pMf,iNode);            
    printf("OK!\n");
    
    FILE *wfp;                        
    if((wfp=fopen("dat.txt","w+"))==NULL)
        printf("Cann't open the file... ");

    for(i=0;i<iNode;i++)
        fprintf(wfp,"%f\n", pMf[i]);
    
    printf("OK!\n");

    fclose(wfp);
}
catch(...)
{    
    printf("Error occured...\n");
}


    printf("Please press Enter to exit...");
    getchar();
}

int Gauss(double a[],double b[],int n)
{
    int *js,l,k,i,j,is,p,q;
    double d,t;
    js=(int*)malloc(n*sizeof(int));
    l=1;
    for(k=0;k<=n-2;k++)
    {
        d=0.0;
        for(i=k;i<=n-1;i++)
            for(j=k;j<=n-1;j++)
            {
                t=fabs(a[i*n+j]);
                if(t>d) { d=t; js[k]=j; is=i;}
            }
            if(d+1.0==1.0) l=0;
            else
            {
                if(js[k]!=k)
                    for(i=0;i<=n-1;i++)
                    {
                        p=i*n+k; q=i*n+js[k];
                        t=a[p]; a[p]=a[q]; a[q]=t;
                    }
                    if(is!=k)
                    {
                        for(j=k;j<=n-1;j++)
                        {
                            p=k*n+j; q=is*n+j;
                            t=a[p]; a[p]=a[q]; a[q]=t;
                        }
                        t=b[k]; b[k]=b[is]; b[is]=t;
                    }
            }
            if(l==0)
            {
                free(js);
                printf("Gauss funtion failed 1...\n");
                return(0);
            }
            d=a[k*n+k];
            for(j=k+1;j<=n-1;j++)
            {
                p=k*n+j; a[p]=a[p]/d;
            }
            b[k]=b[k]/d;
            for(i=k+1;i<=n-1;i++)
            {
                for(j=k+1;j<=n-1;j++)
                {
                    p=i*n+j;
                    a[p]=a[p]-a[i*n+k]*a[k*n+j];
                }
                b[i]=b[i]-a[i*n+k]*b[k];
            }
    }
    d=a[(n-1)*n+n-1];
    if(fabs(d)+1.0==1.0)
    {
        free(js);
        printf("Gauss funtion failed 2...\n");
        return(0);
    }
    b[n-1]=b[n-1]/d;
    for(i=n-2;i>=0;i--)
    {
        t=0.0;
        for(j=i+1;j<=n-1;j++)
            t=t+a[i*n+j]*b[j];
        b[i]=b[i]-t;
    }
    js[n-1]=n-1;
    for(k=n-1;k>=0;k--)
        if(js[k]!=k)
        {
            t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;
        }
    free(js);
    return(1);
}
*/
#endif

文件2（画图类）：MySecondFrameFrm.cpp

///-----------------------------------------------------------------
///
/// @file      MySecondFrameFrm.cpp
/// @author    kangjiming
/// Created:   5/8/2012 1:19:13 PM
/// @section   DESCRIPTION
///            MySecondFrameFrm class implementation
///
///------------------------------------------------------------------

#include "MySecondFrameFrm.h"
#include <stdlib.h>
#include <time.h>


#define Random(x) (rand() % x)
#define Randomize() (srand((unsigned int)time(NULL)))
//Do not add custom headers between
//Header Include Start and Header Include End
//wxDev-C++ designer will remove them
Header Include Start
Header Include End

//----------------------------------------------------------------------------
// MySecondFrameFrm
//----------------------------------------------------------------------------
//Add Custom Events only in the appropriate block.
//Code added in other places will be removed by wxDev-C++
Event Table Start
BEGIN_EVENT_TABLE(MySecondFrameFrm,wxFrame)
    Manual Code Start
    Manual Code End
    EVT_PAINT(MySecondFrameFrm::OnPaint)

    EVT_CLOSE(MySecondFrameFrm::OnClose)
END_EVENT_TABLE()
Event Table End
/*double GaussIntegral(int n,int js[],    //¼ÆËã¶àÖØ»ý·ÖµÄ¸ß˹·½·¨        
             double(*fn)(int n,double x[]),
             void(*fnGaussLimit)(int j, int n,double x[],double y[]));  */

 double* pMatrix;                    
    double* pMf;                        
    ElementTriangle* pE;
    
MySecondFrameFrm::MySecondFrameFrm(wxWindow *parent, wxWindowID id, const wxString &title, const wxPoint &position, const wxSize& size, long style)
: wxFrame(parent, id, title, position, size, style)
{
    CreateGUIControls();
}

MySecondFrameFrm::~MySecondFrameFrm()
{
}

void MySecondFrameFrm::CreateGUIControls()
{
    //Do not add custom code between
    //GUI Items Creation Start and GUI Items Creation End
    //wxDev-C++ designer will remove them.
    //Add the custom code before or after the blocks
    GUI Items Creation Start

    SetTitle(wxT("MySecondFrame"));
    SetIcon(wxNullIcon);
    SetSize(8,8,721,460);
    Center();
    
    GUI Items Creation End
}

void MySecondFrameFrm::OnClose(wxCloseEvent& event)
{
    Destroy();
}

//画出三角形单元（同时迭代剖分确定三角形的位置）
void MySecondFrameFrm::show( ElementTriangle *pE,wxPaintDC &pdc)
{
    static int m=0;
 
    m++;
    if(m>30)return;
    static wxPoint tmp[3];
       static wxPoint Water[3];
    for(int s=0;s<3;s++){
        tmp[s].x= pE->nd[s].x;tmp[s].y= pE->nd[s].y;
        }
     
        Water[0].x=(tmp[0].x+tmp[2].x)/2; Water[0].y=(tmp[0].y+tmp[2].y)/2;
        Water[1].x=(tmp[0].x+tmp[1].x)/2; Water[1].y=(tmp[0].y+tmp[1].y)/2;
        Water[2].x=(tmp[2].x+tmp[1].x)/2; Water[2].y=(tmp[2].y+tmp[1].y)/2;
        
            pE->leftChild=(ElementTriangle*)malloc(sizeof(ElementTriangle));
            pE->leftChild->nd[0].x=Water[0].x;pE->leftChild->nd[0].y=Water[0].y;
            pE->leftChild->nd[1].x=Water[1].x;pE->leftChild->nd[1].y=Water[1].y;
            pE->leftChild->nd[2].x=Water[2].x;pE->leftChild->nd[2].y=Water[2].y;
            
            pE->rightNeighbor=(ElementTriangle*)malloc(sizeof(ElementTriangle));
            pE->rightNeighbor->nd[0].x=tmp[0].x;pE->rightNeighbor->nd[0].y=tmp[0].y;
            pE->rightNeighbor->nd[1].x=Water[1].x;pE->rightNeighbor->nd[1].y=Water[1].y;
            pE->rightNeighbor->nd[2].x=Water[0].x;pE->rightNeighbor->nd[2].y=Water[0].y;
            
            pE->rightNeighbor->rightNeighbor=(ElementTriangle*)malloc(sizeof(ElementTriangle));
            pE->rightNeighbor->rightNeighbor->nd[0].x=Water[1].x;pE->rightNeighbor->rightNeighbor->nd[0].y=Water[1].y;
            pE->rightNeighbor->rightNeighbor->nd[1].x=tmp[1].x;pE->rightNeighbor->rightNeighbor->nd[1].y=tmp[1].y;
            pE->rightNeighbor->rightNeighbor->nd[2].x=Water[2].x;pE->rightNeighbor->rightNeighbor->nd[2].y=Water[2].y;
    
                pE->rightNeighbor->rightNeighbor->rightNeighbor=(ElementTriangle*)malloc(sizeof(ElementTriangle));
            pE->rightNeighbor->rightNeighbor->rightNeighbor->nd[0].x=Water[2].x;pE->rightNeighbor->rightNeighbor->rightNeighbor->nd[0].y=Water[2].y;
            pE->rightNeighbor->rightNeighbor->rightNeighbor->nd[1].x=tmp[2].x;pE->rightNeighbor->rightNeighbor->rightNeighbor->nd[1].y=tmp[2].y;
            pE->rightNeighbor->rightNeighbor->rightNeighbor->nd[2].x=Water[0].x;pE->rightNeighbor->rightNeighbor->rightNeighbor->nd[2].y=Water[0].y;
         pdc.DrawPolygon(3, Water);
    

    //    show(pE->leftChild,pdc);
           
            show(pE->rightNeighbor,pdc);
        pE=    pE->rightNeighbor;
        
            
        
            
}
void MySecondFrameFrm::prepareRoot()
{
       int Left, Top, Right, Bottom;
    GetClientSize(&Right, &Bottom);

  //  Right *= 3; Right /= 4;
  //  Bottom *= 3; Bottom /= 4;
    Left = 0;
     wxPaintDC dc(this);
    PrepareDC(dc);
    dc.SetBackground(*wxLIGHT_GREY_BRUSH);
    dc.Clear();
    wxPoint Water[3];
    
    //pMatrix=(double*)malloc(iNode*iNode*sizeof(double));
    pE=(ElementTriangle*)malloc(sizeof(ElementTriangle));

            Water[0].x=0; Water[0].y=0;
        Water[1].x=Right; Water[1].y=Bottom;
        Water[2].x=0; Water[2].y=Bottom;
        
        
        pE->nd[0].x=Water[0].x;pE->nd[0].y=Water[0].y;
            pE->nd[1].x=Water[1].x;pE->nd[1].y=Water[1].y;
            pE->nd[2].x=Water[2].x;pE->nd[2].y=Water[2].y;
        
         dc.DrawPolygon(3, Water);
         dc.SetBrush(WaterBrush);
 wxPaintDC &pdc=dc;
 ElementTriangle * root=pE;
 
     root->rightNeighbor=(ElementTriangle*)malloc(sizeof(ElementTriangle));
            pE->rightNeighbor->nd[0].x=0;pE->rightNeighbor->nd[0].y=0;
            pE->rightNeighbor->nd[1].x=Right;pE->rightNeighbor->nd[1].y=0;
            pE->rightNeighbor->nd[2].x=Right;pE->rightNeighbor->nd[2].y=Bottom;
show(root,pdc);
//    free(pMf);
        free(pE);
   //         free(pMatrix);    
}
void MySecondFrameFrm::OnPaint(wxPaintEvent& WXUNUSED(event))
{
    
    //
   // wxPaintDC dc(this);
   // PrepareDC(dc);
   // dc.SetBackground(*wxLIGHT_GREY_BRUSH);
  //  dc.Clear();
    
    

    int Left, Top, Right, Bottom;
    GetClientSize(&Right, &Bottom);

    Right *= 3; Right /= 4;
    Bottom *= 3; Bottom /= 4;
    Left = 0;
   // Top = Bottom/8;

    int num=10;
    int step=(Right-Left)/num;
      int num2=10;
    int step2=(Bottom-Top)/num2;
   // wxPoint Water[3];
    
    
  // nx=num;
  // ny=num2;                
//     Lx=(Right-Left);
  //   Ly=Bottom-Top;                
//    double    D;                        
//    iNode=(nx+1)*(nx+1);    
        
    //Water[0].x = Left;            Water[0].y = Top;
  //  Water[1].x = Right;           Water[1].y = Top;
  //  Water[2].x = Right;  Water[2].y = Bottom;
  //  Water[3].x = 0;        Water[3].y = Bottom;
  //  wxString ss = wxString::Format(wxT("%d"),Left);
wxMessageBox(NULL,'hello world!',ss.wc_str() ,1);
   // dc.SetBrush(WaterBrush);
    //dc.DrawPolygon(4, Water);
   
   /*
        pMatrix=(double*)malloc(iNode*iNode*sizeof(double));
    pE=(ElementTriangle*)malloc(nx*ny*2*sizeof(ElementTriangle));
    pMf=(double*)malloc(iNode*sizeof(double));
    for(i=0;i<iNode*iNode;i++)
        pMatrix[i]=0;
    for(i=0;i<iNode;i++)
        pMf[i]=0;
*/

prepareRoot();
    /*
    for(j=0;j<nx;j++)
    for(i=0;i<ny;i++)
    {    
        //for the first triangle in the rectangle
        pE[i*2+j*ny*2].nd[0].x=(Lx/nx)*i;            
        pE[i*2+j*ny*2].nd[0].y=(Ly/ny)*j;
        pE[i*2+j*ny*2].nd[0].number=i+j*(nx+1);            //NO.0
        
        pE[i*2+j*ny*2].nd[1].x=(Lx/nx)*(i+1);
        pE[i*2+j*ny*2].nd[1].y=(Ly/ny)*(j+1);
        pE[i*2+j*ny*2].nd[1].number=i+1+(nx+1)*(j+1);    //NO.1
        
        pE[i*2+j*ny*2].nd[2].x=(Lx/nx)*i;
        pE[i*2+j*ny*2].nd[2].y=(Ly/ny)*(j+1);
        pE[i*2+j*ny*2].nd[2].number=i+(nx+1)*(j+1);        //NO.2
        //for the second triangle in the rectangle
        pE[i*2+j*ny*2+1].nd[0].x=(Lx/nx)*i;    
        pE[i*2+j*ny*2+1].nd[0].y=(Ly/ny)*j;
        pE[i*2+j*ny*2+1].nd[0].number=i+j*(nx+1);        //NO.0
        pE[i*2+j*ny*2+1].nd[1].x=(Lx/nx)*(i+1);
        pE[i*2+j*ny*2+1].nd[1].y=(Ly/ny)*j;
        pE[i*2+j*ny*2+1].nd[1].number=i+j*(nx+1)+1;        //NO.1
        pE[i*2+j*ny*2+1].nd[2].x=(Lx/nx)*(i+1);
        pE[i*2+j*ny*2+1].nd[2].y=(Ly/ny)*(j+1);
        pE[i*2+j*ny*2+1].nd[2].number=i+1+(nx+1)*(j+1);    //NO.2
        
        Water[0].x=pE[i*2+j*ny*2].nd[0].x; Water[0].y=pE[i*2+j*ny*2].nd[0].y;
        Water[1].x=pE[i*2+j*ny*2].nd[1].x; Water[1].y=pE[i*2+j*ny*2].nd[1].y;
        Water[2].x=pE[i*2+j*ny*2].nd[2].x; Water[2].y=pE[i*2+j*ny*2].nd[2].y;
        
         dc.DrawPolygon(3, Water);
             Water[0].x=pE[i*2+j*ny*2+1].nd[0].x; Water[0].y=pE[i*2+j*ny*2+1].nd[0].y;
        Water[1].x=pE[i*2+j*ny*2+1].nd[1].x; Water[1].y=pE[i*2+j*ny*2+1].nd[1].y;
        Water[2].x=pE[i*2+j*ny*2+1].nd[2].x; Water[2].y=pE[i*2+j*ny*2+1].nd[2].y;
         dc.DrawPolygon(3, Water);
    }
    */
  //  for(int i=0;i<num;i++){
    //    for(int j=0;j<num2;j++){
  //     pE[i][j][0].x = Left+i*step+j*step2;            pE[i][j][0].y = Top+i*step+j*step2;
  //  pE[i][j][1].x =  Left+i*step*2 +j*step2;           pE[i][j][1].y = Top+i*step*2+j*step2;
  //  pE[i][j][2].x =Left+i*step*2+2*j*step2;          pE[i][j][2].y =  Top+2*j*step2+i*step*2;
   
   //       dc.DrawPolygon(3, Water[i][j]);
   //     }
   //     }
        //run();
              

}