// This code conforms with the UFC specification version 1.4 // and was automatically generated by FFC version 0.9.2+. // // This code was generated with the option '-l dolfin' and // contains DOLFIN-specific wrappers that depend on DOLFIN. // // This code was generated with the following parameters: // // cache_dir: '' // convert_exceptions_to_warnings: False // cpp_optimize: False // epsilon: 1e-14 // form_postfix: True // format: 'dolfin' // log_level: 10 // log_prefix: '' // optimize: True // output_dir: '.' // precision: 15 // quadrature_degree: 'auto' // quadrature_rule: 'auto' // representation: 'auto' // split: False #ifndef __ELASTICITY_H #define __ELASTICITY_H #include #include #include #include /// This class defines the interface for a finite element. class elasticity_finite_element_0: public ufc::finite_element { public: /// Constructor elasticity_finite_element_0() : ufc::finite_element() { // Do nothing } /// Destructor virtual ~elasticity_finite_element_0() { // Do nothing } /// Return a string identifying the finite element virtual const char* signature() const { return "FiniteElement('Discontinuous Lagrange', Cell('tetrahedron', 1, Space(3)), 0)"; } /// Return the cell shape virtual ufc::shape cell_shape() const { return ufc::tetrahedron; } /// Return the dimension of the finite element function space virtual unsigned int space_dimension() const { return 1; } /// Return the rank of the value space virtual unsigned int value_rank() const { return 0; } /// Return the dimension of the value space for axis i virtual unsigned int value_dimension(unsigned int i) const { return 1; } /// Evaluate basis function i at given point in cell virtual void evaluate_basis(unsigned int i, double* values, const double* coordinates, const ufc::cell& c) const { // Extract vertex coordinates // Compute Jacobian of affine map from reference cell // Compute sub determinants // Compute determinant of Jacobian // Compute inverse of Jacobian // Compute constants // Get coordinates and map to the reference (FIAT) element // Reset values. *values = 0.000000000000000; // Array of basisvalues. double basisvalues[1] = {0.000000000000000}; // Declare helper variables. // Compute basisvalues. basisvalues[0] = 1.000000000000000; // Table(s) of coefficients. static const double coefficients0[1] = \ {1.000000000000000}; // Compute value(s). for (unsigned int r = 0; r < 1; r++) { *values += coefficients0[r]*basisvalues[r]; }// end loop over 'r' } /// Evaluate all basis functions at given point in cell virtual void evaluate_basis_all(double* values, const double* coordinates, const ufc::cell& c) const { // Element is constant, calling evaluate_basis. evaluate_basis(0, values, coordinates, c); } /// 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 * x = c.coordinates; // Compute Jacobian of affine map from reference cell const double J_00 = x[1][0] - x[0][0]; const double J_01 = x[2][0] - x[0][0]; const double J_02 = x[3][0] - x[0][0]; const double J_10 = x[1][1] - x[0][1]; const double J_11 = x[2][1] - x[0][1]; const double J_12 = x[3][1] - x[0][1]; const double J_20 = x[1][2] - x[0][2]; const double J_21 = x[2][2] - x[0][2]; const double J_22 = x[3][2] - x[0][2]; // Compute sub determinants const double d_00 = J_11*J_22 - J_12*J_21; const double d_01 = J_12*J_20 - J_10*J_22; const double d_02 = J_10*J_21 - J_11*J_20; const double d_10 = J_02*J_21 - J_01*J_22; const double d_11 = J_00*J_22 - J_02*J_20; const double d_12 = J_01*J_20 - J_00*J_21; const double d_20 = J_01*J_12 - J_02*J_11; const double d_21 = J_02*J_10 - J_00*J_12; const double d_22 = J_00*J_11 - J_01*J_10; // Compute determinant of Jacobian double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20; // Compute inverse of Jacobian const double K_00 = d_00 / detJ; const double K_01 = d_10 / detJ; const double K_02 = d_20 / detJ; const double K_10 = d_01 / detJ; const double K_11 = d_11 / detJ; const double K_12 = d_21 / detJ; const double K_20 = d_02 / detJ; const double K_21 = d_12 / detJ; const double K_22 = d_22 / detJ; // Compute constants // Get coordinates and map to the reference (FIAT) element // Compute number of derivatives. unsigned int num_derivatives = 1; for (unsigned int r = 0; r < n; r++) { num_derivatives *= 3; }// end loop over 'r' // Declare pointer to two dimensional array that holds combinations of derivatives and initialise unsigned int **combinations = new unsigned int *[num_derivatives]; for (unsigned int row = 0; row < num_derivatives; row++) { combinations[row] = new unsigned int [n]; for (unsigned int col = 0; col < n; col++) combinations[row][col] = 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 > 2) combinations[row][col] = 0; else { combinations[row][col] += 1; break; } } } } // Compute inverse of Jacobian const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}}; // 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. Assuming that values is always an array. for (unsigned int r = 0; r < num_derivatives; r++) { values[r] = 0.000000000000000; }// end loop over 'r' // Array of basisvalues. double basisvalues[1] = {0.000000000000000}; // Declare helper variables. // Compute basisvalues. basisvalues[0] = 1.000000000000000; // Table(s) of coefficients. static const double coefficients0[1] = \ {1.000000000000000}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[1][1] = \ {{0.000000000000000}}; static const double dmats1[1][1] = \ {{0.000000000000000}}; static const double dmats2[1][1] = \ {{0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[1][1] = \ {{1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[1][1] = \ {{1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 1; t++) { for (unsigned int u = 0; u < 1; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 1; t++) { for (unsigned int u = 0; u < 1; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 1; t++) { for (unsigned int u = 0; u < 1; u++) { for (unsigned int tu = 0; tu < 1; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 1; t++) { for (unsigned int u = 0; u < 1; u++) { for (unsigned int tu = 0; tu < 1; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 1; t++) { for (unsigned int u = 0; u < 1; u++) { for (unsigned int tu = 0; tu < 1; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 1; s++) { for (unsigned int t = 0; t < 1; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; } /// Evaluate order n derivatives of all basis functions at given point in cell virtual void evaluate_basis_derivatives_all(unsigned int n, double* values, const double* coordinates, const ufc::cell& c) const { // Element is constant, calling evaluate_basis_derivatives. evaluate_basis_derivatives(0, n, values, coordinates, c); } /// Evaluate linear functional for dof i on the function f virtual double evaluate_dof(unsigned int i, const ufc::function& f, const ufc::cell& c) const { // Declare variables for result of evaluation. double vals[1]; // Declare variable for physical coordinates. double y[3]; const double * const * x = c.coordinates; switch (i) { case 0: { y[0] = 0.250000000000000*x[0][0] + 0.250000000000000*x[1][0] + 0.250000000000000*x[2][0] + 0.250000000000000*x[3][0]; y[1] = 0.250000000000000*x[0][1] + 0.250000000000000*x[1][1] + 0.250000000000000*x[2][1] + 0.250000000000000*x[3][1]; y[2] = 0.250000000000000*x[0][2] + 0.250000000000000*x[1][2] + 0.250000000000000*x[2][2] + 0.250000000000000*x[3][2]; f.evaluate(vals, y, c); return vals[0]; break; } } return 0.000000000000000; } /// Evaluate linear functionals for all dofs on the function f virtual void evaluate_dofs(double* values, const ufc::function& f, const ufc::cell& c) const { // Declare variables for result of evaluation. double vals[1]; // Declare variable for physical coordinates. double y[3]; const double * const * x = c.coordinates; y[0] = 0.250000000000000*x[0][0] + 0.250000000000000*x[1][0] + 0.250000000000000*x[2][0] + 0.250000000000000*x[3][0]; y[1] = 0.250000000000000*x[0][1] + 0.250000000000000*x[1][1] + 0.250000000000000*x[2][1] + 0.250000000000000*x[3][1]; y[2] = 0.250000000000000*x[0][2] + 0.250000000000000*x[1][2] + 0.250000000000000*x[2][2] + 0.250000000000000*x[3][2]; f.evaluate(vals, y, c); values[0] = vals[0]; } /// Interpolate vertex values from dof values virtual void interpolate_vertex_values(double* vertex_values, const double* dof_values, const ufc::cell& c) const { // Evaluate function and change variables vertex_values[0] = dof_values[0]; vertex_values[1] = dof_values[0]; vertex_values[2] = dof_values[0]; vertex_values[3] = dof_values[0]; } /// Return the number of sub elements (for a mixed element) virtual unsigned int num_sub_elements() const { return 0; } /// Create a new finite element for sub element i (for a mixed element) virtual ufc::finite_element* create_sub_element(unsigned int i) const { return 0; } }; /// This class defines the interface for a finite element. class elasticity_finite_element_1: public ufc::finite_element { public: /// Constructor elasticity_finite_element_1() : ufc::finite_element() { // Do nothing } /// Destructor virtual ~elasticity_finite_element_1() { // Do nothing } /// Return a string identifying the finite element virtual const char* signature() const { return "FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1)"; } /// Return the cell shape virtual ufc::shape cell_shape() const { return ufc::tetrahedron; } /// Return the dimension of the finite element function space virtual unsigned int space_dimension() const { return 4; } /// Return the rank of the value space virtual unsigned int value_rank() const { return 0; } /// Return the dimension of the value space for axis i virtual unsigned int value_dimension(unsigned int i) const { return 1; } /// Evaluate basis function i at given point in cell virtual void evaluate_basis(unsigned int i, double* values, const double* coordinates, const ufc::cell& c) const { // Extract vertex coordinates const double * const * x = c.coordinates; // Compute Jacobian of affine map from reference cell const double J_00 = x[1][0] - x[0][0]; const double J_01 = x[2][0] - x[0][0]; const double J_02 = x[3][0] - x[0][0]; const double J_10 = x[1][1] - x[0][1]; const double J_11 = x[2][1] - x[0][1]; const double J_12 = x[3][1] - x[0][1]; const double J_20 = x[1][2] - x[0][2]; const double J_21 = x[2][2] - x[0][2]; const double J_22 = x[3][2] - x[0][2]; // Compute sub determinants const double d_00 = J_11*J_22 - J_12*J_21; const double d_01 = J_12*J_20 - J_10*J_22; const double d_02 = J_10*J_21 - J_11*J_20; const double d_10 = J_02*J_21 - J_01*J_22; const double d_11 = J_00*J_22 - J_02*J_20; const double d_12 = J_01*J_20 - J_00*J_21; const double d_20 = J_01*J_12 - J_02*J_11; const double d_21 = J_02*J_10 - J_00*J_12; const double d_22 = J_00*J_11 - J_01*J_10; // Compute determinant of Jacobian double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20; // Compute inverse of Jacobian // Compute constants const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0]; const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1]; const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2]; // Get coordinates and map to the reference (FIAT) element double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ; double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ; double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ; // Reset values. *values = 0.000000000000000; switch (i) { case 0: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { *values += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 1: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { *values += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 2: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { *values += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 3: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { *values += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } } } /// Evaluate all basis functions at given point in cell virtual void evaluate_basis_all(double* values, const double* coordinates, const ufc::cell& c) const { // Helper variable to hold values of a single dof. double dof_values = 0.000000000000000; // Loop dofs and call evaluate_basis. for (unsigned int r = 0; r < 4; r++) { evaluate_basis(r, &dof_values, coordinates, c); values[r] = dof_values; }// end loop over 'r' } /// 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 * x = c.coordinates; // Compute Jacobian of affine map from reference cell const double J_00 = x[1][0] - x[0][0]; const double J_01 = x[2][0] - x[0][0]; const double J_02 = x[3][0] - x[0][0]; const double J_10 = x[1][1] - x[0][1]; const double J_11 = x[2][1] - x[0][1]; const double J_12 = x[3][1] - x[0][1]; const double J_20 = x[1][2] - x[0][2]; const double J_21 = x[2][2] - x[0][2]; const double J_22 = x[3][2] - x[0][2]; // Compute sub determinants const double d_00 = J_11*J_22 - J_12*J_21; const double d_01 = J_12*J_20 - J_10*J_22; const double d_02 = J_10*J_21 - J_11*J_20; const double d_10 = J_02*J_21 - J_01*J_22; const double d_11 = J_00*J_22 - J_02*J_20; const double d_12 = J_01*J_20 - J_00*J_21; const double d_20 = J_01*J_12 - J_02*J_11; const double d_21 = J_02*J_10 - J_00*J_12; const double d_22 = J_00*J_11 - J_01*J_10; // Compute determinant of Jacobian double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20; // Compute inverse of Jacobian const double K_00 = d_00 / detJ; const double K_01 = d_10 / detJ; const double K_02 = d_20 / detJ; const double K_10 = d_01 / detJ; const double K_11 = d_11 / detJ; const double K_12 = d_21 / detJ; const double K_20 = d_02 / detJ; const double K_21 = d_12 / detJ; const double K_22 = d_22 / detJ; // Compute constants const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0]; const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1]; const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2]; // Get coordinates and map to the reference (FIAT) element double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ; double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ; double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ; // Compute number of derivatives. unsigned int num_derivatives = 1; for (unsigned int r = 0; r < n; r++) { num_derivatives *= 3; }// end loop over 'r' // Declare pointer to two dimensional array that holds combinations of derivatives and initialise unsigned int **combinations = new unsigned int *[num_derivatives]; for (unsigned int row = 0; row < num_derivatives; row++) { combinations[row] = new unsigned int [n]; for (unsigned int col = 0; col < n; col++) combinations[row][col] = 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 > 2) combinations[row][col] = 0; else { combinations[row][col] += 1; break; } } } } // Compute inverse of Jacobian const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}}; // 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. Assuming that values is always an array. for (unsigned int r = 0; r < num_derivatives; r++) { values[r] = 0.000000000000000; }// end loop over 'r' switch (i) { case 0: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 1: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 2: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 3: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } } } /// Evaluate order n derivatives of all basis functions at given point in cell virtual void evaluate_basis_derivatives_all(unsigned int n, double* values, const double* coordinates, const ufc::cell& c) const { // Compute number of derivatives. unsigned int num_derivatives = 1; for (unsigned int r = 0; r < n; r++) { num_derivatives *= 3; }// end loop over 'r' // Helper variable to hold values of a single dof. double *dof_values = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { dof_values[r] = 0.000000000000000; }// end loop over 'r' // Loop dofs and call evaluate_basis_derivatives. for (unsigned int r = 0; r < 4; r++) { evaluate_basis_derivatives(r, n, dof_values, coordinates, c); for (unsigned int s = 0; s < num_derivatives; s++) { values[r*num_derivatives + s] = dof_values[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer. delete [] dof_values; } /// Evaluate linear functional for dof i on the function f virtual double evaluate_dof(unsigned int i, const ufc::function& f, const ufc::cell& c) const { // Declare variables for result of evaluation. double vals[1]; // Declare variable for physical coordinates. double y[3]; const double * const * x = c.coordinates; switch (i) { case 0: { y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); return vals[0]; break; } case 1: { y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); return vals[0]; break; } case 2: { y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); return vals[0]; break; } case 3: { y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); return vals[0]; break; } } return 0.000000000000000; } /// Evaluate linear functionals for all dofs on the function f virtual void evaluate_dofs(double* values, const ufc::function& f, const ufc::cell& c) const { // Declare variables for result of evaluation. double vals[1]; // Declare variable for physical coordinates. double y[3]; const double * const * x = c.coordinates; y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); values[0] = vals[0]; y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); values[1] = vals[0]; y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); values[2] = vals[0]; y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); values[3] = vals[0]; } /// Interpolate vertex values from dof values virtual void interpolate_vertex_values(double* vertex_values, const double* dof_values, const ufc::cell& c) const { // Evaluate function and change variables vertex_values[0] = dof_values[0]; vertex_values[1] = dof_values[1]; vertex_values[2] = dof_values[2]; vertex_values[3] = dof_values[3]; } /// Return the number of sub elements (for a mixed element) virtual unsigned int num_sub_elements() const { return 0; } /// Create a new finite element for sub element i (for a mixed element) virtual ufc::finite_element* create_sub_element(unsigned int i) const { return 0; } }; /// This class defines the interface for a finite element. class elasticity_finite_element_2: public ufc::finite_element { public: /// Constructor elasticity_finite_element_2() : ufc::finite_element() { // Do nothing } /// Destructor virtual ~elasticity_finite_element_2() { // Do nothing } /// Return a string identifying the finite element virtual const char* signature() const { return "VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3)"; } /// Return the cell shape virtual ufc::shape cell_shape() const { return ufc::tetrahedron; } /// Return the dimension of the finite element function space virtual unsigned int space_dimension() const { return 12; } /// Return the rank of the value space virtual unsigned int value_rank() const { return 1; } /// Return the dimension of the value space for axis i virtual unsigned int value_dimension(unsigned int i) const { switch (i) { case 0: { return 3; break; } } return 0; } /// Evaluate basis function i at given point in cell virtual void evaluate_basis(unsigned int i, double* values, const double* coordinates, const ufc::cell& c) const { // Extract vertex coordinates const double * const * x = c.coordinates; // Compute Jacobian of affine map from reference cell const double J_00 = x[1][0] - x[0][0]; const double J_01 = x[2][0] - x[0][0]; const double J_02 = x[3][0] - x[0][0]; const double J_10 = x[1][1] - x[0][1]; const double J_11 = x[2][1] - x[0][1]; const double J_12 = x[3][1] - x[0][1]; const double J_20 = x[1][2] - x[0][2]; const double J_21 = x[2][2] - x[0][2]; const double J_22 = x[3][2] - x[0][2]; // Compute sub determinants const double d_00 = J_11*J_22 - J_12*J_21; const double d_01 = J_12*J_20 - J_10*J_22; const double d_02 = J_10*J_21 - J_11*J_20; const double d_10 = J_02*J_21 - J_01*J_22; const double d_11 = J_00*J_22 - J_02*J_20; const double d_12 = J_01*J_20 - J_00*J_21; const double d_20 = J_01*J_12 - J_02*J_11; const double d_21 = J_02*J_10 - J_00*J_12; const double d_22 = J_00*J_11 - J_01*J_10; // Compute determinant of Jacobian double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20; // Compute inverse of Jacobian // Compute constants const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0]; const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1]; const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2]; // Get coordinates and map to the reference (FIAT) element double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ; double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ; double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ; // Reset values. values[0] = 0.000000000000000; values[1] = 0.000000000000000; values[2] = 0.000000000000000; switch (i) { case 0: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[0] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 1: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[0] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 2: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[0] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 3: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[0] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 4: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[1] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 5: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[1] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 6: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[1] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 7: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[1] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 8: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[2] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 9: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[2] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 10: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[2] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } case 11: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Compute value(s). for (unsigned int r = 0; r < 4; r++) { values[2] += coefficients0[r]*basisvalues[r]; }// end loop over 'r' break; } } } /// Evaluate all basis functions at given point in cell virtual void evaluate_basis_all(double* values, const double* coordinates, const ufc::cell& c) const { // Helper variable to hold values of a single dof. double dof_values[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000}; // Loop dofs and call evaluate_basis. for (unsigned int r = 0; r < 12; r++) { evaluate_basis(r, dof_values, coordinates, c); for (unsigned int s = 0; s < 3; s++) { values[r*3 + s] = dof_values[s]; }// end loop over 's' }// end loop over 'r' } /// 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 * x = c.coordinates; // Compute Jacobian of affine map from reference cell const double J_00 = x[1][0] - x[0][0]; const double J_01 = x[2][0] - x[0][0]; const double J_02 = x[3][0] - x[0][0]; const double J_10 = x[1][1] - x[0][1]; const double J_11 = x[2][1] - x[0][1]; const double J_12 = x[3][1] - x[0][1]; const double J_20 = x[1][2] - x[0][2]; const double J_21 = x[2][2] - x[0][2]; const double J_22 = x[3][2] - x[0][2]; // Compute sub determinants const double d_00 = J_11*J_22 - J_12*J_21; const double d_01 = J_12*J_20 - J_10*J_22; const double d_02 = J_10*J_21 - J_11*J_20; const double d_10 = J_02*J_21 - J_01*J_22; const double d_11 = J_00*J_22 - J_02*J_20; const double d_12 = J_01*J_20 - J_00*J_21; const double d_20 = J_01*J_12 - J_02*J_11; const double d_21 = J_02*J_10 - J_00*J_12; const double d_22 = J_00*J_11 - J_01*J_10; // Compute determinant of Jacobian double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20; // Compute inverse of Jacobian const double K_00 = d_00 / detJ; const double K_01 = d_10 / detJ; const double K_02 = d_20 / detJ; const double K_10 = d_01 / detJ; const double K_11 = d_11 / detJ; const double K_12 = d_21 / detJ; const double K_20 = d_02 / detJ; const double K_21 = d_12 / detJ; const double K_22 = d_22 / detJ; // Compute constants const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0]; const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1]; const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2]; // Get coordinates and map to the reference (FIAT) element double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ; double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ; double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ; // Compute number of derivatives. unsigned int num_derivatives = 1; for (unsigned int r = 0; r < n; r++) { num_derivatives *= 3; }// end loop over 'r' // Declare pointer to two dimensional array that holds combinations of derivatives and initialise unsigned int **combinations = new unsigned int *[num_derivatives]; for (unsigned int row = 0; row < num_derivatives; row++) { combinations[row] = new unsigned int [n]; for (unsigned int col = 0; col < n; col++) combinations[row][col] = 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 > 2) combinations[row][col] = 0; else { combinations[row][col] += 1; break; } } } } // Compute inverse of Jacobian const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}}; // 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. Assuming that values is always an array. for (unsigned int r = 0; r < 3*num_derivatives; r++) { values[r] = 0.000000000000000; }// end loop over 'r' switch (i) { case 0: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 1: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 2: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 3: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 4: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 5: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 6: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 7: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 8: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[2*num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 9: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[2*num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 10: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[2*num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } case 11: { // Array of basisvalues. double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}; // Declare helper variables. unsigned int rr = 0; unsigned int ss = 0; double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X); // Compute basisvalues. basisvalues[0] = 1.000000000000000; basisvalues[1] = tmp0; for (unsigned int r = 0; r < 1; r++) { rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2; ss = r*(r + 1)*(r + 2)/6; basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000); }// end loop over 'r' for (unsigned int r = 0; r < 1; r++) { for (unsigned int s = 0; s < 1 - r; s++) { rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1; ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2; basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s)); }// end loop over 's' }// end loop over 'r' for (unsigned int r = 0; r < 2; r++) { for (unsigned int s = 0; s < 2 - r; s++) { for (unsigned int t = 0; t < 2 - r - s; t++) { rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t; basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t)); }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Table(s) of coefficients. static const double coefficients0[4] = \ {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979}; // Tables of derivatives of the polynomial base (transpose). static const double dmats0[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats1[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; static const double dmats2[4][4] = \ {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}}; // Compute reference derivatives. // Declare pointer to array of derivatives on FIAT element. double *derivatives = new double[num_derivatives]; for (unsigned int r = 0; r < num_derivatives; r++) { derivatives[r] = 0.000000000000000; }// end loop over 'r' // Declare derivative matrix (of polynomial basis). double dmats[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Declare (auxiliary) derivative matrix (of polynomial basis). double dmats_old[4][4] = \ {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000}, {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}}; // Loop possible derivatives. for (unsigned int r = 0; r < num_derivatives; r++) { // Resetting dmats values to compute next derivative. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats[t][u] = 0.000000000000000; if (t == u) { dmats[t][u] = 1.000000000000000; } }// end loop over 'u' }// end loop over 't' // Looping derivative order to generate dmats. for (unsigned int s = 0; s < n; s++) { // Updating dmats_old with new values and resetting dmats. for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { dmats_old[t][u] = dmats[t][u]; dmats[t][u] = 0.000000000000000; }// end loop over 'u' }// end loop over 't' // Update dmats using an inner product. if (combinations[r][s] == 0) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 1) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } if (combinations[r][s] == 2) { for (unsigned int t = 0; t < 4; t++) { for (unsigned int u = 0; u < 4; u++) { for (unsigned int tu = 0; tu < 4; tu++) { dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u]; }// end loop over 'tu' }// end loop over 'u' }// end loop over 't' } }// end loop over 's' for (unsigned int s = 0; s < 4; s++) { for (unsigned int t = 0; t < 4; t++) { derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t]; }// end loop over 't' }// end loop over 's' }// end loop over 'r' // Transform derivatives back to physical element for (unsigned int r = 0; r < num_derivatives; r++) { for (unsigned int s = 0; s < num_derivatives; s++) { values[2*num_derivatives + r] += transform[r][s]*derivatives[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer to array of derivatives on FIAT element delete [] derivatives; // Delete pointer to array of combinations of derivatives and transform for (unsigned int r = 0; r < num_derivatives; r++) { delete [] combinations[r]; }// end loop over 'r' delete [] combinations; for (unsigned int r = 0; r < num_derivatives; r++) { delete [] transform[r]; }// end loop over 'r' delete [] transform; break; } } } /// Evaluate order n derivatives of all basis functions at given point in cell virtual void evaluate_basis_derivatives_all(unsigned int n, double* values, const double* coordinates, const ufc::cell& c) const { // Compute number of derivatives. unsigned int num_derivatives = 1; for (unsigned int r = 0; r < n; r++) { num_derivatives *= 3; }// end loop over 'r' // Helper variable to hold values of a single dof. double *dof_values = new double[3*num_derivatives]; for (unsigned int r = 0; r < 3*num_derivatives; r++) { dof_values[r] = 0.000000000000000; }// end loop over 'r' // Loop dofs and call evaluate_basis_derivatives. for (unsigned int r = 0; r < 12; r++) { evaluate_basis_derivatives(r, n, dof_values, coordinates, c); for (unsigned int s = 0; s < 3*num_derivatives; s++) { values[r*3*num_derivatives + s] = dof_values[s]; }// end loop over 's' }// end loop over 'r' // Delete pointer. delete [] dof_values; } /// Evaluate linear functional for dof i on the function f virtual double evaluate_dof(unsigned int i, const ufc::function& f, const ufc::cell& c) const { // Declare variables for result of evaluation. double vals[3]; // Declare variable for physical coordinates. double y[3]; const double * const * x = c.coordinates; switch (i) { case 0: { y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); return vals[0]; break; } case 1: { y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); return vals[0]; break; } case 2: { y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); return vals[0]; break; } case 3: { y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); return vals[0]; break; } case 4: { y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); return vals[1]; break; } case 5: { y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); return vals[1]; break; } case 6: { y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); return vals[1]; break; } case 7: { y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); return vals[1]; break; } case 8: { y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); return vals[2]; break; } case 9: { y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); return vals[2]; break; } case 10: { y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); return vals[2]; break; } case 11: { y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); return vals[2]; break; } } return 0.000000000000000; } /// Evaluate linear functionals for all dofs on the function f virtual void evaluate_dofs(double* values, const ufc::function& f, const ufc::cell& c) const { // Declare variables for result of evaluation. double vals[3]; // Declare variable for physical coordinates. double y[3]; const double * const * x = c.coordinates; y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); values[0] = vals[0]; y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); values[1] = vals[0]; y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); values[2] = vals[0]; y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); values[3] = vals[0]; y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); values[4] = vals[1]; y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); values[5] = vals[1]; y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); values[6] = vals[1]; y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); values[7] = vals[1]; y[0] = x[0][0]; y[1] = x[0][1]; y[2] = x[0][2]; f.evaluate(vals, y, c); values[8] = vals[2]; y[0] = x[1][0]; y[1] = x[1][1]; y[2] = x[1][2]; f.evaluate(vals, y, c); values[9] = vals[2]; y[0] = x[2][0]; y[1] = x[2][1]; y[2] = x[2][2]; f.evaluate(vals, y, c); values[10] = vals[2]; y[0] = x[3][0]; y[1] = x[3][1]; y[2] = x[3][2]; f.evaluate(vals, y, c); values[11] = vals[2]; } /// Interpolate vertex values from dof values virtual void interpolate_vertex_values(double* vertex_values, const double* dof_values, const ufc::cell& c) const { // Evaluate function and change variables vertex_values[0] = dof_values[0]; vertex_values[3] = dof_values[1]; vertex_values[6] = dof_values[2]; vertex_values[9] = dof_values[3]; // Evaluate function and change variables vertex_values[1] = dof_values[4]; vertex_values[4] = dof_values[5]; vertex_values[7] = dof_values[6]; vertex_values[10] = dof_values[7]; // Evaluate function and change variables vertex_values[2] = dof_values[8]; vertex_values[5] = dof_values[9]; vertex_values[8] = dof_values[10]; vertex_values[11] = dof_values[11]; } /// Return the number of sub elements (for a mixed element) virtual unsigned int num_sub_elements() const { return 3; } /// Create a new finite element for sub element i (for a mixed element) virtual ufc::finite_element* create_sub_element(unsigned int i) const { switch (i) { case 0: { return new elasticity_finite_element_1(); break; } case 1: { return new elasticity_finite_element_1(); break; } case 2: { return new elasticity_finite_element_1(); break; } } return 0; } }; /// This class defines the interface for a local-to-global mapping of /// degrees of freedom (dofs). class elasticity_dof_map_0: public ufc::dof_map { private: unsigned int _global_dimension; public: /// Constructor elasticity_dof_map_0() : ufc::dof_map() { _global_dimension = 0; } /// Destructor virtual ~elasticity_dof_map_0() { // Do nothing } /// Return a string identifying the dof map virtual const char* signature() const { return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Cell('tetrahedron', 1, Space(3)), 0)"; } /// Return true iff mesh entities of topological dimension d are needed virtual bool needs_mesh_entities(unsigned int d) const { switch (d) { case 0: { return false; break; } case 1: { return false; break; } case 2: { return false; break; } case 3: { return true; break; } } return false; } /// Initialize dof map for mesh (return true iff init_cell() is needed) virtual bool init_mesh(const ufc::mesh& m) { _global_dimension = m.num_entities[3]; return false; } /// Initialize dof map for given cell virtual void init_cell(const ufc::mesh& m, const ufc::cell& c) { // Do nothing } /// Finish initialization of dof map for cells virtual void init_cell_finalize() { // Do nothing } /// Return the dimension of the global finite element function space virtual unsigned int global_dimension() const { return _global_dimension; } /// Return the dimension of the local finite element function space for a cell virtual unsigned int local_dimension(const ufc::cell& c) const { return 1; } /// Return the maximum dimension of the local finite element function space virtual unsigned int max_local_dimension() const { return 1; } // Return the geometric dimension of the coordinates this dof map provides virtual unsigned int geometric_dimension() const { return 3; } /// Return the number of dofs on each cell facet virtual unsigned int num_facet_dofs() const { return 0; } /// Return the number of dofs associated with each cell entity of dimension d virtual unsigned int num_entity_dofs(unsigned int d) const { switch (d) { case 0: { return 0; break; } case 1: { return 0; break; } case 2: { return 0; break; } case 3: { return 1; break; } } return 0; } /// Tabulate the local-to-global mapping of dofs on a cell virtual void tabulate_dofs(unsigned int* dofs, const ufc::mesh& m, const ufc::cell& c) const { dofs[0] = c.entity_indices[3][0]; } /// Tabulate the local-to-local mapping from facet dofs to cell dofs virtual void tabulate_facet_dofs(unsigned int* dofs, unsigned int facet) const { switch (facet) { case 0: { break; } case 1: { break; } case 2: { break; } case 3: { break; } } } /// Tabulate the local-to-local mapping of dofs on entity (d, i) virtual void tabulate_entity_dofs(unsigned int* dofs, unsigned int d, unsigned int i) const { if (d > 3) { throw std::runtime_error("d is larger than dimension (3)"); } switch (d) { case 0: { break; } case 1: { break; } case 2: { break; } case 3: { if (i > 0) { throw std::runtime_error("i is larger than number of entities (0)"); } dofs[0] = 0; break; } } } /// Tabulate the coordinates of all dofs on a cell virtual void tabulate_coordinates(double** coordinates, const ufc::cell& c) const { const double * const * x = c.coordinates; coordinates[0][0] = 0.250000000000000*x[0][0] + 0.250000000000000*x[1][0] + 0.250000000000000*x[2][0] + 0.250000000000000*x[3][0]; coordinates[0][1] = 0.250000000000000*x[0][1] + 0.250000000000000*x[1][1] + 0.250000000000000*x[2][1] + 0.250000000000000*x[3][1]; coordinates[0][2] = 0.250000000000000*x[0][2] + 0.250000000000000*x[1][2] + 0.250000000000000*x[2][2] + 0.250000000000000*x[3][2]; } /// Return the number of sub dof maps (for a mixed element) virtual unsigned int num_sub_dof_maps() const { return 0; } /// Create a new dof_map for sub dof map i (for a mixed element) virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const { return 0; } }; /// This class defines the interface for a local-to-global mapping of /// degrees of freedom (dofs). class elasticity_dof_map_1: public ufc::dof_map { private: unsigned int _global_dimension; public: /// Constructor elasticity_dof_map_1() : ufc::dof_map() { _global_dimension = 0; } /// Destructor virtual ~elasticity_dof_map_1() { // Do nothing } /// Return a string identifying the dof map virtual const char* signature() const { return "FFC dofmap for FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1)"; } /// Return true iff mesh entities of topological dimension d are needed virtual bool needs_mesh_entities(unsigned int d) const { switch (d) { case 0: { return true; break; } case 1: { return false; break; } case 2: { return false; break; } case 3: { return false; break; } } return false; } /// Initialize dof map for mesh (return true iff init_cell() is needed) virtual bool init_mesh(const ufc::mesh& m) { _global_dimension = m.num_entities[0]; return false; } /// Initialize dof map for given cell virtual void init_cell(const ufc::mesh& m, const ufc::cell& c) { // Do nothing } /// Finish initialization of dof map for cells virtual void init_cell_finalize() { // Do nothing } /// Return the dimension of the global finite element function space virtual unsigned int global_dimension() const { return _global_dimension; } /// Return the dimension of the local finite element function space for a cell virtual unsigned int local_dimension(const ufc::cell& c) const { return 4; } /// Return the maximum dimension of the local finite element function space virtual unsigned int max_local_dimension() const { return 4; } // Return the geometric dimension of the coordinates this dof map provides virtual unsigned int geometric_dimension() const { return 3; } /// Return the number of dofs on each cell facet virtual unsigned int num_facet_dofs() const { return 3; } /// Return the number of dofs associated with each cell entity of dimension d virtual unsigned int num_entity_dofs(unsigned int d) const { switch (d) { case 0: { return 1; break; } case 1: { return 0; break; } case 2: { return 0; break; } case 3: { return 0; break; } } return 0; } /// Tabulate the local-to-global mapping of dofs on a cell virtual void tabulate_dofs(unsigned int* dofs, const ufc::mesh& m, const ufc::cell& c) const { dofs[0] = c.entity_indices[0][0]; dofs[1] = c.entity_indices[0][1]; dofs[2] = c.entity_indices[0][2]; dofs[3] = c.entity_indices[0][3]; } /// Tabulate the local-to-local mapping from facet dofs to cell dofs virtual void tabulate_facet_dofs(unsigned int* dofs, unsigned int facet) const { switch (facet) { case 0: { dofs[0] = 1; dofs[1] = 2; dofs[2] = 3; break; } case 1: { dofs[0] = 0; dofs[1] = 2; dofs[2] = 3; break; } case 2: { dofs[0] = 0; dofs[1] = 1; dofs[2] = 3; break; } case 3: { dofs[0] = 0; dofs[1] = 1; dofs[2] = 2; break; } } } /// Tabulate the local-to-local mapping of dofs on entity (d, i) virtual void tabulate_entity_dofs(unsigned int* dofs, unsigned int d, unsigned int i) const { if (d > 3) { throw std::runtime_error("d is larger than dimension (3)"); } switch (d) { case 0: { if (i > 3) { throw std::runtime_error("i is larger than number of entities (3)"); } switch (i) { case 0: { dofs[0] = 0; break; } case 1: { dofs[0] = 1; break; } case 2: { dofs[0] = 2; break; } case 3: { dofs[0] = 3; break; } } break; } case 1: { break; } case 2: { break; } case 3: { break; } } } /// Tabulate the coordinates of all dofs on a cell virtual void tabulate_coordinates(double** coordinates, const ufc::cell& c) const { const double * const * x = c.coordinates; coordinates[0][0] = x[0][0]; coordinates[0][1] = x[0][1]; coordinates[0][2] = x[0][2]; coordinates[1][0] = x[1][0]; coordinates[1][1] = x[1][1]; coordinates[1][2] = x[1][2]; coordinates[2][0] = x[2][0]; coordinates[2][1] = x[2][1]; coordinates[2][2] = x[2][2]; coordinates[3][0] = x[3][0]; coordinates[3][1] = x[3][1]; coordinates[3][2] = x[3][2]; } /// Return the number of sub dof maps (for a mixed element) virtual unsigned int num_sub_dof_maps() const { return 0; } /// Create a new dof_map for sub dof map i (for a mixed element) virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const { return 0; } }; /// This class defines the interface for a local-to-global mapping of /// degrees of freedom (dofs). class elasticity_dof_map_2: public ufc::dof_map { private: unsigned int _global_dimension; public: /// Constructor elasticity_dof_map_2() : ufc::dof_map() { _global_dimension = 0; } /// Destructor virtual ~elasticity_dof_map_2() { // Do nothing } /// Return a string identifying the dof map virtual const char* signature() const { return "FFC dofmap for VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3)"; } /// Return true iff mesh entities of topological dimension d are needed virtual bool needs_mesh_entities(unsigned int d) const { switch (d) { case 0: { return true; break; } case 1: { return false; break; } case 2: { return false; break; } case 3: { return false; break; } } return false; } /// Initialize dof map for mesh (return true iff init_cell() is needed) virtual bool init_mesh(const ufc::mesh& m) { _global_dimension = 3.000000000000000*m.num_entities[0]; return false; } /// Initialize dof map for given cell virtual void init_cell(const ufc::mesh& m, const ufc::cell& c) { // Do nothing } /// Finish initialization of dof map for cells virtual void init_cell_finalize() { // Do nothing } /// Return the dimension of the global finite element function space virtual unsigned int global_dimension() const { return _global_dimension; } /// Return the dimension of the local finite element function space for a cell virtual unsigned int local_dimension(const ufc::cell& c) const { return 12; } /// Return the maximum dimension of the local finite element function space virtual unsigned int max_local_dimension() const { return 12; } // Return the geometric dimension of the coordinates this dof map provides virtual unsigned int geometric_dimension() const { return 3; } /// Return the number of dofs on each cell facet virtual unsigned int num_facet_dofs() const { return 9; } /// Return the number of dofs associated with each cell entity of dimension d virtual unsigned int num_entity_dofs(unsigned int d) const { switch (d) { case 0: { return 3; break; } case 1: { return 0; break; } case 2: { return 0; break; } case 3: { return 0; break; } } return 0; } /// Tabulate the local-to-global mapping of dofs on a cell virtual void tabulate_dofs(unsigned int* dofs, const ufc::mesh& m, const ufc::cell& c) const { unsigned int offset = 0; dofs[0] = offset + c.entity_indices[0][0]; dofs[1] = offset + c.entity_indices[0][1]; dofs[2] = offset + c.entity_indices[0][2]; dofs[3] = offset + c.entity_indices[0][3]; offset += m.num_entities[0]; dofs[4] = offset + c.entity_indices[0][0]; dofs[5] = offset + c.entity_indices[0][1]; dofs[6] = offset + c.entity_indices[0][2]; dofs[7] = offset + c.entity_indices[0][3]; offset += m.num_entities[0]; dofs[8] = offset + c.entity_indices[0][0]; dofs[9] = offset + c.entity_indices[0][1]; dofs[10] = offset + c.entity_indices[0][2]; dofs[11] = offset + c.entity_indices[0][3]; offset += m.num_entities[0]; } /// Tabulate the local-to-local mapping from facet dofs to cell dofs virtual void tabulate_facet_dofs(unsigned int* dofs, unsigned int facet) const { switch (facet) { case 0: { dofs[0] = 1; dofs[1] = 2; dofs[2] = 3; dofs[3] = 5; dofs[4] = 6; dofs[5] = 7; dofs[6] = 9; dofs[7] = 10; dofs[8] = 11; break; } case 1: { dofs[0] = 0; dofs[1] = 2; dofs[2] = 3; dofs[3] = 4; dofs[4] = 6; dofs[5] = 7; dofs[6] = 8; dofs[7] = 10; dofs[8] = 11; break; } case 2: { dofs[0] = 0; dofs[1] = 1; dofs[2] = 3; dofs[3] = 4; dofs[4] = 5; dofs[5] = 7; dofs[6] = 8; dofs[7] = 9; dofs[8] = 11; break; } case 3: { dofs[0] = 0; dofs[1] = 1; dofs[2] = 2; dofs[3] = 4; dofs[4] = 5; dofs[5] = 6; dofs[6] = 8; dofs[7] = 9; dofs[8] = 10; break; } } } /// Tabulate the local-to-local mapping of dofs on entity (d, i) virtual void tabulate_entity_dofs(unsigned int* dofs, unsigned int d, unsigned int i) const { if (d > 3) { throw std::runtime_error("d is larger than dimension (3)"); } switch (d) { case 0: { if (i > 3) { throw std::runtime_error("i is larger than number of entities (3)"); } switch (i) { case 0: { dofs[0] = 0; dofs[1] = 4; dofs[2] = 8; break; } case 1: { dofs[0] = 1; dofs[1] = 5; dofs[2] = 9; break; } case 2: { dofs[0] = 2; dofs[1] = 6; dofs[2] = 10; break; } case 3: { dofs[0] = 3; dofs[1] = 7; dofs[2] = 11; break; } } break; } case 1: { break; } case 2: { break; } case 3: { break; } } } /// Tabulate the coordinates of all dofs on a cell virtual void tabulate_coordinates(double** coordinates, const ufc::cell& c) const { const double * const * x = c.coordinates; coordinates[0][0] = x[0][0]; coordinates[0][1] = x[0][1]; coordinates[0][2] = x[0][2]; coordinates[1][0] = x[1][0]; coordinates[1][1] = x[1][1]; coordinates[1][2] = x[1][2]; coordinates[2][0] = x[2][0]; coordinates[2][1] = x[2][1]; coordinates[2][2] = x[2][2]; coordinates[3][0] = x[3][0]; coordinates[3][1] = x[3][1]; coordinates[3][2] = x[3][2]; coordinates[4][0] = x[0][0]; coordinates[4][1] = x[0][1]; coordinates[4][2] = x[0][2]; coordinates[5][0] = x[1][0]; coordinates[5][1] = x[1][1]; coordinates[5][2] = x[1][2]; coordinates[6][0] = x[2][0]; coordinates[6][1] = x[2][1]; coordinates[6][2] = x[2][2]; coordinates[7][0] = x[3][0]; coordinates[7][1] = x[3][1]; coordinates[7][2] = x[3][2]; coordinates[8][0] = x[0][0]; coordinates[8][1] = x[0][1]; coordinates[8][2] = x[0][2]; coordinates[9][0] = x[1][0]; coordinates[9][1] = x[1][1]; coordinates[9][2] = x[1][2]; coordinates[10][0] = x[2][0]; coordinates[10][1] = x[2][1]; coordinates[10][2] = x[2][2]; coordinates[11][0] = x[3][0]; coordinates[11][1] = x[3][1]; coordinates[11][2] = x[3][2]; } /// Return the number of sub dof maps (for a mixed element) virtual unsigned int num_sub_dof_maps() const { return 3; } /// Create a new dof_map for sub dof map i (for a mixed element) virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const { switch (i) { case 0: { return new elasticity_dof_map_1(); break; } case 1: { return new elasticity_dof_map_1(); break; } case 2: { return new elasticity_dof_map_1(); break; } } return 0; } }; /// This class defines the interface for the tabulation of the cell /// tensor corresponding to the local contribution to a form from /// the integral over a cell. class elasticity_cell_integral_0_0: public ufc::cell_integral { public: /// Constructor elasticity_cell_integral_0_0() : ufc::cell_integral() { // Do nothing } /// Destructor virtual ~elasticity_cell_integral_0_0() { // Do nothing } /// Tabulate the tensor for the contribution from a local cell virtual void tabulate_tensor(double* A, const double * const * w, const ufc::cell& c) const { // Extract vertex coordinates const double * const * x = c.coordinates; // Compute Jacobian of affine map from reference cell const double J_00 = x[1][0] - x[0][0]; const double J_01 = x[2][0] - x[0][0]; const double J_02 = x[3][0] - x[0][0]; const double J_10 = x[1][1] - x[0][1]; const double J_11 = x[2][1] - x[0][1]; const double J_12 = x[3][1] - x[0][1]; const double J_20 = x[1][2] - x[0][2]; const double J_21 = x[2][2] - x[0][2]; const double J_22 = x[3][2] - x[0][2]; // Compute sub determinants const double d_00 = J_11*J_22 - J_12*J_21; const double d_01 = J_12*J_20 - J_10*J_22; const double d_02 = J_10*J_21 - J_11*J_20; const double d_10 = J_02*J_21 - J_01*J_22; const double d_11 = J_00*J_22 - J_02*J_20; const double d_12 = J_01*J_20 - J_00*J_21; const double d_20 = J_01*J_12 - J_02*J_11; const double d_21 = J_02*J_10 - J_00*J_12; const double d_22 = J_00*J_11 - J_01*J_10; // Compute determinant of Jacobian double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20; // Compute inverse of Jacobian const double K_00 = d_00 / detJ; const double K_01 = d_10 / detJ; const double K_02 = d_20 / detJ; const double K_10 = d_01 / detJ; const double K_11 = d_11 / detJ; const double K_12 = d_21 / detJ; const double K_20 = d_02 / detJ; const double K_21 = d_12 / detJ; const double K_22 = d_22 / detJ; // Set scale factor const double det = std::abs(detJ); // Array of quadrature weights. static const double W1 = 0.166666666666666; // Quadrature points on the UFC reference element: (0.250000000000000, 0.250000000000000, 0.250000000000000) // Value of basis functions at quadrature points. static const double FE1_C0_D001[1][2] = \ {{-1.000000000000000, 1.000000000000000}}; // Array of non-zero columns static const unsigned int nzc1[2] = {0, 3}; // Array of non-zero columns static const unsigned int nzc3[2] = {0, 1}; // Array of non-zero columns static const unsigned int nzc9[2] = {8, 11}; // Array of non-zero columns static const unsigned int nzc2[2] = {0, 2}; // Array of non-zero columns static const unsigned int nzc6[2] = {4, 6}; // Array of non-zero columns static const unsigned int nzc7[2] = {4, 5}; // Array of non-zero columns static const unsigned int nzc11[2] = {8, 9}; // Array of non-zero columns static const unsigned int nzc10[2] = {8, 10}; // Array of non-zero columns static const unsigned int nzc5[2] = {4, 7}; // Reset values in the element tensor. for (unsigned int r = 0; r < 144; r++) { A[r] = 0.000000000000000; }// end loop over 'r' // Number of operations to compute geometry constants: 385. double G[45]; G[0] = W1*det*(K_01*K_01*w[1][0] + w[0][0]*(2.000000000000000*K_01*K_01 + K_00*K_00 + K_02*K_02)); G[1] = K_00*K_01*W1*det*(w[0][0] + w[1][0]); G[2] = W1*det*(K_11*K_22*w[0][0] + K_12*K_21*w[1][0]); G[3] = W1*det*(K_00*K_00*w[1][0] + w[0][0]*(2.000000000000000*K_00*K_00 + K_01*K_01 + K_02*K_02)); G[4] = W1*det*(K_01*K_21*(w[1][0] + 2.000000000000000*w[0][0]) + w[0][0]*(K_00*K_20 + K_02*K_22)); G[5] = K_10*K_12*W1*det*(w[0][0] + w[1][0]); G[6] = W1*det*(K_00*K_11*w[1][0] + K_01*K_10*w[0][0]); G[7] = W1*det*(K_12*K_22*w[1][0] + w[0][0]*(2.000000000000000*K_12*K_22 + K_10*K_20 + K_11*K_21)); G[8] = W1*det*(K_00*K_22*w[0][0] + K_02*K_20*w[1][0]); G[9] = K_00*K_02*W1*det*(w[0][0] + w[1][0]); G[10] = K_10*K_11*W1*det*(w[0][0] + w[1][0]); G[11] = W1*det*(K_10*K_22*w[0][0] + K_12*K_20*w[1][0]); G[12] = W1*det*(K_02*K_22*w[1][0] + w[0][0]*(2.000000000000000*K_02*K_22 + K_00*K_20 + K_01*K_21)); G[13] = W1*det*(K_00*K_10*(w[1][0] + 2.000000000000000*w[0][0]) + w[0][0]*(K_01*K_11 + K_02*K_12)); G[14] = W1*det*(K_01*K_22*w[0][0] + K_02*K_21*w[1][0]); G[15] = W1*det*(K_02*K_02*w[1][0] + w[0][0]*(2.000000000000000*K_02*K_02 + K_00*K_00 + K_01*K_01)); G[16] = W1*det*(K_00*K_21*w[0][0] + K_01*K_20*w[1][0]); G[17] = W1*det*(K_01*K_12*w[0][0] + K_02*K_11*w[1][0]); G[18] = W1*det*(K_10*K_20*w[1][0] + w[0][0]*(2.000000000000000*K_10*K_20 + K_11*K_21 + K_12*K_22)); G[19] = W1*det*(K_01*K_22*w[1][0] + K_02*K_21*w[0][0]); G[20] = W1*det*(K_11*K_22*w[1][0] + K_12*K_21*w[0][0]); G[21] = W1*det*(K_01*K_11*w[1][0] + w[0][0]*(2.000000000000000*K_01*K_11 + K_00*K_10 + K_02*K_12)); G[22] = W1*det*(K_00*K_12*w[0][0] + K_02*K_10*w[1][0]); G[23] = W1*det*(K_10*K_21*w[1][0] + K_11*K_20*w[0][0]); G[24] = K_20*K_21*W1*det*(w[0][0] + w[1][0]); G[25] = W1*det*(K_00*K_11*w[0][0] + K_01*K_10*w[1][0]); G[26] = K_20*K_22*W1*det*(w[0][0] + w[1][0]); G[27] = W1*det*(K_00*K_20*w[1][0] + w[0][0]*(2.000000000000000*K_00*K_20 + K_01*K_21 + K_02*K_22)); G[28] = K_21*K_22*W1*det*(w[0][0] + w[1][0]); G[29] = W1*det*(K_01*K_12*w[1][0] + K_02*K_11*w[0][0]); G[30] = W1*det*(K_00*K_22*w[1][0] + K_02*K_20*w[0][0]); G[31] = W1*det*(K_12*K_12*w[1][0] + w[0][0]*(2.000000000000000*K_12*K_12 + K_10*K_10 + K_11*K_11)); G[32] = W1*det*(K_21*K_21*w[1][0] + w[0][0]*(2.000000000000000*K_21*K_21 + K_20*K_20 + K_22*K_22)); G[33] = K_01*K_02*W1*det*(w[0][0] + w[1][0]); G[34] = W1*det*(K_00*K_21*w[1][0] + K_01*K_20*w[0][0]); G[35] = W1*det*(K_02*K_12*w[1][0] + w[0][0]*(2.000000000000000*K_02*K_12 + K_00*K_10 + K_01*K_11)); G[36] = K_11*K_12*W1*det*(w[0][0] + w[1][0]); G[37] = W1*det*(K_10*K_22*w[1][0] + K_12*K_20*w[0][0]); G[38] = W1*det*(K_00*K_12*w[1][0] + K_02*K_10*w[0][0]); G[39] = W1*det*(K_22*K_22*w[1][0] + w[0][0]*(2.000000000000000*K_22*K_22 + K_20*K_20 + K_21*K_21)); G[40] = W1*det*(K_10*K_10*w[1][0] + w[0][0]*(2.000000000000000*K_10*K_10 + K_11*K_11 + K_12*K_12)); G[41] = W1*det*(K_10*K_21*w[0][0] + K_11*K_20*w[1][0]); G[42] = W1*det*(K_11*K_11*w[1][0] + w[0][0]*(2.000000000000000*K_11*K_11 + K_10*K_10 + K_12*K_12)); G[43] = W1*det*(K_11*K_21*w[1][0] + w[0][0]*(2.000000000000000*K_11*K_21 + K_10*K_20 + K_12*K_22)); G[44] = W1*det*(K_20*K_20*w[1][0] + w[0][0]*(2.000000000000000*K_20*K_20 + K_21*K_21 + K_22*K_22)); // Compute element tensor using UFL quadrature representation // Optimisations: ('eliminate zeros', True), ('ignore ones', True), ('ignore zero tables', True), ('optimisation', 'simplify_expressions'), ('remove zero terms', True) // Loop quadrature points for integral. // Number of operations to compute element tensor for following IP loop = 972 // Only 1 integration point, omitting IP loop. // Number of operations for primary indices: 972 for (unsigned int j = 0; j < 2; j++) { for (unsigned int k = 0; k < 2; k++) { // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[0]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[1]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[2]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[3]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[4]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[5]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[6]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[7]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[8]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[9]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[10]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[11]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[12]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[13]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[14]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[10]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[6]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[15]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[16]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[17]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[18]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[19]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[20]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[5]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[7]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[21]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[22]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[23]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[21]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[12]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[22]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[24]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[25]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[26]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[13]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[27]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[28]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[2]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[23]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[29]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[25]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[30]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[31]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[26]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[8]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[32]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[33]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[34]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[35]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[14]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[36]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[37]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[36]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc10[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[38]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[17]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[39]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[18]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc2[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[40]; // Number of operations to compute entry: 3 A[nzc2[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[37]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[34]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[41]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[33]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[24]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[11]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[42]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[1]; // Number of operations to compute entry: 3 A[nzc9[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[20]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[16]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[29]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc6[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[43]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[27]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[38]; // Number of operations to compute entry: 3 A[nzc7[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[19]; // Number of operations to compute entry: 3 A[nzc10[j]*12 + nzc11[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[35]; // Number of operations to compute entry: 3 A[nzc1[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[44]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc1[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[41]; // Number of operations to compute entry: 3 A[nzc11[j]*12 + nzc3[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[9]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[28]; // Number of operations to compute entry: 3 A[nzc3[j]*12 + nzc9[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[30]; // Number of operations to compute entry: 3 A[nzc6[j]*12 + nzc5[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[43]; // Number of operations to compute entry: 3 A[nzc5[j]*12 + nzc7[k]] += FE1_C0_D001[0][j]*FE1_C0_D001[0][k]*G[4]; }// end loop over 'k' }// end loop over 'j' } }; /// This class defines the interface for the tabulation of the cell /// tensor corresponding to the local contribution to a form from /// the integral over a cell. class elasticity_cell_integral_1_0: public ufc::cell_integral { public: /// Constructor elasticity_cell_integral_1_0() : ufc::cell_integral() { // Do nothing } /// Destructor virtual ~elasticity_cell_integral_1_0() { // Do nothing } /// Tabulate the tensor for the contribution from a local cell virtual void tabulate_tensor(double* A, const double * const * w, const ufc::cell& c) const { // Number of operations (multiply-add pairs) for Jacobian data: 18 // Number of operations (multiply-add pairs) for geometry tensor: 12 // Number of operations (multiply-add pairs) for tensor contraction: 42 // Total number of operations (multiply-add pairs): 72 // Extract vertex coordinates const double * const * x = c.coordinates; // Compute Jacobian of affine map from reference cell const double J_00 = x[1][0] - x[0][0]; const double J_01 = x[2][0] - x[0][0]; const double J_02 = x[3][0] - x[0][0]; const double J_10 = x[1][1] - x[0][1]; const double J_11 = x[2][1] - x[0][1]; const double J_12 = x[3][1] - x[0][1]; const double J_20 = x[1][2] - x[0][2]; const double J_21 = x[2][2] - x[0][2]; const double J_22 = x[3][2] - x[0][2]; // Compute sub determinants const double d_00 = J_11*J_22 - J_12*J_21; const double d_10 = J_02*J_21 - J_01*J_22; const double d_20 = J_01*J_12 - J_02*J_11; // Compute determinant of Jacobian double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20; // Compute inverse of Jacobian // Set scale factor const double det = std::abs(detJ); // Compute geometry tensor const double G0_0 = det*w[0][0]*(1.0); const double G0_1 = det*w[0][1]*(1.0); const double G0_2 = det*w[0][2]*(1.0); const double G0_3 = det*w[0][3]*(1.0); const double G0_4 = det*w[0][4]*(1.0); const double G0_5 = det*w[0][5]*(1.0); const double G0_6 = det*w[0][6]*(1.0); const double G0_7 = det*w[0][7]*(1.0); const double G0_8 = det*w[0][8]*(1.0); const double G0_9 = det*w[0][9]*(1.0); const double G0_10 = det*w[0][10]*(1.0); const double G0_11 = det*w[0][11]*(1.0); // Compute element tensor A[0] = 0.016666666666667*G0_0 + 0.008333333333333*G0_1 + 0.008333333333333*G0_2 + 0.008333333333333*G0_3; A[1] = 0.008333333333333*G0_0 + 0.016666666666667*G0_1 + 0.008333333333333*G0_2 + 0.008333333333333*G0_3; A[2] = 0.008333333333333*G0_0 + 0.008333333333333*G0_1 + 0.016666666666667*G0_2 + 0.008333333333333*G0_3; A[3] = 0.008333333333333*G0_0 + 0.008333333333333*G0_1 + 0.008333333333333*G0_2 + 0.016666666666667*G0_3; A[4] = 0.016666666666667*G0_4 + 0.008333333333333*G0_5 + 0.008333333333333*G0_6 + 0.008333333333333*G0_7; A[5] = 0.008333333333333*G0_4 + 0.016666666666667*G0_5 + 0.008333333333333*G0_6 + 0.008333333333333*G0_7; A[6] = 0.008333333333333*G0_4 + 0.008333333333333*G0_5 + 0.016666666666667*G0_6 + 0.008333333333333*G0_7; A[7] = 0.008333333333333*G0_4 + 0.008333333333333*G0_5 + 0.008333333333333*G0_6 + 0.016666666666667*G0_7; A[8] = 0.016666666666667*G0_8 + 0.008333333333333*G0_9 + 0.008333333333333*G0_10 + 0.008333333333333*G0_11; A[9] = 0.008333333333333*G0_8 + 0.016666666666667*G0_9 + 0.008333333333333*G0_10 + 0.008333333333333*G0_11; A[10] = 0.008333333333333*G0_8 + 0.008333333333333*G0_9 + 0.016666666666667*G0_10 + 0.008333333333333*G0_11; A[11] = 0.008333333333333*G0_8 + 0.008333333333333*G0_9 + 0.008333333333333*G0_10 + 0.016666666666667*G0_11; } }; /// This class defines the interface for the assembly of the global /// tensor corresponding to a form with r + n arguments, that is, a /// mapping /// /// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R /// /// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r /// global tensor A is defined by /// /// A = a(V1, V2, ..., Vr, w1, w2, ..., wn), /// /// where each argument Vj represents the application to the /// sequence of basis functions of Vj and w1, w2, ..., wn are given /// fixed functions (coefficients). class elasticity_form_0: public ufc::form { public: /// Constructor elasticity_form_0() : ufc::form() { // Do nothing } /// Destructor virtual ~elasticity_form_0() { // Do nothing } /// Return a string identifying the form virtual const char* signature() const { return "Form([Integral(IndexSum(IndexSum(Product(Indexed(ComponentTensor(Product(FloatValue(0.5, (), (), {}), Indexed(Sum(ComponentTensor(Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 0), MultiIndex((Index(0),), {Index(0): 3})), MultiIndex((Index(1),), {Index(1): 3})), MultiIndex((Index(1), Index(0)), {Index(0): 3, Index(1): 3})), MultiIndex((Index(2), Index(3)), {Index(2): 3, Index(3): 3})), MultiIndex((Index(3), Index(2)), {Index(2): 3, Index(3): 3})), ComponentTensor(Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 0), MultiIndex((Index(4),), {Index(4): 3})), MultiIndex((Index(5),), {Index(5): 3})), MultiIndex((Index(5), Index(4)), {Index(4): 3, Index(5): 3}))), MultiIndex((Index(6), Index(7)), {Index(7): 3, Index(6): 3}))), MultiIndex((Index(6), Index(7)), {Index(7): 3, Index(6): 3})), MultiIndex((Index(8), Index(9)), {Index(8): 3, Index(9): 3})), Indexed(Sum(ComponentTensor(Product(Indexed(ComponentTensor(Product(FloatValue(0.5, (), (), {}), Indexed(Sum(ComponentTensor(Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 1), MultiIndex((Index(10),), {Index(10): 3})), MultiIndex((Index(11),), {Index(11): 3})), MultiIndex((Index(11), Index(10)), {Index(11): 3, Index(10): 3})), MultiIndex((Index(12), Index(13)), {Index(13): 3, Index(12): 3})), MultiIndex((Index(13), Index(12)), {Index(13): 3, Index(12): 3})), ComponentTensor(Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 1), MultiIndex((Index(14),), {Index(14): 3})), MultiIndex((Index(15),), {Index(15): 3})), MultiIndex((Index(15), Index(14)), {Index(14): 3, Index(15): 3}))), MultiIndex((Index(16), Index(17)), {Index(17): 3, Index(16): 3}))), MultiIndex((Index(16), Index(17)), {Index(17): 3, Index(16): 3})), MultiIndex((Index(18), Index(19)), {Index(19): 3, Index(18): 3})), Product(FloatValue(2.0, (), (), {}), Constant(Cell('tetrahedron', 1, Space(3)), 0))), MultiIndex((Index(18), Index(19)), {Index(19): 3, Index(18): 3})), ComponentTensor(Product(Indexed(Identity(3), MultiIndex((Index(20), Index(21)), {Index(21): 3, Index(20): 3})), Product(Constant(Cell('tetrahedron', 1, Space(3)), 1), IndexSum(Indexed(ComponentTensor(Product(FloatValue(0.5, (), (), {}), Indexed(Sum(ComponentTensor(Indexed(ComponentTensor(Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 1), MultiIndex((Index(22),), {Index(22): 3})), MultiIndex((Index(23),), {Index(23): 3})), MultiIndex((Index(23), Index(22)), {Index(22): 3, Index(23): 3})), MultiIndex((Index(24), Index(25)), {Index(24): 3, Index(25): 3})), MultiIndex((Index(25), Index(24)), {Index(24): 3, Index(25): 3})), ComponentTensor(Indexed(SpatialDerivative(Argument(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 1), MultiIndex((Index(26),), {Index(26): 3})), MultiIndex((Index(27),), {Index(27): 3})), MultiIndex((Index(27), Index(26)), {Index(26): 3, Index(27): 3}))), MultiIndex((Index(28), Index(29)), {Index(29): 3, Index(28): 3}))), MultiIndex((Index(28), Index(29)), {Index(29): 3, Index(28): 3})), MultiIndex((Index(30), Index(30)), {Index(30): 3})), MultiIndex((Index(30),), {Index(30): 3})))), MultiIndex((Index(20), Index(21)), {Index(21): 3, Index(20): 3}))), MultiIndex((Index(8), Index(9)), {Index(8): 3, Index(9): 3}))), MultiIndex((Index(8),), {Index(8): 3})), MultiIndex((Index(9),), {Index(9): 3})), Measure('cell', 0, None))])"; } /// Return the rank of the global tensor (r) virtual unsigned int rank() const { return 2; } /// Return the number of coefficients (n) virtual unsigned int num_coefficients() const { return 2; } /// Return the number of cell integrals virtual unsigned int num_cell_integrals() const { return 1; } /// Return the number of exterior facet integrals virtual unsigned int num_exterior_facet_integrals() const { return 0; } /// Return the number of interior facet integrals virtual unsigned int num_interior_facet_integrals() const { return 0; } /// Create a new finite element for argument function i virtual ufc::finite_element* create_finite_element(unsigned int i) const { switch (i) { case 0: { return new elasticity_finite_element_2(); break; } case 1: { return new elasticity_finite_element_2(); break; } case 2: { return new elasticity_finite_element_0(); break; } case 3: { return new elasticity_finite_element_0(); break; } } return 0; } /// Create a new dof map for argument function i virtual ufc::dof_map* create_dof_map(unsigned int i) const { switch (i) { case 0: { return new elasticity_dof_map_2(); break; } case 1: { return new elasticity_dof_map_2(); break; } case 2: { return new elasticity_dof_map_0(); break; } case 3: { return new elasticity_dof_map_0(); break; } } return 0; } /// Create a new cell integral on sub domain i virtual ufc::cell_integral* create_cell_integral(unsigned int i) const { switch (i) { case 0: { return new elasticity_cell_integral_0_0(); break; } } return 0; } /// Create a new exterior facet integral on sub domain i virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const { return 0; } /// Create a new interior facet integral on sub domain i virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const { return 0; } }; /// This class defines the interface for the assembly of the global /// tensor corresponding to a form with r + n arguments, that is, a /// mapping /// /// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R /// /// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r /// global tensor A is defined by /// /// A = a(V1, V2, ..., Vr, w1, w2, ..., wn), /// /// where each argument Vj represents the application to the /// sequence of basis functions of Vj and w1, w2, ..., wn are given /// fixed functions (coefficients). class elasticity_form_1: public ufc::form { public: /// Constructor elasticity_form_1() : ufc::form() { // Do nothing } /// Destructor virtual ~elasticity_form_1() { // Do nothing } /// Return a string identifying the form virtual const char* signature() const { return "Form([Integral(IndexSum(Product(Indexed(Argument(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 0), MultiIndex((Index(0),), {Index(0): 3})), Indexed(Coefficient(VectorElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1, 3), 0), MultiIndex((Index(0),), {Index(0): 3}))), MultiIndex((Index(0),), {Index(0): 3})), Measure('cell', 0, None))])"; } /// Return the rank of the global tensor (r) virtual unsigned int rank() const { return 1; } /// Return the number of coefficients (n) virtual unsigned int num_coefficients() const { return 1; } /// Return the number of cell integrals virtual unsigned int num_cell_integrals() const { return 1; } /// Return the number of exterior facet integrals virtual unsigned int num_exterior_facet_integrals() const { return 0; } /// Return the number of interior facet integrals virtual unsigned int num_interior_facet_integrals() const { return 0; } /// Create a new finite element for argument function i virtual ufc::finite_element* create_finite_element(unsigned int i) const { switch (i) { case 0: { return new elasticity_finite_element_2(); break; } case 1: { return new elasticity_finite_element_2(); break; } } return 0; } /// Create a new dof map for argument function i virtual ufc::dof_map* create_dof_map(unsigned int i) const { switch (i) { case 0: { return new elasticity_dof_map_2(); break; } case 1: { return new elasticity_dof_map_2(); break; } } return 0; } /// Create a new cell integral on sub domain i virtual ufc::cell_integral* create_cell_integral(unsigned int i) const { switch (i) { case 0: { return new elasticity_cell_integral_1_0(); break; } } return 0; } /// Create a new exterior facet integral on sub domain i virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const { return 0; } /// Create a new interior facet integral on sub domain i virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const { return 0; } }; // DOLFIN wrappers // Standard library includes #include // DOLFIN includes #include #include #include #include #include #include #include namespace Elasticity { class CoefficientSpace_f: public dolfin::FunctionSpace { public: CoefficientSpace_f(const dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } CoefficientSpace_f(dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } CoefficientSpace_f(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } CoefficientSpace_f(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } ~CoefficientSpace_f() { } }; class CoefficientSpace_lmbda: public dolfin::FunctionSpace { public: CoefficientSpace_lmbda(const dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), mesh))) { // Do nothing } CoefficientSpace_lmbda(dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), mesh))) { // Do nothing } CoefficientSpace_lmbda(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), *mesh))) { // Do nothing } CoefficientSpace_lmbda(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), *mesh))) { // Do nothing } ~CoefficientSpace_lmbda() { } }; class CoefficientSpace_mu: public dolfin::FunctionSpace { public: CoefficientSpace_mu(const dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), mesh))) { // Do nothing } CoefficientSpace_mu(dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), mesh))) { // Do nothing } CoefficientSpace_mu(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), *mesh))) { // Do nothing } CoefficientSpace_mu(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_0()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_0()), *mesh))) { // Do nothing } ~CoefficientSpace_mu() { } }; class Form_0_FunctionSpace_0: public dolfin::FunctionSpace { public: Form_0_FunctionSpace_0(const dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } Form_0_FunctionSpace_0(dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } Form_0_FunctionSpace_0(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } Form_0_FunctionSpace_0(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } ~Form_0_FunctionSpace_0() { } }; class Form_0_FunctionSpace_1: public dolfin::FunctionSpace { public: Form_0_FunctionSpace_1(const dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } Form_0_FunctionSpace_1(dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } Form_0_FunctionSpace_1(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } Form_0_FunctionSpace_1(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } ~Form_0_FunctionSpace_1() { } }; typedef CoefficientSpace_mu Form_0_FunctionSpace_2; typedef CoefficientSpace_lmbda Form_0_FunctionSpace_3; class Form_0: public dolfin::Form { public: // Constructor Form_0(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1): dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1) { _function_spaces[0] = reference_to_no_delete_pointer(V0); _function_spaces[1] = reference_to_no_delete_pointer(V1); _ufc_form = boost::shared_ptr(new elasticity_form_0()); } // Constructor Form_0(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1, const dolfin::GenericFunction& mu, const dolfin::GenericFunction& lmbda): dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1) { _function_spaces[0] = reference_to_no_delete_pointer(V0); _function_spaces[1] = reference_to_no_delete_pointer(V1); this->mu = mu; this->lmbda = lmbda; _ufc_form = boost::shared_ptr(new elasticity_form_0()); } // Constructor Form_0(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1, boost::shared_ptr mu, boost::shared_ptr lmbda): dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1) { _function_spaces[0] = reference_to_no_delete_pointer(V0); _function_spaces[1] = reference_to_no_delete_pointer(V1); this->mu = *mu; this->lmbda = *lmbda; _ufc_form = boost::shared_ptr(new elasticity_form_0()); } // Constructor Form_0(boost::shared_ptr V0, boost::shared_ptr V1): dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1) { _function_spaces[0] = V0; _function_spaces[1] = V1; _ufc_form = boost::shared_ptr(new elasticity_form_0()); } // Constructor Form_0(boost::shared_ptr V0, boost::shared_ptr V1, const dolfin::GenericFunction& mu, const dolfin::GenericFunction& lmbda): dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1) { _function_spaces[0] = V0; _function_spaces[1] = V1; this->mu = mu; this->lmbda = lmbda; _ufc_form = boost::shared_ptr(new elasticity_form_0()); } // Constructor Form_0(boost::shared_ptr V0, boost::shared_ptr V1, boost::shared_ptr mu, boost::shared_ptr lmbda): dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1) { _function_spaces[0] = V0; _function_spaces[1] = V1; this->mu = *mu; this->lmbda = *lmbda; _ufc_form = boost::shared_ptr(new elasticity_form_0()); } // Destructor ~Form_0() {} /// Return the number of the coefficient with this name virtual dolfin::uint coefficient_number(const std::string& name) const { if (name == "mu") return 0; else if (name == "lmbda") return 1; dolfin::error("Invalid coefficient."); return 0; } /// Return the name of the coefficient with this number virtual std::string coefficient_name(dolfin::uint i) const { switch (i) { case 0: return "mu"; case 1: return "lmbda"; } dolfin::error("Invalid coefficient."); return "unnamed"; } // Typedefs typedef Form_0_FunctionSpace_0 TestSpace; typedef Form_0_FunctionSpace_1 TrialSpace; typedef Form_0_FunctionSpace_2 CoefficientSpace_mu; typedef Form_0_FunctionSpace_3 CoefficientSpace_lmbda; // Coefficients dolfin::CoefficientAssigner mu; dolfin::CoefficientAssigner lmbda; }; class Form_1_FunctionSpace_0: public dolfin::FunctionSpace { public: Form_1_FunctionSpace_0(const dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } Form_1_FunctionSpace_0(dolfin::Mesh& mesh): dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh), boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), mesh))) { // Do nothing } Form_1_FunctionSpace_0(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } Form_1_FunctionSpace_0(boost::shared_ptr mesh): dolfin::FunctionSpace(mesh, boost::shared_ptr(new dolfin::FiniteElement(boost::shared_ptr(new elasticity_finite_element_2()))), boost::shared_ptr(new dolfin::DofMap(boost::shared_ptr(new elasticity_dof_map_2()), *mesh))) { // Do nothing } ~Form_1_FunctionSpace_0() { } }; typedef CoefficientSpace_f Form_1_FunctionSpace_1; class Form_1: public dolfin::Form { public: // Constructor Form_1(const dolfin::FunctionSpace& V0): dolfin::Form(1, 1), f(*this, 0) { _function_spaces[0] = reference_to_no_delete_pointer(V0); _ufc_form = boost::shared_ptr(new elasticity_form_1()); } // Constructor Form_1(const dolfin::FunctionSpace& V0, const dolfin::GenericFunction& f): dolfin::Form(1, 1), f(*this, 0) { _function_spaces[0] = reference_to_no_delete_pointer(V0); this->f = f; _ufc_form = boost::shared_ptr(new elasticity_form_1()); } // Constructor Form_1(const dolfin::FunctionSpace& V0, boost::shared_ptr f): dolfin::Form(1, 1), f(*this, 0) { _function_spaces[0] = reference_to_no_delete_pointer(V0); this->f = *f; _ufc_form = boost::shared_ptr(new elasticity_form_1()); } // Constructor Form_1(boost::shared_ptr V0): dolfin::Form(1, 1), f(*this, 0) { _function_spaces[0] = V0; _ufc_form = boost::shared_ptr(new elasticity_form_1()); } // Constructor Form_1(boost::shared_ptr V0, const dolfin::GenericFunction& f): dolfin::Form(1, 1), f(*this, 0) { _function_spaces[0] = V0; this->f = f; _ufc_form = boost::shared_ptr(new elasticity_form_1()); } // Constructor Form_1(boost::shared_ptr V0, boost::shared_ptr f): dolfin::Form(1, 1), f(*this, 0) { _function_spaces[0] = V0; this->f = *f; _ufc_form = boost::shared_ptr(new elasticity_form_1()); } // Destructor ~Form_1() {} /// Return the number of the coefficient with this name virtual dolfin::uint coefficient_number(const std::string& name) const { if (name == "f") return 0; dolfin::error("Invalid coefficient."); return 0; } /// Return the name of the coefficient with this number virtual std::string coefficient_name(dolfin::uint i) const { switch (i) { case 0: return "f"; } dolfin::error("Invalid coefficient."); return "unnamed"; } // Typedefs typedef Form_1_FunctionSpace_0 TestSpace; typedef Form_1_FunctionSpace_1 CoefficientSpace_f; // Coefficients dolfin::CoefficientAssigner f; }; // Class typedefs typedef Form_0 BilinearForm; typedef Form_1 LinearForm; typedef Form_0::TestSpace FunctionSpace; } #endif