~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/elements/ffc_04.h

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
// This code conforms with the UFC specification version 1.0
2
 
// and was automatically generated by FFC version 0.4.3.
3
 
 
4
 
#ifndef __FFC_04_H
5
 
#define __FFC_04_H
6
 
 
7
 
#include <cmath>
8
 
#include <stdexcept>
9
 
#include <ufc.h>
10
 
 
11
 
/// This class defines the interface for a finite element.
12
 
 
13
 
class ffc_04_finite_element_0: public ufc::finite_element
14
 
{
15
 
public:
16
 
 
17
 
  /// Constructor
18
 
  ffc_04_finite_element_0() : ufc::finite_element()
19
 
  {
20
 
    // Do nothing
21
 
  }
22
 
 
23
 
  /// Destructor
24
 
  virtual ~ffc_04_finite_element_0()
25
 
  {
26
 
    // Do nothing
27
 
  }
28
 
 
29
 
  /// Return a string identifying the finite element
30
 
  virtual const char* signature() const
31
 
  {
32
 
    return "Lagrange finite element of degree 2 on a tetrahedron";
33
 
  }
34
 
 
35
 
  /// Return the cell shape
36
 
  virtual ufc::shape cell_shape() const
37
 
  {
38
 
    return ufc::tetrahedron;
39
 
  }
40
 
 
41
 
  /// Return the dimension of the finite element function space
42
 
  virtual unsigned int space_dimension() const
43
 
  {
44
 
    return 10;
45
 
  }
46
 
 
47
 
  /// Return the rank of the value space
48
 
  virtual unsigned int value_rank() const
49
 
  {
50
 
    return 0;
51
 
  }
52
 
 
53
 
  /// Return the dimension of the value space for axis i
54
 
  virtual unsigned int value_dimension(unsigned int i) const
55
 
  {
56
 
    return 1;
57
 
  }
58
 
 
59
 
  /// Evaluate basis function i at given point in cell
60
 
  virtual void evaluate_basis(unsigned int i,
61
 
                              double* values,
62
 
                              const double* coordinates,
63
 
                              const ufc::cell& c) const
64
 
  {
65
 
    // Extract vertex coordinates
66
 
    const double * const * element_coordinates = c.coordinates;
67
 
    
68
 
    // Compute Jacobian of affine map from reference cell
69
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
70
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
71
 
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
72
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
73
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
74
 
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
75
 
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
76
 
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
77
 
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
78
 
      
79
 
    // Compute sub determinants
80
 
    const double d00 = J_11*J_22 - J_12*J_21;
81
 
    const double d01 = J_12*J_20 - J_10*J_22;
82
 
    const double d02 = J_10*J_21 - J_11*J_20;
83
 
    
84
 
    const double d10 = J_02*J_21 - J_01*J_22;
85
 
    const double d11 = J_00*J_22 - J_02*J_20;
86
 
    const double d12 = J_01*J_20 - J_00*J_21;
87
 
    
88
 
    const double d20 = J_01*J_12 - J_02*J_11;
89
 
    const double d21 = J_02*J_10 - J_00*J_12;
90
 
    const double d22 = J_00*J_11 - J_01*J_10;
91
 
      
92
 
    // Compute determinant of Jacobian
93
 
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
94
 
    
95
 
    // Compute inverse of Jacobian
96
 
    
97
 
    // Compute constants
98
 
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
99
 
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
100
 
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
101
 
    
102
 
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
103
 
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
104
 
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
105
 
    
106
 
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
107
 
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
108
 
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
109
 
    
110
 
    // Get coordinates and map to the UFC reference element
111
 
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
112
 
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
113
 
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
114
 
    
115
 
    // Map coordinates to the reference cube
116
 
    if (std::abs(y + z - 1.0) < 1e-14)
117
 
      x = 1.0;
118
 
    else
119
 
      x = -2.0 * x/(y + z - 1.0) - 1.0;
120
 
    if (std::abs(z - 1.0) < 1e-14)
121
 
      y = -1.0;
122
 
    else
123
 
      y = 2.0 * y/(1.0 - z) - 1.0;
124
 
    z = 2.0 * z - 1.0;
125
 
    
126
 
    // Reset values
127
 
    *values = 0;
128
 
    
129
 
    // Map degree of freedom to element degree of freedom
130
 
    const unsigned int dof = i;
131
 
    
132
 
    // Generate scalings
133
 
    const double scalings_y_0 = 1;
134
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
135
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
136
 
    const double scalings_z_0 = 1;
137
 
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
138
 
    const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
139
 
    
140
 
    // Compute psitilde_a
141
 
    const double psitilde_a_0 = 1;
142
 
    const double psitilde_a_1 = x;
143
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
144
 
    
145
 
    // Compute psitilde_bs
146
 
    const double psitilde_bs_0_0 = 1;
147
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
148
 
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
149
 
    const double psitilde_bs_1_0 = 1;
150
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
151
 
    const double psitilde_bs_2_0 = 1;
152
 
    
153
 
    // Compute psitilde_cs
154
 
    const double psitilde_cs_00_0 = 1;
155
 
    const double psitilde_cs_00_1 = 2*z + 1;
156
 
    const double psitilde_cs_00_2 = 0.3125*psitilde_cs_00_1 + 1.875*z*psitilde_cs_00_1 - 0.5625*psitilde_cs_00_0;
157
 
    const double psitilde_cs_01_0 = 1;
158
 
    const double psitilde_cs_01_1 = 3*z + 2;
159
 
    const double psitilde_cs_02_0 = 1;
160
 
    const double psitilde_cs_10_0 = 1;
161
 
    const double psitilde_cs_10_1 = 3*z + 2;
162
 
    const double psitilde_cs_11_0 = 1;
163
 
    const double psitilde_cs_20_0 = 1;
164
 
    
165
 
    // Compute basisvalues
166
 
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
167
 
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
168
 
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
169
 
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
170
 
    const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
171
 
    const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
172
 
    const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
173
 
    const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
174
 
    const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
175
 
    const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
176
 
    
177
 
    // Table(s) of coefficients
178
 
    const static double coefficients0[10][10] = \
179
 
    {{-0.0577350269189626, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
180
 
    {-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
181
 
    {-0.0577350269189625, 0, 0.0702728368926306, -0.0248451997499977, 0, 0, 0.0872871560943969, 0, -0.0475131096733199, 0.0167984210226323},
182
 
    {-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
183
 
    {0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
184
 
    {0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
185
 
    {0.23094010767585, 0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0, 0.100790526135794, -0.0872871560943969, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
186
 
    {0.23094010767585, -0.121716123890037, -0.0702728368926306, 0.0993807989999906, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
187
 
    {0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
188
 
    {0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
189
 
    
190
 
    // Extract relevant coefficients
191
 
    const double coeff0_0 = coefficients0[dof][0];
192
 
    const double coeff0_1 = coefficients0[dof][1];
193
 
    const double coeff0_2 = coefficients0[dof][2];
194
 
    const double coeff0_3 = coefficients0[dof][3];
195
 
    const double coeff0_4 = coefficients0[dof][4];
196
 
    const double coeff0_5 = coefficients0[dof][5];
197
 
    const double coeff0_6 = coefficients0[dof][6];
198
 
    const double coeff0_7 = coefficients0[dof][7];
199
 
    const double coeff0_8 = coefficients0[dof][8];
200
 
    const double coeff0_9 = coefficients0[dof][9];
201
 
    
202
 
    // Compute value(s)
203
 
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9;
204
 
  }
205
 
 
206
 
  /// Evaluate all basis functions at given point in cell
207
 
  virtual void evaluate_basis_all(double* values,
208
 
                                  const double* coordinates,
209
 
                                  const ufc::cell& c) const
210
 
  {
211
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
212
 
  }
213
 
 
214
 
  /// Evaluate order n derivatives of basis function i at given point in cell
215
 
  virtual void evaluate_basis_derivatives(unsigned int i,
216
 
                                          unsigned int n,
217
 
                                          double* values,
218
 
                                          const double* coordinates,
219
 
                                          const ufc::cell& c) const
220
 
  {
221
 
    // Extract vertex coordinates
222
 
    const double * const * element_coordinates = c.coordinates;
223
 
    
224
 
    // Compute Jacobian of affine map from reference cell
225
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
226
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
227
 
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
228
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
229
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
230
 
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
231
 
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
232
 
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
233
 
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
234
 
      
235
 
    // Compute sub determinants
236
 
    const double d00 = J_11*J_22 - J_12*J_21;
237
 
    const double d01 = J_12*J_20 - J_10*J_22;
238
 
    const double d02 = J_10*J_21 - J_11*J_20;
239
 
    
240
 
    const double d10 = J_02*J_21 - J_01*J_22;
241
 
    const double d11 = J_00*J_22 - J_02*J_20;
242
 
    const double d12 = J_01*J_20 - J_00*J_21;
243
 
    
244
 
    const double d20 = J_01*J_12 - J_02*J_11;
245
 
    const double d21 = J_02*J_10 - J_00*J_12;
246
 
    const double d22 = J_00*J_11 - J_01*J_10;
247
 
      
248
 
    // Compute determinant of Jacobian
249
 
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
250
 
    
251
 
    // Compute inverse of Jacobian
252
 
    
253
 
    // Compute constants
254
 
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
255
 
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
256
 
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
257
 
    
258
 
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
259
 
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
260
 
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
261
 
    
262
 
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
263
 
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
264
 
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
265
 
    
266
 
    // Get coordinates and map to the UFC reference element
267
 
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
268
 
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
269
 
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
270
 
    
271
 
    // Map coordinates to the reference cube
272
 
    if (std::abs(y + z - 1.0) < 1e-14)
273
 
      x = 1.0;
274
 
    else
275
 
      x = -2.0 * x/(y + z - 1.0) - 1.0;
276
 
    if (std::abs(z - 1.0) < 1e-14)
277
 
      y = -1.0;
278
 
    else
279
 
      y = 2.0 * y/(1.0 - z) - 1.0;
280
 
    z = 2.0 * z - 1.0;
281
 
    
282
 
    // Compute number of derivatives
283
 
    unsigned int num_derivatives = 1;
284
 
    
285
 
    for (unsigned int j = 0; j < n; j++)
286
 
      num_derivatives *= 3;
287
 
    
288
 
    
289
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
290
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
291
 
        
292
 
    for (unsigned int j = 0; j < num_derivatives; j++)
293
 
    {
294
 
      combinations[j] = new unsigned int [n];
295
 
      for (unsigned int k = 0; k < n; k++)
296
 
        combinations[j][k] = 0;
297
 
    }
298
 
        
299
 
    // Generate combinations of derivatives
300
 
    for (unsigned int row = 1; row < num_derivatives; row++)
301
 
    {
302
 
      for (unsigned int num = 0; num < row; num++)
303
 
      {
304
 
        for (unsigned int col = n-1; col+1 > 0; col--)
305
 
        {
306
 
          if (combinations[row][col] + 1 > 2)
307
 
            combinations[row][col] = 0;
308
 
          else
309
 
          {
310
 
            combinations[row][col] += 1;
311
 
            break;
312
 
          }
313
 
        }
314
 
      }
315
 
    }
316
 
    
317
 
    // Compute inverse of Jacobian
318
 
    const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
319
 
    
320
 
    // Declare transformation matrix
321
 
    // Declare pointer to two dimensional array and initialise
322
 
    double **transform = new double *[num_derivatives];
323
 
        
324
 
    for (unsigned int j = 0; j < num_derivatives; j++)
325
 
    {
326
 
      transform[j] = new double [num_derivatives];
327
 
      for (unsigned int k = 0; k < num_derivatives; k++)
328
 
        transform[j][k] = 1;
329
 
    }
330
 
    
331
 
    // Construct transformation matrix
332
 
    for (unsigned int row = 0; row < num_derivatives; row++)
333
 
    {
334
 
      for (unsigned int col = 0; col < num_derivatives; col++)
335
 
      {
336
 
        for (unsigned int k = 0; k < n; k++)
337
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
338
 
      }
339
 
    }
340
 
    
341
 
    // Reset values
342
 
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
343
 
      values[j] = 0;
344
 
    
345
 
    // Map degree of freedom to element degree of freedom
346
 
    const unsigned int dof = i;
347
 
    
348
 
    // Generate scalings
349
 
    const double scalings_y_0 = 1;
350
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
351
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
352
 
    const double scalings_z_0 = 1;
353
 
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
354
 
    const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
355
 
    
356
 
    // Compute psitilde_a
357
 
    const double psitilde_a_0 = 1;
358
 
    const double psitilde_a_1 = x;
359
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
360
 
    
361
 
    // Compute psitilde_bs
362
 
    const double psitilde_bs_0_0 = 1;
363
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
364
 
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
365
 
    const double psitilde_bs_1_0 = 1;
366
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
367
 
    const double psitilde_bs_2_0 = 1;
368
 
    
369
 
    // Compute psitilde_cs
370
 
    const double psitilde_cs_00_0 = 1;
371
 
    const double psitilde_cs_00_1 = 2*z + 1;
372
 
    const double psitilde_cs_00_2 = 0.3125*psitilde_cs_00_1 + 1.875*z*psitilde_cs_00_1 - 0.5625*psitilde_cs_00_0;
373
 
    const double psitilde_cs_01_0 = 1;
374
 
    const double psitilde_cs_01_1 = 3*z + 2;
375
 
    const double psitilde_cs_02_0 = 1;
376
 
    const double psitilde_cs_10_0 = 1;
377
 
    const double psitilde_cs_10_1 = 3*z + 2;
378
 
    const double psitilde_cs_11_0 = 1;
379
 
    const double psitilde_cs_20_0 = 1;
380
 
    
381
 
    // Compute basisvalues
382
 
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
383
 
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
384
 
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
385
 
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
386
 
    const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
387
 
    const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
388
 
    const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
389
 
    const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
390
 
    const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
391
 
    const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
392
 
    
393
 
    // Table(s) of coefficients
394
 
    const static double coefficients0[10][10] = \
395
 
    {{-0.0577350269189626, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
396
 
    {-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
397
 
    {-0.0577350269189625, 0, 0.0702728368926306, -0.0248451997499977, 0, 0, 0.0872871560943969, 0, -0.0475131096733199, 0.0167984210226323},
398
 
    {-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
399
 
    {0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
400
 
    {0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
401
 
    {0.23094010767585, 0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0, 0.100790526135794, -0.0872871560943969, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
402
 
    {0.23094010767585, -0.121716123890037, -0.0702728368926306, 0.0993807989999906, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
403
 
    {0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
404
 
    {0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
405
 
    
406
 
    // Interesting (new) part
407
 
    // Tables of derivatives of the polynomial base (transpose)
408
 
    const static double dmats0[10][10] = \
409
 
    {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
410
 
    {6.32455532033676, 0, 0, 0, 0, 0, 0, 0, 0, 0},
411
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
412
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
413
 
    {0, 11.2249721603218, 0, 0, 0, 0, 0, 0, 0, 0},
414
 
    {4.58257569495584, 0, 8.36660026534076, -1.18321595661992, 0, 0, 0, 0, 0, 0},
415
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
416
 
    {3.74165738677394, 0, 0, 8.69482604771366, 0, 0, 0, 0, 0, 0},
417
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
418
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
419
 
    
420
 
    const static double dmats1[10][10] = \
421
 
    {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
422
 
    {3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
423
 
    {5.47722557505166, 0, 0, 0, 0, 0, 0, 0, 0, 0},
424
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
425
 
    {2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
426
 
    {2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
427
 
    {-2.64575131106459, 0, 9.66091783079296, 0.683130051063973, 0, 0, 0, 0, 0, 0},
428
 
    {1.87082869338697, 0, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
429
 
    {3.24037034920393, 0, 0, 7.52994023880668, 0, 0, 0, 0, 0, 0},
430
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
431
 
    
432
 
    const static double dmats2[10][10] = \
433
 
    {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
434
 
    {3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
435
 
    {1.82574185835055, 0, 0, 0, 0, 0, 0, 0, 0, 0},
436
 
    {5.16397779494322, 0, 0, 0, 0, 0, 0, 0, 0, 0},
437
 
    {2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
438
 
    {2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
439
 
    {1.3228756555323, 0, 3.86436713231718, -0.341565025531986, 0, 0, 0, 0, 0, 0},
440
 
    {1.87082869338697, 7.09929573971954, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
441
 
    {1.08012344973464, 0, 7.09929573971954, 2.50998007960223, 0, 0, 0, 0, 0, 0},
442
 
    {-3.81881307912987, 0, 0, 8.87411967464942, 0, 0, 0, 0, 0, 0}};
443
 
    
444
 
    // Compute reference derivatives
445
 
    // Declare pointer to array of derivatives on FIAT element
446
 
    double *derivatives = new double [num_derivatives];
447
 
    
448
 
    // Declare coefficients
449
 
    double coeff0_0 = 0;
450
 
    double coeff0_1 = 0;
451
 
    double coeff0_2 = 0;
452
 
    double coeff0_3 = 0;
453
 
    double coeff0_4 = 0;
454
 
    double coeff0_5 = 0;
455
 
    double coeff0_6 = 0;
456
 
    double coeff0_7 = 0;
457
 
    double coeff0_8 = 0;
458
 
    double coeff0_9 = 0;
459
 
    
460
 
    // Declare new coefficients
461
 
    double new_coeff0_0 = 0;
462
 
    double new_coeff0_1 = 0;
463
 
    double new_coeff0_2 = 0;
464
 
    double new_coeff0_3 = 0;
465
 
    double new_coeff0_4 = 0;
466
 
    double new_coeff0_5 = 0;
467
 
    double new_coeff0_6 = 0;
468
 
    double new_coeff0_7 = 0;
469
 
    double new_coeff0_8 = 0;
470
 
    double new_coeff0_9 = 0;
471
 
    
472
 
    // Loop possible derivatives
473
 
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
474
 
    {
475
 
      // Get values from coefficients array
476
 
      new_coeff0_0 = coefficients0[dof][0];
477
 
      new_coeff0_1 = coefficients0[dof][1];
478
 
      new_coeff0_2 = coefficients0[dof][2];
479
 
      new_coeff0_3 = coefficients0[dof][3];
480
 
      new_coeff0_4 = coefficients0[dof][4];
481
 
      new_coeff0_5 = coefficients0[dof][5];
482
 
      new_coeff0_6 = coefficients0[dof][6];
483
 
      new_coeff0_7 = coefficients0[dof][7];
484
 
      new_coeff0_8 = coefficients0[dof][8];
485
 
      new_coeff0_9 = coefficients0[dof][9];
486
 
    
487
 
      // Loop derivative order
488
 
      for (unsigned int j = 0; j < n; j++)
489
 
      {
490
 
        // Update old coefficients
491
 
        coeff0_0 = new_coeff0_0;
492
 
        coeff0_1 = new_coeff0_1;
493
 
        coeff0_2 = new_coeff0_2;
494
 
        coeff0_3 = new_coeff0_3;
495
 
        coeff0_4 = new_coeff0_4;
496
 
        coeff0_5 = new_coeff0_5;
497
 
        coeff0_6 = new_coeff0_6;
498
 
        coeff0_7 = new_coeff0_7;
499
 
        coeff0_8 = new_coeff0_8;
500
 
        coeff0_9 = new_coeff0_9;
501
 
    
502
 
        if(combinations[deriv_num][j] == 0)
503
 
        {
504
 
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0];
505
 
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1];
506
 
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2];
507
 
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3];
508
 
          new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4];
509
 
          new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5];
510
 
          new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6];
511
 
          new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7];
512
 
          new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8];
513
 
          new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9];
514
 
        }
515
 
        if(combinations[deriv_num][j] == 1)
516
 
        {
517
 
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0];
518
 
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1];
519
 
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2];
520
 
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3];
521
 
          new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4];
522
 
          new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5];
523
 
          new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6];
524
 
          new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7];
525
 
          new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8];
526
 
          new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9];
527
 
        }
528
 
        if(combinations[deriv_num][j] == 2)
529
 
        {
530
 
          new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0] + coeff0_4*dmats2[4][0] + coeff0_5*dmats2[5][0] + coeff0_6*dmats2[6][0] + coeff0_7*dmats2[7][0] + coeff0_8*dmats2[8][0] + coeff0_9*dmats2[9][0];
531
 
          new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1] + coeff0_4*dmats2[4][1] + coeff0_5*dmats2[5][1] + coeff0_6*dmats2[6][1] + coeff0_7*dmats2[7][1] + coeff0_8*dmats2[8][1] + coeff0_9*dmats2[9][1];
532
 
          new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2] + coeff0_4*dmats2[4][2] + coeff0_5*dmats2[5][2] + coeff0_6*dmats2[6][2] + coeff0_7*dmats2[7][2] + coeff0_8*dmats2[8][2] + coeff0_9*dmats2[9][2];
533
 
          new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3] + coeff0_4*dmats2[4][3] + coeff0_5*dmats2[5][3] + coeff0_6*dmats2[6][3] + coeff0_7*dmats2[7][3] + coeff0_8*dmats2[8][3] + coeff0_9*dmats2[9][3];
534
 
          new_coeff0_4 = coeff0_0*dmats2[0][4] + coeff0_1*dmats2[1][4] + coeff0_2*dmats2[2][4] + coeff0_3*dmats2[3][4] + coeff0_4*dmats2[4][4] + coeff0_5*dmats2[5][4] + coeff0_6*dmats2[6][4] + coeff0_7*dmats2[7][4] + coeff0_8*dmats2[8][4] + coeff0_9*dmats2[9][4];
535
 
          new_coeff0_5 = coeff0_0*dmats2[0][5] + coeff0_1*dmats2[1][5] + coeff0_2*dmats2[2][5] + coeff0_3*dmats2[3][5] + coeff0_4*dmats2[4][5] + coeff0_5*dmats2[5][5] + coeff0_6*dmats2[6][5] + coeff0_7*dmats2[7][5] + coeff0_8*dmats2[8][5] + coeff0_9*dmats2[9][5];
536
 
          new_coeff0_6 = coeff0_0*dmats2[0][6] + coeff0_1*dmats2[1][6] + coeff0_2*dmats2[2][6] + coeff0_3*dmats2[3][6] + coeff0_4*dmats2[4][6] + coeff0_5*dmats2[5][6] + coeff0_6*dmats2[6][6] + coeff0_7*dmats2[7][6] + coeff0_8*dmats2[8][6] + coeff0_9*dmats2[9][6];
537
 
          new_coeff0_7 = coeff0_0*dmats2[0][7] + coeff0_1*dmats2[1][7] + coeff0_2*dmats2[2][7] + coeff0_3*dmats2[3][7] + coeff0_4*dmats2[4][7] + coeff0_5*dmats2[5][7] + coeff0_6*dmats2[6][7] + coeff0_7*dmats2[7][7] + coeff0_8*dmats2[8][7] + coeff0_9*dmats2[9][7];
538
 
          new_coeff0_8 = coeff0_0*dmats2[0][8] + coeff0_1*dmats2[1][8] + coeff0_2*dmats2[2][8] + coeff0_3*dmats2[3][8] + coeff0_4*dmats2[4][8] + coeff0_5*dmats2[5][8] + coeff0_6*dmats2[6][8] + coeff0_7*dmats2[7][8] + coeff0_8*dmats2[8][8] + coeff0_9*dmats2[9][8];
539
 
          new_coeff0_9 = coeff0_0*dmats2[0][9] + coeff0_1*dmats2[1][9] + coeff0_2*dmats2[2][9] + coeff0_3*dmats2[3][9] + coeff0_4*dmats2[4][9] + coeff0_5*dmats2[5][9] + coeff0_6*dmats2[6][9] + coeff0_7*dmats2[7][9] + coeff0_8*dmats2[8][9] + coeff0_9*dmats2[9][9];
540
 
        }
541
 
    
542
 
      }
543
 
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
544
 
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9;
545
 
    }
546
 
    
547
 
    // Transform derivatives back to physical element
548
 
    for (unsigned int row = 0; row < num_derivatives; row++)
549
 
    {
550
 
      for (unsigned int col = 0; col < num_derivatives; col++)
551
 
      {
552
 
        values[row] += transform[row][col]*derivatives[col];
553
 
      }
554
 
    }
555
 
    // Delete pointer to array of derivatives on FIAT element
556
 
    delete [] derivatives;
557
 
    
558
 
    // Delete pointer to array of combinations of derivatives and transform
559
 
    for (unsigned int row = 0; row < num_derivatives; row++)
560
 
    {
561
 
      delete [] combinations[row];
562
 
      delete [] transform[row];
563
 
    }
564
 
    
565
 
    delete [] combinations;
566
 
    delete [] transform;
567
 
  }
568
 
 
569
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
570
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
571
 
                                              double* values,
572
 
                                              const double* coordinates,
573
 
                                              const ufc::cell& c) const
574
 
  {
575
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
576
 
  }
577
 
 
578
 
  /// Evaluate linear functional for dof i on the function f
579
 
  virtual double evaluate_dof(unsigned int i,
580
 
                              const ufc::function& f,
581
 
                              const ufc::cell& c) const
582
 
  {
583
 
    // The reference points, direction and weights:
584
 
    const static double X[10][1][3] = {{{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}, {{0, 0.5, 0.5}}, {{0.5, 0, 0.5}}, {{0.5, 0.5, 0}}, {{0, 0, 0.5}}, {{0, 0.5, 0}}, {{0.5, 0, 0}}};
585
 
    const static double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
586
 
    const static double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
587
 
    
588
 
    const double * const * x = c.coordinates;
589
 
    double result = 0.0;
590
 
    // Iterate over the points:
591
 
    // Evaluate basis functions for affine mapping
592
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
593
 
    const double w1 = X[i][0][0];
594
 
    const double w2 = X[i][0][1];
595
 
    const double w3 = X[i][0][2];
596
 
    
597
 
    // Compute affine mapping y = F(X)
598
 
    double y[3];
599
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
600
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
601
 
    y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
602
 
    
603
 
    // Evaluate function at physical points
604
 
    double values[1];
605
 
    f.evaluate(values, y, c);
606
 
    
607
 
    // Map function values using appropriate mapping
608
 
    // Affine map: Do nothing
609
 
    
610
 
    // Note that we do not map the weights (yet).
611
 
    
612
 
    // Take directional components
613
 
    for(int k = 0; k < 1; k++)
614
 
      result += values[k]*D[i][0][k];
615
 
    // Multiply by weights 
616
 
    result *= W[i][0];
617
 
    
618
 
    return result;
619
 
  }
620
 
 
621
 
  /// Evaluate linear functionals for all dofs on the function f
622
 
  virtual void evaluate_dofs(double* values,
623
 
                             const ufc::function& f,
624
 
                             const ufc::cell& c) const
625
 
  {
626
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
627
 
  }
628
 
 
629
 
  /// Interpolate vertex values from dof values
630
 
  virtual void interpolate_vertex_values(double* vertex_values,
631
 
                                         const double* dof_values,
632
 
                                         const ufc::cell& c) const
633
 
  {
634
 
    // Evaluate at vertices and use affine mapping
635
 
    vertex_values[0] = dof_values[0];
636
 
    vertex_values[1] = dof_values[1];
637
 
    vertex_values[2] = dof_values[2];
638
 
    vertex_values[3] = dof_values[3];
639
 
  }
640
 
 
641
 
  /// Return the number of sub elements (for a mixed element)
642
 
  virtual unsigned int num_sub_elements() const
643
 
  {
644
 
    return 1;
645
 
  }
646
 
 
647
 
  /// Create a new finite element for sub element i (for a mixed element)
648
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
649
 
  {
650
 
    return new ffc_04_finite_element_0();
651
 
  }
652
 
 
653
 
};
654
 
 
655
 
/// This class defines the interface for a local-to-global mapping of
656
 
/// degrees of freedom (dofs).
657
 
 
658
 
class ffc_04_dof_map_0: public ufc::dof_map
659
 
{
660
 
private:
661
 
 
662
 
  unsigned int __global_dimension;
663
 
 
664
 
public:
665
 
 
666
 
  /// Constructor
667
 
  ffc_04_dof_map_0() : ufc::dof_map()
668
 
  {
669
 
    __global_dimension = 0;
670
 
  }
671
 
 
672
 
  /// Destructor
673
 
  virtual ~ffc_04_dof_map_0()
674
 
  {
675
 
    // Do nothing
676
 
  }
677
 
 
678
 
  /// Return a string identifying the dof map
679
 
  virtual const char* signature() const
680
 
  {
681
 
    return "FFC dof map for Lagrange finite element of degree 2 on a tetrahedron";
682
 
  }
683
 
 
684
 
  /// Return true iff mesh entities of topological dimension d are needed
685
 
  virtual bool needs_mesh_entities(unsigned int d) const
686
 
  {
687
 
    switch ( d )
688
 
    {
689
 
    case 0:
690
 
      return true;
691
 
      break;
692
 
    case 1:
693
 
      return true;
694
 
      break;
695
 
    case 2:
696
 
      return false;
697
 
      break;
698
 
    case 3:
699
 
      return false;
700
 
      break;
701
 
    }
702
 
    return false;
703
 
  }
704
 
 
705
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
706
 
  virtual bool init_mesh(const ufc::mesh& m)
707
 
  {
708
 
    __global_dimension = m.num_entities[0] + m.num_entities[1];
709
 
    return false;
710
 
  }
711
 
 
712
 
  /// Initialize dof map for given cell
713
 
  virtual void init_cell(const ufc::mesh& m,
714
 
                         const ufc::cell& c)
715
 
  {
716
 
    // Do nothing
717
 
  }
718
 
 
719
 
  /// Finish initialization of dof map for cells
720
 
  virtual void init_cell_finalize()
721
 
  {
722
 
    // Do nothing
723
 
  }
724
 
 
725
 
  /// Return the dimension of the global finite element function space
726
 
  virtual unsigned int global_dimension() const
727
 
  {
728
 
    return __global_dimension;
729
 
  }
730
 
 
731
 
  /// Return the dimension of the local finite element function space
732
 
  virtual unsigned int local_dimension() const
733
 
  {
734
 
    return 10;
735
 
  }
736
 
 
737
 
  // Return the geometric dimension of the coordinates this dof map provides
738
 
  virtual unsigned int geometric_dimension() const
739
 
  {
740
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
741
 
  }
742
 
 
743
 
  /// Return the number of dofs on each cell facet
744
 
  virtual unsigned int num_facet_dofs() const
745
 
  {
746
 
    return 6;
747
 
  }
748
 
 
749
 
  /// Return the number of dofs associated with each cell entity of dimension d
750
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
751
 
  {
752
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
753
 
  }
754
 
 
755
 
  /// Tabulate the local-to-global mapping of dofs on a cell
756
 
  virtual void tabulate_dofs(unsigned int* dofs,
757
 
                             const ufc::mesh& m,
758
 
                             const ufc::cell& c) const
759
 
  {
760
 
    dofs[0] = c.entity_indices[0][0];
761
 
    dofs[1] = c.entity_indices[0][1];
762
 
    dofs[2] = c.entity_indices[0][2];
763
 
    dofs[3] = c.entity_indices[0][3];
764
 
    unsigned int offset = m.num_entities[0];
765
 
    dofs[4] = offset + c.entity_indices[1][0];
766
 
    dofs[5] = offset + c.entity_indices[1][1];
767
 
    dofs[6] = offset + c.entity_indices[1][2];
768
 
    dofs[7] = offset + c.entity_indices[1][3];
769
 
    dofs[8] = offset + c.entity_indices[1][4];
770
 
    dofs[9] = offset + c.entity_indices[1][5];
771
 
  }
772
 
 
773
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
774
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
775
 
                                   unsigned int facet) const
776
 
  {
777
 
    switch ( facet )
778
 
    {
779
 
    case 0:
780
 
      dofs[0] = 1;
781
 
      dofs[1] = 2;
782
 
      dofs[2] = 3;
783
 
      dofs[3] = 4;
784
 
      dofs[4] = 5;
785
 
      dofs[5] = 6;
786
 
      break;
787
 
    case 1:
788
 
      dofs[0] = 0;
789
 
      dofs[1] = 2;
790
 
      dofs[2] = 3;
791
 
      dofs[3] = 4;
792
 
      dofs[4] = 7;
793
 
      dofs[5] = 8;
794
 
      break;
795
 
    case 2:
796
 
      dofs[0] = 0;
797
 
      dofs[1] = 1;
798
 
      dofs[2] = 3;
799
 
      dofs[3] = 5;
800
 
      dofs[4] = 7;
801
 
      dofs[5] = 9;
802
 
      break;
803
 
    case 3:
804
 
      dofs[0] = 0;
805
 
      dofs[1] = 1;
806
 
      dofs[2] = 2;
807
 
      dofs[3] = 6;
808
 
      dofs[4] = 8;
809
 
      dofs[5] = 9;
810
 
      break;
811
 
    }
812
 
  }
813
 
 
814
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
815
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
816
 
                                    unsigned int d, unsigned int i) const
817
 
  {
818
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
819
 
  }
820
 
 
821
 
  /// Tabulate the coordinates of all dofs on a cell
822
 
  virtual void tabulate_coordinates(double** coordinates,
823
 
                                    const ufc::cell& c) const
824
 
  {
825
 
    const double * const * x = c.coordinates;
826
 
    coordinates[0][0] = x[0][0];
827
 
    coordinates[0][1] = x[0][1];
828
 
    coordinates[0][2] = x[0][2];
829
 
    coordinates[1][0] = x[1][0];
830
 
    coordinates[1][1] = x[1][1];
831
 
    coordinates[1][2] = x[1][2];
832
 
    coordinates[2][0] = x[2][0];
833
 
    coordinates[2][1] = x[2][1];
834
 
    coordinates[2][2] = x[2][2];
835
 
    coordinates[3][0] = x[3][0];
836
 
    coordinates[3][1] = x[3][1];
837
 
    coordinates[3][2] = x[3][2];
838
 
    coordinates[4][0] = 0.5*x[2][0] + 0.5*x[3][0];
839
 
    coordinates[4][1] = 0.5*x[2][1] + 0.5*x[3][1];
840
 
    coordinates[4][2] = 0.5*x[2][2] + 0.5*x[3][2];
841
 
    coordinates[5][0] = 0.5*x[1][0] + 0.5*x[3][0];
842
 
    coordinates[5][1] = 0.5*x[1][1] + 0.5*x[3][1];
843
 
    coordinates[5][2] = 0.5*x[1][2] + 0.5*x[3][2];
844
 
    coordinates[6][0] = 0.5*x[1][0] + 0.5*x[2][0];
845
 
    coordinates[6][1] = 0.5*x[1][1] + 0.5*x[2][1];
846
 
    coordinates[6][2] = 0.5*x[1][2] + 0.5*x[2][2];
847
 
    coordinates[7][0] = 0.5*x[0][0] + 0.5*x[3][0];
848
 
    coordinates[7][1] = 0.5*x[0][1] + 0.5*x[3][1];
849
 
    coordinates[7][2] = 0.5*x[0][2] + 0.5*x[3][2];
850
 
    coordinates[8][0] = 0.5*x[0][0] + 0.5*x[2][0];
851
 
    coordinates[8][1] = 0.5*x[0][1] + 0.5*x[2][1];
852
 
    coordinates[8][2] = 0.5*x[0][2] + 0.5*x[2][2];
853
 
    coordinates[9][0] = 0.5*x[0][0] + 0.5*x[1][0];
854
 
    coordinates[9][1] = 0.5*x[0][1] + 0.5*x[1][1];
855
 
    coordinates[9][2] = 0.5*x[0][2] + 0.5*x[1][2];
856
 
  }
857
 
 
858
 
  /// Return the number of sub dof maps (for a mixed element)
859
 
  virtual unsigned int num_sub_dof_maps() const
860
 
  {
861
 
    return 1;
862
 
  }
863
 
 
864
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
865
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
866
 
  {
867
 
    return new ffc_04_dof_map_0();
868
 
  }
869
 
 
870
 
};
871
 
 
872
 
#endif