p5tri.h

finite element library for mathematic majored research

    const double basisvalue15 = 5.744562647*psitilde_a_5*scalings_y_5*psitilde_bs_5_0;
    const double basisvalue16 = 5.196152423*psitilde_a_4*scalings_y_4*psitilde_bs_4_1;
    const double basisvalue17 = 4.582575695*psitilde_a_3*scalings_y_3*psitilde_bs_3_2;
    const double basisvalue18 = 3.872983346*psitilde_a_2*scalings_y_2*psitilde_bs_2_3;
    const double basisvalue19 = 3*psitilde_a_1*scalings_y_1*psitilde_bs_1_4;
    const double basisvalue20 = 1.732050808*psitilde_a_0*scalings_y_0*psitilde_bs_0_5;
    
    // Table(s) of coefficients
    const static double coefficients0[21][21] = \
    {{0.0154328861, -0.009450674049, -0.005456349206, 0.03018753073, 0.02338316075, 0.01350027416, -0.02319976058, -0.01960737637, -0.01518780842, -0.008768685283, 0.02138986513, 0.0188640879, 0.01594306415, 0.01234944439, 0.007129955042, -0.01798917576, -0.01627182171, -0.01435039787, -0.01212829982, -0.009394540645, -0.00542394057},
    {0.0154328861, 0.009450674049, -0.005456349206, 0.03018753073, -0.02338316075, 0.01350027416, 0.02319976058, -0.01960737637, 0.01518780842, -0.008768685283, 0.02138986513, -0.0188640879, 0.01594306415, -0.01234944439, 0.007129955042, 0.01798917576, -0.01627182171, 0.01435039787, -0.01212829982, 0.009394540645, -0.00542394057},
    {0.0154328861, 0, 0.01091269841, 0, 0, 0.04050082247, 0, 0, 0, 0.03507474113, 0, 0, 0, 0, 0.03564977521, 0, 0, 0, 0, 0, 0.03254364342},
    {0.03507474113, 0.1260089873, -0.08432539683, 0.1018829162, -0.02864437192, 0.002025041123, 0.07733253528, -0.02614316849, 0, 0.005845790188, 0.03564977521, 0.0125760586, -0.0385290717, 0.04425217572, -0.02970814601, 0, 0.03254364342, -0.05166143235, 0.05821583914, -0.05260942761, 0.03254364342},
    {0.03507474113, -0.03293416714, 0.1215277778, -0.01132032402, 0.07541069343, -0.02227545236, 0, 0.1045726739, -0.1012520562, 0.05845790188, 0, 0.0565922637, -0.009300120754, -0.04013569426, 0.04159140441, 0, 0, 0.05166143235, -0.1018777185, 0.1183712121, -0.08135910855},
    {0.03507474113, 0.08877905925, -0.08928571429, 0.03018753073, 0.04325884739, -0.05940120628, 0, 0.02614316849, 0.1012520562, -0.1169158038, 0, 0, 0.07440096603, -0.03087361097, -0.0118832584, 0, 0, 0, 0.06791847899, -0.131523569, 0.1084788114},
    {0.03507474113, -0.01002344217, 0.1512896825, 0, 0.05027379562, 0.09315189167, 0, 0, 0.05062602808, 0.06430369207, 0, 0, 0, 0.0720384256, -0.02376651681, 0, 0, 0, 0, 0.06576178451, -0.08135910855},
    {0.03507474113, -0.1260089873, -0.08432539683, 0.1018829162, 0.02864437192, 0.002025041123, -0.07733253528, -0.02614316849, 0, 0.005845790188, 0.03564977521, -0.0125760586, -0.0385290717, -0.04425217572, -0.02970814601, 0, 0.03254364342, 0.05166143235, 0.05821583914, 0.05260942761, 0.03254364342},
    {0.03507474113, 0.03293416714, 0.1215277778, -0.01132032402, -0.07541069343, -0.02227545236, 0, 0.1045726739, 0.1012520562, 0.05845790188, 0, -0.0565922637, -0.009300120754, 0.04013569426, 0.04159140441, 0, 0, -0.05166143235, -0.1018777185, -0.1183712121, -0.08135910855},
    {0.03507474113, -0.08877905925, -0.08928571429, 0.03018753073, -0.04325884739, -0.05940120628, 0, 0.02614316849, -0.1012520562, -0.1169158038, 0, 0, 0.07440096603, 0.03087361097, -0.0118832584, 0, 0, 0, 0.06791847899, 0.131523569, 0.1084788114},
    {0.03507474113, 0.01002344217, 0.1512896825, 0, -0.05027379562, 0.09315189167, 0, 0, -0.05062602808, 0.06430369207, 0, 0, 0, -0.0720384256, -0.02376651681, 0, 0, 0, 0, -0.06576178451, -0.08135910855},
    {0.03507474113, -0.1360324295, -0.06696428571, 0.03697972514, 0.07891816754, 0.06007621999, -0.007733253528, -0.04575054485, -0.05568863089, -0.03799763622, -0.06416959538, -0.0125760586, 0.01727165283, 0.02778624987, 0.02020153929, 0.08994587878, 0.04881546513, 0.02009055702, 0.002425659964, -0.005636724387, -0.00542394057},
    {0.03507474113, 0.1217132264, -0.03224206349, -0.06716725587, -0.03215184604, 0.02767556202, 0.1391985635, 0.06535792122, -0.01012520562, -0.02922895094, 0.04277973025, -0.0565922637, -0.03321471698, 0.009262083291, 0.02020153929, -0.1798917576, -0.03254364342, 0.0114803183, 0.009702639856, -0.001878908129, -0.00542394057},
    {0.03507474113, -0.1217132264, -0.03224206349, -0.06716725587, 0.03215184604, 0.02767556202, -0.1391985635, 0.06535792122, 0.01012520562, -0.02922895094, 0.04277973025, 0.0565922637, -0.03321471698, -0.009262083291, 0.02020153929, 0.1798917576, -0.03254364342, -0.0114803183, 0.009702639856, 0.001878908129, -0.00542394057},
    {0.03507474113, 0.1360324295, -0.06696428571, 0.03697972514, -0.07891816754, 0.06007621999, 0.007733253528, -0.04575054485, 0.05568863089, -0.03799763622, -0.06416959538, 0.0125760586, 0.01727165283, -0.02778624987, 0.02020153929, -0.08994587878, 0.04881546513, -0.02009055702, 0.002425659964, 0.005636724387, -0.00542394057},
    {0.280597929, -0.2147880466, -0.1240079365, 0.09056259218, 0.1169158038, 0, 0.1546650706, 0.07842950546, 0.1012520562, 0.08184106264, -0.1425991008, 0.0251521172, -0.0318861283, -0.107028518, -0.09506606722, 0, -0.1301745737, -0.1033228647, -0.02910791957, 0.0263047138, 0.03254364342},
    {0.03507474113, 0, -0.1240079365, -0.3848910168, 0, 0.1012520562, 0, -0.1045726739, 0, 0.04676632151, 0.2138986513, 0, 0.1408304, 0, -0.09506606722, 0, 0.1952618605, 0, -0.05821583914, 0, 0.03254364342},
    {0.280597929, 0.2147880466, -0.1240079365, 0.09056259218, -0.1169158038, 0, -0.1546650706, 0.07842950546, -0.1012520562, 0.08184106264, -0.1425991008, -0.0251521172, -0.0318861283, 0.107028518, -0.09506606722, 0, -0.1301745737, 0.1033228647, -0.02910791957, -0.0263047138, 0.03254364342},
    {0.03507474113, -0.1073940233, 0.06200396825, 0.01132032402, -0.3069039849, -0.2531301404, 0, -0.1045726739, 0, 0.04676632151, 0, 0.1697767911, 0.009300120754, 0.126581805, 0.1723072468, 0, 0, 0.154984297, 0.1018777185, -0.03945707071, -0.08135910855},
    {0.03507474113, 0.1073940233, 0.06200396825, 0.01132032402, 0.3069039849, -0.2531301404, 0, -0.1045726739, 0, 0.04676632151, 0, -0.1697767911, 0.009300120754, -0.126581805, 0.1723072468, 0, 0, -0.154984297, 0.1018777185, 0.03945707071, -0.08135910855},
    {0.280597929, 0, 0.248015873, -0.06037506145, 0, 0.1350027416, 0, -0.05228633697, 0, -0.2104484468, 0, 0, -0.1488019321, 0, -0.1425991008, 0, 0, 0, -0.135836958, 0, 0.1084788114}};
    
    // Extract relevant coefficients
    const double coeff0_0 = coefficients0[dof][0];
    const double coeff0_1 = coefficients0[dof][1];
    const double coeff0_2 = coefficients0[dof][2];
    const double coeff0_3 = coefficients0[dof][3];
    const double coeff0_4 = coefficients0[dof][4];
    const double coeff0_5 = coefficients0[dof][5];
    const double coeff0_6 = coefficients0[dof][6];
    const double coeff0_7 = coefficients0[dof][7];
    const double coeff0_8 = coefficients0[dof][8];
    const double coeff0_9 = coefficients0[dof][9];
    const double coeff0_10 = coefficients0[dof][10];
    const double coeff0_11 = coefficients0[dof][11];
    const double coeff0_12 = coefficients0[dof][12];
    const double coeff0_13 = coefficients0[dof][13];
    const double coeff0_14 = coefficients0[dof][14];
    const double coeff0_15 = coefficients0[dof][15];
    const double coeff0_16 = coefficients0[dof][16];
    const double coeff0_17 = coefficients0[dof][17];
    const double coeff0_18 = coefficients0[dof][18];
    const double coeff0_19 = coefficients0[dof][19];
    const double coeff0_20 = coefficients0[dof][20];
    
    // Compute value(s)
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14 + coeff0_15*basisvalue15 + coeff0_16*basisvalue16 + coeff0_17*basisvalue17 + coeff0_18*basisvalue18 + coeff0_19*basisvalue19 + coeff0_20*basisvalue20;
  }

  /// Evaluate order n derivatives of basis function i at given point in cell
  virtual void evaluate_basis_derivatives(unsigned int i,
                                          unsigned int n,
                                          double* values,
                                          const double* coordinates,
                                          const ufc::cell& c) const
  {
    // Extract vertex coordinates
    const double * const * element_coordinates = c.coordinates;
    
    // Compute Jacobian of affine map from reference cell
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
      
    // Compute determinant of Jacobian
    const double detJ = J_00*J_11 - J_01*J_10;
    
    // Get coordinates and map to the reference (UFC) element
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
                element_coordinates[0][0]*element_coordinates[2][1] +\
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
                element_coordinates[1][0]*element_coordinates[0][1] -\
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
    
    // Map coordinates to the reference square
    if (std::abs(y - 1.0) < 1e-09)
      x = -1.0;
    else
      x = 2.0 *x/(1.0 - y) - 1.0;
    y = 2.0*y - 1.0;
    
    // Compute number of derivatives
    unsigned int num_derivatives = 1;
    
    for (unsigned int j = 0; j < n; j++)
      num_derivatives *= 2;
    
    
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
    unsigned int **combinations = new unsigned int *[num_derivatives];
        
    for (unsigned int j = 0; j < num_derivatives; j++)
    {
      combinations[j] = new unsigned int [n];
      for (unsigned int k = 0; k < n; k++)
        combinations[j][k] = 0;
    }
        
    // Generate combinations of derivatives
    for (unsigned int row = 1; row < num_derivatives; row++)
    {
      for (unsigned int num = 0; num < row; num++)
      {
        for (unsigned int col = n-1; col+1 > 0; col--)
        {
          if (combinations[row][col] + 1 > 1)
            combinations[row][col] = 0;
          else
          {
            combinations[row][col] += 1;
            break;
          }
        }
      }
    }
    
    // Compute inverse of Jacobian
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
    
    // Declare transformation matrix
    // Declare pointer to two dimensional array and initialise
    double **transform = new double *[num_derivatives];
        
    for (unsigned int j = 0; j < num_derivatives; j++)
    {
      transform[j] = new double [num_derivatives];
      for (unsigned int k = 0; k < num_derivatives; k++)
        transform[j][k] = 1;
    }
    
    // Construct transformation matrix
    for (unsigned int row = 0; row < num_derivatives; row++)
    {
      for (unsigned int col = 0; col < num_derivatives; col++)
      {
        for (unsigned int k = 0; k < n; k++)
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
      }
    }
    
    // Reset values