volume_and_surface_area.cpp

算法的一些集合

#include "vs.h"
static const double PI = 3.141592654;
static const double deg = PI/180.0;
void sphere(double* X, double& Z, double phi, double theta) {
	X[0] = sin(phi*deg)*cos(theta*deg);
   X[1] = sin(phi*deg)*sin(theta*deg);
   Z = cos(phi*deg);
}

int main() {
   Quadrature qp(2, 9);  // 2-dimension, 3x3 Gaussian integration
   H1 ZAI(2, (double*)0, qp),                 // Natrual Coordinates
      N = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE( // Shape Functions
          "int, int, Quadrature", 9/*nen*/, 2/*nsd*/, qp),
      Zai, Eta;
	Zai &= ZAI[0]; Eta &= ZAI[1];
   // initial four corner nodes
   N[0] = (1-Zai)*(1-Eta)/4; N[1] = (1+Zai)*(1-Eta)/4;
   N[2] = (1+Zai)*(1+Eta)/4; N[3] = (1-Zai)*(1+Eta)/4;
   // add ceter node
   N[8] = (1-Zai.pow(2))*(1-Eta.pow(2));
   // modification to four corner nodes due to the presence of the center node
   N[0] -= N[8]/4; N[1] -= N[8]/4; N[2] -= N[8]/4; N[3] -= N[8]/4;
   // add four edge nodes
   N[4] = ((1-Zai.pow(2))*(1-Eta)-N[8])/2; N[5] = ((1-Eta.pow(2))*(1+Zai)-N[8])/2;
   N[6] = ((1-Zai.pow(2))*(1+Eta)-N[8])/2; N[7] = ((1-Eta.pow(2))*(1-Zai)-N[8])/2;
   // modification to four corner nodes due to the presence of the four edge nodes
   N[0] -= (N[4]+N[7])/2; N[1] -= (N[4]+N[5])/2;
   N[2] -= (N[5]+N[6])/2; N[3] -= (N[6]+N[7])/2;
#if defined(__THREE_ELEMENTS)
	const int ELEMENT_NO = 3;
   const double theta = 90.0; const double phi = 90.0;
   double X[ELEMENT_NO][9][2], Z[ELEMENT_NO][9];
   sphere(X[0][0], Z[0][0], 0.0, 0.0);                   sphere(X[0][1], Z[0][1], 1.0/2.0*phi, 0.0);
   sphere(X[0][2], Z[0][2], 1.0/2.0*phi, 1.0/2.0*theta); sphere(X[0][3], Z[0][3], 1.0/2.0*phi, theta);
   sphere(X[0][4], Z[0][4], 1.0/4.0*phi, 0.0);           sphere(X[0][5], Z[0][5], 1.0/2.0*phi, 1.0/4.0*theta);
   sphere(X[0][6], Z[0][6], 1.0/2.0*phi, 3.0/4.0*theta); sphere(X[0][7], Z[0][7], 1.0/4.0*phi, theta);
   sphere(X[0][8], Z[0][8], 1.0/4.0*phi, 1.0/2.0*theta);
   sphere(X[1][0], Z[1][0], 1.0/2.0*phi, 0.0);           sphere(X[1][1], Z[1][1], phi, 0.0);
   sphere(X[1][2], Z[1][2], phi, 1.0/2.0*theta);         sphere(X[1][3], Z[1][3], 1.0/2.0*phi, 1.0/2.0*theta);
   sphere(X[1][4], Z[1][4], 3.0/4.0*phi, 0.0);           sphere(X[1][5], Z[1][5], phi, 1.0/4.0*theta);
   sphere(X[1][6], Z[1][6], 3.0/4.0*phi, 1.0/2.0*theta); sphere(X[1][7], Z[1][7], 1.0/2.0*phi, 1.0/4.0*theta);
   sphere(X[1][8], Z[1][8], 3.0/4.0*phi, 1.0/4.0*theta);
   sphere(X[2][0], Z[2][0], 1.0/2.0*phi, theta);         sphere(X[2][1], Z[2][1], 1.0/2.0*phi, 1.0/2.0*theta);
   sphere(X[2][2], Z[2][2], phi, 1.0/2.0*theta);         sphere(X[2][3], Z[2][3], phi, theta);
   sphere(X[2][4], Z[2][4], 1.0/2.0*phi, 3.0/4.0*theta); sphere(X[2][5], Z[2][5], 3.0/4.0*phi, 1.0/2.0*theta);
   sphere(X[2][6], Z[2][6], phi, 3.0/4.0*theta);        sphere(X[2][7], Z[2][7], 3.0/4.0*phi, theta);
   sphere(X[2][8], Z[2][8], 3.0/4.0*phi, 3.0/4.0*theta);
#else
	const int ELEMENT_NO = 1;
   const double theta = 90.0; const double phi = 90.0;
   double X[ELEMENT_NO][9][2], Z[ELEMENT_NO][9];
   sphere(X[0][0], Z[0][0], 0.0, 0.0);                   sphere(X[0][1], Z[0][1], phi, 0.0);
   sphere(X[0][2], Z[0][2], phi, 1.0/2.0*theta);         sphere(X[0][3], Z[0][3], phi, theta);
   sphere(X[0][4], Z[0][4], 1.0/2.0*phi, 0.0);           sphere(X[0][5], Z[0][5], phi, 1.0/4.0*theta);
   sphere(X[0][6], Z[0][6], phi, 3.0/4.0*theta);         sphere(X[0][7], Z[0][7], 1.0/2.0*phi, theta);
   sphere(X[0][8], Z[0][8], 1.0/2.0*phi, 1.0/2.0*theta);
#endif             
   C0 vol(0.0), area(0.0);
   for(int i = 0; i < ELEMENT_NO; i++) {
   	//  coordinate transformation
   	C0 X_(9, 2, X[i][0]);
   	H1 x = N*X_;
      C0 Z_(9, Z[i]);      
#if defined(__SURFACE_AREA)
		H1 z = N*Z_;
      H0 dz_dx = d(z) * d(x).inverse(),
         da = sqrt((dz_dx[0]).pow(2)+(dz_dx[1]).pow(2)+1.0);
      area += da | J(d(x).det());
	}
   cout << (8.0*90.0/theta*area) << endl;
#else
#if defined(__ANALYTICAL_GEOMETRY)
      H0 z = ((H0)N)*Z_;
#else
      H0 z = sqrt(1-((H0)x[0]).pow(2)-((H0)x[1]).pow(2));
#endif
      vol +=  z | J(d(x).det());
   }
	cout << (8.0*90.0/theta*vol) << endl;
#endif
	return 0;
}