dirichlet.cpp

算法的一些集合

#include "vs.h"
#define PI 3.141592654

int main() {
	const int n = 16;
//------------------------------------------------------------------------------
// Example 6-1-1, Dirichlet boundary condition, Reddy, 1986, p. 265
//           -d^2 u
//           -------  = cos PI x, 0 < x < 1, u(0) = u(1) = 0
//            d x^2
//           u_exact = (cos(PI x) + 2x -1) / PI^2
//           u_n = sum(j = 1, N) c_j * sin(j PI x)
//------------------------------------------------------------------------------
	// A. Extended Simpson's Integration Formula
	double w[35] = {1.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0,
                            4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0,
                            4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0,
                            4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0, 4.0/3.0, 2.0/3.0,
                            4.0/3.0, 1.0/3.0};
   Quadrature qp(w, 0.0, 1.0, 35);
   J d_l(1.0/34.0);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      phi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
                                      n/*vector size*/, 1/*spatial dim.*/, qp);
   for(int i = 0; i < n; i++) phi[i] = sin(((double)i+1.0)*PI*x);
   // C. Variational Formulation
   H0 d_phi_2(n, (double*)0, qp);              // LHS bilinear form
   for(i = 0; i < n; i++) d_phi_2[i] = d(phi[i]).pow(2);
   C0 M = d_phi_2 | d_l;
   C0 b = ( ((H0)phi) * cos(PI*((H0)x)) ) | d_l; // RHS forced vector

   // D. Solutions--the Ritz Coefficients
   C0 c(n, (double*)NULL);
   for(i = 0; i < n; i++) c[i] = b[i] / M[i];
   //double semi[8]; // coefficients of semi-analytical Rayleight-Ritz approximation
   //semi[0] = semi[2] = semi[4] = semi[6] = 0.0;
   //for(int i = 1; i < 8; i=i+2)
   //	semi[i] = 4.0/( pow(PI,3.0)*((double)(i+1))*(pow(((double)(i+1)),2.0)-1.0) );
   for(i = 0; i < n; i++) cout << c[i] << ", "; cout << endl;
   //for(int i = 0; i < 8; i++) cout << semi[i] << ", "; cout << endl;
   //for(int i = 0; i < 8; i=i+2)                             // 0, 2, 4, 6
   //	assert(fabs( ((double)c[i]) - semi[i] ) < 1.e-4);
   //for(int i = 1; i < 8; i=i+2)                             // 1, 3, 5, 7
   //	assert(fabs( ( ((double)c[i]) - semi[i] ) / semi[i] ) < 0.01);
	return 0;
}