11
11
#include <stdexcept>
15
/// This class defines the interface for a finite element.
17
class UFC_ProjectionBilinearForm_finite_element_0: public ufc::finite_element
22
UFC_ProjectionBilinearForm_finite_element_0() : ufc::finite_element()
28
virtual ~UFC_ProjectionBilinearForm_finite_element_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "FiniteElement('Lagrange', 'tetrahedron', 2)";
39
/// Return the cell shape
40
virtual ufc::shape cell_shape() const
42
return ufc::tetrahedron;
45
/// Return the dimension of the finite element function space
46
virtual unsigned int space_dimension() const
51
/// Return the rank of the value space
52
virtual unsigned int value_rank() const
57
/// Return the dimension of the value space for axis i
58
virtual unsigned int value_dimension(unsigned int i) const
63
/// Evaluate basis function i at given point in cell
64
virtual void evaluate_basis(unsigned int i,
66
const double* coordinates,
67
const ufc::cell& c) const
69
// Extract vertex coordinates
70
const double * const * element_coordinates = c.coordinates;
72
// Compute Jacobian of affine map from reference cell
73
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
74
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
75
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
76
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
77
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
78
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
79
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
80
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
81
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
83
// Compute sub determinants
84
const double d00 = J_11*J_22 - J_12*J_21;
85
const double d01 = J_12*J_20 - J_10*J_22;
86
const double d02 = J_10*J_21 - J_11*J_20;
88
const double d10 = J_02*J_21 - J_01*J_22;
89
const double d11 = J_00*J_22 - J_02*J_20;
90
const double d12 = J_01*J_20 - J_00*J_21;
92
const double d20 = J_01*J_12 - J_02*J_11;
93
const double d21 = J_02*J_10 - J_00*J_12;
94
const double d22 = J_00*J_11 - J_01*J_10;
96
// Compute determinant of Jacobian
97
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
99
// Compute inverse of Jacobian
102
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
103
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
104
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
106
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
107
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
108
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
110
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
111
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
112
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
114
// Get coordinates and map to the UFC reference element
115
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
116
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
117
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
119
// Map coordinates to the reference cube
120
if (std::abs(y + z - 1.0) < 1e-14)
123
x = -2.0 * x/(y + z - 1.0) - 1.0;
124
if (std::abs(z - 1.0) < 1e-14)
127
y = 2.0 * y/(1.0 - z) - 1.0;
133
// Map degree of freedom to element degree of freedom
134
const unsigned int dof = i;
137
const double scalings_y_0 = 1;
138
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
139
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
140
const double scalings_z_0 = 1;
141
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
142
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
144
// Compute psitilde_a
145
const double psitilde_a_0 = 1;
146
const double psitilde_a_1 = x;
147
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
149
// Compute psitilde_bs
150
const double psitilde_bs_0_0 = 1;
151
const double psitilde_bs_0_1 = 1.5*y + 0.5;
152
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;
153
const double psitilde_bs_1_0 = 1;
154
const double psitilde_bs_1_1 = 2.5*y + 1.5;
155
const double psitilde_bs_2_0 = 1;
157
// Compute psitilde_cs
158
const double psitilde_cs_00_0 = 1;
159
const double psitilde_cs_00_1 = 2*z + 1;
160
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;
161
const double psitilde_cs_01_0 = 1;
162
const double psitilde_cs_01_1 = 3*z + 2;
163
const double psitilde_cs_02_0 = 1;
164
const double psitilde_cs_10_0 = 1;
165
const double psitilde_cs_10_1 = 3*z + 2;
166
const double psitilde_cs_11_0 = 1;
167
const double psitilde_cs_20_0 = 1;
169
// Compute basisvalues
170
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
171
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
172
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
173
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
174
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
175
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
176
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
177
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
178
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
179
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
181
// Table(s) of coefficients
182
const static double coefficients0[10][10] = \
183
{{-0.0577350269189625, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
184
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499976, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
185
{-0.0577350269189626, 0, 0.0702728368926306, -0.0248451997499976, 0, 0, 0.087287156094397, 0, -0.0475131096733199, 0.0167984210226323},
186
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
187
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
188
{0.23094010767585, 0.121716123890037, -0.0702728368926306, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
189
{0.23094010767585, 0.121716123890037, 0.0702728368926307, -0.0993807989999907, 0, 0.100790526135794, -0.087287156094397, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
190
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999906, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
191
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
192
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
194
// Extract relevant coefficients
195
const double coeff0_0 = coefficients0[dof][0];
196
const double coeff0_1 = coefficients0[dof][1];
197
const double coeff0_2 = coefficients0[dof][2];
198
const double coeff0_3 = coefficients0[dof][3];
199
const double coeff0_4 = coefficients0[dof][4];
200
const double coeff0_5 = coefficients0[dof][5];
201
const double coeff0_6 = coefficients0[dof][6];
202
const double coeff0_7 = coefficients0[dof][7];
203
const double coeff0_8 = coefficients0[dof][8];
204
const double coeff0_9 = coefficients0[dof][9];
207
*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;
210
/// Evaluate all basis functions at given point in cell
211
virtual void evaluate_basis_all(double* values,
212
const double* coordinates,
213
const ufc::cell& c) const
215
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
218
/// Evaluate order n derivatives of basis function i at given point in cell
219
virtual void evaluate_basis_derivatives(unsigned int i,
222
const double* coordinates,
223
const ufc::cell& c) const
225
// Extract vertex coordinates
226
const double * const * element_coordinates = c.coordinates;
228
// Compute Jacobian of affine map from reference cell
229
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
230
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
231
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
232
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
233
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
234
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
235
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
236
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
237
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
239
// Compute sub determinants
240
const double d00 = J_11*J_22 - J_12*J_21;
241
const double d01 = J_12*J_20 - J_10*J_22;
242
const double d02 = J_10*J_21 - J_11*J_20;
244
const double d10 = J_02*J_21 - J_01*J_22;
245
const double d11 = J_00*J_22 - J_02*J_20;
246
const double d12 = J_01*J_20 - J_00*J_21;
248
const double d20 = J_01*J_12 - J_02*J_11;
249
const double d21 = J_02*J_10 - J_00*J_12;
250
const double d22 = J_00*J_11 - J_01*J_10;
252
// Compute determinant of Jacobian
253
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
255
// Compute inverse of Jacobian
258
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
259
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
260
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
262
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
263
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
264
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
266
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
267
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
268
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
270
// Get coordinates and map to the UFC reference element
271
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
272
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
273
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
275
// Map coordinates to the reference cube
276
if (std::abs(y + z - 1.0) < 1e-14)
279
x = -2.0 * x/(y + z - 1.0) - 1.0;
280
if (std::abs(z - 1.0) < 1e-14)
283
y = 2.0 * y/(1.0 - z) - 1.0;
286
// Compute number of derivatives
287
unsigned int num_derivatives = 1;
289
for (unsigned int j = 0; j < n; j++)
290
num_derivatives *= 3;
293
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
294
unsigned int **combinations = new unsigned int *[num_derivatives];
296
for (unsigned int j = 0; j < num_derivatives; j++)
298
combinations[j] = new unsigned int [n];
299
for (unsigned int k = 0; k < n; k++)
300
combinations[j][k] = 0;
303
// Generate combinations of derivatives
304
for (unsigned int row = 1; row < num_derivatives; row++)
306
for (unsigned int num = 0; num < row; num++)
308
for (unsigned int col = n-1; col+1 > 0; col--)
310
if (combinations[row][col] + 1 > 2)
311
combinations[row][col] = 0;
314
combinations[row][col] += 1;
321
// Compute inverse of Jacobian
322
const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
324
// Declare transformation matrix
325
// Declare pointer to two dimensional array and initialise
326
double **transform = new double *[num_derivatives];
328
for (unsigned int j = 0; j < num_derivatives; j++)
330
transform[j] = new double [num_derivatives];
331
for (unsigned int k = 0; k < num_derivatives; k++)
335
// Construct transformation matrix
336
for (unsigned int row = 0; row < num_derivatives; row++)
338
for (unsigned int col = 0; col < num_derivatives; col++)
340
for (unsigned int k = 0; k < n; k++)
341
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
346
for (unsigned int j = 0; j < 1*num_derivatives; j++)
349
// Map degree of freedom to element degree of freedom
350
const unsigned int dof = i;
353
const double scalings_y_0 = 1;
354
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
355
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
356
const double scalings_z_0 = 1;
357
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
358
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
360
// Compute psitilde_a
361
const double psitilde_a_0 = 1;
362
const double psitilde_a_1 = x;
363
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
365
// Compute psitilde_bs
366
const double psitilde_bs_0_0 = 1;
367
const double psitilde_bs_0_1 = 1.5*y + 0.5;
368
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;
369
const double psitilde_bs_1_0 = 1;
370
const double psitilde_bs_1_1 = 2.5*y + 1.5;
371
const double psitilde_bs_2_0 = 1;
373
// Compute psitilde_cs
374
const double psitilde_cs_00_0 = 1;
375
const double psitilde_cs_00_1 = 2*z + 1;
376
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;
377
const double psitilde_cs_01_0 = 1;
378
const double psitilde_cs_01_1 = 3*z + 2;
379
const double psitilde_cs_02_0 = 1;
380
const double psitilde_cs_10_0 = 1;
381
const double psitilde_cs_10_1 = 3*z + 2;
382
const double psitilde_cs_11_0 = 1;
383
const double psitilde_cs_20_0 = 1;
385
// Compute basisvalues
386
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
387
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
388
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
389
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
390
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
391
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
392
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
393
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
394
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
395
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
397
// Table(s) of coefficients
398
const static double coefficients0[10][10] = \
399
{{-0.0577350269189625, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
400
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499976, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
401
{-0.0577350269189626, 0, 0.0702728368926306, -0.0248451997499976, 0, 0, 0.087287156094397, 0, -0.0475131096733199, 0.0167984210226323},
402
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
403
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
404
{0.23094010767585, 0.121716123890037, -0.0702728368926306, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
405
{0.23094010767585, 0.121716123890037, 0.0702728368926307, -0.0993807989999907, 0, 0.100790526135794, -0.087287156094397, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
406
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999906, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
407
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
408
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
410
// Interesting (new) part
411
// Tables of derivatives of the polynomial base (transpose)
412
const static double dmats0[10][10] = \
413
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
414
{6.32455532033676, 0, 0, 0, 0, 0, 0, 0, 0, 0},
415
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
416
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
417
{0, 11.2249721603218, 0, 0, 0, 0, 0, 0, 0, 0},
418
{4.58257569495584, 0, 8.36660026534076, -1.18321595661992, 0, 0, 0, 0, 0, 0},
419
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
420
{3.74165738677394, 0, 0, 8.69482604771366, 0, 0, 0, 0, 0, 0},
421
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
422
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
424
const static double dmats1[10][10] = \
425
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
426
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
427
{5.47722557505166, 0, 0, 0, 0, 0, 0, 0, 0, 0},
428
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
429
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825974, 0, 0, 0, 0, 0, 0},
430
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
431
{-2.64575131106459, 0, 9.66091783079296, 0.683130051063974, 0, 0, 0, 0, 0, 0},
432
{1.87082869338697, 0, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
433
{3.24037034920393, 0, 0, 7.52994023880668, 0, 0, 0, 0, 0, 0},
434
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
436
const static double dmats2[10][10] = \
437
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
438
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
439
{1.82574185835055, 0, 0, 0, 0, 0, 0, 0, 0, 0},
440
{5.16397779494322, 0, 0, 0, 0, 0, 0, 0, 0, 0},
441
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
442
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
443
{1.32287565553229, 0, 3.86436713231718, -0.341565025531986, 0, 0, 0, 0, 0, 0},
444
{1.87082869338697, 7.09929573971954, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
445
{1.08012344973464, 0, 7.09929573971954, 2.50998007960223, 0, 0, 0, 0, 0, 0},
446
{-3.81881307912986, 0, 0, 8.87411967464942, 0, 0, 0, 0, 0, 0}};
448
// Compute reference derivatives
449
// Declare pointer to array of derivatives on FIAT element
450
double *derivatives = new double [num_derivatives];
452
// Declare coefficients
464
// Declare new coefficients
465
double new_coeff0_0 = 0;
466
double new_coeff0_1 = 0;
467
double new_coeff0_2 = 0;
468
double new_coeff0_3 = 0;
469
double new_coeff0_4 = 0;
470
double new_coeff0_5 = 0;
471
double new_coeff0_6 = 0;
472
double new_coeff0_7 = 0;
473
double new_coeff0_8 = 0;
474
double new_coeff0_9 = 0;
476
// Loop possible derivatives
477
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
479
// Get values from coefficients array
480
new_coeff0_0 = coefficients0[dof][0];
481
new_coeff0_1 = coefficients0[dof][1];
482
new_coeff0_2 = coefficients0[dof][2];
483
new_coeff0_3 = coefficients0[dof][3];
484
new_coeff0_4 = coefficients0[dof][4];
485
new_coeff0_5 = coefficients0[dof][5];
486
new_coeff0_6 = coefficients0[dof][6];
487
new_coeff0_7 = coefficients0[dof][7];
488
new_coeff0_8 = coefficients0[dof][8];
489
new_coeff0_9 = coefficients0[dof][9];
491
// Loop derivative order
492
for (unsigned int j = 0; j < n; j++)
494
// Update old coefficients
495
coeff0_0 = new_coeff0_0;
496
coeff0_1 = new_coeff0_1;
497
coeff0_2 = new_coeff0_2;
498
coeff0_3 = new_coeff0_3;
499
coeff0_4 = new_coeff0_4;
500
coeff0_5 = new_coeff0_5;
501
coeff0_6 = new_coeff0_6;
502
coeff0_7 = new_coeff0_7;
503
coeff0_8 = new_coeff0_8;
504
coeff0_9 = new_coeff0_9;
506
if(combinations[deriv_num][j] == 0)
508
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];
509
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];
510
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];
511
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];
512
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];
513
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];
514
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];
515
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];
516
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];
517
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];
519
if(combinations[deriv_num][j] == 1)
521
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];
522
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];
523
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];
524
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];
525
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];
526
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];
527
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];
528
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];
529
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];
530
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];
532
if(combinations[deriv_num][j] == 2)
534
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];
535
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];
536
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];
537
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];
538
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];
539
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];
540
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];
541
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];
542
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];
543
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];
547
// Compute derivatives on reference element as dot product of coefficients and basisvalues
548
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;
551
// Transform derivatives back to physical element
552
for (unsigned int row = 0; row < num_derivatives; row++)
554
for (unsigned int col = 0; col < num_derivatives; col++)
556
values[row] += transform[row][col]*derivatives[col];
559
// Delete pointer to array of derivatives on FIAT element
560
delete [] derivatives;
562
// Delete pointer to array of combinations of derivatives and transform
563
for (unsigned int row = 0; row < num_derivatives; row++)
565
delete [] combinations[row];
566
delete [] transform[row];
569
delete [] combinations;
573
/// Evaluate order n derivatives of all basis functions at given point in cell
574
virtual void evaluate_basis_derivatives_all(unsigned int n,
576
const double* coordinates,
577
const ufc::cell& c) const
579
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
582
/// Evaluate linear functional for dof i on the function f
583
virtual double evaluate_dof(unsigned int i,
584
const ufc::function& f,
585
const ufc::cell& c) const
587
// The reference points, direction and weights:
588
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}}};
589
const static double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
590
const static double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
592
const double * const * x = c.coordinates;
594
// Iterate over the points:
595
// Evaluate basis functions for affine mapping
596
const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
597
const double w1 = X[i][0][0];
598
const double w2 = X[i][0][1];
599
const double w3 = X[i][0][2];
601
// Compute affine mapping y = F(X)
603
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
604
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
605
y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
607
// Evaluate function at physical points
609
f.evaluate(values, y, c);
611
// Map function values using appropriate mapping
612
// Affine map: Do nothing
614
// Note that we do not map the weights (yet).
616
// Take directional components
617
for(int k = 0; k < 1; k++)
618
result += values[k]*D[i][0][k];
619
// Multiply by weights
625
/// Evaluate linear functionals for all dofs on the function f
626
virtual void evaluate_dofs(double* values,
627
const ufc::function& f,
628
const ufc::cell& c) const
630
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
633
/// Interpolate vertex values from dof values
634
virtual void interpolate_vertex_values(double* vertex_values,
635
const double* dof_values,
636
const ufc::cell& c) const
638
// Evaluate at vertices and use affine mapping
639
vertex_values[0] = dof_values[0];
640
vertex_values[1] = dof_values[1];
641
vertex_values[2] = dof_values[2];
642
vertex_values[3] = dof_values[3];
645
/// Return the number of sub elements (for a mixed element)
646
virtual unsigned int num_sub_elements() const
651
/// Create a new finite element for sub element i (for a mixed element)
652
virtual ufc::finite_element* create_sub_element(unsigned int i) const
654
return new UFC_ProjectionBilinearForm_finite_element_0();
659
/// This class defines the interface for a finite element.
661
class UFC_ProjectionBilinearForm_finite_element_1: public ufc::finite_element
666
UFC_ProjectionBilinearForm_finite_element_1() : ufc::finite_element()
672
virtual ~UFC_ProjectionBilinearForm_finite_element_1()
677
/// Return a string identifying the finite element
678
virtual const char* signature() const
680
return "FiniteElement('Lagrange', 'tetrahedron', 2)";
683
/// Return the cell shape
684
virtual ufc::shape cell_shape() const
686
return ufc::tetrahedron;
689
/// Return the dimension of the finite element function space
690
virtual unsigned int space_dimension() const
695
/// Return the rank of the value space
696
virtual unsigned int value_rank() const
701
/// Return the dimension of the value space for axis i
702
virtual unsigned int value_dimension(unsigned int i) const
707
/// Evaluate basis function i at given point in cell
708
virtual void evaluate_basis(unsigned int i,
710
const double* coordinates,
711
const ufc::cell& c) const
713
// Extract vertex coordinates
714
const double * const * element_coordinates = c.coordinates;
716
// Compute Jacobian of affine map from reference cell
717
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
718
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
719
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
720
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
721
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
722
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
723
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
724
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
725
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
727
// Compute sub determinants
728
const double d00 = J_11*J_22 - J_12*J_21;
729
const double d01 = J_12*J_20 - J_10*J_22;
730
const double d02 = J_10*J_21 - J_11*J_20;
732
const double d10 = J_02*J_21 - J_01*J_22;
733
const double d11 = J_00*J_22 - J_02*J_20;
734
const double d12 = J_01*J_20 - J_00*J_21;
736
const double d20 = J_01*J_12 - J_02*J_11;
737
const double d21 = J_02*J_10 - J_00*J_12;
738
const double d22 = J_00*J_11 - J_01*J_10;
740
// Compute determinant of Jacobian
741
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
743
// Compute inverse of Jacobian
746
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
747
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
748
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
750
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
751
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
752
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
754
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
755
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
756
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
758
// Get coordinates and map to the UFC reference element
759
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
760
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
761
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
763
// Map coordinates to the reference cube
764
if (std::abs(y + z - 1.0) < 1e-14)
767
x = -2.0 * x/(y + z - 1.0) - 1.0;
768
if (std::abs(z - 1.0) < 1e-14)
771
y = 2.0 * y/(1.0 - z) - 1.0;
777
// Map degree of freedom to element degree of freedom
778
const unsigned int dof = i;
781
const double scalings_y_0 = 1;
782
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
783
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
784
const double scalings_z_0 = 1;
785
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
786
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
788
// Compute psitilde_a
789
const double psitilde_a_0 = 1;
790
const double psitilde_a_1 = x;
791
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
793
// Compute psitilde_bs
794
const double psitilde_bs_0_0 = 1;
795
const double psitilde_bs_0_1 = 1.5*y + 0.5;
796
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;
797
const double psitilde_bs_1_0 = 1;
798
const double psitilde_bs_1_1 = 2.5*y + 1.5;
799
const double psitilde_bs_2_0 = 1;
801
// Compute psitilde_cs
802
const double psitilde_cs_00_0 = 1;
803
const double psitilde_cs_00_1 = 2*z + 1;
804
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;
805
const double psitilde_cs_01_0 = 1;
806
const double psitilde_cs_01_1 = 3*z + 2;
807
const double psitilde_cs_02_0 = 1;
808
const double psitilde_cs_10_0 = 1;
809
const double psitilde_cs_10_1 = 3*z + 2;
810
const double psitilde_cs_11_0 = 1;
811
const double psitilde_cs_20_0 = 1;
813
// Compute basisvalues
814
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
815
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
816
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
817
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
818
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
819
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
820
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
821
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
822
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
823
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
825
// Table(s) of coefficients
826
const static double coefficients0[10][10] = \
827
{{-0.0577350269189625, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
828
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499976, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
829
{-0.0577350269189626, 0, 0.0702728368926306, -0.0248451997499976, 0, 0, 0.087287156094397, 0, -0.0475131096733199, 0.0167984210226323},
830
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
831
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
832
{0.23094010767585, 0.121716123890037, -0.0702728368926306, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
833
{0.23094010767585, 0.121716123890037, 0.0702728368926307, -0.0993807989999907, 0, 0.100790526135794, -0.087287156094397, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
834
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999906, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
835
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
836
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
838
// Extract relevant coefficients
839
const double coeff0_0 = coefficients0[dof][0];
840
const double coeff0_1 = coefficients0[dof][1];
841
const double coeff0_2 = coefficients0[dof][2];
842
const double coeff0_3 = coefficients0[dof][3];
843
const double coeff0_4 = coefficients0[dof][4];
844
const double coeff0_5 = coefficients0[dof][5];
845
const double coeff0_6 = coefficients0[dof][6];
846
const double coeff0_7 = coefficients0[dof][7];
847
const double coeff0_8 = coefficients0[dof][8];
848
const double coeff0_9 = coefficients0[dof][9];
851
*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;
854
/// Evaluate all basis functions at given point in cell
855
virtual void evaluate_basis_all(double* values,
856
const double* coordinates,
857
const ufc::cell& c) const
859
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
862
/// Evaluate order n derivatives of basis function i at given point in cell
863
virtual void evaluate_basis_derivatives(unsigned int i,
866
const double* coordinates,
867
const ufc::cell& c) const
869
// Extract vertex coordinates
870
const double * const * element_coordinates = c.coordinates;
872
// Compute Jacobian of affine map from reference cell
873
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
874
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
875
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
876
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
877
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
878
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
879
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
880
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
881
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
883
// Compute sub determinants
884
const double d00 = J_11*J_22 - J_12*J_21;
885
const double d01 = J_12*J_20 - J_10*J_22;
886
const double d02 = J_10*J_21 - J_11*J_20;
888
const double d10 = J_02*J_21 - J_01*J_22;
889
const double d11 = J_00*J_22 - J_02*J_20;
890
const double d12 = J_01*J_20 - J_00*J_21;
892
const double d20 = J_01*J_12 - J_02*J_11;
893
const double d21 = J_02*J_10 - J_00*J_12;
894
const double d22 = J_00*J_11 - J_01*J_10;
896
// Compute determinant of Jacobian
897
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
899
// Compute inverse of Jacobian
902
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
903
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
904
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
906
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
907
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
908
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
910
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
911
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
912
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
914
// Get coordinates and map to the UFC reference element
915
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
916
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
917
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
919
// Map coordinates to the reference cube
920
if (std::abs(y + z - 1.0) < 1e-14)
923
x = -2.0 * x/(y + z - 1.0) - 1.0;
924
if (std::abs(z - 1.0) < 1e-14)
927
y = 2.0 * y/(1.0 - z) - 1.0;
930
// Compute number of derivatives
931
unsigned int num_derivatives = 1;
933
for (unsigned int j = 0; j < n; j++)
934
num_derivatives *= 3;
937
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
938
unsigned int **combinations = new unsigned int *[num_derivatives];
940
for (unsigned int j = 0; j < num_derivatives; j++)
942
combinations[j] = new unsigned int [n];
943
for (unsigned int k = 0; k < n; k++)
944
combinations[j][k] = 0;
947
// Generate combinations of derivatives
948
for (unsigned int row = 1; row < num_derivatives; row++)
950
for (unsigned int num = 0; num < row; num++)
952
for (unsigned int col = n-1; col+1 > 0; col--)
954
if (combinations[row][col] + 1 > 2)
955
combinations[row][col] = 0;
958
combinations[row][col] += 1;
965
// Compute inverse of Jacobian
966
const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
968
// Declare transformation matrix
969
// Declare pointer to two dimensional array and initialise
970
double **transform = new double *[num_derivatives];
972
for (unsigned int j = 0; j < num_derivatives; j++)
974
transform[j] = new double [num_derivatives];
975
for (unsigned int k = 0; k < num_derivatives; k++)
979
// Construct transformation matrix
980
for (unsigned int row = 0; row < num_derivatives; row++)
982
for (unsigned int col = 0; col < num_derivatives; col++)
984
for (unsigned int k = 0; k < n; k++)
985
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
990
for (unsigned int j = 0; j < 1*num_derivatives; j++)
993
// Map degree of freedom to element degree of freedom
994
const unsigned int dof = i;
997
const double scalings_y_0 = 1;
998
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
999
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
1000
const double scalings_z_0 = 1;
1001
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
1002
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
1004
// Compute psitilde_a
1005
const double psitilde_a_0 = 1;
1006
const double psitilde_a_1 = x;
1007
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
1009
// Compute psitilde_bs
1010
const double psitilde_bs_0_0 = 1;
1011
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1012
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;
1013
const double psitilde_bs_1_0 = 1;
1014
const double psitilde_bs_1_1 = 2.5*y + 1.5;
1015
const double psitilde_bs_2_0 = 1;
1017
// Compute psitilde_cs
1018
const double psitilde_cs_00_0 = 1;
1019
const double psitilde_cs_00_1 = 2*z + 1;
1020
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;
1021
const double psitilde_cs_01_0 = 1;
1022
const double psitilde_cs_01_1 = 3*z + 2;
1023
const double psitilde_cs_02_0 = 1;
1024
const double psitilde_cs_10_0 = 1;
1025
const double psitilde_cs_10_1 = 3*z + 2;
1026
const double psitilde_cs_11_0 = 1;
1027
const double psitilde_cs_20_0 = 1;
1029
// Compute basisvalues
1030
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
1031
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
1032
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
1033
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
1034
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
1035
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
1036
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
1037
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
1038
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
1039
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
1041
// Table(s) of coefficients
1042
const static double coefficients0[10][10] = \
1043
{{-0.0577350269189625, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
1044
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499976, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
1045
{-0.0577350269189626, 0, 0.0702728368926306, -0.0248451997499976, 0, 0, 0.087287156094397, 0, -0.0475131096733199, 0.0167984210226323},
1046
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
1047
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
1048
{0.23094010767585, 0.121716123890037, -0.0702728368926306, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
1049
{0.23094010767585, 0.121716123890037, 0.0702728368926307, -0.0993807989999907, 0, 0.100790526135794, -0.087287156094397, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
1050
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999906, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
1051
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
1052
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
1054
// Interesting (new) part
1055
// Tables of derivatives of the polynomial base (transpose)
1056
const static double dmats0[10][10] = \
1057
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1058
{6.32455532033676, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1059
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1060
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1061
{0, 11.2249721603218, 0, 0, 0, 0, 0, 0, 0, 0},
1062
{4.58257569495584, 0, 8.36660026534076, -1.18321595661992, 0, 0, 0, 0, 0, 0},
1063
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1064
{3.74165738677394, 0, 0, 8.69482604771366, 0, 0, 0, 0, 0, 0},
1065
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1066
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1068
const static double dmats1[10][10] = \
1069
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1070
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1071
{5.47722557505166, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1072
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1073
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825974, 0, 0, 0, 0, 0, 0},
1074
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
1075
{-2.64575131106459, 0, 9.66091783079296, 0.683130051063974, 0, 0, 0, 0, 0, 0},
1076
{1.87082869338697, 0, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
1077
{3.24037034920393, 0, 0, 7.52994023880668, 0, 0, 0, 0, 0, 0},
1078
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1080
const static double dmats2[10][10] = \
1081
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1082
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1083
{1.82574185835055, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1084
{5.16397779494322, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1085
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
1086
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
1087
{1.32287565553229, 0, 3.86436713231718, -0.341565025531986, 0, 0, 0, 0, 0, 0},
1088
{1.87082869338697, 7.09929573971954, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
1089
{1.08012344973464, 0, 7.09929573971954, 2.50998007960223, 0, 0, 0, 0, 0, 0},
1090
{-3.81881307912986, 0, 0, 8.87411967464942, 0, 0, 0, 0, 0, 0}};
1092
// Compute reference derivatives
1093
// Declare pointer to array of derivatives on FIAT element
1094
double *derivatives = new double [num_derivatives];
1096
// Declare coefficients
1097
double coeff0_0 = 0;
1098
double coeff0_1 = 0;
1099
double coeff0_2 = 0;
1100
double coeff0_3 = 0;
1101
double coeff0_4 = 0;
1102
double coeff0_5 = 0;
1103
double coeff0_6 = 0;
1104
double coeff0_7 = 0;
1105
double coeff0_8 = 0;
1106
double coeff0_9 = 0;
1108
// Declare new coefficients
1109
double new_coeff0_0 = 0;
1110
double new_coeff0_1 = 0;
1111
double new_coeff0_2 = 0;
1112
double new_coeff0_3 = 0;
1113
double new_coeff0_4 = 0;
1114
double new_coeff0_5 = 0;
1115
double new_coeff0_6 = 0;
1116
double new_coeff0_7 = 0;
1117
double new_coeff0_8 = 0;
1118
double new_coeff0_9 = 0;
1120
// Loop possible derivatives
1121
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1123
// Get values from coefficients array
1124
new_coeff0_0 = coefficients0[dof][0];
1125
new_coeff0_1 = coefficients0[dof][1];
1126
new_coeff0_2 = coefficients0[dof][2];
1127
new_coeff0_3 = coefficients0[dof][3];
1128
new_coeff0_4 = coefficients0[dof][4];
1129
new_coeff0_5 = coefficients0[dof][5];
1130
new_coeff0_6 = coefficients0[dof][6];
1131
new_coeff0_7 = coefficients0[dof][7];
1132
new_coeff0_8 = coefficients0[dof][8];
1133
new_coeff0_9 = coefficients0[dof][9];
1135
// Loop derivative order
1136
for (unsigned int j = 0; j < n; j++)
1138
// Update old coefficients
1139
coeff0_0 = new_coeff0_0;
1140
coeff0_1 = new_coeff0_1;
1141
coeff0_2 = new_coeff0_2;
1142
coeff0_3 = new_coeff0_3;
1143
coeff0_4 = new_coeff0_4;
1144
coeff0_5 = new_coeff0_5;
1145
coeff0_6 = new_coeff0_6;
1146
coeff0_7 = new_coeff0_7;
1147
coeff0_8 = new_coeff0_8;
1148
coeff0_9 = new_coeff0_9;
1150
if(combinations[deriv_num][j] == 0)
1152
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];
1153
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];
1154
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];
1155
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];
1156
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];
1157
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];
1158
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];
1159
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];
1160
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];
1161
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];
1163
if(combinations[deriv_num][j] == 1)
1165
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];
1166
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];
1167
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];
1168
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];
1169
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];
1170
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];
1171
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];
1172
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];
1173
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];
1174
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];
1176
if(combinations[deriv_num][j] == 2)
1178
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];
1179
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];
1180
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];
1181
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];
1182
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];
1183
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];
1184
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];
1185
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];
1186
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];
1187
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];
1191
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1192
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;
1195
// Transform derivatives back to physical element
1196
for (unsigned int row = 0; row < num_derivatives; row++)
1198
for (unsigned int col = 0; col < num_derivatives; col++)
1200
values[row] += transform[row][col]*derivatives[col];
1203
// Delete pointer to array of derivatives on FIAT element
1204
delete [] derivatives;
1206
// Delete pointer to array of combinations of derivatives and transform
1207
for (unsigned int row = 0; row < num_derivatives; row++)
1209
delete [] combinations[row];
1210
delete [] transform[row];
1213
delete [] combinations;
1214
delete [] transform;
1217
/// Evaluate order n derivatives of all basis functions at given point in cell
1218
virtual void evaluate_basis_derivatives_all(unsigned int n,
1220
const double* coordinates,
1221
const ufc::cell& c) const
1223
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1226
/// Evaluate linear functional for dof i on the function f
1227
virtual double evaluate_dof(unsigned int i,
1228
const ufc::function& f,
1229
const ufc::cell& c) const
1231
// The reference points, direction and weights:
1232
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}}};
1233
const static double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
1234
const static double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
1236
const double * const * x = c.coordinates;
1237
double result = 0.0;
1238
// Iterate over the points:
1239
// Evaluate basis functions for affine mapping
1240
const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
1241
const double w1 = X[i][0][0];
1242
const double w2 = X[i][0][1];
1243
const double w3 = X[i][0][2];
1245
// Compute affine mapping y = F(X)
1247
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
1248
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
1249
y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
1251
// Evaluate function at physical points
1253
f.evaluate(values, y, c);
1255
// Map function values using appropriate mapping
1256
// Affine map: Do nothing
1258
// Note that we do not map the weights (yet).
1260
// Take directional components
1261
for(int k = 0; k < 1; k++)
1262
result += values[k]*D[i][0][k];
1263
// Multiply by weights
1269
/// Evaluate linear functionals for all dofs on the function f
1270
virtual void evaluate_dofs(double* values,
1271
const ufc::function& f,
1272
const ufc::cell& c) const
1274
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1277
/// Interpolate vertex values from dof values
1278
virtual void interpolate_vertex_values(double* vertex_values,
1279
const double* dof_values,
1280
const ufc::cell& c) const
1282
// Evaluate at vertices and use affine mapping
1283
vertex_values[0] = dof_values[0];
1284
vertex_values[1] = dof_values[1];
1285
vertex_values[2] = dof_values[2];
1286
vertex_values[3] = dof_values[3];
1289
/// Return the number of sub elements (for a mixed element)
1290
virtual unsigned int num_sub_elements() const
1295
/// Create a new finite element for sub element i (for a mixed element)
1296
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1298
return new UFC_ProjectionBilinearForm_finite_element_1();
1303
/// This class defines the interface for a local-to-global mapping of
1304
/// degrees of freedom (dofs).
1306
class UFC_ProjectionBilinearForm_dof_map_0: public ufc::dof_map
1310
unsigned int __global_dimension;
1315
UFC_ProjectionBilinearForm_dof_map_0() : ufc::dof_map()
1317
__global_dimension = 0;
1321
virtual ~UFC_ProjectionBilinearForm_dof_map_0()
1326
/// Return a string identifying the dof map
1327
virtual const char* signature() const
1329
return "FFC dof map for FiniteElement('Lagrange', 'tetrahedron', 2)";
1332
/// Return true iff mesh entities of topological dimension d are needed
1333
virtual bool needs_mesh_entities(unsigned int d) const
1353
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1354
virtual bool init_mesh(const ufc::mesh& m)
1356
__global_dimension = m.num_entities[0] + m.num_entities[1];
1360
/// Initialize dof map for given cell
1361
virtual void init_cell(const ufc::mesh& m,
1367
/// Finish initialization of dof map for cells
1368
virtual void init_cell_finalize()
1373
/// Return the dimension of the global finite element function space
1374
virtual unsigned int global_dimension() const
1376
return __global_dimension;
1379
/// Return the dimension of the local finite element function space
1380
virtual unsigned int local_dimension() const
1385
// Return the geometric dimension of the coordinates this dof map provides
1386
virtual unsigned int geometric_dimension() const
1391
/// Return the number of dofs on each cell facet
1392
virtual unsigned int num_facet_dofs() const
1397
/// Return the number of dofs associated with each cell entity of dimension d
1398
virtual unsigned int num_entity_dofs(unsigned int d) const
1400
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1403
/// Tabulate the local-to-global mapping of dofs on a cell
1404
virtual void tabulate_dofs(unsigned int* dofs,
1406
const ufc::cell& c) const
1408
dofs[0] = c.entity_indices[0][0];
1409
dofs[1] = c.entity_indices[0][1];
1410
dofs[2] = c.entity_indices[0][2];
1411
dofs[3] = c.entity_indices[0][3];
1412
unsigned int offset = m.num_entities[0];
1413
dofs[4] = offset + c.entity_indices[1][0];
1414
dofs[5] = offset + c.entity_indices[1][1];
1415
dofs[6] = offset + c.entity_indices[1][2];
1416
dofs[7] = offset + c.entity_indices[1][3];
1417
dofs[8] = offset + c.entity_indices[1][4];
1418
dofs[9] = offset + c.entity_indices[1][5];
1421
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1422
virtual void tabulate_facet_dofs(unsigned int* dofs,
1423
unsigned int facet) const
1462
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1463
virtual void tabulate_entity_dofs(unsigned int* dofs,
1464
unsigned int d, unsigned int i) const
1466
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1469
/// Tabulate the coordinates of all dofs on a cell
1470
virtual void tabulate_coordinates(double** coordinates,
1471
const ufc::cell& c) const
1473
const double * const * x = c.coordinates;
1474
coordinates[0][0] = x[0][0];
1475
coordinates[0][1] = x[0][1];
1476
coordinates[0][2] = x[0][2];
1477
coordinates[1][0] = x[1][0];
1478
coordinates[1][1] = x[1][1];
1479
coordinates[1][2] = x[1][2];
1480
coordinates[2][0] = x[2][0];
1481
coordinates[2][1] = x[2][1];
1482
coordinates[2][2] = x[2][2];
1483
coordinates[3][0] = x[3][0];
1484
coordinates[3][1] = x[3][1];
1485
coordinates[3][2] = x[3][2];
1486
coordinates[4][0] = 0.5*x[2][0] + 0.5*x[3][0];
1487
coordinates[4][1] = 0.5*x[2][1] + 0.5*x[3][1];
1488
coordinates[4][2] = 0.5*x[2][2] + 0.5*x[3][2];
1489
coordinates[5][0] = 0.5*x[1][0] + 0.5*x[3][0];
1490
coordinates[5][1] = 0.5*x[1][1] + 0.5*x[3][1];
1491
coordinates[5][2] = 0.5*x[1][2] + 0.5*x[3][2];
1492
coordinates[6][0] = 0.5*x[1][0] + 0.5*x[2][0];
1493
coordinates[6][1] = 0.5*x[1][1] + 0.5*x[2][1];
1494
coordinates[6][2] = 0.5*x[1][2] + 0.5*x[2][2];
1495
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[3][0];
1496
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[3][1];
1497
coordinates[7][2] = 0.5*x[0][2] + 0.5*x[3][2];
1498
coordinates[8][0] = 0.5*x[0][0] + 0.5*x[2][0];
1499
coordinates[8][1] = 0.5*x[0][1] + 0.5*x[2][1];
1500
coordinates[8][2] = 0.5*x[0][2] + 0.5*x[2][2];
1501
coordinates[9][0] = 0.5*x[0][0] + 0.5*x[1][0];
1502
coordinates[9][1] = 0.5*x[0][1] + 0.5*x[1][1];
1503
coordinates[9][2] = 0.5*x[0][2] + 0.5*x[1][2];
1506
/// Return the number of sub dof maps (for a mixed element)
1507
virtual unsigned int num_sub_dof_maps() const
1512
/// Create a new dof_map for sub dof map i (for a mixed element)
1513
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1515
return new UFC_ProjectionBilinearForm_dof_map_0();
1520
/// This class defines the interface for a local-to-global mapping of
1521
/// degrees of freedom (dofs).
1523
class UFC_ProjectionBilinearForm_dof_map_1: public ufc::dof_map
1527
unsigned int __global_dimension;
1532
UFC_ProjectionBilinearForm_dof_map_1() : ufc::dof_map()
1534
__global_dimension = 0;
1538
virtual ~UFC_ProjectionBilinearForm_dof_map_1()
1543
/// Return a string identifying the dof map
1544
virtual const char* signature() const
1546
return "FFC dof map for FiniteElement('Lagrange', 'tetrahedron', 2)";
1549
/// Return true iff mesh entities of topological dimension d are needed
1550
virtual bool needs_mesh_entities(unsigned int d) const
1570
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1571
virtual bool init_mesh(const ufc::mesh& m)
1573
__global_dimension = m.num_entities[0] + m.num_entities[1];
1577
/// Initialize dof map for given cell
1578
virtual void init_cell(const ufc::mesh& m,
1584
/// Finish initialization of dof map for cells
1585
virtual void init_cell_finalize()
1590
/// Return the dimension of the global finite element function space
1591
virtual unsigned int global_dimension() const
1593
return __global_dimension;
1596
/// Return the dimension of the local finite element function space
1597
virtual unsigned int local_dimension() const
1602
// Return the geometric dimension of the coordinates this dof map provides
1603
virtual unsigned int geometric_dimension() const
1608
/// Return the number of dofs on each cell facet
1609
virtual unsigned int num_facet_dofs() const
1614
/// Return the number of dofs associated with each cell entity of dimension d
1615
virtual unsigned int num_entity_dofs(unsigned int d) const
1617
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1620
/// Tabulate the local-to-global mapping of dofs on a cell
1621
virtual void tabulate_dofs(unsigned int* dofs,
1623
const ufc::cell& c) const
1625
dofs[0] = c.entity_indices[0][0];
1626
dofs[1] = c.entity_indices[0][1];
1627
dofs[2] = c.entity_indices[0][2];
1628
dofs[3] = c.entity_indices[0][3];
1629
unsigned int offset = m.num_entities[0];
1630
dofs[4] = offset + c.entity_indices[1][0];
1631
dofs[5] = offset + c.entity_indices[1][1];
1632
dofs[6] = offset + c.entity_indices[1][2];
1633
dofs[7] = offset + c.entity_indices[1][3];
1634
dofs[8] = offset + c.entity_indices[1][4];
1635
dofs[9] = offset + c.entity_indices[1][5];
1638
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1639
virtual void tabulate_facet_dofs(unsigned int* dofs,
1640
unsigned int facet) const
1679
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1680
virtual void tabulate_entity_dofs(unsigned int* dofs,
1681
unsigned int d, unsigned int i) const
1683
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1686
/// Tabulate the coordinates of all dofs on a cell
1687
virtual void tabulate_coordinates(double** coordinates,
1688
const ufc::cell& c) const
1690
const double * const * x = c.coordinates;
1691
coordinates[0][0] = x[0][0];
1692
coordinates[0][1] = x[0][1];
1693
coordinates[0][2] = x[0][2];
1694
coordinates[1][0] = x[1][0];
1695
coordinates[1][1] = x[1][1];
1696
coordinates[1][2] = x[1][2];
1697
coordinates[2][0] = x[2][0];
1698
coordinates[2][1] = x[2][1];
1699
coordinates[2][2] = x[2][2];
1700
coordinates[3][0] = x[3][0];
1701
coordinates[3][1] = x[3][1];
1702
coordinates[3][2] = x[3][2];
1703
coordinates[4][0] = 0.5*x[2][0] + 0.5*x[3][0];
1704
coordinates[4][1] = 0.5*x[2][1] + 0.5*x[3][1];
1705
coordinates[4][2] = 0.5*x[2][2] + 0.5*x[3][2];
1706
coordinates[5][0] = 0.5*x[1][0] + 0.5*x[3][0];
1707
coordinates[5][1] = 0.5*x[1][1] + 0.5*x[3][1];
1708
coordinates[5][2] = 0.5*x[1][2] + 0.5*x[3][2];
1709
coordinates[6][0] = 0.5*x[1][0] + 0.5*x[2][0];
1710
coordinates[6][1] = 0.5*x[1][1] + 0.5*x[2][1];
1711
coordinates[6][2] = 0.5*x[1][2] + 0.5*x[2][2];
1712
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[3][0];
1713
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[3][1];
1714
coordinates[7][2] = 0.5*x[0][2] + 0.5*x[3][2];
1715
coordinates[8][0] = 0.5*x[0][0] + 0.5*x[2][0];
1716
coordinates[8][1] = 0.5*x[0][1] + 0.5*x[2][1];
1717
coordinates[8][2] = 0.5*x[0][2] + 0.5*x[2][2];
1718
coordinates[9][0] = 0.5*x[0][0] + 0.5*x[1][0];
1719
coordinates[9][1] = 0.5*x[0][1] + 0.5*x[1][1];
1720
coordinates[9][2] = 0.5*x[0][2] + 0.5*x[1][2];
1723
/// Return the number of sub dof maps (for a mixed element)
1724
virtual unsigned int num_sub_dof_maps() const
1729
/// Create a new dof_map for sub dof map i (for a mixed element)
1730
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1732
return new UFC_ProjectionBilinearForm_dof_map_1();
15
/// This class defines the interface for a finite element.
17
class projection_0_finite_element_0: public ufc::finite_element
22
projection_0_finite_element_0() : ufc::finite_element()
28
virtual ~projection_0_finite_element_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "FiniteElement('Lagrange', 'tetrahedron', 2)";
39
/// Return the cell shape
40
virtual ufc::shape cell_shape() const
42
return ufc::tetrahedron;
45
/// Return the dimension of the finite element function space
46
virtual unsigned int space_dimension() const
51
/// Return the rank of the value space
52
virtual unsigned int value_rank() const
57
/// Return the dimension of the value space for axis i
58
virtual unsigned int value_dimension(unsigned int i) const
63
/// Evaluate basis function i at given point in cell
64
virtual void evaluate_basis(unsigned int i,
66
const double* coordinates,
67
const ufc::cell& c) const
69
// Extract vertex coordinates
70
const double * const * element_coordinates = c.coordinates;
72
// Compute Jacobian of affine map from reference cell
73
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
74
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
75
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
76
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
77
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
78
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
79
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
80
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
81
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
83
// Compute sub determinants
84
const double d00 = J_11*J_22 - J_12*J_21;
85
const double d01 = J_12*J_20 - J_10*J_22;
86
const double d02 = J_10*J_21 - J_11*J_20;
88
const double d10 = J_02*J_21 - J_01*J_22;
89
const double d11 = J_00*J_22 - J_02*J_20;
90
const double d12 = J_01*J_20 - J_00*J_21;
92
const double d20 = J_01*J_12 - J_02*J_11;
93
const double d21 = J_02*J_10 - J_00*J_12;
94
const double d22 = J_00*J_11 - J_01*J_10;
96
// Compute determinant of Jacobian
97
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
99
// Compute inverse of Jacobian
102
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
103
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
104
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
106
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
107
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
108
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
110
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
111
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
112
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
114
// Get coordinates and map to the UFC reference element
115
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
116
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
117
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
119
// Map coordinates to the reference cube
120
if (std::abs(y + z - 1.0) < 1e-14)
123
x = -2.0 * x/(y + z - 1.0) - 1.0;
124
if (std::abs(z - 1.0) < 1e-14)
127
y = 2.0 * y/(1.0 - z) - 1.0;
133
// Map degree of freedom to element degree of freedom
134
const unsigned int dof = i;
137
const double scalings_y_0 = 1;
138
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
139
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
140
const double scalings_z_0 = 1;
141
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
142
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
144
// Compute psitilde_a
145
const double psitilde_a_0 = 1;
146
const double psitilde_a_1 = x;
147
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
149
// Compute psitilde_bs
150
const double psitilde_bs_0_0 = 1;
151
const double psitilde_bs_0_1 = 1.5*y + 0.5;
152
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;
153
const double psitilde_bs_1_0 = 1;
154
const double psitilde_bs_1_1 = 2.5*y + 1.5;
155
const double psitilde_bs_2_0 = 1;
157
// Compute psitilde_cs
158
const double psitilde_cs_00_0 = 1;
159
const double psitilde_cs_00_1 = 2*z + 1;
160
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;
161
const double psitilde_cs_01_0 = 1;
162
const double psitilde_cs_01_1 = 3*z + 2;
163
const double psitilde_cs_02_0 = 1;
164
const double psitilde_cs_10_0 = 1;
165
const double psitilde_cs_10_1 = 3*z + 2;
166
const double psitilde_cs_11_0 = 1;
167
const double psitilde_cs_20_0 = 1;
169
// Compute basisvalues
170
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
171
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
172
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
173
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
174
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
175
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
176
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
177
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
178
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
179
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
181
// Table(s) of coefficients
182
static const double coefficients0[10][10] = \
183
{{-0.0577350269189626, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.0503952630678969, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
184
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
185
{-0.0577350269189626, 0, 0.0702728368926307, -0.0248451997499977, 0, 0, 0.0872871560943969, 0, -0.0475131096733199, 0.0167984210226323},
186
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
187
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
188
{0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
189
{0.23094010767585, 0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0, 0.100790526135794, -0.0872871560943969, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
190
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
191
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
192
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
194
// Extract relevant coefficients
195
const double coeff0_0 = coefficients0[dof][0];
196
const double coeff0_1 = coefficients0[dof][1];
197
const double coeff0_2 = coefficients0[dof][2];
198
const double coeff0_3 = coefficients0[dof][3];
199
const double coeff0_4 = coefficients0[dof][4];
200
const double coeff0_5 = coefficients0[dof][5];
201
const double coeff0_6 = coefficients0[dof][6];
202
const double coeff0_7 = coefficients0[dof][7];
203
const double coeff0_8 = coefficients0[dof][8];
204
const double coeff0_9 = coefficients0[dof][9];
207
*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;
210
/// Evaluate all basis functions at given point in cell
211
virtual void evaluate_basis_all(double* values,
212
const double* coordinates,
213
const ufc::cell& c) const
215
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
218
/// Evaluate order n derivatives of basis function i at given point in cell
219
virtual void evaluate_basis_derivatives(unsigned int i,
222
const double* coordinates,
223
const ufc::cell& c) const
225
// Extract vertex coordinates
226
const double * const * element_coordinates = c.coordinates;
228
// Compute Jacobian of affine map from reference cell
229
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
230
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
231
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
232
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
233
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
234
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
235
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
236
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
237
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
239
// Compute sub determinants
240
const double d00 = J_11*J_22 - J_12*J_21;
241
const double d01 = J_12*J_20 - J_10*J_22;
242
const double d02 = J_10*J_21 - J_11*J_20;
244
const double d10 = J_02*J_21 - J_01*J_22;
245
const double d11 = J_00*J_22 - J_02*J_20;
246
const double d12 = J_01*J_20 - J_00*J_21;
248
const double d20 = J_01*J_12 - J_02*J_11;
249
const double d21 = J_02*J_10 - J_00*J_12;
250
const double d22 = J_00*J_11 - J_01*J_10;
252
// Compute determinant of Jacobian
253
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
255
// Compute inverse of Jacobian
258
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
259
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
260
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
262
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
263
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
264
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
266
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
267
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
268
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
270
// Get coordinates and map to the UFC reference element
271
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
272
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
273
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
275
// Map coordinates to the reference cube
276
if (std::abs(y + z - 1.0) < 1e-14)
279
x = -2.0 * x/(y + z - 1.0) - 1.0;
280
if (std::abs(z - 1.0) < 1e-14)
283
y = 2.0 * y/(1.0 - z) - 1.0;
286
// Compute number of derivatives
287
unsigned int num_derivatives = 1;
289
for (unsigned int j = 0; j < n; j++)
290
num_derivatives *= 3;
293
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
294
unsigned int **combinations = new unsigned int *[num_derivatives];
296
for (unsigned int j = 0; j < num_derivatives; j++)
298
combinations[j] = new unsigned int [n];
299
for (unsigned int k = 0; k < n; k++)
300
combinations[j][k] = 0;
303
// Generate combinations of derivatives
304
for (unsigned int row = 1; row < num_derivatives; row++)
306
for (unsigned int num = 0; num < row; num++)
308
for (unsigned int col = n-1; col+1 > 0; col--)
310
if (combinations[row][col] + 1 > 2)
311
combinations[row][col] = 0;
314
combinations[row][col] += 1;
321
// Compute inverse of Jacobian
322
const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
324
// Declare transformation matrix
325
// Declare pointer to two dimensional array and initialise
326
double **transform = new double *[num_derivatives];
328
for (unsigned int j = 0; j < num_derivatives; j++)
330
transform[j] = new double [num_derivatives];
331
for (unsigned int k = 0; k < num_derivatives; k++)
335
// Construct transformation matrix
336
for (unsigned int row = 0; row < num_derivatives; row++)
338
for (unsigned int col = 0; col < num_derivatives; col++)
340
for (unsigned int k = 0; k < n; k++)
341
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
346
for (unsigned int j = 0; j < 1*num_derivatives; j++)
349
// Map degree of freedom to element degree of freedom
350
const unsigned int dof = i;
353
const double scalings_y_0 = 1;
354
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
355
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
356
const double scalings_z_0 = 1;
357
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
358
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
360
// Compute psitilde_a
361
const double psitilde_a_0 = 1;
362
const double psitilde_a_1 = x;
363
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
365
// Compute psitilde_bs
366
const double psitilde_bs_0_0 = 1;
367
const double psitilde_bs_0_1 = 1.5*y + 0.5;
368
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;
369
const double psitilde_bs_1_0 = 1;
370
const double psitilde_bs_1_1 = 2.5*y + 1.5;
371
const double psitilde_bs_2_0 = 1;
373
// Compute psitilde_cs
374
const double psitilde_cs_00_0 = 1;
375
const double psitilde_cs_00_1 = 2*z + 1;
376
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;
377
const double psitilde_cs_01_0 = 1;
378
const double psitilde_cs_01_1 = 3*z + 2;
379
const double psitilde_cs_02_0 = 1;
380
const double psitilde_cs_10_0 = 1;
381
const double psitilde_cs_10_1 = 3*z + 2;
382
const double psitilde_cs_11_0 = 1;
383
const double psitilde_cs_20_0 = 1;
385
// Compute basisvalues
386
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
387
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
388
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
389
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
390
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
391
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
392
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
393
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
394
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
395
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
397
// Table(s) of coefficients
398
static const double coefficients0[10][10] = \
399
{{-0.0577350269189626, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.0503952630678969, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
400
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
401
{-0.0577350269189626, 0, 0.0702728368926307, -0.0248451997499977, 0, 0, 0.0872871560943969, 0, -0.0475131096733199, 0.0167984210226323},
402
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
403
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
404
{0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
405
{0.23094010767585, 0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0, 0.100790526135794, -0.0872871560943969, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
406
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
407
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
408
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
410
// Interesting (new) part
411
// Tables of derivatives of the polynomial base (transpose)
412
static const double dmats0[10][10] = \
413
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
414
{6.32455532033676, 0, 0, 0, 0, 0, 0, 0, 0, 0},
415
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
416
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
417
{0, 11.2249721603218, 0, 0, 0, 0, 0, 0, 0, 0},
418
{4.58257569495584, 0, 8.36660026534076, -1.18321595661992, 0, 0, 0, 0, 0, 0},
419
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
420
{3.74165738677394, 0, 0, 8.69482604771366, 0, 0, 0, 0, 0, 0},
421
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
422
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
424
static const double dmats1[10][10] = \
425
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
426
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
427
{5.47722557505166, 0, 0, 0, 0, 0, 0, 0, 0, 0},
428
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
429
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
430
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
431
{-2.64575131106459, 0, 9.66091783079296, 0.683130051063973, 0, 0, 0, 0, 0, 0},
432
{1.87082869338697, 0, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
433
{3.24037034920393, 0, 0, 7.52994023880668, 0, 0, 0, 0, 0, 0},
434
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
436
static const double dmats2[10][10] = \
437
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
438
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
439
{1.82574185835055, 0, 0, 0, 0, 0, 0, 0, 0, 0},
440
{5.16397779494322, 0, 0, 0, 0, 0, 0, 0, 0, 0},
441
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
442
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
443
{1.3228756555323, 0, 3.86436713231718, -0.341565025531986, 0, 0, 0, 0, 0, 0},
444
{1.87082869338697, 7.09929573971954, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
445
{1.08012344973464, 0, 7.09929573971954, 2.50998007960223, 0, 0, 0, 0, 0, 0},
446
{-3.81881307912987, 0, 0, 8.87411967464942, 0, 0, 0, 0, 0, 0}};
448
// Compute reference derivatives
449
// Declare pointer to array of derivatives on FIAT element
450
double *derivatives = new double [num_derivatives];
452
// Declare coefficients
464
// Declare new coefficients
465
double new_coeff0_0 = 0;
466
double new_coeff0_1 = 0;
467
double new_coeff0_2 = 0;
468
double new_coeff0_3 = 0;
469
double new_coeff0_4 = 0;
470
double new_coeff0_5 = 0;
471
double new_coeff0_6 = 0;
472
double new_coeff0_7 = 0;
473
double new_coeff0_8 = 0;
474
double new_coeff0_9 = 0;
476
// Loop possible derivatives
477
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
479
// Get values from coefficients array
480
new_coeff0_0 = coefficients0[dof][0];
481
new_coeff0_1 = coefficients0[dof][1];
482
new_coeff0_2 = coefficients0[dof][2];
483
new_coeff0_3 = coefficients0[dof][3];
484
new_coeff0_4 = coefficients0[dof][4];
485
new_coeff0_5 = coefficients0[dof][5];
486
new_coeff0_6 = coefficients0[dof][6];
487
new_coeff0_7 = coefficients0[dof][7];
488
new_coeff0_8 = coefficients0[dof][8];
489
new_coeff0_9 = coefficients0[dof][9];
491
// Loop derivative order
492
for (unsigned int j = 0; j < n; j++)
494
// Update old coefficients
495
coeff0_0 = new_coeff0_0;
496
coeff0_1 = new_coeff0_1;
497
coeff0_2 = new_coeff0_2;
498
coeff0_3 = new_coeff0_3;
499
coeff0_4 = new_coeff0_4;
500
coeff0_5 = new_coeff0_5;
501
coeff0_6 = new_coeff0_6;
502
coeff0_7 = new_coeff0_7;
503
coeff0_8 = new_coeff0_8;
504
coeff0_9 = new_coeff0_9;
506
if(combinations[deriv_num][j] == 0)
508
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];
509
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];
510
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];
511
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];
512
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];
513
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];
514
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];
515
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];
516
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];
517
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];
519
if(combinations[deriv_num][j] == 1)
521
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];
522
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];
523
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];
524
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];
525
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];
526
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];
527
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];
528
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];
529
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];
530
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];
532
if(combinations[deriv_num][j] == 2)
534
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];
535
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];
536
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];
537
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];
538
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];
539
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];
540
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];
541
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];
542
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];
543
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];
547
// Compute derivatives on reference element as dot product of coefficients and basisvalues
548
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;
551
// Transform derivatives back to physical element
552
for (unsigned int row = 0; row < num_derivatives; row++)
554
for (unsigned int col = 0; col < num_derivatives; col++)
556
values[row] += transform[row][col]*derivatives[col];
559
// Delete pointer to array of derivatives on FIAT element
560
delete [] derivatives;
562
// Delete pointer to array of combinations of derivatives and transform
563
for (unsigned int row = 0; row < num_derivatives; row++)
565
delete [] combinations[row];
566
delete [] transform[row];
569
delete [] combinations;
573
/// Evaluate order n derivatives of all basis functions at given point in cell
574
virtual void evaluate_basis_derivatives_all(unsigned int n,
576
const double* coordinates,
577
const ufc::cell& c) const
579
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
582
/// Evaluate linear functional for dof i on the function f
583
virtual double evaluate_dof(unsigned int i,
584
const ufc::function& f,
585
const ufc::cell& c) const
587
// The reference points, direction and weights:
588
static const 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}}};
589
static const double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
590
static const double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
592
const double * const * x = c.coordinates;
594
// Iterate over the points:
595
// Evaluate basis functions for affine mapping
596
const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
597
const double w1 = X[i][0][0];
598
const double w2 = X[i][0][1];
599
const double w3 = X[i][0][2];
601
// Compute affine mapping y = F(X)
603
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
604
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
605
y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
607
// Evaluate function at physical points
609
f.evaluate(values, y, c);
611
// Map function values using appropriate mapping
612
// Affine map: Do nothing
614
// Note that we do not map the weights (yet).
616
// Take directional components
617
for(int k = 0; k < 1; k++)
618
result += values[k]*D[i][0][k];
619
// Multiply by weights
625
/// Evaluate linear functionals for all dofs on the function f
626
virtual void evaluate_dofs(double* values,
627
const ufc::function& f,
628
const ufc::cell& c) const
630
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
633
/// Interpolate vertex values from dof values
634
virtual void interpolate_vertex_values(double* vertex_values,
635
const double* dof_values,
636
const ufc::cell& c) const
638
// Evaluate at vertices and use affine mapping
639
vertex_values[0] = dof_values[0];
640
vertex_values[1] = dof_values[1];
641
vertex_values[2] = dof_values[2];
642
vertex_values[3] = dof_values[3];
645
/// Return the number of sub elements (for a mixed element)
646
virtual unsigned int num_sub_elements() const
651
/// Create a new finite element for sub element i (for a mixed element)
652
virtual ufc::finite_element* create_sub_element(unsigned int i) const
654
return new projection_0_finite_element_0();
659
/// This class defines the interface for a finite element.
661
class projection_0_finite_element_1: public ufc::finite_element
666
projection_0_finite_element_1() : ufc::finite_element()
672
virtual ~projection_0_finite_element_1()
677
/// Return a string identifying the finite element
678
virtual const char* signature() const
680
return "FiniteElement('Lagrange', 'tetrahedron', 2)";
683
/// Return the cell shape
684
virtual ufc::shape cell_shape() const
686
return ufc::tetrahedron;
689
/// Return the dimension of the finite element function space
690
virtual unsigned int space_dimension() const
695
/// Return the rank of the value space
696
virtual unsigned int value_rank() const
701
/// Return the dimension of the value space for axis i
702
virtual unsigned int value_dimension(unsigned int i) const
707
/// Evaluate basis function i at given point in cell
708
virtual void evaluate_basis(unsigned int i,
710
const double* coordinates,
711
const ufc::cell& c) const
713
// Extract vertex coordinates
714
const double * const * element_coordinates = c.coordinates;
716
// Compute Jacobian of affine map from reference cell
717
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
718
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
719
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
720
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
721
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
722
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
723
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
724
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
725
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
727
// Compute sub determinants
728
const double d00 = J_11*J_22 - J_12*J_21;
729
const double d01 = J_12*J_20 - J_10*J_22;
730
const double d02 = J_10*J_21 - J_11*J_20;
732
const double d10 = J_02*J_21 - J_01*J_22;
733
const double d11 = J_00*J_22 - J_02*J_20;
734
const double d12 = J_01*J_20 - J_00*J_21;
736
const double d20 = J_01*J_12 - J_02*J_11;
737
const double d21 = J_02*J_10 - J_00*J_12;
738
const double d22 = J_00*J_11 - J_01*J_10;
740
// Compute determinant of Jacobian
741
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
743
// Compute inverse of Jacobian
746
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
747
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
748
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
750
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
751
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
752
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
754
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
755
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
756
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
758
// Get coordinates and map to the UFC reference element
759
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
760
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
761
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
763
// Map coordinates to the reference cube
764
if (std::abs(y + z - 1.0) < 1e-14)
767
x = -2.0 * x/(y + z - 1.0) - 1.0;
768
if (std::abs(z - 1.0) < 1e-14)
771
y = 2.0 * y/(1.0 - z) - 1.0;
777
// Map degree of freedom to element degree of freedom
778
const unsigned int dof = i;
781
const double scalings_y_0 = 1;
782
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
783
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
784
const double scalings_z_0 = 1;
785
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
786
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
788
// Compute psitilde_a
789
const double psitilde_a_0 = 1;
790
const double psitilde_a_1 = x;
791
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
793
// Compute psitilde_bs
794
const double psitilde_bs_0_0 = 1;
795
const double psitilde_bs_0_1 = 1.5*y + 0.5;
796
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;
797
const double psitilde_bs_1_0 = 1;
798
const double psitilde_bs_1_1 = 2.5*y + 1.5;
799
const double psitilde_bs_2_0 = 1;
801
// Compute psitilde_cs
802
const double psitilde_cs_00_0 = 1;
803
const double psitilde_cs_00_1 = 2*z + 1;
804
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;
805
const double psitilde_cs_01_0 = 1;
806
const double psitilde_cs_01_1 = 3*z + 2;
807
const double psitilde_cs_02_0 = 1;
808
const double psitilde_cs_10_0 = 1;
809
const double psitilde_cs_10_1 = 3*z + 2;
810
const double psitilde_cs_11_0 = 1;
811
const double psitilde_cs_20_0 = 1;
813
// Compute basisvalues
814
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
815
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
816
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
817
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
818
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
819
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
820
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
821
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
822
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
823
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
825
// Table(s) of coefficients
826
static const double coefficients0[10][10] = \
827
{{-0.0577350269189626, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.0503952630678969, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
828
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
829
{-0.0577350269189626, 0, 0.0702728368926307, -0.0248451997499977, 0, 0, 0.0872871560943969, 0, -0.0475131096733199, 0.0167984210226323},
830
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
831
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
832
{0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
833
{0.23094010767585, 0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0, 0.100790526135794, -0.0872871560943969, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
834
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
835
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
836
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
838
// Extract relevant coefficients
839
const double coeff0_0 = coefficients0[dof][0];
840
const double coeff0_1 = coefficients0[dof][1];
841
const double coeff0_2 = coefficients0[dof][2];
842
const double coeff0_3 = coefficients0[dof][3];
843
const double coeff0_4 = coefficients0[dof][4];
844
const double coeff0_5 = coefficients0[dof][5];
845
const double coeff0_6 = coefficients0[dof][6];
846
const double coeff0_7 = coefficients0[dof][7];
847
const double coeff0_8 = coefficients0[dof][8];
848
const double coeff0_9 = coefficients0[dof][9];
851
*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;
854
/// Evaluate all basis functions at given point in cell
855
virtual void evaluate_basis_all(double* values,
856
const double* coordinates,
857
const ufc::cell& c) const
859
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
862
/// Evaluate order n derivatives of basis function i at given point in cell
863
virtual void evaluate_basis_derivatives(unsigned int i,
866
const double* coordinates,
867
const ufc::cell& c) const
869
// Extract vertex coordinates
870
const double * const * element_coordinates = c.coordinates;
872
// Compute Jacobian of affine map from reference cell
873
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
874
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
875
const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
876
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
877
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
878
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
879
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
880
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
881
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
883
// Compute sub determinants
884
const double d00 = J_11*J_22 - J_12*J_21;
885
const double d01 = J_12*J_20 - J_10*J_22;
886
const double d02 = J_10*J_21 - J_11*J_20;
888
const double d10 = J_02*J_21 - J_01*J_22;
889
const double d11 = J_00*J_22 - J_02*J_20;
890
const double d12 = J_01*J_20 - J_00*J_21;
892
const double d20 = J_01*J_12 - J_02*J_11;
893
const double d21 = J_02*J_10 - J_00*J_12;
894
const double d22 = J_00*J_11 - J_01*J_10;
896
// Compute determinant of Jacobian
897
double detJ = J_00*d00 + J_10*d10 + J_20*d20;
899
// Compute inverse of Jacobian
902
const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
903
+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
904
+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
906
const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
907
+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
908
+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
910
const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
911
+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
912
+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
914
// Get coordinates and map to the UFC reference element
915
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
916
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
917
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
919
// Map coordinates to the reference cube
920
if (std::abs(y + z - 1.0) < 1e-14)
923
x = -2.0 * x/(y + z - 1.0) - 1.0;
924
if (std::abs(z - 1.0) < 1e-14)
927
y = 2.0 * y/(1.0 - z) - 1.0;
930
// Compute number of derivatives
931
unsigned int num_derivatives = 1;
933
for (unsigned int j = 0; j < n; j++)
934
num_derivatives *= 3;
937
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
938
unsigned int **combinations = new unsigned int *[num_derivatives];
940
for (unsigned int j = 0; j < num_derivatives; j++)
942
combinations[j] = new unsigned int [n];
943
for (unsigned int k = 0; k < n; k++)
944
combinations[j][k] = 0;
947
// Generate combinations of derivatives
948
for (unsigned int row = 1; row < num_derivatives; row++)
950
for (unsigned int num = 0; num < row; num++)
952
for (unsigned int col = n-1; col+1 > 0; col--)
954
if (combinations[row][col] + 1 > 2)
955
combinations[row][col] = 0;
958
combinations[row][col] += 1;
965
// Compute inverse of Jacobian
966
const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
968
// Declare transformation matrix
969
// Declare pointer to two dimensional array and initialise
970
double **transform = new double *[num_derivatives];
972
for (unsigned int j = 0; j < num_derivatives; j++)
974
transform[j] = new double [num_derivatives];
975
for (unsigned int k = 0; k < num_derivatives; k++)
979
// Construct transformation matrix
980
for (unsigned int row = 0; row < num_derivatives; row++)
982
for (unsigned int col = 0; col < num_derivatives; col++)
984
for (unsigned int k = 0; k < n; k++)
985
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
990
for (unsigned int j = 0; j < 1*num_derivatives; j++)
993
// Map degree of freedom to element degree of freedom
994
const unsigned int dof = i;
997
const double scalings_y_0 = 1;
998
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
999
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
1000
const double scalings_z_0 = 1;
1001
const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
1002
const double scalings_z_2 = scalings_z_1*(0.5 - 0.5*z);
1004
// Compute psitilde_a
1005
const double psitilde_a_0 = 1;
1006
const double psitilde_a_1 = x;
1007
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
1009
// Compute psitilde_bs
1010
const double psitilde_bs_0_0 = 1;
1011
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1012
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;
1013
const double psitilde_bs_1_0 = 1;
1014
const double psitilde_bs_1_1 = 2.5*y + 1.5;
1015
const double psitilde_bs_2_0 = 1;
1017
// Compute psitilde_cs
1018
const double psitilde_cs_00_0 = 1;
1019
const double psitilde_cs_00_1 = 2*z + 1;
1020
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;
1021
const double psitilde_cs_01_0 = 1;
1022
const double psitilde_cs_01_1 = 3*z + 2;
1023
const double psitilde_cs_02_0 = 1;
1024
const double psitilde_cs_10_0 = 1;
1025
const double psitilde_cs_10_1 = 3*z + 2;
1026
const double psitilde_cs_11_0 = 1;
1027
const double psitilde_cs_20_0 = 1;
1029
// Compute basisvalues
1030
const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
1031
const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
1032
const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
1033
const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
1034
const double basisvalue4 = 5.1234753829798*psitilde_a_2*scalings_y_2*psitilde_bs_2_0*scalings_z_2*psitilde_cs_20_0;
1035
const double basisvalue5 = 3.96862696659689*psitilde_a_1*scalings_y_1*psitilde_bs_1_1*scalings_z_2*psitilde_cs_11_0;
1036
const double basisvalue6 = 2.29128784747792*psitilde_a_0*scalings_y_0*psitilde_bs_0_2*scalings_z_2*psitilde_cs_02_0;
1037
const double basisvalue7 = 3.24037034920393*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_1;
1038
const double basisvalue8 = 1.87082869338697*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_1;
1039
const double basisvalue9 = 1.3228756555323*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_2;
1041
// Table(s) of coefficients
1042
static const double coefficients0[10][10] = \
1043
{{-0.0577350269189626, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.0503952630678969, 0.0290957186981323, 0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
1044
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, 0.0290957186981323, -0.0411475599898912, 0.0237565548366599, 0.0167984210226323},
1045
{-0.0577350269189626, 0, 0.0702728368926307, -0.0248451997499977, 0, 0, 0.0872871560943969, 0, -0.0475131096733199, 0.0167984210226323},
1046
{-0.0577350269189626, 0, 0, 0.074535599249993, 0, 0, 0, 0, 0, 0.100790526135794},
1047
{0.23094010767585, 0, 0.140545673785261, 0.0993807989999906, 0, 0, 0, 0, 0.1187827741833, -0.0671936840905293},
1048
{0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, 0.102868899974728, -0.0593913870916499, -0.0671936840905293},
1049
{0.23094010767585, 0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0, 0.100790526135794, -0.0872871560943969, -0.0205737799949456, -0.01187827741833, 0.0167984210226323},
1050
{0.23094010767585, -0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0, 0, 0, -0.102868899974728, -0.0593913870916499, -0.0671936840905293},
1051
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999907, 0, -0.100790526135794, -0.0872871560943969, 0.0205737799949456, -0.01187827741833, 0.0167984210226323},
1052
{0.23094010767585, 0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0, 0.0290957186981323, 0, 0.02375655483666, 0.0167984210226323}};
1054
// Interesting (new) part
1055
// Tables of derivatives of the polynomial base (transpose)
1056
static const double dmats0[10][10] = \
1057
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1058
{6.32455532033676, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1059
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1060
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1061
{0, 11.2249721603218, 0, 0, 0, 0, 0, 0, 0, 0},
1062
{4.58257569495584, 0, 8.36660026534076, -1.18321595661992, 0, 0, 0, 0, 0, 0},
1063
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1064
{3.74165738677394, 0, 0, 8.69482604771366, 0, 0, 0, 0, 0, 0},
1065
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1066
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1068
static const double dmats1[10][10] = \
1069
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1070
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1071
{5.47722557505166, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1072
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1073
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
1074
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
1075
{-2.64575131106459, 0, 9.66091783079296, 0.683130051063973, 0, 0, 0, 0, 0, 0},
1076
{1.87082869338697, 0, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
1077
{3.24037034920393, 0, 0, 7.52994023880668, 0, 0, 0, 0, 0, 0},
1078
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1080
static const double dmats2[10][10] = \
1081
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1082
{3.16227766016838, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1083
{1.82574185835055, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1084
{5.16397779494322, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1085
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0, 0, 0, 0, 0, 0},
1086
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0, 0, 0, 0, 0, 0},
1087
{1.3228756555323, 0, 3.86436713231718, -0.341565025531986, 0, 0, 0, 0, 0, 0},
1088
{1.87082869338697, 7.09929573971954, 0, 4.34741302385683, 0, 0, 0, 0, 0, 0},
1089
{1.08012344973464, 0, 7.09929573971954, 2.50998007960223, 0, 0, 0, 0, 0, 0},
1090
{-3.81881307912987, 0, 0, 8.87411967464942, 0, 0, 0, 0, 0, 0}};
1092
// Compute reference derivatives
1093
// Declare pointer to array of derivatives on FIAT element
1094
double *derivatives = new double [num_derivatives];
1096
// Declare coefficients
1097
double coeff0_0 = 0;
1098
double coeff0_1 = 0;
1099
double coeff0_2 = 0;
1100
double coeff0_3 = 0;
1101
double coeff0_4 = 0;
1102
double coeff0_5 = 0;
1103
double coeff0_6 = 0;
1104
double coeff0_7 = 0;
1105
double coeff0_8 = 0;
1106
double coeff0_9 = 0;
1108
// Declare new coefficients
1109
double new_coeff0_0 = 0;
1110
double new_coeff0_1 = 0;
1111
double new_coeff0_2 = 0;
1112
double new_coeff0_3 = 0;
1113
double new_coeff0_4 = 0;
1114
double new_coeff0_5 = 0;
1115
double new_coeff0_6 = 0;
1116
double new_coeff0_7 = 0;
1117
double new_coeff0_8 = 0;
1118
double new_coeff0_9 = 0;
1120
// Loop possible derivatives
1121
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1123
// Get values from coefficients array
1124
new_coeff0_0 = coefficients0[dof][0];
1125
new_coeff0_1 = coefficients0[dof][1];
1126
new_coeff0_2 = coefficients0[dof][2];
1127
new_coeff0_3 = coefficients0[dof][3];
1128
new_coeff0_4 = coefficients0[dof][4];
1129
new_coeff0_5 = coefficients0[dof][5];
1130
new_coeff0_6 = coefficients0[dof][6];
1131
new_coeff0_7 = coefficients0[dof][7];
1132
new_coeff0_8 = coefficients0[dof][8];
1133
new_coeff0_9 = coefficients0[dof][9];
1135
// Loop derivative order
1136
for (unsigned int j = 0; j < n; j++)
1138
// Update old coefficients
1139
coeff0_0 = new_coeff0_0;
1140
coeff0_1 = new_coeff0_1;
1141
coeff0_2 = new_coeff0_2;
1142
coeff0_3 = new_coeff0_3;
1143
coeff0_4 = new_coeff0_4;
1144
coeff0_5 = new_coeff0_5;
1145
coeff0_6 = new_coeff0_6;
1146
coeff0_7 = new_coeff0_7;
1147
coeff0_8 = new_coeff0_8;
1148
coeff0_9 = new_coeff0_9;
1150
if(combinations[deriv_num][j] == 0)
1152
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];
1153
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];
1154
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];
1155
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];
1156
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];
1157
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];
1158
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];
1159
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];
1160
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];
1161
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];
1163
if(combinations[deriv_num][j] == 1)
1165
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];
1166
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];
1167
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];
1168
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];
1169
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];
1170
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];
1171
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];
1172
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];
1173
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];
1174
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];
1176
if(combinations[deriv_num][j] == 2)
1178
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];
1179
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];
1180
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];
1181
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];
1182
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];
1183
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];
1184
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];
1185
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];
1186
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];
1187
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];
1191
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1192
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;
1195
// Transform derivatives back to physical element
1196
for (unsigned int row = 0; row < num_derivatives; row++)
1198
for (unsigned int col = 0; col < num_derivatives; col++)
1200
values[row] += transform[row][col]*derivatives[col];
1203
// Delete pointer to array of derivatives on FIAT element
1204
delete [] derivatives;
1206
// Delete pointer to array of combinations of derivatives and transform
1207
for (unsigned int row = 0; row < num_derivatives; row++)
1209
delete [] combinations[row];
1210
delete [] transform[row];
1213
delete [] combinations;
1214
delete [] transform;
1217
/// Evaluate order n derivatives of all basis functions at given point in cell
1218
virtual void evaluate_basis_derivatives_all(unsigned int n,
1220
const double* coordinates,
1221
const ufc::cell& c) const
1223
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1226
/// Evaluate linear functional for dof i on the function f
1227
virtual double evaluate_dof(unsigned int i,
1228
const ufc::function& f,
1229
const ufc::cell& c) const
1231
// The reference points, direction and weights:
1232
static const 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}}};
1233
static const double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
1234
static const double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
1236
const double * const * x = c.coordinates;
1237
double result = 0.0;
1238
// Iterate over the points:
1239
// Evaluate basis functions for affine mapping
1240
const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
1241
const double w1 = X[i][0][0];
1242
const double w2 = X[i][0][1];
1243
const double w3 = X[i][0][2];
1245
// Compute affine mapping y = F(X)
1247
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
1248
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
1249
y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
1251
// Evaluate function at physical points
1253
f.evaluate(values, y, c);
1255
// Map function values using appropriate mapping
1256
// Affine map: Do nothing
1258
// Note that we do not map the weights (yet).
1260
// Take directional components
1261
for(int k = 0; k < 1; k++)
1262
result += values[k]*D[i][0][k];
1263
// Multiply by weights
1269
/// Evaluate linear functionals for all dofs on the function f
1270
virtual void evaluate_dofs(double* values,
1271
const ufc::function& f,
1272
const ufc::cell& c) const
1274
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1277
/// Interpolate vertex values from dof values
1278
virtual void interpolate_vertex_values(double* vertex_values,
1279
const double* dof_values,
1280
const ufc::cell& c) const
1282
// Evaluate at vertices and use affine mapping
1283
vertex_values[0] = dof_values[0];
1284
vertex_values[1] = dof_values[1];
1285
vertex_values[2] = dof_values[2];
1286
vertex_values[3] = dof_values[3];
1289
/// Return the number of sub elements (for a mixed element)
1290
virtual unsigned int num_sub_elements() const
1295
/// Create a new finite element for sub element i (for a mixed element)
1296
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1298
return new projection_0_finite_element_1();
1303
/// This class defines the interface for a local-to-global mapping of
1304
/// degrees of freedom (dofs).
1306
class projection_0_dof_map_0: public ufc::dof_map
1310
unsigned int __global_dimension;
1315
projection_0_dof_map_0() : ufc::dof_map()
1317
__global_dimension = 0;
1321
virtual ~projection_0_dof_map_0()
1326
/// Return a string identifying the dof map
1327
virtual const char* signature() const
1329
return "FFC dof map for FiniteElement('Lagrange', 'tetrahedron', 2)";
1332
/// Return true iff mesh entities of topological dimension d are needed
1333
virtual bool needs_mesh_entities(unsigned int d) const
1353
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1354
virtual bool init_mesh(const ufc::mesh& m)
1356
__global_dimension = m.num_entities[0] + m.num_entities[1];
1360
/// Initialize dof map for given cell
1361
virtual void init_cell(const ufc::mesh& m,
1367
/// Finish initialization of dof map for cells
1368
virtual void init_cell_finalize()
1373
/// Return the dimension of the global finite element function space
1374
virtual unsigned int global_dimension() const
1376
return __global_dimension;
1379
/// Return the dimension of the local finite element function space for a cell
1380
virtual unsigned int local_dimension(const ufc::cell& c) const
1385
/// Return the maximum dimension of the local finite element function space
1386
virtual unsigned int max_local_dimension() const
1391
// Return the geometric dimension of the coordinates this dof map provides
1392
virtual unsigned int geometric_dimension() const
1397
/// Return the number of dofs on each cell facet
1398
virtual unsigned int num_facet_dofs() const
1403
/// Return the number of dofs associated with each cell entity of dimension d
1404
virtual unsigned int num_entity_dofs(unsigned int d) const
1406
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1409
/// Tabulate the local-to-global mapping of dofs on a cell
1410
virtual void tabulate_dofs(unsigned int* dofs,
1412
const ufc::cell& c) const
1414
dofs[0] = c.entity_indices[0][0];
1415
dofs[1] = c.entity_indices[0][1];
1416
dofs[2] = c.entity_indices[0][2];
1417
dofs[3] = c.entity_indices[0][3];
1418
unsigned int offset = m.num_entities[0];
1419
dofs[4] = offset + c.entity_indices[1][0];
1420
dofs[5] = offset + c.entity_indices[1][1];
1421
dofs[6] = offset + c.entity_indices[1][2];
1422
dofs[7] = offset + c.entity_indices[1][3];
1423
dofs[8] = offset + c.entity_indices[1][4];
1424
dofs[9] = offset + c.entity_indices[1][5];
1427
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1428
virtual void tabulate_facet_dofs(unsigned int* dofs,
1429
unsigned int facet) const
1468
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1469
virtual void tabulate_entity_dofs(unsigned int* dofs,
1470
unsigned int d, unsigned int i) const
1472
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1475
/// Tabulate the coordinates of all dofs on a cell
1476
virtual void tabulate_coordinates(double** coordinates,
1477
const ufc::cell& c) const
1479
const double * const * x = c.coordinates;
1480
coordinates[0][0] = x[0][0];
1481
coordinates[0][1] = x[0][1];
1482
coordinates[0][2] = x[0][2];
1483
coordinates[1][0] = x[1][0];
1484
coordinates[1][1] = x[1][1];
1485
coordinates[1][2] = x[1][2];
1486
coordinates[2][0] = x[2][0];
1487
coordinates[2][1] = x[2][1];
1488
coordinates[2][2] = x[2][2];
1489
coordinates[3][0] = x[3][0];
1490
coordinates[3][1] = x[3][1];
1491
coordinates[3][2] = x[3][2];
1492
coordinates[4][0] = 0.5*x[2][0] + 0.5*x[3][0];
1493
coordinates[4][1] = 0.5*x[2][1] + 0.5*x[3][1];
1494
coordinates[4][2] = 0.5*x[2][2] + 0.5*x[3][2];
1495
coordinates[5][0] = 0.5*x[1][0] + 0.5*x[3][0];
1496
coordinates[5][1] = 0.5*x[1][1] + 0.5*x[3][1];
1497
coordinates[5][2] = 0.5*x[1][2] + 0.5*x[3][2];
1498
coordinates[6][0] = 0.5*x[1][0] + 0.5*x[2][0];
1499
coordinates[6][1] = 0.5*x[1][1] + 0.5*x[2][1];
1500
coordinates[6][2] = 0.5*x[1][2] + 0.5*x[2][2];
1501
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[3][0];
1502
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[3][1];
1503
coordinates[7][2] = 0.5*x[0][2] + 0.5*x[3][2];
1504
coordinates[8][0] = 0.5*x[0][0] + 0.5*x[2][0];
1505
coordinates[8][1] = 0.5*x[0][1] + 0.5*x[2][1];
1506
coordinates[8][2] = 0.5*x[0][2] + 0.5*x[2][2];
1507
coordinates[9][0] = 0.5*x[0][0] + 0.5*x[1][0];
1508
coordinates[9][1] = 0.5*x[0][1] + 0.5*x[1][1];
1509
coordinates[9][2] = 0.5*x[0][2] + 0.5*x[1][2];
1512
/// Return the number of sub dof maps (for a mixed element)
1513
virtual unsigned int num_sub_dof_maps() const
1518
/// Create a new dof_map for sub dof map i (for a mixed element)
1519
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1521
return new projection_0_dof_map_0();
1526
/// This class defines the interface for a local-to-global mapping of
1527
/// degrees of freedom (dofs).
1529
class projection_0_dof_map_1: public ufc::dof_map
1533
unsigned int __global_dimension;
1538
projection_0_dof_map_1() : ufc::dof_map()
1540
__global_dimension = 0;
1544
virtual ~projection_0_dof_map_1()
1549
/// Return a string identifying the dof map
1550
virtual const char* signature() const
1552
return "FFC dof map for FiniteElement('Lagrange', 'tetrahedron', 2)";
1555
/// Return true iff mesh entities of topological dimension d are needed
1556
virtual bool needs_mesh_entities(unsigned int d) const
1576
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1577
virtual bool init_mesh(const ufc::mesh& m)
1579
__global_dimension = m.num_entities[0] + m.num_entities[1];
1583
/// Initialize dof map for given cell
1584
virtual void init_cell(const ufc::mesh& m,
1590
/// Finish initialization of dof map for cells
1591
virtual void init_cell_finalize()
1596
/// Return the dimension of the global finite element function space
1597
virtual unsigned int global_dimension() const
1599
return __global_dimension;
1602
/// Return the dimension of the local finite element function space for a cell
1603
virtual unsigned int local_dimension(const ufc::cell& c) const
1608
/// Return the maximum dimension of the local finite element function space
1609
virtual unsigned int max_local_dimension() const
1614
// Return the geometric dimension of the coordinates this dof map provides
1615
virtual unsigned int geometric_dimension() const
1620
/// Return the number of dofs on each cell facet
1621
virtual unsigned int num_facet_dofs() const
1626
/// Return the number of dofs associated with each cell entity of dimension d
1627
virtual unsigned int num_entity_dofs(unsigned int d) const
1629
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1632
/// Tabulate the local-to-global mapping of dofs on a cell
1633
virtual void tabulate_dofs(unsigned int* dofs,
1635
const ufc::cell& c) const
1637
dofs[0] = c.entity_indices[0][0];
1638
dofs[1] = c.entity_indices[0][1];
1639
dofs[2] = c.entity_indices[0][2];
1640
dofs[3] = c.entity_indices[0][3];
1641
unsigned int offset = m.num_entities[0];
1642
dofs[4] = offset + c.entity_indices[1][0];
1643
dofs[5] = offset + c.entity_indices[1][1];
1644
dofs[6] = offset + c.entity_indices[1][2];
1645
dofs[7] = offset + c.entity_indices[1][3];
1646
dofs[8] = offset + c.entity_indices[1][4];
1647
dofs[9] = offset + c.entity_indices[1][5];
1650
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1651
virtual void tabulate_facet_dofs(unsigned int* dofs,
1652
unsigned int facet) const
1691
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1692
virtual void tabulate_entity_dofs(unsigned int* dofs,
1693
unsigned int d, unsigned int i) const
1695
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1698
/// Tabulate the coordinates of all dofs on a cell
1699
virtual void tabulate_coordinates(double** coordinates,
1700
const ufc::cell& c) const
1702
const double * const * x = c.coordinates;
1703
coordinates[0][0] = x[0][0];
1704
coordinates[0][1] = x[0][1];
1705
coordinates[0][2] = x[0][2];
1706
coordinates[1][0] = x[1][0];
1707
coordinates[1][1] = x[1][1];
1708
coordinates[1][2] = x[1][2];
1709
coordinates[2][0] = x[2][0];
1710
coordinates[2][1] = x[2][1];
1711
coordinates[2][2] = x[2][2];
1712
coordinates[3][0] = x[3][0];
1713
coordinates[3][1] = x[3][1];
1714
coordinates[3][2] = x[3][2];
1715
coordinates[4][0] = 0.5*x[2][0] + 0.5*x[3][0];
1716
coordinates[4][1] = 0.5*x[2][1] + 0.5*x[3][1];
1717
coordinates[4][2] = 0.5*x[2][2] + 0.5*x[3][2];
1718
coordinates[5][0] = 0.5*x[1][0] + 0.5*x[3][0];
1719
coordinates[5][1] = 0.5*x[1][1] + 0.5*x[3][1];
1720
coordinates[5][2] = 0.5*x[1][2] + 0.5*x[3][2];
1721
coordinates[6][0] = 0.5*x[1][0] + 0.5*x[2][0];
1722
coordinates[6][1] = 0.5*x[1][1] + 0.5*x[2][1];
1723
coordinates[6][2] = 0.5*x[1][2] + 0.5*x[2][2];
1724
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[3][0];
1725
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[3][1];
1726
coordinates[7][2] = 0.5*x[0][2] + 0.5*x[3][2];
1727
coordinates[8][0] = 0.5*x[0][0] + 0.5*x[2][0];
1728
coordinates[8][1] = 0.5*x[0][1] + 0.5*x[2][1];
1729
coordinates[8][2] = 0.5*x[0][2] + 0.5*x[2][2];
1730
coordinates[9][0] = 0.5*x[0][0] + 0.5*x[1][0];
1731
coordinates[9][1] = 0.5*x[0][1] + 0.5*x[1][1];
1732
coordinates[9][2] = 0.5*x[0][2] + 0.5*x[1][2];
1735
/// Return the number of sub dof maps (for a mixed element)
1736
virtual unsigned int num_sub_dof_maps() const
1741
/// Create a new dof_map for sub dof map i (for a mixed element)
1742
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1744
return new projection_0_dof_map_1();