rkpmdlg.cpp

无网格法：RKPM,计算平板弯曲

// RKPMDlg.cpp : implementation file
//

#include "stdafx.h"
#include "RKPM.h"
#include "RKPMDlg.h"
#include "fstream.h"
#include "iostream.h"
#include "iomanip.h"
#include "matlab.h"
#include "math.h"
#include "libmatlb.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CAboutDlg dialog used for App About

class CAboutDlg : public CDialog
{
public:
	CAboutDlg();

// Dialog Data
	//{{AFX_DATA(CAboutDlg)
	enum { IDD = IDD_ABOUTBOX };
	//}}AFX_DATA

	// ClassWizard generated virtual function overrides
	//{{AFX_VIRTUAL(CAboutDlg)
	protected:
	virtual void DoDataExchange(CDataExchange* pDX);    // DDX/DDV support
	//}}AFX_VIRTUAL

// Implementation
protected:
	//{{AFX_MSG(CAboutDlg)
	//}}AFX_MSG
	DECLARE_MESSAGE_MAP()
};

CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD)
{
	//{{AFX_DATA_INIT(CAboutDlg)
	//}}AFX_DATA_INIT
}

void CAboutDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CAboutDlg)
	//}}AFX_DATA_MAP
}

BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)
	//{{AFX_MSG_MAP(CAboutDlg)
		// No message handlers
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CRKPMDlg dialog

CRKPMDlg::CRKPMDlg(CWnd* pParent /*=NULL*/)
	: CDialog(CRKPMDlg::IDD, pParent)
{
	//{{AFX_DATA_INIT(CRKPMDlg)
	m_banlength = 1.0f;
	m_banwidth = 0.5f;
	m_hdotnumber = 10;
	m_vdotnumber = 5;
	m_ax = 0.2;
	m_ay = 0.2;
	m_banheigth = 0.01;
	m_e =2.1e11;
	m_u = 0.3;
	m_alpha = 1e20;
	m_f =1000.0;
	//}}AFX_DATA_INIT
	// Note that LoadIcon does not require a subsequent DestroyIcon in Win32
	HI=m_banlength/(2*(m_hdotnumber-1));
	HJ=m_banwidth/(2*(m_vdotnumber-1));
	m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}

void CRKPMDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CRKPMDlg)
	DDX_Text(pDX, IDC_banlength, m_banlength);
	DDX_Text(pDX, IDC_banwidth, m_banwidth);
	DDX_Text(pDX, IDC_hdotnumber, m_hdotnumber);
	DDX_Text(pDX, IDC_vdotnumber, m_vdotnumber);
	DDX_Text(pDX, IDC_x_expand, m_ax);
	DDX_Text(pDX, IDC_y_expand, m_ay);
	DDX_Text(pDX, IDC_heigth, m_banheigth);
	DDX_Text(pDX, IDC_E, m_e);
	DDX_Text(pDX, IDC_u, m_u);
	DDX_Text(pDX, IDC_alpha, m_alpha);
	DDX_Text(pDX, IDC_FORCE, m_f);
	//}}AFX_DATA_MAP
	NP=m_hdotnumber*m_vdotnumber;
}

BEGIN_MESSAGE_MAP(CRKPMDlg, CDialog)
	//{{AFX_MSG_MAP(CRKPMDlg)
	ON_WM_SYSCOMMAND()
	ON_WM_PAINT()
	ON_WM_QUERYDRAGICON()
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CRKPMDlg message handlers

BOOL CRKPMDlg::OnInitDialog()
{
	CDialog::OnInitDialog();

	// Add "About..." menu item to system menu.

	// IDM_ABOUTBOX must be in the system command range.
	ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);
	ASSERT(IDM_ABOUTBOX < 0xF000);

	CMenu* pSysMenu = GetSystemMenu(FALSE);
	if (pSysMenu != NULL)
	{
		CString strAboutMenu;
		strAboutMenu.LoadString(IDS_ABOUTBOX);
		if (!strAboutMenu.IsEmpty())
		{
			pSysMenu->AppendMenu(MF_SEPARATOR);
			pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
		}
	}

	// Set the icon for this dialog.  The framework does this automatically
	//  when the application's main window is not a dialog
	SetIcon(m_hIcon, TRUE);			// Set big icon
	SetIcon(m_hIcon, FALSE);		// Set small icon
	
	// TODO: Add extra initialization here
	
	return TRUE;  // return TRUE  unless you set the focus to a control
}

void CRKPMDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
	if ((nID & 0xFFF0) == IDM_ABOUTBOX)
	{
		CAboutDlg dlgAbout;
		dlgAbout.DoModal();
	}
	else
	{
		CDialog::OnSysCommand(nID, lParam);
	}
}

// If you add a minimize button to your dialog, you will need the code below
//  to draw the icon.  For MFC applications using the document/view model,
//  this is automatically done for you by the framework.

void CRKPMDlg::OnPaint() 
{
	if (IsIconic())
	{
		CPaintDC dc(this); // device context for painting

		SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);

		// Center icon in client rectangle
		int cxIcon = GetSystemMetrics(SM_CXICON);
		int cyIcon = GetSystemMetrics(SM_CYICON);
		CRect rect;
		GetClientRect(&rect);
		int x = (rect.Width() - cxIcon + 1) / 2;
		int y = (rect.Height() - cyIcon + 1) / 2;

		// Draw the icon
		dc.DrawIcon(x, y, m_hIcon);
	}
	else
	{
		CDialog::OnPaint();
	}
}

// The system calls this to obtain the cursor to display while the user drags
//  the minimized window.
HCURSOR CRKPMDlg::OnQueryDragIcon()
{
	return (HCURSOR) m_hIcon;
}

void CRKPMDlg::setcordinate()
{
   for(int i=0;i<=m_vdotnumber-1;i++)
  {
	 for(int j=0;j<=m_hdotnumber-1;j++)
	 {
         x[i*m_hdotnumber+j]=j*(m_banlength/(m_hdotnumber-1));  
	 }
  }
  for(int k=0;k<=m_vdotnumber;k++)
  {
	  for(int p=0;p<=m_hdotnumber;p++)
	  {
		  y[k*m_hdotnumber+p]=k*(m_banwidth/(m_vdotnumber-1));
	  }
  }
}

void CRKPMDlg::outputcordinate()
{
     ofstream outfile;
     outfile.open("D:\\Paper\\Program\\RKPM\\Result\\cordinate.txt");
	 outfile<<setiosflags(ios::left);
	 outfile<<setw(10)<<"节点号"<<setw(10)<<"x坐标"<<setw(10)<<"y坐标"<<endl;
	 for(int i=0;i<NP;i++)
	 {
		 outfile<<setw(10)<<i<<setw(10)<<x[i]<<setw(10)<<y[i]<<endl;
	 }
	 outfile.close();

	 outfile.open("D:\\Paper\\Program\\RKPM\\Result\\x.txt");
     outfile<<setiosflags(ios::left);
     for(int j=0;j<NP;j++)
	 {
		 outfile<<x[j]<<endl;
	 }
	 outfile.close();

	 outfile.open("D:\\Paper\\Program\\RKPM\\Result\\y.txt");
     outfile<<setiosflags(ios::left);
     for(int k=0;k<NP;k++)
	 {
		 outfile<<y[k]<<endl;
	 }
	 outfile.close();
}




void CRKPMDlg::setinsertfun()
{
   static double a[4]={-0.577350269189626,0.577350269189626,
	             0.577350269189626,-0.577350269189626};//积分点坐标
   static double b[4]={-0.577350269189626,-0.577350269189626,
	             0.577350269189626,0.577350269189626};//积分点坐标 
   for(int i=0;i<4;i++)
	{
		N1[i]=(1-a[i])*(1-b[i])/4;
	    N2[i]=(1+a[i])*(1-b[i])/4;
        N3[i]=(1+a[i])*(1+b[i])/4;
        N4[i]=(1-a[i])*(1+b[i])/4;
	}
   ofstream N;
		N.open("D:\\Paper\\Program\\RKPM\\Result\\n.txt");
		for(int n=0;n<4;n++)
			N<<setw(10)<<N1[n]<<setw(10)<<N2[n]<<setw(10)<<N3[n]<<setw(10)<<N4[n]<<endl;
		N.close();
}

void CRKPMDlg::setkernal(double xi,double yi,double xj,double yj)
{
   double s1,s2;
   s1=fabs(xi-xj)/m_ax;
   s2=fabs(yi-yj)/m_ay;
   if(s1<=1 && s2<=1)
	   ker=(1-6*pow(s1,2)+8*pow(s1,3)-3*pow(s1,4))*(1-6*pow(s2,2)+8*pow(s2,3)-3*pow(s2,4));
   else
	   ker=0;
} 

void CRKPMDlg::setkernalofdx(double xi,double yi,double xj,double yj)
{
   double s1,s2;
   s1=fabs(xi-xj)/m_ax;
   s2=fabs(yi-yj)/m_ay;
   if(s1<=1 && s2<=1)
	   kerofdx=(-12*s1+24*pow(s1,2)-12*pow(s1,3))*(1-6*pow(s2,2)+8*pow(s2,3)-3*pow(s2,4))/(m_ax*m_ax*m_ax);
   else
	   kerofdx=0;
} 

void CRKPMDlg::setkernalofd2x(double xi,double yi,double xj,double yj)
{
   double s1,s2;
   s1=fabs(xi-xj)/m_ax;
   s2=fabs(yi-yj)/m_ay;
   if(s1<=1 && s2<=1)
	   kerofd2x=(-12+48*s1-36*pow(s1,2))*(1-6*pow(s2,2)+8*pow(s2,3)-3*pow(s2,4))/(pow(m_ax,6));
   else
	   kerofd2x=0;
}

void CRKPMDlg::setkernalofdy(double xi,double yi,double xj,double yj)
{
   double s1,s2;
   s1=fabs(xi-xj)/m_ax;
   s2=fabs(yi-yj)/m_ay;
   if(s1<=1 && s2<=1)
	   kerofdy=(1-6*pow(s1,2)+8*pow(s1,3)-3*pow(s1,4))*(-12*s2+24*pow(s2,2)-12*pow(s2,3))/(pow(m_ay,3));
   else
	   kerofdy=0;
}

void CRKPMDlg::setkernalofd2y(double xi, double yi, double xj, double yj)
{
   double s1,s2;
   s1=fabs(xi-xj)/m_ax;
   s2=fabs(yi-yj)/m_ay;
   if(s1<=1 && s2<=1)
	   kerofd2y=(1-6*pow(s1,2)+8*pow(s1,3)-3*pow(s1,4))*(-12+48*s2-36*pow(s2,2))/(pow(m_ay,6));
   else
	   kerofd2y=0;
}


void CRKPMDlg::setkernalofdxdy(double xi, double yi, double xj, double yj)
{
   double s1,s2;
   s1=fabs(xi-xj)/m_ax;
   s2=fabs(yi-yj)/m_ay;
   if(s1<=1 && s2<=1)
	   kerofdxdy=(-12*s1+24*pow(s1,2)-12*pow(s1,3))*(-12*s2+24*pow(s2,2)-12*pow(s2,3))/(pow(m_ax,3)*pow(m_ay,3));
   else
	   kerofdxdy=0; 
}

void CRKPMDlg::setgaosicordinate()
{
  setinsertfun();
  setcordinate();
  //参照公式推导（5.5）
	for(int ni=0;ni<(m_vdotnumber-1);ni++)
		for(int nj=0;nj<(m_hdotnumber-1);nj++)
			for(int nk=0;nk<4;nk++)
			{
                gx[4*(m_hdotnumber-1)*ni+4*nj+nk]=N1[nk]*x[ni*m_hdotnumber+nj]+
					                              N2[nk]*x[ni*m_hdotnumber+nj+1]+
											      N3[nk]*x[(ni+1)*m_hdotnumber+nj+1]+
											      N4[nk]*x[(ni+1)*m_hdotnumber+nj];
                gy[4*(m_hdotnumber-1)*ni+4*nj+nk]=N1[nk]*y[ni*m_hdotnumber+nj]+
					                              N2[nk]*y[ni*m_hdotnumber+nj+1]+
											      N3[nk]*y[(ni+1)*m_hdotnumber+nj+1]+
											      N4[nk]*y[(ni+1)*m_hdotnumber+nj];
			}
	ofstream outfile1;
	outfile1.open("D:\\Paper\\Program\\RKPM\\Result\\gaosi.txt");
    outfile1<<setiosflags(ios::left);
	for(int pi=0;pi<(4*(m_hdotnumber-1)*(m_vdotnumber-1));pi++)
		outfile1<<setw(20)<<pi<<setw(20)<<gx[pi]<<setw(20)<<gy[pi]<<endl;

	outfile1.close();
}


//计算雅克比矩阵，参照公式计算（5.13）
void CRKPMDlg::setya()
{
    setgaosicordinate();
	static double a[4]={-0.57735,0.57735,0.57735,-0.57735};//积分点坐标
    static double b[4]={-0.57735,-0.57735,0.57735,0.57735};//积分点坐标
	for(int ni=0;ni<(m_vdotnumber-1);ni++)
		for(int nj=0;nj<(m_hdotnumber-1);nj++)
			for(int nk=0;nk<4;nk++)
			{
		       ya[4*(m_hdotnumber-1)*ni+4*nj+nk]=(((-1+b[nk])*x[ni*m_hdotnumber+nj]+(1-b[nk])*x[ni*m_hdotnumber+nj+1]+(1+b[nk])*x[(ni+1)*m_hdotnumber+nj+1]+(-1-b[nk])*x[(ni+1)*m_hdotnumber+nj])*
			                                      ((-1+a[nk])*y[ni*m_hdotnumber+nj]+(-1-a[nk])*y[ni*m_hdotnumber+nj+1]+(1+a[nk])*y[(ni+1)*m_hdotnumber+nj+1]+(1-a[nk])*y[(ni+1)*m_hdotnumber+nj])-