boundary_element_method.cpp

算法的一些集合

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <iomanip.h>
#include "vs.h"
ofstream ofs("temp.cpp", ios::out | ios::trunc);
static const double PI = 3.141592654;
static C0 A(32, 32, (double*)0);
static C0 H(32, 16, A, 0, 0);
static C0 mG(32, 16, A, 0, 16);
static C0 f(32, (double*)0);
double w[9] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
               		 64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0},
				  x[8][2] = {{-1.0, -0.75}, {-0.75, -0.5}, {-0.5, -0.25}, {-0.25, -0.0},
           					 {0.0, 0.25}, {0.25, 0.5}, {0.5, 0.75}, {0.75, 1.0}},
				  y[8][2] = {{-1.0, -0.75}, {-0.75, -0.5}, {-0.5, -0.25}, {-0.25, -0.0},
           					 {0.0, 0.25}, {0.25, 0.5}, {0.5, 0.75}, {0.75, 1.0}},
				  zai[8] = {-0.875, -0.625, -0.375, -0.125, 0.125, 0.375, 0.625, 0.875},
				  eta[8] = {-0.875, -0.625, -0.375, -0.125, 0.125, 0.375, 0.625, 0.875};
              
void LHS_0_15() {
   // A. LHS
   // A.1 top and bottom nodes (node # 0-15); q = 0
   for(int i = 0; i < 16; i++) H[i][i] = -1.0/2.0;
   for(i = 0; i < 8; i++) {
   	for(int j = 0; j < 8; j++) {
   		Quadrature qp(w, x[j][0], x[j][1], 9);
      	J d_l(0.25/8.0);
      	H0 X(qp),
         	R_pow_2 = 4+(X-zai[i]).pow(2);
   		H[i][j+8] = H[i+8][j] = 1.0/(PI*R_pow_2) | d_l;
      }
   }
   for(i = 0; i < 8; i++) {
   	for(int j = 0; j < 8; j++) {
      	 Quadrature qp(w, y[j][0], y[j][1], 9);
          J d_l(0.25/8.0);
          H0 Y(qp),
            R0 = sqrt(pow(1.0-zai[i], 2)+(Y-1).pow(2)), R1 = sqrt(pow(1.0+zai[i], 2)+(Y-1).pow(2)),
            R2 = sqrt(pow(1.0-zai[i], 2)+(Y+1).pow(2)), R3 = sqrt(pow(1.0+zai[i], 2)+(Y+1).pow(2));
          mG[i][j] = 1.0/(2.0*PI)*log(R0) | d_l;   mG[i][j+8] = 1.0/(2.0*PI)*log(R1) | d_l;
          mG[i+8][j] = 1.0/(2.0*PI)*log(R2) | d_l; mG[i+8][j+8] = 1.0/(2.0*PI)*log(R3) | d_l;
      }
   }
}
void LHS_16_31() {
   // A.2: right-and-left nodes (node # 16-31); u(16-23) = 100; u(24-31) = 0
   double mG_diag = 0.125 / PI * (log(0.125)-1.0);
   for(int i = 0; i < 16; i++) mG[i+16][i] = mG_diag;
   for(i = 0; i < 8; i++) {
   	for(int j = 0; j < 8; j++) {
   		Quadrature qp(w, y[j][0], y[j][1], 9);
      	J d_l(0.25/8.0);
      	H0 Y(qp);
      	if(i != j) {
         	if(i<j) mG[i+16][j] = mG[i+24][j+8] = 1.0/(2.0*PI) * log(Y-eta[i]) | d_l;
            else mG[i+16][j] = mG[i+24][j+8] = 1.0/(2.0*PI) * log(eta[i]-Y) | d_l;
         }
      }
   }
   for(i = 0; i < 8; i++)
   for(int j = 0; j < 8; j++) {
   	Quadrature qp(w, y[j][0], y[j][1], 9);
      J d_l(0.25/8.0);
      H0 Y(qp),
      	R = sqrt(4+(Y-eta[i]).pow(2));
   	mG[i+16][j+8] = mG[i+24][j] = 1.0/(2.0*PI) * log(R) | d_l;
   }
   for(i = 0; i < 8; i++)
   for(int j = 0; j < 8; j++) {
   		Quadrature qp(w, x[j][0], x[j][1], 9);
      	J d_l(0.25/8.0);
      	H0 X(qp),
            R0_2 = pow(1-eta[i], 2) +(X-1).pow(2), R1_2 = pow(1+eta[i], 2) +(X-1).pow(2),
            R2_2 = pow(1-eta[i], 2) +(X+1).pow(2), R3_2 = pow(1+eta[i], 2) +(X+1).pow(2);
   		H[i+16][j] = (1.0-eta[i])/(2.0*PI*R0_2) | d_l; H[i+16][j+8] = (1.0+eta[i])/(2.0*PI*R1_2) | d_l;
         H[i+24][j] = (1.0-eta[i])/(2.0*PI*R2_2) | d_l; H[i+24][j+8] = (1.0+eta[i])/(2.0*PI*R3_2) | d_l;
   }
}
void RHS() {
	// B. Boundary conditions(RHS)
   for(int i = 0; i < 8; i++) {
   	for(int j = 0; j < 8; j++) {
   		Quadrature qp(w, y[j][0], y[j][1], 9);
         J d_l(0.25/8.0);
         H0 Y(qp),
            R0_2 = pow(1-zai[i], 2)+(Y-1).pow(2),
            R1_2 = pow(1-zai[i], 2)+(Y+1).pow(2);
         f[i] += (-100.0*(1-zai[i]))/(2.0*PI*R0_2) | d_l;
         f[i+8] += (-100.0*(1-zai[i]))/(2.0*PI*R1_2) | d_l;
      }
   }
   for(i = 16; i < 24; i++) f[i] = 100.0 / 2.0;
   for(i = 24; i < 32; i++) {
   	for(int j = 0; j < 8; j++) {
   		Quadrature qp(w, y[j][0], y[j][1], 9);
      	J d_l(0.25/8.0);
      	H0 Y(qp), R_pow_2 = 4+(Y-eta[i-24]).pow(2);
   		f[i] += - 100.0 / (PI*R_pow_2) | d_l;
      }
   }
}
void T_recovery(C0& T_gamma, C0& q_gamma, C0& T) {
	for(int i = 0; i < 8; i++)
   for(int j = 0; j < 8; j++) {
   	for(int k = 0; k < 8; k++) {
      	Quadrature qpx(w, x[k][0], x[k][1], 9);
      	J d_l(0.25/8.0);
      	H0 X(qpx),
      		R0_2 = pow((1-eta[i]), 2) + (X-zai[j]).pow(2),
         	R1_2 = pow((1+eta[i]), 2) + (X-zai[j]).pow(2);
      	T[i][j] += T_gamma[k] /(2.0*PI)*(1-eta[i])/R0_2 | d_l;
      	T[i][j] += T_gamma[k+8] /(2.0*PI)*(1+eta[i])/R1_2 | d_l;
      	Quadrature qpy(w, y[k][0], y[k][1], 9);
      	H0 Y(qpy),
         	R2_2 = pow((1-zai[j]), 2) + (Y-eta[i]).pow(2),
         	R3_2 = pow((zai[j]+1), 2) + (Y-eta[i]).pow(2);
      	T[i][j] += 100.0 / (2.0*PI)*(1-zai[j])/R2_2 | d_l;
      	T[i][j] += q_gamma[k] / (2.0*PI)*log(sqrt(R2_2)) | d_l;
      	T[i][j] += q_gamma[k+8] / (2.0*PI)*log(sqrt(R3_2)) | d_l;
      }
   }
}
void q_recovery(C0& T_gamma, C0& q_gamma, C0& q_x, C0& q_y) {
	double nx, ny;
   nx = 1.0;
	for(int i = 0; i < 8; i++)
   for(int j = 0; j < 8; j++) {
   	for(int k = 0; k < 8; k++) {
      	Quadrature qpx(w, x[k][0], x[k][1], 9);
      	J d_l(0.25/8.0);
      	H0 X(qpx),
      		R0_2 = pow((1-eta[i]), 2) + (X-zai[j]).pow(2),
         	R1_2 = pow((1+eta[i]), 2) + (X-zai[j]).pow(2);
         ny = 1.0;
      	q_x[i][j] -= (ny*T_gamma[k] /(2.0*PI*R0_2)*(2.0*(X-zai[j])*(1-eta[i])/R0_2)) | d_l;
      	q_y[i][j] -= (ny*T_gamma[k] /(2.0*PI*R0_2)*(2.0*(1-eta[i])*(1-eta[i])/R0_2-1.0)) | d_l;
         ny = -1.0;
      	q_x[i][j] -= (ny*T_gamma[k+8] /(2.0*PI*R1_2)*(2.0*(X-zai[j])*(-1-eta[i])/R1_2)) | d_l;
	     	q_y[i][j] -= (ny*T_gamma[k+8] /(2.0*PI*R1_2)*(2.0*(-1-eta[i])*(-1-eta[i])/R1_2-1.0)) | d_l;

      	Quadrature qpy(w, y[k][0], y[k][1], 9);
      	H0 Y(qpy),
         	R2_2 = pow((1-zai[j]), 2) + (Y-eta[i]).pow(2),
         	R3_2 = pow((zai[j]+1), 2) + (Y-eta[i]).pow(2);
      	q_x[i][j] -= (nx*100.0 / (2.0*PI*R2_2)*(2.0*(1-zai[j])*(1-zai[j])/R2_2 - 1.0)) | d_l;
      	q_y[i][j] -= (nx*100.0 / (2.0*PI*R2_2)*(2.0*(Y-eta[i])*(1-zai[j])/R2_2 )) | d_l;
      	q_x[i][j] += (q_gamma[k] / (2.0*PI*R2_2)*(1-zai[j])) | d_l;
      	q_y[i][j] += (q_gamma[k] / (2.0*PI*R2_2)*(Y-eta[i])) | d_l;
      	q_x[i][j] += (q_gamma[k+8] / (2.0*PI*R3_2)*(-1-zai[j])) | d_l;
      	q_y[i][j] += (q_gamma[k+8] / (2.0*PI*R3_2)*(Y-eta[i])) | d_l;
      }
   }
}
int main() {
	LHS_0_15();
	LHS_16_31();
	RHS();
	// C. Matrix solver
	C0 Y = f / A, T_gamma(16, Y, 0), q_gamma(16, Y, 16),
      T(8, 8, (double*)0), q_x(8, 8, (double*)0), q_y(8, 8, (double*)0);
   ofs << T_gamma << endl;
   ofs << q_gamma << endl;
	T_recovery(T_gamma, q_gamma, T);
   ofs << T << endl;
	/*C0 T_gamma(16, (double*)0), q_gamma(16, (double*)0), q_x(8, 8, (double*)0), q_y(8, 8, (double*)0);
   for(int i = 0; i < 8; i++) {
   	T_gamma[i+8] = T_gamma[i] = 100.0/16.0*((double)(i*2+1));
      q_gamma[i+8] = -(q_gamma[i] = -50.0);
   }
   ofs << T_gamma << endl;
   ofs << q_gamma << endl;*/
   q_recovery(T_gamma, q_gamma, q_x, q_y);
   ofs << q_x << endl;
   ofs << q_y << endl;
   ofs.close();
	return 0;
}