lagrange_multiplier.cpp

算法的一些集合

#include "vs.h"
#define L_ 1.0
#define E_ 1.0
#define I_ 1.0
#define f_ 1.0
#define M_ 1.0

int main() {
	{ // One parameter approximation--1. define the matrix element by element
	// A. two-point Trapezoidal Integration Formula
	double weight[2] = {0.5, 0.5};
   Quadrature qp(weight, 0.0, L_, 2);
   J d_l(L_);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      psi = x, w = x;
   double lambda = 1.0;
   // C. Variational Formulation
   H0 b = INTEGRABLE_MATRIX("int, int, Quadrature", 3, 3, qp);
   b[0][0] = (E_*I_)*d(psi)*d(psi);
   b[0][2] = b[2][0] = ((H0)psi)*lambda;
   b[1][2] = b[2][1] = d(w)*lambda;
   b[0][1] = b[1][0] =  b[1][1] = b[2][2] = 0.0;
   C0 B = b | d_l;                // LSH bilinear form
   C0 l(3, (double*)NULL);          // RHS vector
   l[0] = -L_*M_; l[1] = ( ((H0)w)*f_ ) | d_l; l[2] = 0.0;

   // D. Solutions--the Ritz Coefficients
   C0 c = l / B;
   for(int i = 0; i < 3; i++) cout << c[i] << ", "; cout << endl;
   }
	{ // One parameter approximation--2. use matrix concatenation
	// A. two-point Trapezoidal Integration Formula
	double weight[2] = {0.5, 0.5};
   Quadrature qp(weight, 0.0, L_, 2);
   J d_l(L_);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      psi = x, w = x;
   H0 null(qp);
   null = 0.0;
   double lambda = 1.0;
   // C. Variational Formulation
   C0 B = ( ( (E_*I_*(d(psi)*d(psi)))|     null      | (((H0)psi)*lambda) )&
            (     null               |     null      | (d(w)*lambda)      )&
            ( (((H0)psi)*lambda)     | (d(w)*lambda) |     null           )
          ) | d_l;                                 // LSH bilinear form
   C0 l = ( C0(-L_*M_) &
            ( (((H0)w)*f_) | d_l ) &
            C0(0.0)
          );    // RHS vector

   // D. Solutions--the Ritz Coefficients
   C0 c = l / B;
   for(int i = 0; i < 3; i++) cout << c[i] << ", "; cout << endl;
   }
	{ // One parameter approximation--3. use Cartesian basis expression
	// A. two-point Trapezoidal Integration Formula
	double weight[2] = {0.5, 0.5};
   Quadrature qp(weight, 0.0, L_, 2);
   J d_l(L_);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      psi = x, w = x;
   double lambda = 1.0;
   // C. Variational Formulation
   C0 e(3);
   C0 B = ( (E_*I_*(d(psi)*d(psi)))*(e[0]%e[0])+
            0.0                    *(e[0]%e[1])+
            (((H0)psi)*lambda)     *(e[0]%e[2])+
            0.0                    *(e[1]%e[0])+
            0.0                    *(e[1]%e[1])+
            (d(w)*lambda)          *(e[1]%e[2])+
            (((H0)psi)*lambda)     *(e[2]%e[0])+
            (d(w)*lambda)          *(e[2]%e[1])+
            0.0                    *(e[2]%e[2])
          ) | d_l;                                     // LSH bilinear form
   C0 l = ((-L_*M_)*e[0] + ((((H0)w)*f_) | d_l)*e[1] + 0.0*e[2]); // RHS vector

   // D. Solutions--the Ritz Coefficients
   C0 c = l / B;
   for(int i = 0; i < 3; i++) cout << c[i] << ", "; cout << endl;
   }
   { // Two parameter approximation--1. use matrix concatenation for the formulation
	// A. Simpson's Integration Formula
	double weight[3] = {1.0/3.0, 4.0/3.0, 1.0/3.0};
   Quadrature qp(weight, 0.0, L_, 3);
   J d_l(L_/2.0);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      psi    = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      w      = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      lambda = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp);
   psi[0]    = x;   psi[1]    = x.pow(2);  // psi
   w[0]      = x;   w[1]      = x.pow(2);  // w
   lambda[0] = 1.0; lambda[1] = x;         // lambda
   
   // C. Weak Formulation
   H0 null = INTEGRABLE_MATRIX("int, int, Quadrature", 2, 2, qp); null = 0.0;
   C0 B =(( (E_*I_*(d(psi)*(~d(psi))))|        null            | (((H0)psi)%((H0)lambda)) )&
          (      null                 |        null            | (d(w)(0)%((H0)lambda))   )&
          ( (((H0)lambda)%((H0)psi))  | (((H0)lambda)%d(w)(0)) |          null            )
         ) | d_l;
   C0 M_delta_psi(2, (double*)NULL), zero(2, (double*)NULL);
   M_delta_psi[0] = -M_*L_; M_delta_psi[1] = -M_*L_*L_; zero = 0.0;
   C0 l = ( M_delta_psi & ( (((H0)w)*f_) | d_l) & zero);

   // D. Solutions--the Ritz Coefficients
   C0 c = l / B;
   for(int i = 0; i < 6; i++) cout << c[i] << ", "; cout << endl;
   }
   { // Two parameter approximation--2. use referenced vector/matrix for the formulation
	// A. Simpson's Integration Formula
	double weight[3] = {1.0/3.0, 4.0/3.0, 1.0/3.0};
   Quadrature qp(weight, 0.0, L_, 3);
   J d_l(L_/2.0);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      psi    = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      w      = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      lambda = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp);
   psi[0]    = x;   psi[1]    = x.pow(2);  // psi
   w[0]      = x;   w[1]      = x.pow(2);  // w
   lambda[0] = 1.0; lambda[1] = x;         // lambda
   
   // C. Weak Formulation
   H0 b = INTEGRABLE_MATRIX("int, int, Quadrature", 6, 6, qp),
      b00 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 0, 0, qp), // referenced submatrix
      b01 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 0, 2, qp),
      b02 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 0, 4, qp),
      b10 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 2, 0, qp),
      b11 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 2, 2, qp),
      b12 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 2, 4, qp),
      b20 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 4, 0, qp),
      b21 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 4, 2, qp),
      b22 = INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, b, 4, 4, qp);
   b00 = (E_*I_)*d(psi)*(~d(psi));b01 = 0.0;                  b02 = ((H0)psi)%((H0)lambda);
   b10 = 0.0;                     b11 = 0.0;                  b12 = d(w)(0)%((H0)lambda);
   b20 = ((H0)lambda)%((H0)psi);  b21 = ((H0)lambda)%d(w)(0); b22 = 0.0;
   C0 B = b | d_l;                    // LSH bilinear form
   H0 F = INTEGRABLE_VECTOR("int, Quadrature", 6, qp),
      F0  = INTEGRABLE_VECTOR("int, H0&, int, const Quadrature&", 2, F, 0, qp),                  // referenced vector
      F1  = INTEGRABLE_VECTOR("int, H0&, int, const Quadrature&", 2, F, 2, qp),
      F2  = INTEGRABLE_VECTOR("int, H0&, int, const Quadrature&", 2, F, 4, qp);
   F1 = (((H0)w)*f_); F0 = F2 = 0.0;
   C0 l = F |d_l; // f vector size = 6
   l[0] = -M_*L_; l[1] = -M_*L_*L_;

   // D. Solutions--the Ritz Coefficients
   C0 c = l / B;
   for(int i = 0; i < 6; i++) cout << c[i] << ", "; cout << endl;
   }
   { // Two parameter approximation--3. use equal-partition submatrix/subvector for the formulation
	// A. Simpson's Integration Formula
	double weight[3] = {1.0/3.0, 4.0/3.0, 1.0/3.0};
   Quadrature qp(weight, 0.0, L_, 3);
   J d_l(L_/2.0);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      psi    = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      w      = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      lambda = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp);
   psi[0]    = x;   psi[1]    = x.pow(2);  // psi
   w[0]      = x;   w[1]      = x.pow(2);  // w
   lambda[0] = 1.0; lambda[1] = x;         // lambda
   
   // C. Weak Formulation
   H0 b = INTEGRABLE_MATRIX("int, int, Quadrature", 6, 6, qp),
      bs = INTEGRABLE_SUBMATRIX("int, int, H0&", 2, 2, b);
   bs(0,0) = E_*I_*(d(psi)*(~d(psi)));bs(0,1) = 0.0;                  bs(0,2) = ((H0)psi)%((H0)lambda);
   bs(1,0) = 0.0;                     bs(1,1) = 0.0;                  bs(1,2) = d(w)(0)%((H0)lambda);
   bs(2,0) = ((H0)lambda)%((H0)psi);  bs(2,1) = ((H0)lambda)%d(w)(0); bs(2,2) = 0.0;
   H0 F = INTEGRABLE_VECTOR("int, Quadrature", 6, qp),
      Fs = INTEGRABLE_SUBVECTOR("int, H0&", 2, F);
   Fs(1) = (((H0)w)*f_); Fs(0) = Fs(2) = 0.0;
   C0 B = b | d_l;                    // LSH bilinear form
   C0 l = F | d_l; // f vector size = 6
   l[0] = -M_*L_; l[1] = -M_*L_*L_;

   // D. Solutions--the Ritz Coefficients
   C0 c = l / B;
   for(int i = 0; i < 6; i++) cout << c[i] << ", "; cout << endl;
   }
   {  // Two parameter approximation--4. use Cartesian basis expression for the formulation of the LHS bilinear form
	// A. Simpson's Integration Formula
	double weight[3] = {1.0/3.0, 4.0/3.0, 1.0/3.0};
   Quadrature qp(weight, 0.0, L_, 3);
   J d_l(L_/2.0);      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      psi    = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      w      = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      lambda = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp);
   psi[0]    = x;   psi[1]    = x.pow(2);  // psi
   w[0]      = x;   w[1]      = x.pow(2);  // w
   lambda[0] = 1.0; lambda[1] = x;         // lambda

   // C. Weak Formulation
   C0 e(2), E(3);  // Cartesian basis
   C0 B =+( ((E_*I_)*(d(psi)*(~d(psi))))*((e%e)*(E[0]%E[0]))+
            0.0                         *((e%e)*(E[0]%E[1]))+
            (((H0)psi)%((H0)lambda))    *((e%e)*(E[0]%E[2]))+
            0.0                         *((e%e)*(E[1]%E[0]))+
            0.0                         *((e%e)*(E[1]%E[1]))+
            (d(w)(0)%((H0)lambda))      *((e%e)*(E[1]%E[2]))+
            (((H0)lambda)%((H0)psi))    *((e%e)*(E[2]%E[0]))+
            (((H0)lambda)%d(w)(0))      *((e%e)*(E[2]%E[1]))+
            0.0                         *((e%e)*(E[2]%E[2]))
           ) | d_l; // LSH bilinear form; use operator+() to convert to integrable_matrix_type
   C0 M_delta_psi(2, (double*)NULL); M_delta_psi[0] = -M_*L_; M_delta_psi[1] = -M_*L_*L_;
   C0 l = +( M_delta_psi             *(e*E[0]) +
             ( (((H0)w)*f_) | d_l )*(e*E[1]) +
             0.0                     *(e*E[2])
           );

   // D. Solutions--the Ritz Coefficients
   C0 c = l / B;
   for(int i = 0; i < 6; i++) cout << c[i] << ", "; cout << endl;
   }
   return 0;
}