11
11
#include <stdexcept>
15
/// This class defines the interface for a finite element.
17
class UFC_Velocity_finite_element_0_0: public ufc::finite_element
22
UFC_Velocity_finite_element_0_0() : ufc::finite_element()
28
virtual ~UFC_Velocity_finite_element_0_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "FiniteElement('Lagrange', 'triangle', 2)";
39
/// Return the cell shape
40
virtual ufc::shape cell_shape() const
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_10 = element_coordinates[1][1] - element_coordinates[0][1];
76
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
78
// Compute determinant of Jacobian
79
const double detJ = J_00*J_11 - J_01*J_10;
81
// Compute inverse of Jacobian
83
// Get coordinates and map to the reference (UFC) element
84
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
85
element_coordinates[0][0]*element_coordinates[2][1] +\
86
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
87
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
88
element_coordinates[1][0]*element_coordinates[0][1] -\
89
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
91
// Map coordinates to the reference square
92
if (std::abs(y - 1.0) < 1e-14)
95
x = 2.0 *x/(1.0 - y) - 1.0;
101
// Map degree of freedom to element degree of freedom
102
const unsigned int dof = i;
105
const double scalings_y_0 = 1;
106
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
107
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
109
// Compute psitilde_a
110
const double psitilde_a_0 = 1;
111
const double psitilde_a_1 = x;
112
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
114
// Compute psitilde_bs
115
const double psitilde_bs_0_0 = 1;
116
const double psitilde_bs_0_1 = 1.5*y + 0.5;
117
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;
118
const double psitilde_bs_1_0 = 1;
119
const double psitilde_bs_1_1 = 2.5*y + 1.5;
120
const double psitilde_bs_2_0 = 1;
122
// Compute basisvalues
123
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
124
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
125
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
126
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
127
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
128
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
130
// Table(s) of coefficients
131
const static double coefficients0[6][6] = \
132
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
133
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
134
{0, 0, 0.2, 0, 0, 0.163299316185545},
135
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
136
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
137
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
139
// Extract relevant coefficients
140
const double coeff0_0 = coefficients0[dof][0];
141
const double coeff0_1 = coefficients0[dof][1];
142
const double coeff0_2 = coefficients0[dof][2];
143
const double coeff0_3 = coefficients0[dof][3];
144
const double coeff0_4 = coefficients0[dof][4];
145
const double coeff0_5 = coefficients0[dof][5];
148
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
151
/// Evaluate all basis functions at given point in cell
152
virtual void evaluate_basis_all(double* values,
153
const double* coordinates,
154
const ufc::cell& c) const
156
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
159
/// Evaluate order n derivatives of basis function i at given point in cell
160
virtual void evaluate_basis_derivatives(unsigned int i,
163
const double* coordinates,
164
const ufc::cell& c) const
166
// Extract vertex coordinates
167
const double * const * element_coordinates = c.coordinates;
169
// Compute Jacobian of affine map from reference cell
170
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
171
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
172
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
173
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
175
// Compute determinant of Jacobian
176
const double detJ = J_00*J_11 - J_01*J_10;
178
// Compute inverse of Jacobian
180
// Get coordinates and map to the reference (UFC) element
181
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
182
element_coordinates[0][0]*element_coordinates[2][1] +\
183
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
184
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
185
element_coordinates[1][0]*element_coordinates[0][1] -\
186
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
188
// Map coordinates to the reference square
189
if (std::abs(y - 1.0) < 1e-14)
192
x = 2.0 *x/(1.0 - y) - 1.0;
195
// Compute number of derivatives
196
unsigned int num_derivatives = 1;
198
for (unsigned int j = 0; j < n; j++)
199
num_derivatives *= 2;
202
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
203
unsigned int **combinations = new unsigned int *[num_derivatives];
205
for (unsigned int j = 0; j < num_derivatives; j++)
207
combinations[j] = new unsigned int [n];
208
for (unsigned int k = 0; k < n; k++)
209
combinations[j][k] = 0;
212
// Generate combinations of derivatives
213
for (unsigned int row = 1; row < num_derivatives; row++)
215
for (unsigned int num = 0; num < row; num++)
217
for (unsigned int col = n-1; col+1 > 0; col--)
219
if (combinations[row][col] + 1 > 1)
220
combinations[row][col] = 0;
223
combinations[row][col] += 1;
230
// Compute inverse of Jacobian
231
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
233
// Declare transformation matrix
234
// Declare pointer to two dimensional array and initialise
235
double **transform = new double *[num_derivatives];
237
for (unsigned int j = 0; j < num_derivatives; j++)
239
transform[j] = new double [num_derivatives];
240
for (unsigned int k = 0; k < num_derivatives; k++)
244
// Construct transformation matrix
245
for (unsigned int row = 0; row < num_derivatives; row++)
247
for (unsigned int col = 0; col < num_derivatives; col++)
249
for (unsigned int k = 0; k < n; k++)
250
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
255
for (unsigned int j = 0; j < 1*num_derivatives; j++)
258
// Map degree of freedom to element degree of freedom
259
const unsigned int dof = i;
262
const double scalings_y_0 = 1;
263
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
264
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
266
// Compute psitilde_a
267
const double psitilde_a_0 = 1;
268
const double psitilde_a_1 = x;
269
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
271
// Compute psitilde_bs
272
const double psitilde_bs_0_0 = 1;
273
const double psitilde_bs_0_1 = 1.5*y + 0.5;
274
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;
275
const double psitilde_bs_1_0 = 1;
276
const double psitilde_bs_1_1 = 2.5*y + 1.5;
277
const double psitilde_bs_2_0 = 1;
279
// Compute basisvalues
280
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
281
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
282
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
283
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
284
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
285
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
287
// Table(s) of coefficients
288
const static double coefficients0[6][6] = \
289
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
290
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
291
{0, 0, 0.2, 0, 0, 0.163299316185545},
292
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
293
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
294
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
296
// Interesting (new) part
297
// Tables of derivatives of the polynomial base (transpose)
298
const static double dmats0[6][6] = \
300
{4.89897948556635, 0, 0, 0, 0, 0},
302
{0, 9.48683298050514, 0, 0, 0, 0},
303
{4, 0, 7.07106781186548, 0, 0, 0},
306
const static double dmats1[6][6] = \
308
{2.44948974278318, 0, 0, 0, 0, 0},
309
{4.24264068711928, 0, 0, 0, 0, 0},
310
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
311
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
312
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
314
// Compute reference derivatives
315
// Declare pointer to array of derivatives on FIAT element
316
double *derivatives = new double [num_derivatives];
318
// Declare coefficients
326
// Declare new coefficients
327
double new_coeff0_0 = 0;
328
double new_coeff0_1 = 0;
329
double new_coeff0_2 = 0;
330
double new_coeff0_3 = 0;
331
double new_coeff0_4 = 0;
332
double new_coeff0_5 = 0;
334
// Loop possible derivatives
335
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
337
// Get values from coefficients array
338
new_coeff0_0 = coefficients0[dof][0];
339
new_coeff0_1 = coefficients0[dof][1];
340
new_coeff0_2 = coefficients0[dof][2];
341
new_coeff0_3 = coefficients0[dof][3];
342
new_coeff0_4 = coefficients0[dof][4];
343
new_coeff0_5 = coefficients0[dof][5];
345
// Loop derivative order
346
for (unsigned int j = 0; j < n; j++)
348
// Update old coefficients
349
coeff0_0 = new_coeff0_0;
350
coeff0_1 = new_coeff0_1;
351
coeff0_2 = new_coeff0_2;
352
coeff0_3 = new_coeff0_3;
353
coeff0_4 = new_coeff0_4;
354
coeff0_5 = new_coeff0_5;
356
if(combinations[deriv_num][j] == 0)
358
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];
359
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];
360
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];
361
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];
362
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];
363
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];
365
if(combinations[deriv_num][j] == 1)
367
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];
368
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];
369
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];
370
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];
371
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];
372
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];
376
// Compute derivatives on reference element as dot product of coefficients and basisvalues
377
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;
380
// Transform derivatives back to physical element
381
for (unsigned int row = 0; row < num_derivatives; row++)
383
for (unsigned int col = 0; col < num_derivatives; col++)
385
values[row] += transform[row][col]*derivatives[col];
388
// Delete pointer to array of derivatives on FIAT element
389
delete [] derivatives;
391
// Delete pointer to array of combinations of derivatives and transform
392
for (unsigned int row = 0; row < num_derivatives; row++)
394
delete [] combinations[row];
395
delete [] transform[row];
398
delete [] combinations;
402
/// Evaluate order n derivatives of all basis functions at given point in cell
403
virtual void evaluate_basis_derivatives_all(unsigned int n,
405
const double* coordinates,
406
const ufc::cell& c) const
408
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
411
/// Evaluate linear functional for dof i on the function f
412
virtual double evaluate_dof(unsigned int i,
413
const ufc::function& f,
414
const ufc::cell& c) const
416
// The reference points, direction and weights:
417
const static double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
418
const static double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
419
const static double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
421
const double * const * x = c.coordinates;
423
// Iterate over the points:
424
// Evaluate basis functions for affine mapping
425
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
426
const double w1 = X[i][0][0];
427
const double w2 = X[i][0][1];
429
// Compute affine mapping y = F(X)
431
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
432
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
434
// Evaluate function at physical points
436
f.evaluate(values, y, c);
438
// Map function values using appropriate mapping
439
// Affine map: Do nothing
441
// Note that we do not map the weights (yet).
443
// Take directional components
444
for(int k = 0; k < 1; k++)
445
result += values[k]*D[i][0][k];
446
// Multiply by weights
452
/// Evaluate linear functionals for all dofs on the function f
453
virtual void evaluate_dofs(double* values,
454
const ufc::function& f,
455
const ufc::cell& c) const
457
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
460
/// Interpolate vertex values from dof values
461
virtual void interpolate_vertex_values(double* vertex_values,
462
const double* dof_values,
463
const ufc::cell& c) const
465
// Evaluate at vertices and use affine mapping
466
vertex_values[0] = dof_values[0];
467
vertex_values[1] = dof_values[1];
468
vertex_values[2] = dof_values[2];
471
/// Return the number of sub elements (for a mixed element)
472
virtual unsigned int num_sub_elements() const
477
/// Create a new finite element for sub element i (for a mixed element)
478
virtual ufc::finite_element* create_sub_element(unsigned int i) const
480
return new UFC_Velocity_finite_element_0_0();
485
/// This class defines the interface for a finite element.
487
class UFC_Velocity_finite_element_0_1: public ufc::finite_element
492
UFC_Velocity_finite_element_0_1() : ufc::finite_element()
498
virtual ~UFC_Velocity_finite_element_0_1()
503
/// Return a string identifying the finite element
504
virtual const char* signature() const
506
return "FiniteElement('Lagrange', 'triangle', 2)";
509
/// Return the cell shape
510
virtual ufc::shape cell_shape() const
512
return ufc::triangle;
515
/// Return the dimension of the finite element function space
516
virtual unsigned int space_dimension() const
521
/// Return the rank of the value space
522
virtual unsigned int value_rank() const
527
/// Return the dimension of the value space for axis i
528
virtual unsigned int value_dimension(unsigned int i) const
533
/// Evaluate basis function i at given point in cell
534
virtual void evaluate_basis(unsigned int i,
536
const double* coordinates,
537
const ufc::cell& c) const
539
// Extract vertex coordinates
540
const double * const * element_coordinates = c.coordinates;
542
// Compute Jacobian of affine map from reference cell
543
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
544
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
545
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
546
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
548
// Compute determinant of Jacobian
549
const double detJ = J_00*J_11 - J_01*J_10;
551
// Compute inverse of Jacobian
553
// Get coordinates and map to the reference (UFC) element
554
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
555
element_coordinates[0][0]*element_coordinates[2][1] +\
556
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
557
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
558
element_coordinates[1][0]*element_coordinates[0][1] -\
559
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
561
// Map coordinates to the reference square
562
if (std::abs(y - 1.0) < 1e-14)
565
x = 2.0 *x/(1.0 - y) - 1.0;
571
// Map degree of freedom to element degree of freedom
572
const unsigned int dof = i;
575
const double scalings_y_0 = 1;
576
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
577
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
579
// Compute psitilde_a
580
const double psitilde_a_0 = 1;
581
const double psitilde_a_1 = x;
582
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
584
// Compute psitilde_bs
585
const double psitilde_bs_0_0 = 1;
586
const double psitilde_bs_0_1 = 1.5*y + 0.5;
587
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;
588
const double psitilde_bs_1_0 = 1;
589
const double psitilde_bs_1_1 = 2.5*y + 1.5;
590
const double psitilde_bs_2_0 = 1;
592
// Compute basisvalues
593
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
594
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
595
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
596
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
597
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
598
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
600
// Table(s) of coefficients
601
const static double coefficients0[6][6] = \
602
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
603
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
604
{0, 0, 0.2, 0, 0, 0.163299316185545},
605
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
606
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
607
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
609
// Extract relevant coefficients
610
const double coeff0_0 = coefficients0[dof][0];
611
const double coeff0_1 = coefficients0[dof][1];
612
const double coeff0_2 = coefficients0[dof][2];
613
const double coeff0_3 = coefficients0[dof][3];
614
const double coeff0_4 = coefficients0[dof][4];
615
const double coeff0_5 = coefficients0[dof][5];
618
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
621
/// Evaluate all basis functions at given point in cell
622
virtual void evaluate_basis_all(double* values,
623
const double* coordinates,
624
const ufc::cell& c) const
626
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
629
/// Evaluate order n derivatives of basis function i at given point in cell
630
virtual void evaluate_basis_derivatives(unsigned int i,
633
const double* coordinates,
634
const ufc::cell& c) const
636
// Extract vertex coordinates
637
const double * const * element_coordinates = c.coordinates;
639
// Compute Jacobian of affine map from reference cell
640
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
641
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
642
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
643
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
645
// Compute determinant of Jacobian
646
const double detJ = J_00*J_11 - J_01*J_10;
648
// Compute inverse of Jacobian
650
// Get coordinates and map to the reference (UFC) element
651
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
652
element_coordinates[0][0]*element_coordinates[2][1] +\
653
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
654
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
655
element_coordinates[1][0]*element_coordinates[0][1] -\
656
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
658
// Map coordinates to the reference square
659
if (std::abs(y - 1.0) < 1e-14)
662
x = 2.0 *x/(1.0 - y) - 1.0;
665
// Compute number of derivatives
666
unsigned int num_derivatives = 1;
668
for (unsigned int j = 0; j < n; j++)
669
num_derivatives *= 2;
672
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
673
unsigned int **combinations = new unsigned int *[num_derivatives];
675
for (unsigned int j = 0; j < num_derivatives; j++)
677
combinations[j] = new unsigned int [n];
678
for (unsigned int k = 0; k < n; k++)
679
combinations[j][k] = 0;
682
// Generate combinations of derivatives
683
for (unsigned int row = 1; row < num_derivatives; row++)
685
for (unsigned int num = 0; num < row; num++)
687
for (unsigned int col = n-1; col+1 > 0; col--)
689
if (combinations[row][col] + 1 > 1)
690
combinations[row][col] = 0;
693
combinations[row][col] += 1;
700
// Compute inverse of Jacobian
701
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
703
// Declare transformation matrix
704
// Declare pointer to two dimensional array and initialise
705
double **transform = new double *[num_derivatives];
707
for (unsigned int j = 0; j < num_derivatives; j++)
709
transform[j] = new double [num_derivatives];
710
for (unsigned int k = 0; k < num_derivatives; k++)
714
// Construct transformation matrix
715
for (unsigned int row = 0; row < num_derivatives; row++)
717
for (unsigned int col = 0; col < num_derivatives; col++)
719
for (unsigned int k = 0; k < n; k++)
720
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
725
for (unsigned int j = 0; j < 1*num_derivatives; j++)
728
// Map degree of freedom to element degree of freedom
729
const unsigned int dof = i;
732
const double scalings_y_0 = 1;
733
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
734
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
736
// Compute psitilde_a
737
const double psitilde_a_0 = 1;
738
const double psitilde_a_1 = x;
739
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
741
// Compute psitilde_bs
742
const double psitilde_bs_0_0 = 1;
743
const double psitilde_bs_0_1 = 1.5*y + 0.5;
744
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;
745
const double psitilde_bs_1_0 = 1;
746
const double psitilde_bs_1_1 = 2.5*y + 1.5;
747
const double psitilde_bs_2_0 = 1;
749
// Compute basisvalues
750
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
751
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
752
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
753
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
754
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
755
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
757
// Table(s) of coefficients
758
const static double coefficients0[6][6] = \
759
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
760
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
761
{0, 0, 0.2, 0, 0, 0.163299316185545},
762
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
763
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
764
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
766
// Interesting (new) part
767
// Tables of derivatives of the polynomial base (transpose)
768
const static double dmats0[6][6] = \
770
{4.89897948556635, 0, 0, 0, 0, 0},
772
{0, 9.48683298050514, 0, 0, 0, 0},
773
{4, 0, 7.07106781186548, 0, 0, 0},
776
const static double dmats1[6][6] = \
778
{2.44948974278318, 0, 0, 0, 0, 0},
779
{4.24264068711928, 0, 0, 0, 0, 0},
780
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
781
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
782
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
784
// Compute reference derivatives
785
// Declare pointer to array of derivatives on FIAT element
786
double *derivatives = new double [num_derivatives];
788
// Declare coefficients
796
// Declare new coefficients
797
double new_coeff0_0 = 0;
798
double new_coeff0_1 = 0;
799
double new_coeff0_2 = 0;
800
double new_coeff0_3 = 0;
801
double new_coeff0_4 = 0;
802
double new_coeff0_5 = 0;
804
// Loop possible derivatives
805
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
807
// Get values from coefficients array
808
new_coeff0_0 = coefficients0[dof][0];
809
new_coeff0_1 = coefficients0[dof][1];
810
new_coeff0_2 = coefficients0[dof][2];
811
new_coeff0_3 = coefficients0[dof][3];
812
new_coeff0_4 = coefficients0[dof][4];
813
new_coeff0_5 = coefficients0[dof][5];
815
// Loop derivative order
816
for (unsigned int j = 0; j < n; j++)
818
// Update old coefficients
819
coeff0_0 = new_coeff0_0;
820
coeff0_1 = new_coeff0_1;
821
coeff0_2 = new_coeff0_2;
822
coeff0_3 = new_coeff0_3;
823
coeff0_4 = new_coeff0_4;
824
coeff0_5 = new_coeff0_5;
826
if(combinations[deriv_num][j] == 0)
828
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];
829
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];
830
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];
831
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];
832
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];
833
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];
835
if(combinations[deriv_num][j] == 1)
837
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];
838
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];
839
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];
840
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];
841
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];
842
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];
846
// Compute derivatives on reference element as dot product of coefficients and basisvalues
847
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;
850
// Transform derivatives back to physical element
851
for (unsigned int row = 0; row < num_derivatives; row++)
853
for (unsigned int col = 0; col < num_derivatives; col++)
855
values[row] += transform[row][col]*derivatives[col];
858
// Delete pointer to array of derivatives on FIAT element
859
delete [] derivatives;
861
// Delete pointer to array of combinations of derivatives and transform
862
for (unsigned int row = 0; row < num_derivatives; row++)
864
delete [] combinations[row];
865
delete [] transform[row];
868
delete [] combinations;
872
/// Evaluate order n derivatives of all basis functions at given point in cell
873
virtual void evaluate_basis_derivatives_all(unsigned int n,
875
const double* coordinates,
876
const ufc::cell& c) const
878
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
881
/// Evaluate linear functional for dof i on the function f
882
virtual double evaluate_dof(unsigned int i,
883
const ufc::function& f,
884
const ufc::cell& c) const
886
// The reference points, direction and weights:
887
const static double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
888
const static double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
889
const static double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
891
const double * const * x = c.coordinates;
893
// Iterate over the points:
894
// Evaluate basis functions for affine mapping
895
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
896
const double w1 = X[i][0][0];
897
const double w2 = X[i][0][1];
899
// Compute affine mapping y = F(X)
901
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
902
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
904
// Evaluate function at physical points
906
f.evaluate(values, y, c);
908
// Map function values using appropriate mapping
909
// Affine map: Do nothing
911
// Note that we do not map the weights (yet).
913
// Take directional components
914
for(int k = 0; k < 1; k++)
915
result += values[k]*D[i][0][k];
916
// Multiply by weights
922
/// Evaluate linear functionals for all dofs on the function f
923
virtual void evaluate_dofs(double* values,
924
const ufc::function& f,
925
const ufc::cell& c) const
927
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
930
/// Interpolate vertex values from dof values
931
virtual void interpolate_vertex_values(double* vertex_values,
932
const double* dof_values,
933
const ufc::cell& c) const
935
// Evaluate at vertices and use affine mapping
936
vertex_values[0] = dof_values[0];
937
vertex_values[1] = dof_values[1];
938
vertex_values[2] = dof_values[2];
941
/// Return the number of sub elements (for a mixed element)
942
virtual unsigned int num_sub_elements() const
947
/// Create a new finite element for sub element i (for a mixed element)
948
virtual ufc::finite_element* create_sub_element(unsigned int i) const
950
return new UFC_Velocity_finite_element_0_1();
955
/// This class defines the interface for a finite element.
957
class UFC_Velocity_finite_element_0: public ufc::finite_element
962
UFC_Velocity_finite_element_0() : ufc::finite_element()
968
virtual ~UFC_Velocity_finite_element_0()
15
/// This class defines the interface for a finite element.
17
class velocity_0_finite_element_0_0: public ufc::finite_element
22
velocity_0_finite_element_0_0() : ufc::finite_element()
28
virtual ~velocity_0_finite_element_0_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "FiniteElement('Lagrange', 'triangle', 2)";
39
/// Return the cell shape
40
virtual ufc::shape cell_shape() const
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_10 = element_coordinates[1][1] - element_coordinates[0][1];
76
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
78
// Compute determinant of Jacobian
79
const double detJ = J_00*J_11 - J_01*J_10;
81
// Compute inverse of Jacobian
83
// Get coordinates and map to the reference (UFC) element
84
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
85
element_coordinates[0][0]*element_coordinates[2][1] +\
86
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
87
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
88
element_coordinates[1][0]*element_coordinates[0][1] -\
89
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
91
// Map coordinates to the reference square
92
if (std::abs(y - 1.0) < 1e-14)
95
x = 2.0 *x/(1.0 - y) - 1.0;
101
// Map degree of freedom to element degree of freedom
102
const unsigned int dof = i;
105
const double scalings_y_0 = 1;
106
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
107
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
109
// Compute psitilde_a
110
const double psitilde_a_0 = 1;
111
const double psitilde_a_1 = x;
112
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
114
// Compute psitilde_bs
115
const double psitilde_bs_0_0 = 1;
116
const double psitilde_bs_0_1 = 1.5*y + 0.5;
117
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;
118
const double psitilde_bs_1_0 = 1;
119
const double psitilde_bs_1_1 = 2.5*y + 1.5;
120
const double psitilde_bs_2_0 = 1;
122
// Compute basisvalues
123
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
124
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
125
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
126
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
127
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
128
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
130
// Table(s) of coefficients
131
static const double coefficients0[6][6] = \
132
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
133
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
134
{0, 0, 0.2, 0, 0, 0.163299316185545},
135
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
136
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
137
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
139
// Extract relevant coefficients
140
const double coeff0_0 = coefficients0[dof][0];
141
const double coeff0_1 = coefficients0[dof][1];
142
const double coeff0_2 = coefficients0[dof][2];
143
const double coeff0_3 = coefficients0[dof][3];
144
const double coeff0_4 = coefficients0[dof][4];
145
const double coeff0_5 = coefficients0[dof][5];
148
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
151
/// Evaluate all basis functions at given point in cell
152
virtual void evaluate_basis_all(double* values,
153
const double* coordinates,
154
const ufc::cell& c) const
156
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
159
/// Evaluate order n derivatives of basis function i at given point in cell
160
virtual void evaluate_basis_derivatives(unsigned int i,
163
const double* coordinates,
164
const ufc::cell& c) const
166
// Extract vertex coordinates
167
const double * const * element_coordinates = c.coordinates;
169
// Compute Jacobian of affine map from reference cell
170
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
171
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
172
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
173
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
175
// Compute determinant of Jacobian
176
const double detJ = J_00*J_11 - J_01*J_10;
178
// Compute inverse of Jacobian
180
// Get coordinates and map to the reference (UFC) element
181
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
182
element_coordinates[0][0]*element_coordinates[2][1] +\
183
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
184
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
185
element_coordinates[1][0]*element_coordinates[0][1] -\
186
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
188
// Map coordinates to the reference square
189
if (std::abs(y - 1.0) < 1e-14)
192
x = 2.0 *x/(1.0 - y) - 1.0;
195
// Compute number of derivatives
196
unsigned int num_derivatives = 1;
198
for (unsigned int j = 0; j < n; j++)
199
num_derivatives *= 2;
202
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
203
unsigned int **combinations = new unsigned int *[num_derivatives];
205
for (unsigned int j = 0; j < num_derivatives; j++)
207
combinations[j] = new unsigned int [n];
208
for (unsigned int k = 0; k < n; k++)
209
combinations[j][k] = 0;
212
// Generate combinations of derivatives
213
for (unsigned int row = 1; row < num_derivatives; row++)
215
for (unsigned int num = 0; num < row; num++)
217
for (unsigned int col = n-1; col+1 > 0; col--)
219
if (combinations[row][col] + 1 > 1)
220
combinations[row][col] = 0;
223
combinations[row][col] += 1;
230
// Compute inverse of Jacobian
231
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
233
// Declare transformation matrix
234
// Declare pointer to two dimensional array and initialise
235
double **transform = new double *[num_derivatives];
237
for (unsigned int j = 0; j < num_derivatives; j++)
239
transform[j] = new double [num_derivatives];
240
for (unsigned int k = 0; k < num_derivatives; k++)
244
// Construct transformation matrix
245
for (unsigned int row = 0; row < num_derivatives; row++)
247
for (unsigned int col = 0; col < num_derivatives; col++)
249
for (unsigned int k = 0; k < n; k++)
250
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
255
for (unsigned int j = 0; j < 1*num_derivatives; j++)
258
// Map degree of freedom to element degree of freedom
259
const unsigned int dof = i;
262
const double scalings_y_0 = 1;
263
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
264
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
266
// Compute psitilde_a
267
const double psitilde_a_0 = 1;
268
const double psitilde_a_1 = x;
269
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
271
// Compute psitilde_bs
272
const double psitilde_bs_0_0 = 1;
273
const double psitilde_bs_0_1 = 1.5*y + 0.5;
274
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;
275
const double psitilde_bs_1_0 = 1;
276
const double psitilde_bs_1_1 = 2.5*y + 1.5;
277
const double psitilde_bs_2_0 = 1;
279
// Compute basisvalues
280
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
281
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
282
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
283
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
284
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
285
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
287
// Table(s) of coefficients
288
static const double coefficients0[6][6] = \
289
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
290
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
291
{0, 0, 0.2, 0, 0, 0.163299316185545},
292
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
293
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
294
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
296
// Interesting (new) part
297
// Tables of derivatives of the polynomial base (transpose)
298
static const double dmats0[6][6] = \
300
{4.89897948556635, 0, 0, 0, 0, 0},
302
{0, 9.48683298050514, 0, 0, 0, 0},
303
{4, 0, 7.07106781186548, 0, 0, 0},
306
static const double dmats1[6][6] = \
308
{2.44948974278318, 0, 0, 0, 0, 0},
309
{4.24264068711928, 0, 0, 0, 0, 0},
310
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
311
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
312
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
314
// Compute reference derivatives
315
// Declare pointer to array of derivatives on FIAT element
316
double *derivatives = new double [num_derivatives];
318
// Declare coefficients
326
// Declare new coefficients
327
double new_coeff0_0 = 0;
328
double new_coeff0_1 = 0;
329
double new_coeff0_2 = 0;
330
double new_coeff0_3 = 0;
331
double new_coeff0_4 = 0;
332
double new_coeff0_5 = 0;
334
// Loop possible derivatives
335
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
337
// Get values from coefficients array
338
new_coeff0_0 = coefficients0[dof][0];
339
new_coeff0_1 = coefficients0[dof][1];
340
new_coeff0_2 = coefficients0[dof][2];
341
new_coeff0_3 = coefficients0[dof][3];
342
new_coeff0_4 = coefficients0[dof][4];
343
new_coeff0_5 = coefficients0[dof][5];
345
// Loop derivative order
346
for (unsigned int j = 0; j < n; j++)
348
// Update old coefficients
349
coeff0_0 = new_coeff0_0;
350
coeff0_1 = new_coeff0_1;
351
coeff0_2 = new_coeff0_2;
352
coeff0_3 = new_coeff0_3;
353
coeff0_4 = new_coeff0_4;
354
coeff0_5 = new_coeff0_5;
356
if(combinations[deriv_num][j] == 0)
358
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];
359
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];
360
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];
361
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];
362
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];
363
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];
365
if(combinations[deriv_num][j] == 1)
367
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];
368
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];
369
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];
370
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];
371
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];
372
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];
376
// Compute derivatives on reference element as dot product of coefficients and basisvalues
377
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;
380
// Transform derivatives back to physical element
381
for (unsigned int row = 0; row < num_derivatives; row++)
383
for (unsigned int col = 0; col < num_derivatives; col++)
385
values[row] += transform[row][col]*derivatives[col];
388
// Delete pointer to array of derivatives on FIAT element
389
delete [] derivatives;
391
// Delete pointer to array of combinations of derivatives and transform
392
for (unsigned int row = 0; row < num_derivatives; row++)
394
delete [] combinations[row];
395
delete [] transform[row];
398
delete [] combinations;
402
/// Evaluate order n derivatives of all basis functions at given point in cell
403
virtual void evaluate_basis_derivatives_all(unsigned int n,
405
const double* coordinates,
406
const ufc::cell& c) const
408
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
411
/// Evaluate linear functional for dof i on the function f
412
virtual double evaluate_dof(unsigned int i,
413
const ufc::function& f,
414
const ufc::cell& c) const
416
// The reference points, direction and weights:
417
static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
418
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
419
static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
421
const double * const * x = c.coordinates;
423
// Iterate over the points:
424
// Evaluate basis functions for affine mapping
425
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
426
const double w1 = X[i][0][0];
427
const double w2 = X[i][0][1];
429
// Compute affine mapping y = F(X)
431
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
432
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
434
// Evaluate function at physical points
436
f.evaluate(values, y, c);
438
// Map function values using appropriate mapping
439
// Affine map: Do nothing
441
// Note that we do not map the weights (yet).
443
// Take directional components
444
for(int k = 0; k < 1; k++)
445
result += values[k]*D[i][0][k];
446
// Multiply by weights
452
/// Evaluate linear functionals for all dofs on the function f
453
virtual void evaluate_dofs(double* values,
454
const ufc::function& f,
455
const ufc::cell& c) const
457
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
460
/// Interpolate vertex values from dof values
461
virtual void interpolate_vertex_values(double* vertex_values,
462
const double* dof_values,
463
const ufc::cell& c) const
465
// Evaluate at vertices and use affine mapping
466
vertex_values[0] = dof_values[0];
467
vertex_values[1] = dof_values[1];
468
vertex_values[2] = dof_values[2];
471
/// Return the number of sub elements (for a mixed element)
472
virtual unsigned int num_sub_elements() const
477
/// Create a new finite element for sub element i (for a mixed element)
478
virtual ufc::finite_element* create_sub_element(unsigned int i) const
480
return new velocity_0_finite_element_0_0();
485
/// This class defines the interface for a finite element.
487
class velocity_0_finite_element_0_1: public ufc::finite_element
492
velocity_0_finite_element_0_1() : ufc::finite_element()
498
virtual ~velocity_0_finite_element_0_1()
503
/// Return a string identifying the finite element
504
virtual const char* signature() const
506
return "FiniteElement('Lagrange', 'triangle', 2)";
509
/// Return the cell shape
510
virtual ufc::shape cell_shape() const
512
return ufc::triangle;
515
/// Return the dimension of the finite element function space
516
virtual unsigned int space_dimension() const
521
/// Return the rank of the value space
522
virtual unsigned int value_rank() const
527
/// Return the dimension of the value space for axis i
528
virtual unsigned int value_dimension(unsigned int i) const
533
/// Evaluate basis function i at given point in cell
534
virtual void evaluate_basis(unsigned int i,
536
const double* coordinates,
537
const ufc::cell& c) const
539
// Extract vertex coordinates
540
const double * const * element_coordinates = c.coordinates;
542
// Compute Jacobian of affine map from reference cell
543
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
544
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
545
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
546
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
548
// Compute determinant of Jacobian
549
const double detJ = J_00*J_11 - J_01*J_10;
551
// Compute inverse of Jacobian
553
// Get coordinates and map to the reference (UFC) element
554
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
555
element_coordinates[0][0]*element_coordinates[2][1] +\
556
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
557
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
558
element_coordinates[1][0]*element_coordinates[0][1] -\
559
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
561
// Map coordinates to the reference square
562
if (std::abs(y - 1.0) < 1e-14)
565
x = 2.0 *x/(1.0 - y) - 1.0;
571
// Map degree of freedom to element degree of freedom
572
const unsigned int dof = i;
575
const double scalings_y_0 = 1;
576
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
577
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
579
// Compute psitilde_a
580
const double psitilde_a_0 = 1;
581
const double psitilde_a_1 = x;
582
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
584
// Compute psitilde_bs
585
const double psitilde_bs_0_0 = 1;
586
const double psitilde_bs_0_1 = 1.5*y + 0.5;
587
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;
588
const double psitilde_bs_1_0 = 1;
589
const double psitilde_bs_1_1 = 2.5*y + 1.5;
590
const double psitilde_bs_2_0 = 1;
592
// Compute basisvalues
593
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
594
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
595
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
596
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
597
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
598
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
600
// Table(s) of coefficients
601
static const double coefficients0[6][6] = \
602
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
603
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
604
{0, 0, 0.2, 0, 0, 0.163299316185545},
605
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
606
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
607
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
609
// Extract relevant coefficients
610
const double coeff0_0 = coefficients0[dof][0];
611
const double coeff0_1 = coefficients0[dof][1];
612
const double coeff0_2 = coefficients0[dof][2];
613
const double coeff0_3 = coefficients0[dof][3];
614
const double coeff0_4 = coefficients0[dof][4];
615
const double coeff0_5 = coefficients0[dof][5];
618
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
621
/// Evaluate all basis functions at given point in cell
622
virtual void evaluate_basis_all(double* values,
623
const double* coordinates,
624
const ufc::cell& c) const
626
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
629
/// Evaluate order n derivatives of basis function i at given point in cell
630
virtual void evaluate_basis_derivatives(unsigned int i,
633
const double* coordinates,
634
const ufc::cell& c) const
636
// Extract vertex coordinates
637
const double * const * element_coordinates = c.coordinates;
639
// Compute Jacobian of affine map from reference cell
640
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
641
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
642
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
643
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
645
// Compute determinant of Jacobian
646
const double detJ = J_00*J_11 - J_01*J_10;
648
// Compute inverse of Jacobian
650
// Get coordinates and map to the reference (UFC) element
651
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
652
element_coordinates[0][0]*element_coordinates[2][1] +\
653
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
654
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
655
element_coordinates[1][0]*element_coordinates[0][1] -\
656
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
658
// Map coordinates to the reference square
659
if (std::abs(y - 1.0) < 1e-14)
662
x = 2.0 *x/(1.0 - y) - 1.0;
665
// Compute number of derivatives
666
unsigned int num_derivatives = 1;
668
for (unsigned int j = 0; j < n; j++)
669
num_derivatives *= 2;
672
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
673
unsigned int **combinations = new unsigned int *[num_derivatives];
675
for (unsigned int j = 0; j < num_derivatives; j++)
677
combinations[j] = new unsigned int [n];
678
for (unsigned int k = 0; k < n; k++)
679
combinations[j][k] = 0;
682
// Generate combinations of derivatives
683
for (unsigned int row = 1; row < num_derivatives; row++)
685
for (unsigned int num = 0; num < row; num++)
687
for (unsigned int col = n-1; col+1 > 0; col--)
689
if (combinations[row][col] + 1 > 1)
690
combinations[row][col] = 0;
693
combinations[row][col] += 1;
700
// Compute inverse of Jacobian
701
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
703
// Declare transformation matrix
704
// Declare pointer to two dimensional array and initialise
705
double **transform = new double *[num_derivatives];
707
for (unsigned int j = 0; j < num_derivatives; j++)
709
transform[j] = new double [num_derivatives];
710
for (unsigned int k = 0; k < num_derivatives; k++)
714
// Construct transformation matrix
715
for (unsigned int row = 0; row < num_derivatives; row++)
717
for (unsigned int col = 0; col < num_derivatives; col++)
719
for (unsigned int k = 0; k < n; k++)
720
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
725
for (unsigned int j = 0; j < 1*num_derivatives; j++)
728
// Map degree of freedom to element degree of freedom
729
const unsigned int dof = i;
732
const double scalings_y_0 = 1;
733
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
734
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
736
// Compute psitilde_a
737
const double psitilde_a_0 = 1;
738
const double psitilde_a_1 = x;
739
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
741
// Compute psitilde_bs
742
const double psitilde_bs_0_0 = 1;
743
const double psitilde_bs_0_1 = 1.5*y + 0.5;
744
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;
745
const double psitilde_bs_1_0 = 1;
746
const double psitilde_bs_1_1 = 2.5*y + 1.5;
747
const double psitilde_bs_2_0 = 1;
749
// Compute basisvalues
750
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
751
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
752
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
753
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
754
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
755
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
757
// Table(s) of coefficients
758
static const double coefficients0[6][6] = \
759
{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
760
{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
761
{0, 0, 0.2, 0, 0, 0.163299316185545},
762
{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
763
{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
764
{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
766
// Interesting (new) part
767
// Tables of derivatives of the polynomial base (transpose)
768
static const double dmats0[6][6] = \
770
{4.89897948556635, 0, 0, 0, 0, 0},
772
{0, 9.48683298050514, 0, 0, 0, 0},
773
{4, 0, 7.07106781186548, 0, 0, 0},
776
static const double dmats1[6][6] = \
778
{2.44948974278318, 0, 0, 0, 0, 0},
779
{4.24264068711928, 0, 0, 0, 0, 0},
780
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
781
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
782
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
784
// Compute reference derivatives
785
// Declare pointer to array of derivatives on FIAT element
786
double *derivatives = new double [num_derivatives];
788
// Declare coefficients
796
// Declare new coefficients
797
double new_coeff0_0 = 0;
798
double new_coeff0_1 = 0;
799
double new_coeff0_2 = 0;
800
double new_coeff0_3 = 0;
801
double new_coeff0_4 = 0;
802
double new_coeff0_5 = 0;
804
// Loop possible derivatives
805
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
807
// Get values from coefficients array
808
new_coeff0_0 = coefficients0[dof][0];
809
new_coeff0_1 = coefficients0[dof][1];
810
new_coeff0_2 = coefficients0[dof][2];
811
new_coeff0_3 = coefficients0[dof][3];
812
new_coeff0_4 = coefficients0[dof][4];
813
new_coeff0_5 = coefficients0[dof][5];
815
// Loop derivative order
816
for (unsigned int j = 0; j < n; j++)
818
// Update old coefficients
819
coeff0_0 = new_coeff0_0;
820
coeff0_1 = new_coeff0_1;
821
coeff0_2 = new_coeff0_2;
822
coeff0_3 = new_coeff0_3;
823
coeff0_4 = new_coeff0_4;
824
coeff0_5 = new_coeff0_5;
826
if(combinations[deriv_num][j] == 0)
828
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];
829
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];
830
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];
831
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];
832
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];
833
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];
835
if(combinations[deriv_num][j] == 1)
837
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];
838
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];
839
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];
840
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];
841
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];
842
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];
846
// Compute derivatives on reference element as dot product of coefficients and basisvalues
847
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;
850
// Transform derivatives back to physical element
851
for (unsigned int row = 0; row < num_derivatives; row++)
853
for (unsigned int col = 0; col < num_derivatives; col++)
855
values[row] += transform[row][col]*derivatives[col];
858
// Delete pointer to array of derivatives on FIAT element
859
delete [] derivatives;
861
// Delete pointer to array of combinations of derivatives and transform
862
for (unsigned int row = 0; row < num_derivatives; row++)
864
delete [] combinations[row];
865
delete [] transform[row];
868
delete [] combinations;
872
/// Evaluate order n derivatives of all basis functions at given point in cell
873
virtual void evaluate_basis_derivatives_all(unsigned int n,
875
const double* coordinates,
876
const ufc::cell& c) const
878
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
881
/// Evaluate linear functional for dof i on the function f
882
virtual double evaluate_dof(unsigned int i,
883
const ufc::function& f,
884
const ufc::cell& c) const
886
// The reference points, direction and weights:
887
static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
888
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
889
static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
891
const double * const * x = c.coordinates;
893
// Iterate over the points:
894
// Evaluate basis functions for affine mapping
895
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
896
const double w1 = X[i][0][0];
897
const double w2 = X[i][0][1];
899
// Compute affine mapping y = F(X)
901
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
902
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
904
// Evaluate function at physical points
906
f.evaluate(values, y, c);
908
// Map function values using appropriate mapping
909
// Affine map: Do nothing
911
// Note that we do not map the weights (yet).
913
// Take directional components
914
for(int k = 0; k < 1; k++)
915
result += values[k]*D[i][0][k];
916
// Multiply by weights
922
/// Evaluate linear functionals for all dofs on the function f
923
virtual void evaluate_dofs(double* values,
924
const ufc::function& f,
925
const ufc::cell& c) const
927
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
930
/// Interpolate vertex values from dof values
931
virtual void interpolate_vertex_values(double* vertex_values,
932
const double* dof_values,
933
const ufc::cell& c) const
935
// Evaluate at vertices and use affine mapping
936
vertex_values[0] = dof_values[0];
937
vertex_values[1] = dof_values[1];
938
vertex_values[2] = dof_values[2];
941
/// Return the number of sub elements (for a mixed element)
942
virtual unsigned int num_sub_elements() const
947
/// Create a new finite element for sub element i (for a mixed element)
948
virtual ufc::finite_element* create_sub_element(unsigned int i) const
950
return new velocity_0_finite_element_0_1();
955
/// This class defines the interface for a finite element.
957
class velocity_0_finite_element_0: public ufc::finite_element
962
velocity_0_finite_element_0() : ufc::finite_element()
968
virtual ~velocity_0_finite_element_0()
1645
1645
/// This class defines the interface for a local-to-global mapping of
1646
1646
/// degrees of freedom (dofs).
1648
class UFC_Velocity_dof_map_0_0: public ufc::dof_map
1652
unsigned int __global_dimension;
1657
UFC_Velocity_dof_map_0_0() : ufc::dof_map()
1659
__global_dimension = 0;
1663
virtual ~UFC_Velocity_dof_map_0_0()
1668
/// Return a string identifying the dof map
1669
virtual const char* signature() const
1671
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
1674
/// Return true iff mesh entities of topological dimension d are needed
1675
virtual bool needs_mesh_entities(unsigned int d) const
1692
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1693
virtual bool init_mesh(const ufc::mesh& m)
1695
__global_dimension = m.num_entities[0] + m.num_entities[1];
1699
/// Initialize dof map for given cell
1700
virtual void init_cell(const ufc::mesh& m,
1706
/// Finish initialization of dof map for cells
1707
virtual void init_cell_finalize()
1712
/// Return the dimension of the global finite element function space
1713
virtual unsigned int global_dimension() const
1715
return __global_dimension;
1718
/// Return the dimension of the local finite element function space
1719
virtual unsigned int local_dimension() const
1724
// Return the geometric dimension of the coordinates this dof map provides
1725
virtual unsigned int geometric_dimension() const
1730
/// Return the number of dofs on each cell facet
1731
virtual unsigned int num_facet_dofs() const
1736
/// Return the number of dofs associated with each cell entity of dimension d
1737
virtual unsigned int num_entity_dofs(unsigned int d) const
1739
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1742
/// Tabulate the local-to-global mapping of dofs on a cell
1743
virtual void tabulate_dofs(unsigned int* dofs,
1745
const ufc::cell& c) const
1747
dofs[0] = c.entity_indices[0][0];
1748
dofs[1] = c.entity_indices[0][1];
1749
dofs[2] = c.entity_indices[0][2];
1750
unsigned int offset = m.num_entities[0];
1751
dofs[3] = offset + c.entity_indices[1][0];
1752
dofs[4] = offset + c.entity_indices[1][1];
1753
dofs[5] = offset + c.entity_indices[1][2];
1756
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1757
virtual void tabulate_facet_dofs(unsigned int* dofs,
1758
unsigned int facet) const
1780
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1781
virtual void tabulate_entity_dofs(unsigned int* dofs,
1782
unsigned int d, unsigned int i) const
1784
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1787
/// Tabulate the coordinates of all dofs on a cell
1788
virtual void tabulate_coordinates(double** coordinates,
1789
const ufc::cell& c) const
1791
const double * const * x = c.coordinates;
1792
coordinates[0][0] = x[0][0];
1793
coordinates[0][1] = x[0][1];
1794
coordinates[1][0] = x[1][0];
1795
coordinates[1][1] = x[1][1];
1796
coordinates[2][0] = x[2][0];
1797
coordinates[2][1] = x[2][1];
1798
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1799
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1800
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1801
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1802
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1803
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
1806
/// Return the number of sub dof maps (for a mixed element)
1807
virtual unsigned int num_sub_dof_maps() const
1812
/// Create a new dof_map for sub dof map i (for a mixed element)
1813
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1815
return new UFC_Velocity_dof_map_0_0();
1820
/// This class defines the interface for a local-to-global mapping of
1821
/// degrees of freedom (dofs).
1823
class UFC_Velocity_dof_map_0_1: public ufc::dof_map
1827
unsigned int __global_dimension;
1832
UFC_Velocity_dof_map_0_1() : ufc::dof_map()
1834
__global_dimension = 0;
1838
virtual ~UFC_Velocity_dof_map_0_1()
1843
/// Return a string identifying the dof map
1844
virtual const char* signature() const
1846
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
1849
/// Return true iff mesh entities of topological dimension d are needed
1850
virtual bool needs_mesh_entities(unsigned int d) const
1867
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1868
virtual bool init_mesh(const ufc::mesh& m)
1870
__global_dimension = m.num_entities[0] + m.num_entities[1];
1874
/// Initialize dof map for given cell
1875
virtual void init_cell(const ufc::mesh& m,
1881
/// Finish initialization of dof map for cells
1882
virtual void init_cell_finalize()
1887
/// Return the dimension of the global finite element function space
1888
virtual unsigned int global_dimension() const
1890
return __global_dimension;
1893
/// Return the dimension of the local finite element function space
1894
virtual unsigned int local_dimension() const
1899
// Return the geometric dimension of the coordinates this dof map provides
1900
virtual unsigned int geometric_dimension() const
1905
/// Return the number of dofs on each cell facet
1906
virtual unsigned int num_facet_dofs() const
1911
/// Return the number of dofs associated with each cell entity of dimension d
1912
virtual unsigned int num_entity_dofs(unsigned int d) const
1914
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1917
/// Tabulate the local-to-global mapping of dofs on a cell
1918
virtual void tabulate_dofs(unsigned int* dofs,
1920
const ufc::cell& c) const
1922
dofs[0] = c.entity_indices[0][0];
1923
dofs[1] = c.entity_indices[0][1];
1924
dofs[2] = c.entity_indices[0][2];
1925
unsigned int offset = m.num_entities[0];
1926
dofs[3] = offset + c.entity_indices[1][0];
1927
dofs[4] = offset + c.entity_indices[1][1];
1928
dofs[5] = offset + c.entity_indices[1][2];
1931
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1932
virtual void tabulate_facet_dofs(unsigned int* dofs,
1933
unsigned int facet) const
1955
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1956
virtual void tabulate_entity_dofs(unsigned int* dofs,
1957
unsigned int d, unsigned int i) const
1959
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1962
/// Tabulate the coordinates of all dofs on a cell
1963
virtual void tabulate_coordinates(double** coordinates,
1964
const ufc::cell& c) const
1966
const double * const * x = c.coordinates;
1967
coordinates[0][0] = x[0][0];
1968
coordinates[0][1] = x[0][1];
1969
coordinates[1][0] = x[1][0];
1970
coordinates[1][1] = x[1][1];
1971
coordinates[2][0] = x[2][0];
1972
coordinates[2][1] = x[2][1];
1973
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1974
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1975
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1976
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1977
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1978
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
1981
/// Return the number of sub dof maps (for a mixed element)
1982
virtual unsigned int num_sub_dof_maps() const
1987
/// Create a new dof_map for sub dof map i (for a mixed element)
1988
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1990
return new UFC_Velocity_dof_map_0_1();
1995
/// This class defines the interface for a local-to-global mapping of
1996
/// degrees of freedom (dofs).
1998
class UFC_Velocity_dof_map_0: public ufc::dof_map
2002
unsigned int __global_dimension;
2007
UFC_Velocity_dof_map_0() : ufc::dof_map()
2009
__global_dimension = 0;
2013
virtual ~UFC_Velocity_dof_map_0()
1648
class velocity_0_dof_map_0_0: public ufc::dof_map
1652
unsigned int __global_dimension;
1657
velocity_0_dof_map_0_0() : ufc::dof_map()
1659
__global_dimension = 0;
1663
virtual ~velocity_0_dof_map_0_0()
1668
/// Return a string identifying the dof map
1669
virtual const char* signature() const
1671
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
1674
/// Return true iff mesh entities of topological dimension d are needed
1675
virtual bool needs_mesh_entities(unsigned int d) const
1692
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1693
virtual bool init_mesh(const ufc::mesh& m)
1695
__global_dimension = m.num_entities[0] + m.num_entities[1];
1699
/// Initialize dof map for given cell
1700
virtual void init_cell(const ufc::mesh& m,
1706
/// Finish initialization of dof map for cells
1707
virtual void init_cell_finalize()
1712
/// Return the dimension of the global finite element function space
1713
virtual unsigned int global_dimension() const
1715
return __global_dimension;
1718
/// Return the dimension of the local finite element function space for a cell
1719
virtual unsigned int local_dimension(const ufc::cell& c) const
1724
/// Return the maximum dimension of the local finite element function space
1725
virtual unsigned int max_local_dimension() const
1730
// Return the geometric dimension of the coordinates this dof map provides
1731
virtual unsigned int geometric_dimension() const
1736
/// Return the number of dofs on each cell facet
1737
virtual unsigned int num_facet_dofs() const
1742
/// Return the number of dofs associated with each cell entity of dimension d
1743
virtual unsigned int num_entity_dofs(unsigned int d) const
1745
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1748
/// Tabulate the local-to-global mapping of dofs on a cell
1749
virtual void tabulate_dofs(unsigned int* dofs,
1751
const ufc::cell& c) const
1753
dofs[0] = c.entity_indices[0][0];
1754
dofs[1] = c.entity_indices[0][1];
1755
dofs[2] = c.entity_indices[0][2];
1756
unsigned int offset = m.num_entities[0];
1757
dofs[3] = offset + c.entity_indices[1][0];
1758
dofs[4] = offset + c.entity_indices[1][1];
1759
dofs[5] = offset + c.entity_indices[1][2];
1762
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1763
virtual void tabulate_facet_dofs(unsigned int* dofs,
1764
unsigned int facet) const
1786
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1787
virtual void tabulate_entity_dofs(unsigned int* dofs,
1788
unsigned int d, unsigned int i) const
1790
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1793
/// Tabulate the coordinates of all dofs on a cell
1794
virtual void tabulate_coordinates(double** coordinates,
1795
const ufc::cell& c) const
1797
const double * const * x = c.coordinates;
1798
coordinates[0][0] = x[0][0];
1799
coordinates[0][1] = x[0][1];
1800
coordinates[1][0] = x[1][0];
1801
coordinates[1][1] = x[1][1];
1802
coordinates[2][0] = x[2][0];
1803
coordinates[2][1] = x[2][1];
1804
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1805
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1806
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1807
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1808
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1809
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
1812
/// Return the number of sub dof maps (for a mixed element)
1813
virtual unsigned int num_sub_dof_maps() const
1818
/// Create a new dof_map for sub dof map i (for a mixed element)
1819
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1821
return new velocity_0_dof_map_0_0();
1826
/// This class defines the interface for a local-to-global mapping of
1827
/// degrees of freedom (dofs).
1829
class velocity_0_dof_map_0_1: public ufc::dof_map
1833
unsigned int __global_dimension;
1838
velocity_0_dof_map_0_1() : ufc::dof_map()
1840
__global_dimension = 0;
1844
virtual ~velocity_0_dof_map_0_1()
1849
/// Return a string identifying the dof map
1850
virtual const char* signature() const
1852
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
1855
/// Return true iff mesh entities of topological dimension d are needed
1856
virtual bool needs_mesh_entities(unsigned int d) const
1873
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1874
virtual bool init_mesh(const ufc::mesh& m)
1876
__global_dimension = m.num_entities[0] + m.num_entities[1];
1880
/// Initialize dof map for given cell
1881
virtual void init_cell(const ufc::mesh& m,
1887
/// Finish initialization of dof map for cells
1888
virtual void init_cell_finalize()
1893
/// Return the dimension of the global finite element function space
1894
virtual unsigned int global_dimension() const
1896
return __global_dimension;
1899
/// Return the dimension of the local finite element function space for a cell
1900
virtual unsigned int local_dimension(const ufc::cell& c) const
1905
/// Return the maximum dimension of the local finite element function space
1906
virtual unsigned int max_local_dimension() const
1911
// Return the geometric dimension of the coordinates this dof map provides
1912
virtual unsigned int geometric_dimension() const
1917
/// Return the number of dofs on each cell facet
1918
virtual unsigned int num_facet_dofs() const
1923
/// Return the number of dofs associated with each cell entity of dimension d
1924
virtual unsigned int num_entity_dofs(unsigned int d) const
1926
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1929
/// Tabulate the local-to-global mapping of dofs on a cell
1930
virtual void tabulate_dofs(unsigned int* dofs,
1932
const ufc::cell& c) const
1934
dofs[0] = c.entity_indices[0][0];
1935
dofs[1] = c.entity_indices[0][1];
1936
dofs[2] = c.entity_indices[0][2];
1937
unsigned int offset = m.num_entities[0];
1938
dofs[3] = offset + c.entity_indices[1][0];
1939
dofs[4] = offset + c.entity_indices[1][1];
1940
dofs[5] = offset + c.entity_indices[1][2];
1943
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1944
virtual void tabulate_facet_dofs(unsigned int* dofs,
1945
unsigned int facet) const
1967
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1968
virtual void tabulate_entity_dofs(unsigned int* dofs,
1969
unsigned int d, unsigned int i) const
1971
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1974
/// Tabulate the coordinates of all dofs on a cell
1975
virtual void tabulate_coordinates(double** coordinates,
1976
const ufc::cell& c) const
1978
const double * const * x = c.coordinates;
1979
coordinates[0][0] = x[0][0];
1980
coordinates[0][1] = x[0][1];
1981
coordinates[1][0] = x[1][0];
1982
coordinates[1][1] = x[1][1];
1983
coordinates[2][0] = x[2][0];
1984
coordinates[2][1] = x[2][1];
1985
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1986
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1987
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1988
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1989
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1990
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
1993
/// Return the number of sub dof maps (for a mixed element)
1994
virtual unsigned int num_sub_dof_maps() const
1999
/// Create a new dof_map for sub dof map i (for a mixed element)
2000
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2002
return new velocity_0_dof_map_0_1();
2007
/// This class defines the interface for a local-to-global mapping of
2008
/// degrees of freedom (dofs).
2010
class velocity_0_dof_map_0: public ufc::dof_map
2014
unsigned int __global_dimension;
2019
velocity_0_dof_map_0() : ufc::dof_map()
2021
__global_dimension = 0;
2025
virtual ~velocity_0_dof_map_0()
added
removed
demo/pde/dg/advection-diffusion/cpp/Velocity.h
