1
// This code conforms with the UFC specification version 1.0
2
// and was automatically generated by FFC version 0.6.0.
4
// Warning: This code was generated with the option '-l dolfin'
5
// and contains DOLFIN-specific wrappers that depend on DOLFIN.
7
#ifndef __POISSON2D_3_H
8
#define __POISSON2D_3_H
15
/// This class defines the interface for a finite element.
17
class UFC_Poisson2D_3BilinearForm_finite_element_0: public ufc::finite_element
22
UFC_Poisson2D_3BilinearForm_finite_element_0() : ufc::finite_element()
28
virtual ~UFC_Poisson2D_3BilinearForm_finite_element_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "FiniteElement('Lagrange', 'triangle', 3)";
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);
108
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
110
// Compute psitilde_a
111
const double psitilde_a_0 = 1;
112
const double psitilde_a_1 = x;
113
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
114
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
116
// Compute psitilde_bs
117
const double psitilde_bs_0_0 = 1;
118
const double psitilde_bs_0_1 = 1.5*y + 0.5;
119
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;
120
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
121
const double psitilde_bs_1_0 = 1;
122
const double psitilde_bs_1_1 = 2.5*y + 1.5;
123
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
124
const double psitilde_bs_2_0 = 1;
125
const double psitilde_bs_2_1 = 3.5*y + 2.5;
126
const double psitilde_bs_3_0 = 1;
128
// Compute basisvalues
129
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
130
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
131
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
132
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
133
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
134
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
135
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
136
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
137
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
138
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
140
// Table(s) of coefficients
141
const static double coefficients0[10][10] = \
142
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
143
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
144
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
145
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
146
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
147
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
148
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
149
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
150
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
151
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
153
// Extract relevant coefficients
154
const double coeff0_0 = coefficients0[dof][0];
155
const double coeff0_1 = coefficients0[dof][1];
156
const double coeff0_2 = coefficients0[dof][2];
157
const double coeff0_3 = coefficients0[dof][3];
158
const double coeff0_4 = coefficients0[dof][4];
159
const double coeff0_5 = coefficients0[dof][5];
160
const double coeff0_6 = coefficients0[dof][6];
161
const double coeff0_7 = coefficients0[dof][7];
162
const double coeff0_8 = coefficients0[dof][8];
163
const double coeff0_9 = coefficients0[dof][9];
166
*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;
169
/// Evaluate all basis functions at given point in cell
170
virtual void evaluate_basis_all(double* values,
171
const double* coordinates,
172
const ufc::cell& c) const
174
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
177
/// Evaluate order n derivatives of basis function i at given point in cell
178
virtual void evaluate_basis_derivatives(unsigned int i,
181
const double* coordinates,
182
const ufc::cell& c) const
184
// Extract vertex coordinates
185
const double * const * element_coordinates = c.coordinates;
187
// Compute Jacobian of affine map from reference cell
188
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
189
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
190
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
191
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
193
// Compute determinant of Jacobian
194
const double detJ = J_00*J_11 - J_01*J_10;
196
// Compute inverse of Jacobian
198
// Get coordinates and map to the reference (UFC) element
199
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
200
element_coordinates[0][0]*element_coordinates[2][1] +\
201
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
202
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
203
element_coordinates[1][0]*element_coordinates[0][1] -\
204
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
206
// Map coordinates to the reference square
207
if (std::abs(y - 1.0) < 1e-14)
210
x = 2.0 *x/(1.0 - y) - 1.0;
213
// Compute number of derivatives
214
unsigned int num_derivatives = 1;
216
for (unsigned int j = 0; j < n; j++)
217
num_derivatives *= 2;
220
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
221
unsigned int **combinations = new unsigned int *[num_derivatives];
223
for (unsigned int j = 0; j < num_derivatives; j++)
225
combinations[j] = new unsigned int [n];
226
for (unsigned int k = 0; k < n; k++)
227
combinations[j][k] = 0;
230
// Generate combinations of derivatives
231
for (unsigned int row = 1; row < num_derivatives; row++)
233
for (unsigned int num = 0; num < row; num++)
235
for (unsigned int col = n-1; col+1 > 0; col--)
237
if (combinations[row][col] + 1 > 1)
238
combinations[row][col] = 0;
241
combinations[row][col] += 1;
248
// Compute inverse of Jacobian
249
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
251
// Declare transformation matrix
252
// Declare pointer to two dimensional array and initialise
253
double **transform = new double *[num_derivatives];
255
for (unsigned int j = 0; j < num_derivatives; j++)
257
transform[j] = new double [num_derivatives];
258
for (unsigned int k = 0; k < num_derivatives; k++)
262
// Construct transformation matrix
263
for (unsigned int row = 0; row < num_derivatives; row++)
265
for (unsigned int col = 0; col < num_derivatives; col++)
267
for (unsigned int k = 0; k < n; k++)
268
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
273
for (unsigned int j = 0; j < 1*num_derivatives; j++)
276
// Map degree of freedom to element degree of freedom
277
const unsigned int dof = i;
280
const double scalings_y_0 = 1;
281
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
282
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
283
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
285
// Compute psitilde_a
286
const double psitilde_a_0 = 1;
287
const double psitilde_a_1 = x;
288
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
289
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
291
// Compute psitilde_bs
292
const double psitilde_bs_0_0 = 1;
293
const double psitilde_bs_0_1 = 1.5*y + 0.5;
294
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;
295
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
296
const double psitilde_bs_1_0 = 1;
297
const double psitilde_bs_1_1 = 2.5*y + 1.5;
298
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
299
const double psitilde_bs_2_0 = 1;
300
const double psitilde_bs_2_1 = 3.5*y + 2.5;
301
const double psitilde_bs_3_0 = 1;
303
// Compute basisvalues
304
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
305
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
306
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
307
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
308
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
309
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
310
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
311
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
312
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
313
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
315
// Table(s) of coefficients
316
const static double coefficients0[10][10] = \
317
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
318
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
319
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
320
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
321
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
322
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
323
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
324
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
325
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
326
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
328
// Interesting (new) part
329
// Tables of derivatives of the polynomial base (transpose)
330
const static double dmats0[10][10] = \
331
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
332
{4.89897948556636, 0, 0, 0, 0, 0, 0, 0, 0, 0},
333
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
334
{0, 9.48683298050514, 0, 0, 0, 0, 0, 0, 0, 0},
335
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},
336
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
337
{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},
338
{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},
339
{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},
340
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
342
const static double dmats1[10][10] = \
343
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
344
{2.44948974278318, 0, 0, 0, 0, 0, 0, 0, 0, 0},
345
{4.24264068711929, 0, 0, 0, 0, 0, 0, 0, 0, 0},
346
{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},
347
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},
348
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},
349
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},
350
{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},
351
{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},
352
{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};
354
// Compute reference derivatives
355
// Declare pointer to array of derivatives on FIAT element
356
double *derivatives = new double [num_derivatives];
358
// Declare coefficients
370
// Declare new coefficients
371
double new_coeff0_0 = 0;
372
double new_coeff0_1 = 0;
373
double new_coeff0_2 = 0;
374
double new_coeff0_3 = 0;
375
double new_coeff0_4 = 0;
376
double new_coeff0_5 = 0;
377
double new_coeff0_6 = 0;
378
double new_coeff0_7 = 0;
379
double new_coeff0_8 = 0;
380
double new_coeff0_9 = 0;
382
// Loop possible derivatives
383
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
385
// Get values from coefficients array
386
new_coeff0_0 = coefficients0[dof][0];
387
new_coeff0_1 = coefficients0[dof][1];
388
new_coeff0_2 = coefficients0[dof][2];
389
new_coeff0_3 = coefficients0[dof][3];
390
new_coeff0_4 = coefficients0[dof][4];
391
new_coeff0_5 = coefficients0[dof][5];
392
new_coeff0_6 = coefficients0[dof][6];
393
new_coeff0_7 = coefficients0[dof][7];
394
new_coeff0_8 = coefficients0[dof][8];
395
new_coeff0_9 = coefficients0[dof][9];
397
// Loop derivative order
398
for (unsigned int j = 0; j < n; j++)
400
// Update old coefficients
401
coeff0_0 = new_coeff0_0;
402
coeff0_1 = new_coeff0_1;
403
coeff0_2 = new_coeff0_2;
404
coeff0_3 = new_coeff0_3;
405
coeff0_4 = new_coeff0_4;
406
coeff0_5 = new_coeff0_5;
407
coeff0_6 = new_coeff0_6;
408
coeff0_7 = new_coeff0_7;
409
coeff0_8 = new_coeff0_8;
410
coeff0_9 = new_coeff0_9;
412
if(combinations[deriv_num][j] == 0)
414
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];
415
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];
416
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];
417
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];
418
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];
419
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];
420
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];
421
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];
422
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];
423
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];
425
if(combinations[deriv_num][j] == 1)
427
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];
428
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];
429
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];
430
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];
431
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];
432
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];
433
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];
434
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];
435
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];
436
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];
440
// Compute derivatives on reference element as dot product of coefficients and basisvalues
441
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;
444
// Transform derivatives back to physical element
445
for (unsigned int row = 0; row < num_derivatives; row++)
447
for (unsigned int col = 0; col < num_derivatives; col++)
449
values[row] += transform[row][col]*derivatives[col];
452
// Delete pointer to array of derivatives on FIAT element
453
delete [] derivatives;
455
// Delete pointer to array of combinations of derivatives and transform
456
for (unsigned int row = 0; row < num_derivatives; row++)
458
delete [] combinations[row];
459
delete [] transform[row];
462
delete [] combinations;
466
/// Evaluate order n derivatives of all basis functions at given point in cell
467
virtual void evaluate_basis_derivatives_all(unsigned int n,
469
const double* coordinates,
470
const ufc::cell& c) const
472
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
475
/// Evaluate linear functional for dof i on the function f
476
virtual double evaluate_dof(unsigned int i,
477
const ufc::function& f,
478
const ufc::cell& c) const
480
// The reference points, direction and weights:
481
const static double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666666666667, 0.333333333333333}}, {{0.333333333333333, 0.666666666666667}}, {{0, 0.333333333333333}}, {{0, 0.666666666666667}}, {{0.333333333333333, 0}}, {{0.666666666666667, 0}}, {{0.333333333333333, 0.333333333333333}}};
482
const static double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
483
const static double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
485
const double * const * x = c.coordinates;
487
// Iterate over the points:
488
// Evaluate basis functions for affine mapping
489
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
490
const double w1 = X[i][0][0];
491
const double w2 = X[i][0][1];
493
// Compute affine mapping y = F(X)
495
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
496
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
498
// Evaluate function at physical points
500
f.evaluate(values, y, c);
502
// Map function values using appropriate mapping
503
// Affine map: Do nothing
505
// Note that we do not map the weights (yet).
507
// Take directional components
508
for(int k = 0; k < 1; k++)
509
result += values[k]*D[i][0][k];
510
// Multiply by weights
516
/// Evaluate linear functionals for all dofs on the function f
517
virtual void evaluate_dofs(double* values,
518
const ufc::function& f,
519
const ufc::cell& c) const
521
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
524
/// Interpolate vertex values from dof values
525
virtual void interpolate_vertex_values(double* vertex_values,
526
const double* dof_values,
527
const ufc::cell& c) const
529
// Evaluate at vertices and use affine mapping
530
vertex_values[0] = dof_values[0];
531
vertex_values[1] = dof_values[1];
532
vertex_values[2] = dof_values[2];
535
/// Return the number of sub elements (for a mixed element)
536
virtual unsigned int num_sub_elements() const
541
/// Create a new finite element for sub element i (for a mixed element)
542
virtual ufc::finite_element* create_sub_element(unsigned int i) const
544
return new UFC_Poisson2D_3BilinearForm_finite_element_0();
549
/// This class defines the interface for a finite element.
551
class UFC_Poisson2D_3BilinearForm_finite_element_1: public ufc::finite_element
556
UFC_Poisson2D_3BilinearForm_finite_element_1() : ufc::finite_element()
562
virtual ~UFC_Poisson2D_3BilinearForm_finite_element_1()
567
/// Return a string identifying the finite element
568
virtual const char* signature() const
570
return "FiniteElement('Lagrange', 'triangle', 3)";
573
/// Return the cell shape
574
virtual ufc::shape cell_shape() const
576
return ufc::triangle;
579
/// Return the dimension of the finite element function space
580
virtual unsigned int space_dimension() const
585
/// Return the rank of the value space
586
virtual unsigned int value_rank() const
591
/// Return the dimension of the value space for axis i
592
virtual unsigned int value_dimension(unsigned int i) const
597
/// Evaluate basis function i at given point in cell
598
virtual void evaluate_basis(unsigned int i,
600
const double* coordinates,
601
const ufc::cell& c) const
603
// Extract vertex coordinates
604
const double * const * element_coordinates = c.coordinates;
606
// Compute Jacobian of affine map from reference cell
607
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
608
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
609
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
610
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
612
// Compute determinant of Jacobian
613
const double detJ = J_00*J_11 - J_01*J_10;
615
// Compute inverse of Jacobian
617
// Get coordinates and map to the reference (UFC) element
618
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
619
element_coordinates[0][0]*element_coordinates[2][1] +\
620
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
621
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
622
element_coordinates[1][0]*element_coordinates[0][1] -\
623
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
625
// Map coordinates to the reference square
626
if (std::abs(y - 1.0) < 1e-14)
629
x = 2.0 *x/(1.0 - y) - 1.0;
635
// Map degree of freedom to element degree of freedom
636
const unsigned int dof = i;
639
const double scalings_y_0 = 1;
640
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
641
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
642
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
644
// Compute psitilde_a
645
const double psitilde_a_0 = 1;
646
const double psitilde_a_1 = x;
647
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
648
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
650
// Compute psitilde_bs
651
const double psitilde_bs_0_0 = 1;
652
const double psitilde_bs_0_1 = 1.5*y + 0.5;
653
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;
654
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
655
const double psitilde_bs_1_0 = 1;
656
const double psitilde_bs_1_1 = 2.5*y + 1.5;
657
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
658
const double psitilde_bs_2_0 = 1;
659
const double psitilde_bs_2_1 = 3.5*y + 2.5;
660
const double psitilde_bs_3_0 = 1;
662
// Compute basisvalues
663
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
664
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
665
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
666
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
667
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
668
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
669
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
670
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
671
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
672
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
674
// Table(s) of coefficients
675
const static double coefficients0[10][10] = \
676
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
677
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
678
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
679
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
680
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
681
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
682
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
683
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
684
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
685
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
687
// Extract relevant coefficients
688
const double coeff0_0 = coefficients0[dof][0];
689
const double coeff0_1 = coefficients0[dof][1];
690
const double coeff0_2 = coefficients0[dof][2];
691
const double coeff0_3 = coefficients0[dof][3];
692
const double coeff0_4 = coefficients0[dof][4];
693
const double coeff0_5 = coefficients0[dof][5];
694
const double coeff0_6 = coefficients0[dof][6];
695
const double coeff0_7 = coefficients0[dof][7];
696
const double coeff0_8 = coefficients0[dof][8];
697
const double coeff0_9 = coefficients0[dof][9];
700
*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;
703
/// Evaluate all basis functions at given point in cell
704
virtual void evaluate_basis_all(double* values,
705
const double* coordinates,
706
const ufc::cell& c) const
708
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
711
/// Evaluate order n derivatives of basis function i at given point in cell
712
virtual void evaluate_basis_derivatives(unsigned int i,
715
const double* coordinates,
716
const ufc::cell& c) const
718
// Extract vertex coordinates
719
const double * const * element_coordinates = c.coordinates;
721
// Compute Jacobian of affine map from reference cell
722
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
723
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
724
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
725
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
727
// Compute determinant of Jacobian
728
const double detJ = J_00*J_11 - J_01*J_10;
730
// Compute inverse of Jacobian
732
// Get coordinates and map to the reference (UFC) element
733
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
734
element_coordinates[0][0]*element_coordinates[2][1] +\
735
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
736
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
737
element_coordinates[1][0]*element_coordinates[0][1] -\
738
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
740
// Map coordinates to the reference square
741
if (std::abs(y - 1.0) < 1e-14)
744
x = 2.0 *x/(1.0 - y) - 1.0;
747
// Compute number of derivatives
748
unsigned int num_derivatives = 1;
750
for (unsigned int j = 0; j < n; j++)
751
num_derivatives *= 2;
754
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
755
unsigned int **combinations = new unsigned int *[num_derivatives];
757
for (unsigned int j = 0; j < num_derivatives; j++)
759
combinations[j] = new unsigned int [n];
760
for (unsigned int k = 0; k < n; k++)
761
combinations[j][k] = 0;
764
// Generate combinations of derivatives
765
for (unsigned int row = 1; row < num_derivatives; row++)
767
for (unsigned int num = 0; num < row; num++)
769
for (unsigned int col = n-1; col+1 > 0; col--)
771
if (combinations[row][col] + 1 > 1)
772
combinations[row][col] = 0;
775
combinations[row][col] += 1;
782
// Compute inverse of Jacobian
783
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
785
// Declare transformation matrix
786
// Declare pointer to two dimensional array and initialise
787
double **transform = new double *[num_derivatives];
789
for (unsigned int j = 0; j < num_derivatives; j++)
791
transform[j] = new double [num_derivatives];
792
for (unsigned int k = 0; k < num_derivatives; k++)
796
// Construct transformation matrix
797
for (unsigned int row = 0; row < num_derivatives; row++)
799
for (unsigned int col = 0; col < num_derivatives; col++)
801
for (unsigned int k = 0; k < n; k++)
802
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
807
for (unsigned int j = 0; j < 1*num_derivatives; j++)
810
// Map degree of freedom to element degree of freedom
811
const unsigned int dof = i;
814
const double scalings_y_0 = 1;
815
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
816
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
817
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
819
// Compute psitilde_a
820
const double psitilde_a_0 = 1;
821
const double psitilde_a_1 = x;
822
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
823
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
825
// Compute psitilde_bs
826
const double psitilde_bs_0_0 = 1;
827
const double psitilde_bs_0_1 = 1.5*y + 0.5;
828
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;
829
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
830
const double psitilde_bs_1_0 = 1;
831
const double psitilde_bs_1_1 = 2.5*y + 1.5;
832
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
833
const double psitilde_bs_2_0 = 1;
834
const double psitilde_bs_2_1 = 3.5*y + 2.5;
835
const double psitilde_bs_3_0 = 1;
837
// Compute basisvalues
838
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
839
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
840
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
841
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
842
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
843
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
844
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
845
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
846
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
847
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
849
// Table(s) of coefficients
850
const static double coefficients0[10][10] = \
851
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
852
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
853
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
854
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
855
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
856
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
857
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
858
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
859
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
860
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
862
// Interesting (new) part
863
// Tables of derivatives of the polynomial base (transpose)
864
const static double dmats0[10][10] = \
865
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
866
{4.89897948556636, 0, 0, 0, 0, 0, 0, 0, 0, 0},
867
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
868
{0, 9.48683298050514, 0, 0, 0, 0, 0, 0, 0, 0},
869
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},
870
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
871
{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},
872
{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},
873
{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},
874
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
876
const static double dmats1[10][10] = \
877
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
878
{2.44948974278318, 0, 0, 0, 0, 0, 0, 0, 0, 0},
879
{4.24264068711929, 0, 0, 0, 0, 0, 0, 0, 0, 0},
880
{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},
881
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},
882
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},
883
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},
884
{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},
885
{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},
886
{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};
888
// Compute reference derivatives
889
// Declare pointer to array of derivatives on FIAT element
890
double *derivatives = new double [num_derivatives];
892
// Declare coefficients
904
// Declare new coefficients
905
double new_coeff0_0 = 0;
906
double new_coeff0_1 = 0;
907
double new_coeff0_2 = 0;
908
double new_coeff0_3 = 0;
909
double new_coeff0_4 = 0;
910
double new_coeff0_5 = 0;
911
double new_coeff0_6 = 0;
912
double new_coeff0_7 = 0;
913
double new_coeff0_8 = 0;
914
double new_coeff0_9 = 0;
916
// Loop possible derivatives
917
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
919
// Get values from coefficients array
920
new_coeff0_0 = coefficients0[dof][0];
921
new_coeff0_1 = coefficients0[dof][1];
922
new_coeff0_2 = coefficients0[dof][2];
923
new_coeff0_3 = coefficients0[dof][3];
924
new_coeff0_4 = coefficients0[dof][4];
925
new_coeff0_5 = coefficients0[dof][5];
926
new_coeff0_6 = coefficients0[dof][6];
927
new_coeff0_7 = coefficients0[dof][7];
928
new_coeff0_8 = coefficients0[dof][8];
929
new_coeff0_9 = coefficients0[dof][9];
931
// Loop derivative order
932
for (unsigned int j = 0; j < n; j++)
934
// Update old coefficients
935
coeff0_0 = new_coeff0_0;
936
coeff0_1 = new_coeff0_1;
937
coeff0_2 = new_coeff0_2;
938
coeff0_3 = new_coeff0_3;
939
coeff0_4 = new_coeff0_4;
940
coeff0_5 = new_coeff0_5;
941
coeff0_6 = new_coeff0_6;
942
coeff0_7 = new_coeff0_7;
943
coeff0_8 = new_coeff0_8;
944
coeff0_9 = new_coeff0_9;
946
if(combinations[deriv_num][j] == 0)
948
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];
949
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];
950
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];
951
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];
952
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];
953
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];
954
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];
955
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];
956
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];
957
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];
959
if(combinations[deriv_num][j] == 1)
961
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];
962
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];
963
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];
964
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];
965
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];
966
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];
967
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];
968
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];
969
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];
970
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];
974
// Compute derivatives on reference element as dot product of coefficients and basisvalues
975
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;
978
// Transform derivatives back to physical element
979
for (unsigned int row = 0; row < num_derivatives; row++)
981
for (unsigned int col = 0; col < num_derivatives; col++)
983
values[row] += transform[row][col]*derivatives[col];
986
// Delete pointer to array of derivatives on FIAT element
987
delete [] derivatives;
989
// Delete pointer to array of combinations of derivatives and transform
990
for (unsigned int row = 0; row < num_derivatives; row++)
992
delete [] combinations[row];
993
delete [] transform[row];
996
delete [] combinations;
1000
/// Evaluate order n derivatives of all basis functions at given point in cell
1001
virtual void evaluate_basis_derivatives_all(unsigned int n,
1003
const double* coordinates,
1004
const ufc::cell& c) const
1006
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1009
/// Evaluate linear functional for dof i on the function f
1010
virtual double evaluate_dof(unsigned int i,
1011
const ufc::function& f,
1012
const ufc::cell& c) const
1014
// The reference points, direction and weights:
1015
const static double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666666666667, 0.333333333333333}}, {{0.333333333333333, 0.666666666666667}}, {{0, 0.333333333333333}}, {{0, 0.666666666666667}}, {{0.333333333333333, 0}}, {{0.666666666666667, 0}}, {{0.333333333333333, 0.333333333333333}}};
1016
const static double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
1017
const static double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
1019
const double * const * x = c.coordinates;
1020
double result = 0.0;
1021
// Iterate over the points:
1022
// Evaluate basis functions for affine mapping
1023
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1024
const double w1 = X[i][0][0];
1025
const double w2 = X[i][0][1];
1027
// Compute affine mapping y = F(X)
1029
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1030
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1032
// Evaluate function at physical points
1034
f.evaluate(values, y, c);
1036
// Map function values using appropriate mapping
1037
// Affine map: Do nothing
1039
// Note that we do not map the weights (yet).
1041
// Take directional components
1042
for(int k = 0; k < 1; k++)
1043
result += values[k]*D[i][0][k];
1044
// Multiply by weights
1050
/// Evaluate linear functionals for all dofs on the function f
1051
virtual void evaluate_dofs(double* values,
1052
const ufc::function& f,
1053
const ufc::cell& c) const
1055
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1058
/// Interpolate vertex values from dof values
1059
virtual void interpolate_vertex_values(double* vertex_values,
1060
const double* dof_values,
1061
const ufc::cell& c) const
1063
// Evaluate at vertices and use affine mapping
1064
vertex_values[0] = dof_values[0];
1065
vertex_values[1] = dof_values[1];
1066
vertex_values[2] = dof_values[2];
1069
/// Return the number of sub elements (for a mixed element)
1070
virtual unsigned int num_sub_elements() const
1075
/// Create a new finite element for sub element i (for a mixed element)
1076
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1078
return new UFC_Poisson2D_3BilinearForm_finite_element_1();
1083
/// This class defines the interface for a local-to-global mapping of
1084
/// degrees of freedom (dofs).
1086
class UFC_Poisson2D_3BilinearForm_dof_map_0: public ufc::dof_map
1090
unsigned int __global_dimension;
1095
UFC_Poisson2D_3BilinearForm_dof_map_0() : ufc::dof_map()
1097
__global_dimension = 0;
1101
virtual ~UFC_Poisson2D_3BilinearForm_dof_map_0()
1106
/// Return a string identifying the dof map
1107
virtual const char* signature() const
1109
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";
1112
/// Return true iff mesh entities of topological dimension d are needed
1113
virtual bool needs_mesh_entities(unsigned int d) const
1130
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1131
virtual bool init_mesh(const ufc::mesh& m)
1133
__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];
1137
/// Initialize dof map for given cell
1138
virtual void init_cell(const ufc::mesh& m,
1144
/// Finish initialization of dof map for cells
1145
virtual void init_cell_finalize()
1150
/// Return the dimension of the global finite element function space
1151
virtual unsigned int global_dimension() const
1153
return __global_dimension;
1156
/// Return the dimension of the local finite element function space
1157
virtual unsigned int local_dimension() const
1162
// Return the geometric dimension of the coordinates this dof map provides
1163
virtual unsigned int geometric_dimension() const
1168
/// Return the number of dofs on each cell facet
1169
virtual unsigned int num_facet_dofs() const
1174
/// Return the number of dofs associated with each cell entity of dimension d
1175
virtual unsigned int num_entity_dofs(unsigned int d) const
1177
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1180
/// Tabulate the local-to-global mapping of dofs on a cell
1181
virtual void tabulate_dofs(unsigned int* dofs,
1183
const ufc::cell& c) const
1185
dofs[0] = c.entity_indices[0][0];
1186
dofs[1] = c.entity_indices[0][1];
1187
dofs[2] = c.entity_indices[0][2];
1188
unsigned int offset = m.num_entities[0];
1189
dofs[3] = offset + 2*c.entity_indices[1][0];
1190
dofs[4] = offset + 2*c.entity_indices[1][0] + 1;
1191
dofs[5] = offset + 2*c.entity_indices[1][1];
1192
dofs[6] = offset + 2*c.entity_indices[1][1] + 1;
1193
dofs[7] = offset + 2*c.entity_indices[1][2];
1194
dofs[8] = offset + 2*c.entity_indices[1][2] + 1;
1195
offset = offset + 2*m.num_entities[1];
1196
dofs[9] = offset + c.entity_indices[2][0];
1199
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1200
virtual void tabulate_facet_dofs(unsigned int* dofs,
1201
unsigned int facet) const
1226
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1227
virtual void tabulate_entity_dofs(unsigned int* dofs,
1228
unsigned int d, unsigned int i) const
1230
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1233
/// Tabulate the coordinates of all dofs on a cell
1234
virtual void tabulate_coordinates(double** coordinates,
1235
const ufc::cell& c) const
1237
const double * const * x = c.coordinates;
1238
coordinates[0][0] = x[0][0];
1239
coordinates[0][1] = x[0][1];
1240
coordinates[1][0] = x[1][0];
1241
coordinates[1][1] = x[1][1];
1242
coordinates[2][0] = x[2][0];
1243
coordinates[2][1] = x[2][1];
1244
coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];
1245
coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];
1246
coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];
1247
coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];
1248
coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];
1249
coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];
1250
coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];
1251
coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];
1252
coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];
1253
coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];
1254
coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];
1255
coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];
1256
coordinates[9][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
1257
coordinates[9][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
1260
/// Return the number of sub dof maps (for a mixed element)
1261
virtual unsigned int num_sub_dof_maps() const
1266
/// Create a new dof_map for sub dof map i (for a mixed element)
1267
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1269
return new UFC_Poisson2D_3BilinearForm_dof_map_0();
1274
/// This class defines the interface for a local-to-global mapping of
1275
/// degrees of freedom (dofs).
1277
class UFC_Poisson2D_3BilinearForm_dof_map_1: public ufc::dof_map
1281
unsigned int __global_dimension;
1286
UFC_Poisson2D_3BilinearForm_dof_map_1() : ufc::dof_map()
1288
__global_dimension = 0;
1292
virtual ~UFC_Poisson2D_3BilinearForm_dof_map_1()
1297
/// Return a string identifying the dof map
1298
virtual const char* signature() const
1300
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";
1303
/// Return true iff mesh entities of topological dimension d are needed
1304
virtual bool needs_mesh_entities(unsigned int d) const
1321
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1322
virtual bool init_mesh(const ufc::mesh& m)
1324
__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];
1328
/// Initialize dof map for given cell
1329
virtual void init_cell(const ufc::mesh& m,
1335
/// Finish initialization of dof map for cells
1336
virtual void init_cell_finalize()
1341
/// Return the dimension of the global finite element function space
1342
virtual unsigned int global_dimension() const
1344
return __global_dimension;
1347
/// Return the dimension of the local finite element function space
1348
virtual unsigned int local_dimension() const
1353
// Return the geometric dimension of the coordinates this dof map provides
1354
virtual unsigned int geometric_dimension() const
1359
/// Return the number of dofs on each cell facet
1360
virtual unsigned int num_facet_dofs() const
1365
/// Return the number of dofs associated with each cell entity of dimension d
1366
virtual unsigned int num_entity_dofs(unsigned int d) const
1368
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1371
/// Tabulate the local-to-global mapping of dofs on a cell
1372
virtual void tabulate_dofs(unsigned int* dofs,
1374
const ufc::cell& c) const
1376
dofs[0] = c.entity_indices[0][0];
1377
dofs[1] = c.entity_indices[0][1];
1378
dofs[2] = c.entity_indices[0][2];
1379
unsigned int offset = m.num_entities[0];
1380
dofs[3] = offset + 2*c.entity_indices[1][0];
1381
dofs[4] = offset + 2*c.entity_indices[1][0] + 1;
1382
dofs[5] = offset + 2*c.entity_indices[1][1];
1383
dofs[6] = offset + 2*c.entity_indices[1][1] + 1;
1384
dofs[7] = offset + 2*c.entity_indices[1][2];
1385
dofs[8] = offset + 2*c.entity_indices[1][2] + 1;
1386
offset = offset + 2*m.num_entities[1];
1387
dofs[9] = offset + c.entity_indices[2][0];
1390
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1391
virtual void tabulate_facet_dofs(unsigned int* dofs,
1392
unsigned int facet) const
1417
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1418
virtual void tabulate_entity_dofs(unsigned int* dofs,
1419
unsigned int d, unsigned int i) const
1421
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1424
/// Tabulate the coordinates of all dofs on a cell
1425
virtual void tabulate_coordinates(double** coordinates,
1426
const ufc::cell& c) const
1428
const double * const * x = c.coordinates;
1429
coordinates[0][0] = x[0][0];
1430
coordinates[0][1] = x[0][1];
1431
coordinates[1][0] = x[1][0];
1432
coordinates[1][1] = x[1][1];
1433
coordinates[2][0] = x[2][0];
1434
coordinates[2][1] = x[2][1];
1435
coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];
1436
coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];
1437
coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];
1438
coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];
1439
coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];
1440
coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];
1441
coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];
1442
coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];
1443
coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];
1444
coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];
1445
coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];
1446
coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];
1447
coordinates[9][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
1448
coordinates[9][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
1451
/// Return the number of sub dof maps (for a mixed element)
1452
virtual unsigned int num_sub_dof_maps() const
1457
/// Create a new dof_map for sub dof map i (for a mixed element)
1458
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1460
return new UFC_Poisson2D_3BilinearForm_dof_map_1();
1465
/// This class defines the interface for the tabulation of the cell
1466
/// tensor corresponding to the local contribution to a form from
1467
/// the integral over a cell.
1469
class UFC_Poisson2D_3BilinearForm_cell_integral_0: public ufc::cell_integral
1474
UFC_Poisson2D_3BilinearForm_cell_integral_0() : ufc::cell_integral()
1480
virtual ~UFC_Poisson2D_3BilinearForm_cell_integral_0()
1485
/// Tabulate the tensor for the contribution from a local cell
1486
virtual void tabulate_tensor(double* A,
1487
const double * const * w,
1488
const ufc::cell& c) const
1490
// Extract vertex coordinates
1491
const double * const * x = c.coordinates;
1493
// Compute Jacobian of affine map from reference cell
1494
const double J_00 = x[1][0] - x[0][0];
1495
const double J_01 = x[2][0] - x[0][0];
1496
const double J_10 = x[1][1] - x[0][1];
1497
const double J_11 = x[2][1] - x[0][1];
1499
// Compute determinant of Jacobian
1500
double detJ = J_00*J_11 - J_01*J_10;
1502
// Compute inverse of Jacobian
1503
const double Jinv_00 = J_11 / detJ;
1504
const double Jinv_01 = -J_01 / detJ;
1505
const double Jinv_10 = -J_10 / detJ;
1506
const double Jinv_11 = J_00 / detJ;
1509
const double det = std::abs(detJ);
1511
// Number of operations to compute element tensor = 426
1512
// Compute geometry tensors
1513
// Number of operations to compute decalrations = 16
1514
const double G0_0_0 = det*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
1515
const double G0_0_1 = det*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
1516
const double G0_1_0 = det*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
1517
const double G0_1_1 = det*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
1519
// Compute element tensor
1520
// Number of operations to compute tensor = 410
1521
A[0] = 0.425*G0_0_0 + 0.425*G0_0_1 + 0.425*G0_1_0 + 0.425*G0_1_1;
1522
A[1] = -0.0875*G0_0_0 - 0.0875*G0_1_0;
1523
A[2] = -0.0874999999999998*G0_0_1 - 0.0874999999999998*G0_1_1;
1524
A[3] = -0.0374999999999998*G0_0_0 - 0.0374999999999999*G0_0_1 - 0.0374999999999999*G0_1_0 - 0.0375000000000002*G0_1_1;
1525
A[4] = -0.0374999999999997*G0_0_0 - 0.0374999999999996*G0_0_1 - 0.0374999999999997*G0_1_0 - 0.0374999999999995*G0_1_1;
1526
A[5] = 0.0374999999999999*G0_0_0 - 0.674999999999999*G0_0_1 + 0.0374999999999997*G0_1_0 - 0.674999999999999*G0_1_1;
1527
A[6] = 0.0374999999999997*G0_0_0 + 0.3375*G0_0_1 + 0.0374999999999997*G0_1_0 + 0.3375*G0_1_1;
1528
A[7] = -0.674999999999999*G0_0_0 + 0.0374999999999999*G0_0_1 - 0.674999999999999*G0_1_0 + 0.0374999999999998*G0_1_1;
1529
A[8] = 0.3375*G0_0_0 + 0.0375000000000001*G0_0_1 + 0.3375*G0_1_0 + 0.0375000000000003*G0_1_1;
1531
A[10] = -0.0875*G0_0_0 - 0.0875*G0_0_1;
1532
A[11] = 0.425*G0_0_0;
1533
A[12] = 0.0874999999999999*G0_0_1;
1534
A[13] = 0.0375000000000002*G0_0_0 + 0.712499999999999*G0_0_1;
1535
A[14] = 0.0374999999999999*G0_0_0 - 0.3*G0_0_1;
1536
A[15] = -0.0374999999999997*G0_0_0;
1537
A[16] = -0.0375*G0_0_0;
1538
A[17] = 0.3375*G0_0_0 + 0.3*G0_0_1;
1539
A[18] = -0.674999999999999*G0_0_0 - 0.7125*G0_0_1;
1541
A[20] = -0.0874999999999997*G0_1_0 - 0.0874999999999998*G0_1_1;
1542
A[21] = 0.0874999999999999*G0_1_0;
1543
A[22] = 0.424999999999999*G0_1_1;
1544
A[23] = -0.3*G0_1_0 + 0.0375000000000001*G0_1_1;
1545
A[24] = 0.712499999999999*G0_1_0 + 0.0375000000000001*G0_1_1;
1546
A[25] = 0.299999999999999*G0_1_0 + 0.337499999999999*G0_1_1;
1547
A[26] = -0.712499999999999*G0_1_0 - 0.674999999999999*G0_1_1;
1548
A[27] = -0.0375*G0_1_1;
1549
A[28] = -0.0375*G0_1_1;
1551
A[30] = -0.0374999999999998*G0_0_0 - 0.0374999999999999*G0_0_1 - 0.0375*G0_1_0 - 0.0375000000000002*G0_1_1;
1552
A[31] = 0.0375000000000002*G0_0_0 + 0.712499999999999*G0_1_0;
1553
A[32] = -0.3*G0_0_1 + 0.0375000000000001*G0_1_1;
1554
A[33] = 1.6875*G0_0_0 + 0.84375*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;
1555
A[34] = -0.3375*G0_0_0 + 0.843749999999998*G0_0_1 - 0.16875*G0_1_0 - 0.3375*G0_1_1;
1556
A[35] = 0.337500000000001*G0_0_0 + 0.16875*G0_0_1 + 0.168750000000001*G0_1_0;
1557
A[36] = 0.337499999999999*G0_0_0 + 0.168749999999999*G0_0_1 + 0.16875*G0_1_0;
1558
A[37] = 0.168750000000001*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;
1559
A[38] = -0.84375*G0_0_1 - 0.843749999999999*G0_1_0 - 1.6875*G0_1_1;
1560
A[39] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;
1561
A[40] = -0.0374999999999997*G0_0_0 - 0.0374999999999997*G0_0_1 - 0.0374999999999996*G0_1_0 - 0.0374999999999995*G0_1_1;
1562
A[41] = 0.0374999999999999*G0_0_0 - 0.3*G0_1_0;
1563
A[42] = 0.712499999999999*G0_0_1 + 0.0375000000000002*G0_1_1;
1564
A[43] = -0.3375*G0_0_0 - 0.16875*G0_0_1 + 0.843749999999998*G0_1_0 - 0.3375*G0_1_1;
1565
A[44] = 1.6875*G0_0_0 + 0.84375*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;
1566
A[45] = 0.337499999999999*G0_0_0 + 0.168749999999999*G0_0_1 + 0.16875*G0_1_0;
1567
A[46] = -1.6875*G0_0_0 - 0.843749999999998*G0_0_1 - 0.84375*G0_1_0;
1568
A[47] = 0.16875*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;
1569
A[48] = 0.16875*G0_0_1 + 0.168749999999999*G0_1_0 + 0.3375*G0_1_1;
1570
A[49] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;
1571
A[50] = 0.0374999999999998*G0_0_0 + 0.0374999999999997*G0_0_1 - 0.674999999999999*G0_1_0 - 0.674999999999999*G0_1_1;
1572
A[51] = -0.0374999999999997*G0_0_0;
1573
A[52] = 0.299999999999999*G0_0_1 + 0.337499999999999*G0_1_1;
1574
A[53] = 0.337500000000001*G0_0_0 + 0.168750000000001*G0_0_1 + 0.16875*G0_1_0;
1575
A[54] = 0.337499999999999*G0_0_0 + 0.16875*G0_0_1 + 0.168749999999999*G0_1_0;
1576
A[55] = 1.6875*G0_0_0 + 0.843749999999999*G0_0_1 + 0.843749999999999*G0_1_0 + 1.6875*G0_1_1;
1577
A[56] = -0.337499999999999*G0_0_0 - 1.18125*G0_0_1 - 0.168749999999999*G0_1_0 - 1.35*G0_1_1;
1578
A[57] = 0.84375*G0_0_1 + 0.843749999999999*G0_1_0;
1579
A[58] = -0.168750000000001*G0_0_1 - 0.16875*G0_1_0;
1580
A[59] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;
1581
A[60] = 0.0374999999999997*G0_0_0 + 0.0374999999999997*G0_0_1 + 0.3375*G0_1_0 + 0.3375*G0_1_1;
1582
A[61] = -0.0375*G0_0_0;
1583
A[62] = -0.712499999999999*G0_0_1 - 0.674999999999999*G0_1_1;
1584
A[63] = 0.337499999999999*G0_0_0 + 0.16875*G0_0_1 + 0.168749999999999*G0_1_0;
1585
A[64] = -1.6875*G0_0_0 - 0.84375*G0_0_1 - 0.843749999999998*G0_1_0;
1586
A[65] = -0.337499999999999*G0_0_0 - 0.168749999999999*G0_0_1 - 1.18125*G0_1_0 - 1.35*G0_1_1;
1587
A[66] = 1.6875*G0_0_0 + 0.843749999999998*G0_0_1 + 0.843749999999998*G0_1_0 + 1.6875*G0_1_1;
1588
A[67] = -0.16875*G0_0_1 - 0.16875*G0_1_0;
1589
A[68] = -0.168749999999999*G0_0_1 - 0.168749999999999*G0_1_0;
1590
A[69] = 1.0125*G0_0_1 + 1.0125*G0_1_0;
1591
A[70] = -0.674999999999999*G0_0_0 - 0.674999999999999*G0_0_1 + 0.0374999999999999*G0_1_0 + 0.0374999999999998*G0_1_1;
1592
A[71] = 0.3375*G0_0_0 + 0.3*G0_1_0;
1593
A[72] = -0.0375*G0_1_1;
1594
A[73] = 0.16875*G0_0_1 + 0.168750000000001*G0_1_0 + 0.3375*G0_1_1;
1595
A[74] = 0.16875*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;
1596
A[75] = 0.843749999999999*G0_0_1 + 0.84375*G0_1_0;
1597
A[76] = -0.16875*G0_0_1 - 0.16875*G0_1_0;
1598
A[77] = 1.6875*G0_0_0 + 0.843749999999999*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;
1599
A[78] = -1.35*G0_0_0 - 0.16875*G0_0_1 - 1.18125*G0_1_0 - 0.3375*G0_1_1;
1600
A[79] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;
1601
A[80] = 0.3375*G0_0_0 + 0.3375*G0_0_1 + 0.0375000000000001*G0_1_0 + 0.0375000000000003*G0_1_1;
1602
A[81] = -0.674999999999999*G0_0_0 - 0.712499999999999*G0_1_0;
1603
A[82] = -0.0375*G0_1_1;
1604
A[83] = -0.843749999999999*G0_0_1 - 0.84375*G0_1_0 - 1.6875*G0_1_1;
1605
A[84] = 0.168749999999999*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;
1606
A[85] = -0.16875*G0_0_1 - 0.168750000000001*G0_1_0;
1607
A[86] = -0.168749999999999*G0_0_1 - 0.168749999999999*G0_1_0;
1608
A[87] = -1.35*G0_0_0 - 1.18125*G0_0_1 - 0.16875*G0_1_0 - 0.3375*G0_1_1;
1609
A[88] = 1.6875*G0_0_0 + 0.84375*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;
1610
A[89] = 1.0125*G0_0_1 + 1.0125*G0_1_0;
1614
A[93] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;
1615
A[94] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;
1616
A[95] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;
1617
A[96] = 1.0125*G0_0_1 + 1.0125*G0_1_0;
1618
A[97] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;
1619
A[98] = 1.0125*G0_0_1 + 1.0125*G0_1_0;
1620
A[99] = 4.05*G0_0_0 + 2.025*G0_0_1 + 2.025*G0_1_0 + 4.05*G0_1_1;
1625
/// This class defines the interface for the assembly of the global
1626
/// tensor corresponding to a form with r + n arguments, that is, a
1629
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
1631
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
1632
/// global tensor A is defined by
1634
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
1636
/// where each argument Vj represents the application to the
1637
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
1638
/// fixed functions (coefficients).
1640
class UFC_Poisson2D_3BilinearForm: public ufc::form
1645
UFC_Poisson2D_3BilinearForm() : ufc::form()
1651
virtual ~UFC_Poisson2D_3BilinearForm()
1656
/// Return a string identifying the form
1657
virtual const char* signature() const
1659
return "(dXa0[0, 1]/dxb0[0, 1])(dXa1[0, 1]/dxb0[0, 1]) | ((d/dXa0[0, 1])vi0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])*((d/dXa1[0, 1])vi1[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])*dX(0)";
1662
/// Return the rank of the global tensor (r)
1663
virtual unsigned int rank() const
1668
/// Return the number of coefficients (n)
1669
virtual unsigned int num_coefficients() const
1674
/// Return the number of cell integrals
1675
virtual unsigned int num_cell_integrals() const
1680
/// Return the number of exterior facet integrals
1681
virtual unsigned int num_exterior_facet_integrals() const
1686
/// Return the number of interior facet integrals
1687
virtual unsigned int num_interior_facet_integrals() const
1692
/// Create a new finite element for argument function i
1693
virtual ufc::finite_element* create_finite_element(unsigned int i) const
1698
return new UFC_Poisson2D_3BilinearForm_finite_element_0();
1701
return new UFC_Poisson2D_3BilinearForm_finite_element_1();
1707
/// Create a new dof map for argument function i
1708
virtual ufc::dof_map* create_dof_map(unsigned int i) const
1713
return new UFC_Poisson2D_3BilinearForm_dof_map_0();
1716
return new UFC_Poisson2D_3BilinearForm_dof_map_1();
1722
/// Create a new cell integral on sub domain i
1723
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
1725
return new UFC_Poisson2D_3BilinearForm_cell_integral_0();
1728
/// Create a new exterior facet integral on sub domain i
1729
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
1734
/// Create a new interior facet integral on sub domain i
1735
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
1742
/// This class defines the interface for a finite element.
1744
class UFC_Poisson2D_3LinearForm_finite_element_0: public ufc::finite_element
1749
UFC_Poisson2D_3LinearForm_finite_element_0() : ufc::finite_element()
1755
virtual ~UFC_Poisson2D_3LinearForm_finite_element_0()
1760
/// Return a string identifying the finite element
1761
virtual const char* signature() const
1763
return "FiniteElement('Lagrange', 'triangle', 3)";
1766
/// Return the cell shape
1767
virtual ufc::shape cell_shape() const
1769
return ufc::triangle;
1772
/// Return the dimension of the finite element function space
1773
virtual unsigned int space_dimension() const
1778
/// Return the rank of the value space
1779
virtual unsigned int value_rank() const
1784
/// Return the dimension of the value space for axis i
1785
virtual unsigned int value_dimension(unsigned int i) const
1790
/// Evaluate basis function i at given point in cell
1791
virtual void evaluate_basis(unsigned int i,
1793
const double* coordinates,
1794
const ufc::cell& c) const
1796
// Extract vertex coordinates
1797
const double * const * element_coordinates = c.coordinates;
1799
// Compute Jacobian of affine map from reference cell
1800
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1801
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1802
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1803
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1805
// Compute determinant of Jacobian
1806
const double detJ = J_00*J_11 - J_01*J_10;
1808
// Compute inverse of Jacobian
1810
// Get coordinates and map to the reference (UFC) element
1811
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1812
element_coordinates[0][0]*element_coordinates[2][1] +\
1813
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1814
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1815
element_coordinates[1][0]*element_coordinates[0][1] -\
1816
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1818
// Map coordinates to the reference square
1819
if (std::abs(y - 1.0) < 1e-14)
1822
x = 2.0 *x/(1.0 - y) - 1.0;
1828
// Map degree of freedom to element degree of freedom
1829
const unsigned int dof = i;
1831
// Generate scalings
1832
const double scalings_y_0 = 1;
1833
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1834
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
1835
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
1837
// Compute psitilde_a
1838
const double psitilde_a_0 = 1;
1839
const double psitilde_a_1 = x;
1840
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
1841
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
1843
// Compute psitilde_bs
1844
const double psitilde_bs_0_0 = 1;
1845
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1846
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;
1847
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
1848
const double psitilde_bs_1_0 = 1;
1849
const double psitilde_bs_1_1 = 2.5*y + 1.5;
1850
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
1851
const double psitilde_bs_2_0 = 1;
1852
const double psitilde_bs_2_1 = 3.5*y + 2.5;
1853
const double psitilde_bs_3_0 = 1;
1855
// Compute basisvalues
1856
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1857
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1858
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1859
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
1860
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
1861
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
1862
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
1863
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
1864
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
1865
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
1867
// Table(s) of coefficients
1868
const static double coefficients0[10][10] = \
1869
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
1870
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
1871
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
1872
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
1873
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
1874
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
1875
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
1876
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
1877
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
1878
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
1880
// Extract relevant coefficients
1881
const double coeff0_0 = coefficients0[dof][0];
1882
const double coeff0_1 = coefficients0[dof][1];
1883
const double coeff0_2 = coefficients0[dof][2];
1884
const double coeff0_3 = coefficients0[dof][3];
1885
const double coeff0_4 = coefficients0[dof][4];
1886
const double coeff0_5 = coefficients0[dof][5];
1887
const double coeff0_6 = coefficients0[dof][6];
1888
const double coeff0_7 = coefficients0[dof][7];
1889
const double coeff0_8 = coefficients0[dof][8];
1890
const double coeff0_9 = coefficients0[dof][9];
1893
*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;
1896
/// Evaluate all basis functions at given point in cell
1897
virtual void evaluate_basis_all(double* values,
1898
const double* coordinates,
1899
const ufc::cell& c) const
1901
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1904
/// Evaluate order n derivatives of basis function i at given point in cell
1905
virtual void evaluate_basis_derivatives(unsigned int i,
1908
const double* coordinates,
1909
const ufc::cell& c) const
1911
// Extract vertex coordinates
1912
const double * const * element_coordinates = c.coordinates;
1914
// Compute Jacobian of affine map from reference cell
1915
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1916
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1917
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1918
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1920
// Compute determinant of Jacobian
1921
const double detJ = J_00*J_11 - J_01*J_10;
1923
// Compute inverse of Jacobian
1925
// Get coordinates and map to the reference (UFC) element
1926
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1927
element_coordinates[0][0]*element_coordinates[2][1] +\
1928
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1929
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1930
element_coordinates[1][0]*element_coordinates[0][1] -\
1931
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1933
// Map coordinates to the reference square
1934
if (std::abs(y - 1.0) < 1e-14)
1937
x = 2.0 *x/(1.0 - y) - 1.0;
1940
// Compute number of derivatives
1941
unsigned int num_derivatives = 1;
1943
for (unsigned int j = 0; j < n; j++)
1944
num_derivatives *= 2;
1947
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1948
unsigned int **combinations = new unsigned int *[num_derivatives];
1950
for (unsigned int j = 0; j < num_derivatives; j++)
1952
combinations[j] = new unsigned int [n];
1953
for (unsigned int k = 0; k < n; k++)
1954
combinations[j][k] = 0;
1957
// Generate combinations of derivatives
1958
for (unsigned int row = 1; row < num_derivatives; row++)
1960
for (unsigned int num = 0; num < row; num++)
1962
for (unsigned int col = n-1; col+1 > 0; col--)
1964
if (combinations[row][col] + 1 > 1)
1965
combinations[row][col] = 0;
1968
combinations[row][col] += 1;
1975
// Compute inverse of Jacobian
1976
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1978
// Declare transformation matrix
1979
// Declare pointer to two dimensional array and initialise
1980
double **transform = new double *[num_derivatives];
1982
for (unsigned int j = 0; j < num_derivatives; j++)
1984
transform[j] = new double [num_derivatives];
1985
for (unsigned int k = 0; k < num_derivatives; k++)
1986
transform[j][k] = 1;
1989
// Construct transformation matrix
1990
for (unsigned int row = 0; row < num_derivatives; row++)
1992
for (unsigned int col = 0; col < num_derivatives; col++)
1994
for (unsigned int k = 0; k < n; k++)
1995
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2000
for (unsigned int j = 0; j < 1*num_derivatives; j++)
2003
// Map degree of freedom to element degree of freedom
2004
const unsigned int dof = i;
2006
// Generate scalings
2007
const double scalings_y_0 = 1;
2008
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
2009
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
2010
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
2012
// Compute psitilde_a
2013
const double psitilde_a_0 = 1;
2014
const double psitilde_a_1 = x;
2015
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
2016
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
2018
// Compute psitilde_bs
2019
const double psitilde_bs_0_0 = 1;
2020
const double psitilde_bs_0_1 = 1.5*y + 0.5;
2021
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;
2022
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
2023
const double psitilde_bs_1_0 = 1;
2024
const double psitilde_bs_1_1 = 2.5*y + 1.5;
2025
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
2026
const double psitilde_bs_2_0 = 1;
2027
const double psitilde_bs_2_1 = 3.5*y + 2.5;
2028
const double psitilde_bs_3_0 = 1;
2030
// Compute basisvalues
2031
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2032
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
2033
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
2034
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
2035
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
2036
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
2037
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
2038
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
2039
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
2040
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
2042
// Table(s) of coefficients
2043
const static double coefficients0[10][10] = \
2044
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
2045
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
2046
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
2047
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
2048
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
2049
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
2050
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
2051
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
2052
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
2053
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
2055
// Interesting (new) part
2056
// Tables of derivatives of the polynomial base (transpose)
2057
const static double dmats0[10][10] = \
2058
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2059
{4.89897948556636, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2060
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2061
{0, 9.48683298050514, 0, 0, 0, 0, 0, 0, 0, 0},
2062
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},
2063
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2064
{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},
2065
{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},
2066
{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},
2067
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
2069
const static double dmats1[10][10] = \
2070
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2071
{2.44948974278318, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2072
{4.24264068711929, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2073
{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},
2074
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},
2075
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},
2076
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},
2077
{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},
2078
{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},
2079
{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};
2081
// Compute reference derivatives
2082
// Declare pointer to array of derivatives on FIAT element
2083
double *derivatives = new double [num_derivatives];
2085
// Declare coefficients
2086
double coeff0_0 = 0;
2087
double coeff0_1 = 0;
2088
double coeff0_2 = 0;
2089
double coeff0_3 = 0;
2090
double coeff0_4 = 0;
2091
double coeff0_5 = 0;
2092
double coeff0_6 = 0;
2093
double coeff0_7 = 0;
2094
double coeff0_8 = 0;
2095
double coeff0_9 = 0;
2097
// Declare new coefficients
2098
double new_coeff0_0 = 0;
2099
double new_coeff0_1 = 0;
2100
double new_coeff0_2 = 0;
2101
double new_coeff0_3 = 0;
2102
double new_coeff0_4 = 0;
2103
double new_coeff0_5 = 0;
2104
double new_coeff0_6 = 0;
2105
double new_coeff0_7 = 0;
2106
double new_coeff0_8 = 0;
2107
double new_coeff0_9 = 0;
2109
// Loop possible derivatives
2110
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
2112
// Get values from coefficients array
2113
new_coeff0_0 = coefficients0[dof][0];
2114
new_coeff0_1 = coefficients0[dof][1];
2115
new_coeff0_2 = coefficients0[dof][2];
2116
new_coeff0_3 = coefficients0[dof][3];
2117
new_coeff0_4 = coefficients0[dof][4];
2118
new_coeff0_5 = coefficients0[dof][5];
2119
new_coeff0_6 = coefficients0[dof][6];
2120
new_coeff0_7 = coefficients0[dof][7];
2121
new_coeff0_8 = coefficients0[dof][8];
2122
new_coeff0_9 = coefficients0[dof][9];
2124
// Loop derivative order
2125
for (unsigned int j = 0; j < n; j++)
2127
// Update old coefficients
2128
coeff0_0 = new_coeff0_0;
2129
coeff0_1 = new_coeff0_1;
2130
coeff0_2 = new_coeff0_2;
2131
coeff0_3 = new_coeff0_3;
2132
coeff0_4 = new_coeff0_4;
2133
coeff0_5 = new_coeff0_5;
2134
coeff0_6 = new_coeff0_6;
2135
coeff0_7 = new_coeff0_7;
2136
coeff0_8 = new_coeff0_8;
2137
coeff0_9 = new_coeff0_9;
2139
if(combinations[deriv_num][j] == 0)
2141
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];
2142
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];
2143
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];
2144
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];
2145
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];
2146
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];
2147
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];
2148
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];
2149
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];
2150
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];
2152
if(combinations[deriv_num][j] == 1)
2154
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];
2155
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];
2156
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];
2157
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];
2158
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];
2159
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];
2160
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];
2161
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];
2162
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];
2163
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];
2167
// Compute derivatives on reference element as dot product of coefficients and basisvalues
2168
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;
2171
// Transform derivatives back to physical element
2172
for (unsigned int row = 0; row < num_derivatives; row++)
2174
for (unsigned int col = 0; col < num_derivatives; col++)
2176
values[row] += transform[row][col]*derivatives[col];
2179
// Delete pointer to array of derivatives on FIAT element
2180
delete [] derivatives;
2182
// Delete pointer to array of combinations of derivatives and transform
2183
for (unsigned int row = 0; row < num_derivatives; row++)
2185
delete [] combinations[row];
2186
delete [] transform[row];
2189
delete [] combinations;
2190
delete [] transform;
2193
/// Evaluate order n derivatives of all basis functions at given point in cell
2194
virtual void evaluate_basis_derivatives_all(unsigned int n,
2196
const double* coordinates,
2197
const ufc::cell& c) const
2199
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
2202
/// Evaluate linear functional for dof i on the function f
2203
virtual double evaluate_dof(unsigned int i,
2204
const ufc::function& f,
2205
const ufc::cell& c) const
2207
// The reference points, direction and weights:
2208
const static double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666666666667, 0.333333333333333}}, {{0.333333333333333, 0.666666666666667}}, {{0, 0.333333333333333}}, {{0, 0.666666666666667}}, {{0.333333333333333, 0}}, {{0.666666666666667, 0}}, {{0.333333333333333, 0.333333333333333}}};
2209
const static double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
2210
const static double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
2212
const double * const * x = c.coordinates;
2213
double result = 0.0;
2214
// Iterate over the points:
2215
// Evaluate basis functions for affine mapping
2216
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
2217
const double w1 = X[i][0][0];
2218
const double w2 = X[i][0][1];
2220
// Compute affine mapping y = F(X)
2222
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
2223
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
2225
// Evaluate function at physical points
2227
f.evaluate(values, y, c);
2229
// Map function values using appropriate mapping
2230
// Affine map: Do nothing
2232
// Note that we do not map the weights (yet).
2234
// Take directional components
2235
for(int k = 0; k < 1; k++)
2236
result += values[k]*D[i][0][k];
2237
// Multiply by weights
2243
/// Evaluate linear functionals for all dofs on the function f
2244
virtual void evaluate_dofs(double* values,
2245
const ufc::function& f,
2246
const ufc::cell& c) const
2248
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2251
/// Interpolate vertex values from dof values
2252
virtual void interpolate_vertex_values(double* vertex_values,
2253
const double* dof_values,
2254
const ufc::cell& c) const
2256
// Evaluate at vertices and use affine mapping
2257
vertex_values[0] = dof_values[0];
2258
vertex_values[1] = dof_values[1];
2259
vertex_values[2] = dof_values[2];
2262
/// Return the number of sub elements (for a mixed element)
2263
virtual unsigned int num_sub_elements() const
2268
/// Create a new finite element for sub element i (for a mixed element)
2269
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2271
return new UFC_Poisson2D_3LinearForm_finite_element_0();
2276
/// This class defines the interface for a finite element.
2278
class UFC_Poisson2D_3LinearForm_finite_element_1: public ufc::finite_element
2283
UFC_Poisson2D_3LinearForm_finite_element_1() : ufc::finite_element()
2289
virtual ~UFC_Poisson2D_3LinearForm_finite_element_1()
2294
/// Return a string identifying the finite element
2295
virtual const char* signature() const
2297
return "FiniteElement('Lagrange', 'triangle', 3)";
2300
/// Return the cell shape
2301
virtual ufc::shape cell_shape() const
2303
return ufc::triangle;
2306
/// Return the dimension of the finite element function space
2307
virtual unsigned int space_dimension() const
2312
/// Return the rank of the value space
2313
virtual unsigned int value_rank() const
2318
/// Return the dimension of the value space for axis i
2319
virtual unsigned int value_dimension(unsigned int i) const
2324
/// Evaluate basis function i at given point in cell
2325
virtual void evaluate_basis(unsigned int i,
2327
const double* coordinates,
2328
const ufc::cell& c) const
2330
// Extract vertex coordinates
2331
const double * const * element_coordinates = c.coordinates;
2333
// Compute Jacobian of affine map from reference cell
2334
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
2335
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
2336
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
2337
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
2339
// Compute determinant of Jacobian
2340
const double detJ = J_00*J_11 - J_01*J_10;
2342
// Compute inverse of Jacobian
2344
// Get coordinates and map to the reference (UFC) element
2345
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
2346
element_coordinates[0][0]*element_coordinates[2][1] +\
2347
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
2348
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
2349
element_coordinates[1][0]*element_coordinates[0][1] -\
2350
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
2352
// Map coordinates to the reference square
2353
if (std::abs(y - 1.0) < 1e-14)
2356
x = 2.0 *x/(1.0 - y) - 1.0;
2362
// Map degree of freedom to element degree of freedom
2363
const unsigned int dof = i;
2365
// Generate scalings
2366
const double scalings_y_0 = 1;
2367
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
2368
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
2369
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
2371
// Compute psitilde_a
2372
const double psitilde_a_0 = 1;
2373
const double psitilde_a_1 = x;
2374
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
2375
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
2377
// Compute psitilde_bs
2378
const double psitilde_bs_0_0 = 1;
2379
const double psitilde_bs_0_1 = 1.5*y + 0.5;
2380
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;
2381
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
2382
const double psitilde_bs_1_0 = 1;
2383
const double psitilde_bs_1_1 = 2.5*y + 1.5;
2384
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
2385
const double psitilde_bs_2_0 = 1;
2386
const double psitilde_bs_2_1 = 3.5*y + 2.5;
2387
const double psitilde_bs_3_0 = 1;
2389
// Compute basisvalues
2390
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2391
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
2392
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
2393
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
2394
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
2395
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
2396
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
2397
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
2398
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
2399
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
2401
// Table(s) of coefficients
2402
const static double coefficients0[10][10] = \
2403
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
2404
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
2405
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
2406
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
2407
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
2408
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
2409
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
2410
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
2411
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
2412
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
2414
// Extract relevant coefficients
2415
const double coeff0_0 = coefficients0[dof][0];
2416
const double coeff0_1 = coefficients0[dof][1];
2417
const double coeff0_2 = coefficients0[dof][2];
2418
const double coeff0_3 = coefficients0[dof][3];
2419
const double coeff0_4 = coefficients0[dof][4];
2420
const double coeff0_5 = coefficients0[dof][5];
2421
const double coeff0_6 = coefficients0[dof][6];
2422
const double coeff0_7 = coefficients0[dof][7];
2423
const double coeff0_8 = coefficients0[dof][8];
2424
const double coeff0_9 = coefficients0[dof][9];
2427
*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;
2430
/// Evaluate all basis functions at given point in cell
2431
virtual void evaluate_basis_all(double* values,
2432
const double* coordinates,
2433
const ufc::cell& c) const
2435
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
2438
/// Evaluate order n derivatives of basis function i at given point in cell
2439
virtual void evaluate_basis_derivatives(unsigned int i,
2442
const double* coordinates,
2443
const ufc::cell& c) const
2445
// Extract vertex coordinates
2446
const double * const * element_coordinates = c.coordinates;
2448
// Compute Jacobian of affine map from reference cell
2449
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
2450
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
2451
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
2452
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
2454
// Compute determinant of Jacobian
2455
const double detJ = J_00*J_11 - J_01*J_10;
2457
// Compute inverse of Jacobian
2459
// Get coordinates and map to the reference (UFC) element
2460
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
2461
element_coordinates[0][0]*element_coordinates[2][1] +\
2462
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
2463
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
2464
element_coordinates[1][0]*element_coordinates[0][1] -\
2465
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
2467
// Map coordinates to the reference square
2468
if (std::abs(y - 1.0) < 1e-14)
2471
x = 2.0 *x/(1.0 - y) - 1.0;
2474
// Compute number of derivatives
2475
unsigned int num_derivatives = 1;
2477
for (unsigned int j = 0; j < n; j++)
2478
num_derivatives *= 2;
2481
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
2482
unsigned int **combinations = new unsigned int *[num_derivatives];
2484
for (unsigned int j = 0; j < num_derivatives; j++)
2486
combinations[j] = new unsigned int [n];
2487
for (unsigned int k = 0; k < n; k++)
2488
combinations[j][k] = 0;
2491
// Generate combinations of derivatives
2492
for (unsigned int row = 1; row < num_derivatives; row++)
2494
for (unsigned int num = 0; num < row; num++)
2496
for (unsigned int col = n-1; col+1 > 0; col--)
2498
if (combinations[row][col] + 1 > 1)
2499
combinations[row][col] = 0;
2502
combinations[row][col] += 1;
2509
// Compute inverse of Jacobian
2510
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
2512
// Declare transformation matrix
2513
// Declare pointer to two dimensional array and initialise
2514
double **transform = new double *[num_derivatives];
2516
for (unsigned int j = 0; j < num_derivatives; j++)
2518
transform[j] = new double [num_derivatives];
2519
for (unsigned int k = 0; k < num_derivatives; k++)
2520
transform[j][k] = 1;
2523
// Construct transformation matrix
2524
for (unsigned int row = 0; row < num_derivatives; row++)
2526
for (unsigned int col = 0; col < num_derivatives; col++)
2528
for (unsigned int k = 0; k < n; k++)
2529
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2534
for (unsigned int j = 0; j < 1*num_derivatives; j++)
2537
// Map degree of freedom to element degree of freedom
2538
const unsigned int dof = i;
2540
// Generate scalings
2541
const double scalings_y_0 = 1;
2542
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
2543
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
2544
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
2546
// Compute psitilde_a
2547
const double psitilde_a_0 = 1;
2548
const double psitilde_a_1 = x;
2549
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
2550
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
2552
// Compute psitilde_bs
2553
const double psitilde_bs_0_0 = 1;
2554
const double psitilde_bs_0_1 = 1.5*y + 0.5;
2555
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;
2556
const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
2557
const double psitilde_bs_1_0 = 1;
2558
const double psitilde_bs_1_1 = 2.5*y + 1.5;
2559
const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
2560
const double psitilde_bs_2_0 = 1;
2561
const double psitilde_bs_2_1 = 3.5*y + 2.5;
2562
const double psitilde_bs_3_0 = 1;
2564
// Compute basisvalues
2565
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2566
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
2567
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
2568
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
2569
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
2570
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
2571
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
2572
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
2573
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
2574
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
2576
// Table(s) of coefficients
2577
const static double coefficients0[10][10] = \
2578
{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
2579
{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
2580
{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
2581
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
2582
{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
2583
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
2584
{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
2585
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},
2586
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},
2587
{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};
2589
// Interesting (new) part
2590
// Tables of derivatives of the polynomial base (transpose)
2591
const static double dmats0[10][10] = \
2592
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2593
{4.89897948556636, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2594
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2595
{0, 9.48683298050514, 0, 0, 0, 0, 0, 0, 0, 0},
2596
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},
2597
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2598
{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},
2599
{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},
2600
{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},
2601
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
2603
const static double dmats1[10][10] = \
2604
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2605
{2.44948974278318, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2606
{4.24264068711929, 0, 0, 0, 0, 0, 0, 0, 0, 0},
2607
{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},
2608
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},
2609
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},
2610
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},
2611
{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},
2612
{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},
2613
{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};
2615
// Compute reference derivatives
2616
// Declare pointer to array of derivatives on FIAT element
2617
double *derivatives = new double [num_derivatives];
2619
// Declare coefficients
2620
double coeff0_0 = 0;
2621
double coeff0_1 = 0;
2622
double coeff0_2 = 0;
2623
double coeff0_3 = 0;
2624
double coeff0_4 = 0;
2625
double coeff0_5 = 0;
2626
double coeff0_6 = 0;
2627
double coeff0_7 = 0;
2628
double coeff0_8 = 0;
2629
double coeff0_9 = 0;
2631
// Declare new coefficients
2632
double new_coeff0_0 = 0;
2633
double new_coeff0_1 = 0;
2634
double new_coeff0_2 = 0;
2635
double new_coeff0_3 = 0;
2636
double new_coeff0_4 = 0;
2637
double new_coeff0_5 = 0;
2638
double new_coeff0_6 = 0;
2639
double new_coeff0_7 = 0;
2640
double new_coeff0_8 = 0;
2641
double new_coeff0_9 = 0;
2643
// Loop possible derivatives
2644
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
2646
// Get values from coefficients array
2647
new_coeff0_0 = coefficients0[dof][0];
2648
new_coeff0_1 = coefficients0[dof][1];
2649
new_coeff0_2 = coefficients0[dof][2];
2650
new_coeff0_3 = coefficients0[dof][3];
2651
new_coeff0_4 = coefficients0[dof][4];
2652
new_coeff0_5 = coefficients0[dof][5];
2653
new_coeff0_6 = coefficients0[dof][6];
2654
new_coeff0_7 = coefficients0[dof][7];
2655
new_coeff0_8 = coefficients0[dof][8];
2656
new_coeff0_9 = coefficients0[dof][9];
2658
// Loop derivative order
2659
for (unsigned int j = 0; j < n; j++)
2661
// Update old coefficients
2662
coeff0_0 = new_coeff0_0;
2663
coeff0_1 = new_coeff0_1;
2664
coeff0_2 = new_coeff0_2;
2665
coeff0_3 = new_coeff0_3;
2666
coeff0_4 = new_coeff0_4;
2667
coeff0_5 = new_coeff0_5;
2668
coeff0_6 = new_coeff0_6;
2669
coeff0_7 = new_coeff0_7;
2670
coeff0_8 = new_coeff0_8;
2671
coeff0_9 = new_coeff0_9;
2673
if(combinations[deriv_num][j] == 0)
2675
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];
2676
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];
2677
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];
2678
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];
2679
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];
2680
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];
2681
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];
2682
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];
2683
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];
2684
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];
2686
if(combinations[deriv_num][j] == 1)
2688
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];
2689
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];
2690
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];
2691
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];
2692
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];
2693
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];
2694
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];
2695
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];
2696
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];
2697
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];
2701
// Compute derivatives on reference element as dot product of coefficients and basisvalues
2702
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;
2705
// Transform derivatives back to physical element
2706
for (unsigned int row = 0; row < num_derivatives; row++)
2708
for (unsigned int col = 0; col < num_derivatives; col++)
2710
values[row] += transform[row][col]*derivatives[col];
2713
// Delete pointer to array of derivatives on FIAT element
2714
delete [] derivatives;
2716
// Delete pointer to array of combinations of derivatives and transform
2717
for (unsigned int row = 0; row < num_derivatives; row++)
2719
delete [] combinations[row];
2720
delete [] transform[row];
2723
delete [] combinations;
2724
delete [] transform;
2727
/// Evaluate order n derivatives of all basis functions at given point in cell
2728
virtual void evaluate_basis_derivatives_all(unsigned int n,
2730
const double* coordinates,
2731
const ufc::cell& c) const
2733
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
2736
/// Evaluate linear functional for dof i on the function f
2737
virtual double evaluate_dof(unsigned int i,
2738
const ufc::function& f,
2739
const ufc::cell& c) const
2741
// The reference points, direction and weights:
2742
const static double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666666666667, 0.333333333333333}}, {{0.333333333333333, 0.666666666666667}}, {{0, 0.333333333333333}}, {{0, 0.666666666666667}}, {{0.333333333333333, 0}}, {{0.666666666666667, 0}}, {{0.333333333333333, 0.333333333333333}}};
2743
const static double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
2744
const static double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
2746
const double * const * x = c.coordinates;
2747
double result = 0.0;
2748
// Iterate over the points:
2749
// Evaluate basis functions for affine mapping
2750
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
2751
const double w1 = X[i][0][0];
2752
const double w2 = X[i][0][1];
2754
// Compute affine mapping y = F(X)
2756
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
2757
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
2759
// Evaluate function at physical points
2761
f.evaluate(values, y, c);
2763
// Map function values using appropriate mapping
2764
// Affine map: Do nothing
2766
// Note that we do not map the weights (yet).
2768
// Take directional components
2769
for(int k = 0; k < 1; k++)
2770
result += values[k]*D[i][0][k];
2771
// Multiply by weights
2777
/// Evaluate linear functionals for all dofs on the function f
2778
virtual void evaluate_dofs(double* values,
2779
const ufc::function& f,
2780
const ufc::cell& c) const
2782
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2785
/// Interpolate vertex values from dof values
2786
virtual void interpolate_vertex_values(double* vertex_values,
2787
const double* dof_values,
2788
const ufc::cell& c) const
2790
// Evaluate at vertices and use affine mapping
2791
vertex_values[0] = dof_values[0];
2792
vertex_values[1] = dof_values[1];
2793
vertex_values[2] = dof_values[2];
2796
/// Return the number of sub elements (for a mixed element)
2797
virtual unsigned int num_sub_elements() const
2802
/// Create a new finite element for sub element i (for a mixed element)
2803
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2805
return new UFC_Poisson2D_3LinearForm_finite_element_1();
2810
/// This class defines the interface for a local-to-global mapping of
2811
/// degrees of freedom (dofs).
2813
class UFC_Poisson2D_3LinearForm_dof_map_0: public ufc::dof_map
2817
unsigned int __global_dimension;
2822
UFC_Poisson2D_3LinearForm_dof_map_0() : ufc::dof_map()
2824
__global_dimension = 0;
2828
virtual ~UFC_Poisson2D_3LinearForm_dof_map_0()
2833
/// Return a string identifying the dof map
2834
virtual const char* signature() const
2836
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";
2839
/// Return true iff mesh entities of topological dimension d are needed
2840
virtual bool needs_mesh_entities(unsigned int d) const
2857
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2858
virtual bool init_mesh(const ufc::mesh& m)
2860
__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];
2864
/// Initialize dof map for given cell
2865
virtual void init_cell(const ufc::mesh& m,
2871
/// Finish initialization of dof map for cells
2872
virtual void init_cell_finalize()
2877
/// Return the dimension of the global finite element function space
2878
virtual unsigned int global_dimension() const
2880
return __global_dimension;
2883
/// Return the dimension of the local finite element function space
2884
virtual unsigned int local_dimension() const
2889
// Return the geometric dimension of the coordinates this dof map provides
2890
virtual unsigned int geometric_dimension() const
2895
/// Return the number of dofs on each cell facet
2896
virtual unsigned int num_facet_dofs() const
2901
/// Return the number of dofs associated with each cell entity of dimension d
2902
virtual unsigned int num_entity_dofs(unsigned int d) const
2904
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2907
/// Tabulate the local-to-global mapping of dofs on a cell
2908
virtual void tabulate_dofs(unsigned int* dofs,
2910
const ufc::cell& c) const
2912
dofs[0] = c.entity_indices[0][0];
2913
dofs[1] = c.entity_indices[0][1];
2914
dofs[2] = c.entity_indices[0][2];
2915
unsigned int offset = m.num_entities[0];
2916
dofs[3] = offset + 2*c.entity_indices[1][0];
2917
dofs[4] = offset + 2*c.entity_indices[1][0] + 1;
2918
dofs[5] = offset + 2*c.entity_indices[1][1];
2919
dofs[6] = offset + 2*c.entity_indices[1][1] + 1;
2920
dofs[7] = offset + 2*c.entity_indices[1][2];
2921
dofs[8] = offset + 2*c.entity_indices[1][2] + 1;
2922
offset = offset + 2*m.num_entities[1];
2923
dofs[9] = offset + c.entity_indices[2][0];
2926
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2927
virtual void tabulate_facet_dofs(unsigned int* dofs,
2928
unsigned int facet) const
2953
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2954
virtual void tabulate_entity_dofs(unsigned int* dofs,
2955
unsigned int d, unsigned int i) const
2957
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2960
/// Tabulate the coordinates of all dofs on a cell
2961
virtual void tabulate_coordinates(double** coordinates,
2962
const ufc::cell& c) const
2964
const double * const * x = c.coordinates;
2965
coordinates[0][0] = x[0][0];
2966
coordinates[0][1] = x[0][1];
2967
coordinates[1][0] = x[1][0];
2968
coordinates[1][1] = x[1][1];
2969
coordinates[2][0] = x[2][0];
2970
coordinates[2][1] = x[2][1];
2971
coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];
2972
coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];
2973
coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];
2974
coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];
2975
coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];
2976
coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];
2977
coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];
2978
coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];
2979
coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];
2980
coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];
2981
coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];
2982
coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];
2983
coordinates[9][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
2984
coordinates[9][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
2987
/// Return the number of sub dof maps (for a mixed element)
2988
virtual unsigned int num_sub_dof_maps() const
2993
/// Create a new dof_map for sub dof map i (for a mixed element)
2994
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2996
return new UFC_Poisson2D_3LinearForm_dof_map_0();
3001
/// This class defines the interface for a local-to-global mapping of
3002
/// degrees of freedom (dofs).
3004
class UFC_Poisson2D_3LinearForm_dof_map_1: public ufc::dof_map
3008
unsigned int __global_dimension;
3013
UFC_Poisson2D_3LinearForm_dof_map_1() : ufc::dof_map()
3015
__global_dimension = 0;
3019
virtual ~UFC_Poisson2D_3LinearForm_dof_map_1()
3024
/// Return a string identifying the dof map
3025
virtual const char* signature() const
3027
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";
3030
/// Return true iff mesh entities of topological dimension d are needed
3031
virtual bool needs_mesh_entities(unsigned int d) const
3048
/// Initialize dof map for mesh (return true iff init_cell() is needed)
3049
virtual bool init_mesh(const ufc::mesh& m)
3051
__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];
3055
/// Initialize dof map for given cell
3056
virtual void init_cell(const ufc::mesh& m,
3062
/// Finish initialization of dof map for cells
3063
virtual void init_cell_finalize()
3068
/// Return the dimension of the global finite element function space
3069
virtual unsigned int global_dimension() const
3071
return __global_dimension;
3074
/// Return the dimension of the local finite element function space
3075
virtual unsigned int local_dimension() const
3080
// Return the geometric dimension of the coordinates this dof map provides
3081
virtual unsigned int geometric_dimension() const
3086
/// Return the number of dofs on each cell facet
3087
virtual unsigned int num_facet_dofs() const
3092
/// Return the number of dofs associated with each cell entity of dimension d
3093
virtual unsigned int num_entity_dofs(unsigned int d) const
3095
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3098
/// Tabulate the local-to-global mapping of dofs on a cell
3099
virtual void tabulate_dofs(unsigned int* dofs,
3101
const ufc::cell& c) const
3103
dofs[0] = c.entity_indices[0][0];
3104
dofs[1] = c.entity_indices[0][1];
3105
dofs[2] = c.entity_indices[0][2];
3106
unsigned int offset = m.num_entities[0];
3107
dofs[3] = offset + 2*c.entity_indices[1][0];
3108
dofs[4] = offset + 2*c.entity_indices[1][0] + 1;
3109
dofs[5] = offset + 2*c.entity_indices[1][1];
3110
dofs[6] = offset + 2*c.entity_indices[1][1] + 1;
3111
dofs[7] = offset + 2*c.entity_indices[1][2];
3112
dofs[8] = offset + 2*c.entity_indices[1][2] + 1;
3113
offset = offset + 2*m.num_entities[1];
3114
dofs[9] = offset + c.entity_indices[2][0];
3117
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3118
virtual void tabulate_facet_dofs(unsigned int* dofs,
3119
unsigned int facet) const
3144
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3145
virtual void tabulate_entity_dofs(unsigned int* dofs,
3146
unsigned int d, unsigned int i) const
3148
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3151
/// Tabulate the coordinates of all dofs on a cell
3152
virtual void tabulate_coordinates(double** coordinates,
3153
const ufc::cell& c) const
3155
const double * const * x = c.coordinates;
3156
coordinates[0][0] = x[0][0];
3157
coordinates[0][1] = x[0][1];
3158
coordinates[1][0] = x[1][0];
3159
coordinates[1][1] = x[1][1];
3160
coordinates[2][0] = x[2][0];
3161
coordinates[2][1] = x[2][1];
3162
coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];
3163
coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];
3164
coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];
3165
coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];
3166
coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];
3167
coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];
3168
coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];
3169
coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];
3170
coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];
3171
coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];
3172
coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];
3173
coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];
3174
coordinates[9][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
3175
coordinates[9][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
3178
/// Return the number of sub dof maps (for a mixed element)
3179
virtual unsigned int num_sub_dof_maps() const
3184
/// Create a new dof_map for sub dof map i (for a mixed element)
3185
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
3187
return new UFC_Poisson2D_3LinearForm_dof_map_1();
3192
/// This class defines the interface for the tabulation of the cell
3193
/// tensor corresponding to the local contribution to a form from
3194
/// the integral over a cell.
3196
class UFC_Poisson2D_3LinearForm_cell_integral_0: public ufc::cell_integral
3201
UFC_Poisson2D_3LinearForm_cell_integral_0() : ufc::cell_integral()
3207
virtual ~UFC_Poisson2D_3LinearForm_cell_integral_0()
3212
/// Tabulate the tensor for the contribution from a local cell
3213
virtual void tabulate_tensor(double* A,
3214
const double * const * w,
3215
const ufc::cell& c) const
3217
// Extract vertex coordinates
3218
const double * const * x = c.coordinates;
3220
// Compute Jacobian of affine map from reference cell
3221
const double J_00 = x[1][0] - x[0][0];
3222
const double J_01 = x[2][0] - x[0][0];
3223
const double J_10 = x[1][1] - x[0][1];
3224
const double J_11 = x[2][1] - x[0][1];
3226
// Compute determinant of Jacobian
3227
double detJ = J_00*J_11 - J_01*J_10;
3229
// Compute inverse of Jacobian
3232
const double det = std::abs(detJ);
3234
// Number of operations to compute element tensor = 176
3235
// Compute coefficients
3236
const double c0_0_0_0 = w[0][0];
3237
const double c0_0_0_1 = w[0][1];
3238
const double c0_0_0_2 = w[0][2];
3239
const double c0_0_0_3 = w[0][3];
3240
const double c0_0_0_4 = w[0][4];
3241
const double c0_0_0_5 = w[0][5];
3242
const double c0_0_0_6 = w[0][6];
3243
const double c0_0_0_7 = w[0][7];
3244
const double c0_0_0_8 = w[0][8];
3245
const double c0_0_0_9 = w[0][9];
3247
// Compute geometry tensors
3248
// Number of operations to compute decalrations = 10
3249
const double G0_0 = det*c0_0_0_0;
3250
const double G0_1 = det*c0_0_0_1;
3251
const double G0_2 = det*c0_0_0_2;
3252
const double G0_3 = det*c0_0_0_3;
3253
const double G0_4 = det*c0_0_0_4;
3254
const double G0_5 = det*c0_0_0_5;
3255
const double G0_6 = det*c0_0_0_6;
3256
const double G0_7 = det*c0_0_0_7;
3257
const double G0_8 = det*c0_0_0_8;
3258
const double G0_9 = det*c0_0_0_9;
3260
// Compute element tensor
3261
// Number of operations to compute tensor = 166
3262
A[0] = 0.0056547619047619*G0_0 + 0.000818452380952379*G0_1 + 0.000818452380952379*G0_2 + 0.00200892857142857*G0_3 + 0.00200892857142857*G0_4 + 0.0013392857142857*G0_5 + 0.0013392857142857*G0_7 + 0.00267857142857142*G0_9;
3263
A[1] = 0.000818452380952379*G0_0 + 0.0056547619047619*G0_1 + 0.00081845238095238*G0_2 + 0.0013392857142857*G0_3 + 0.00200892857142857*G0_5 + 0.00200892857142858*G0_6 + 0.0013392857142857*G0_8 + 0.00267857142857143*G0_9;
3264
A[2] = 0.000818452380952379*G0_0 + 0.00081845238095238*G0_1 + 0.0056547619047619*G0_2 + 0.00133928571428572*G0_4 + 0.00133928571428571*G0_6 + 0.00200892857142857*G0_7 + 0.00200892857142857*G0_8 + 0.00267857142857142*G0_9;
3265
A[3] = 0.00200892857142857*G0_0 + 0.0013392857142857*G0_1 + 0.0401785714285714*G0_3 - 0.0140625*G0_4 - 0.00401785714285713*G0_5 - 0.0100446428571429*G0_6 - 0.0100446428571428*G0_7 + 0.0200892857142857*G0_8 + 0.0120535714285714*G0_9;
3266
A[4] = 0.00200892857142857*G0_0 + 0.00133928571428572*G0_2 - 0.0140625*G0_3 + 0.0401785714285714*G0_4 - 0.0100446428571429*G0_5 + 0.0200892857142857*G0_6 - 0.00401785714285714*G0_7 - 0.0100446428571428*G0_8 + 0.0120535714285714*G0_9;
3267
A[5] = 0.0013392857142857*G0_0 + 0.00200892857142857*G0_1 - 0.00401785714285713*G0_3 - 0.0100446428571429*G0_4 + 0.0401785714285714*G0_5 - 0.0140625*G0_6 + 0.0200892857142857*G0_7 - 0.0100446428571428*G0_8 + 0.0120535714285714*G0_9;
3268
A[6] = 0.00200892857142858*G0_1 + 0.00133928571428571*G0_2 - 0.0100446428571429*G0_3 + 0.0200892857142857*G0_4 - 0.0140625*G0_5 + 0.0401785714285714*G0_6 - 0.0100446428571429*G0_7 - 0.00401785714285714*G0_8 + 0.0120535714285714*G0_9;
3269
A[7] = 0.0013392857142857*G0_0 + 0.00200892857142857*G0_2 - 0.0100446428571428*G0_3 - 0.00401785714285714*G0_4 + 0.0200892857142857*G0_5 - 0.0100446428571429*G0_6 + 0.0401785714285714*G0_7 - 0.0140625*G0_8 + 0.0120535714285714*G0_9;
3270
A[8] = 0.0013392857142857*G0_1 + 0.00200892857142857*G0_2 + 0.0200892857142857*G0_3 - 0.0100446428571428*G0_4 - 0.0100446428571428*G0_5 - 0.00401785714285714*G0_6 - 0.0140625*G0_7 + 0.0401785714285714*G0_8 + 0.0120535714285714*G0_9;
3271
A[9] = 0.00267857142857142*G0_0 + 0.00267857142857143*G0_1 + 0.00267857142857142*G0_2 + 0.0120535714285714*G0_3 + 0.0120535714285714*G0_4 + 0.0120535714285714*G0_5 + 0.0120535714285714*G0_6 + 0.0120535714285714*G0_7 + 0.0120535714285714*G0_8 + 0.144642857142857*G0_9;
3276
/// This class defines the interface for the assembly of the global
3277
/// tensor corresponding to a form with r + n arguments, that is, a
3280
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3282
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3283
/// global tensor A is defined by
3285
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3287
/// where each argument Vj represents the application to the
3288
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3289
/// fixed functions (coefficients).
3291
class UFC_Poisson2D_3LinearForm: public ufc::form
3296
UFC_Poisson2D_3LinearForm() : ufc::form()
3302
virtual ~UFC_Poisson2D_3LinearForm()
3307
/// Return a string identifying the form
3308
virtual const char* signature() const
3310
return "w0_a0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9] | vi0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]*va0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]*dX(0)";
3313
/// Return the rank of the global tensor (r)
3314
virtual unsigned int rank() const
3319
/// Return the number of coefficients (n)
3320
virtual unsigned int num_coefficients() const
3325
/// Return the number of cell integrals
3326
virtual unsigned int num_cell_integrals() const
3331
/// Return the number of exterior facet integrals
3332
virtual unsigned int num_exterior_facet_integrals() const
3337
/// Return the number of interior facet integrals
3338
virtual unsigned int num_interior_facet_integrals() const
3343
/// Create a new finite element for argument function i
3344
virtual ufc::finite_element* create_finite_element(unsigned int i) const
3349
return new UFC_Poisson2D_3LinearForm_finite_element_0();
3352
return new UFC_Poisson2D_3LinearForm_finite_element_1();
3358
/// Create a new dof map for argument function i
3359
virtual ufc::dof_map* create_dof_map(unsigned int i) const
3364
return new UFC_Poisson2D_3LinearForm_dof_map_0();
3367
return new UFC_Poisson2D_3LinearForm_dof_map_1();
3373
/// Create a new cell integral on sub domain i
3374
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
3376
return new UFC_Poisson2D_3LinearForm_cell_integral_0();
3379
/// Create a new exterior facet integral on sub domain i
3380
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
3385
/// Create a new interior facet integral on sub domain i
3386
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
3395
#include <dolfin/fem/Form.h>
3396
#include <dolfin/fem/FiniteElement.h>
3397
#include <dolfin/fem/DofMap.h>
3398
#include <dolfin/function/Coefficient.h>
3399
#include <dolfin/function/Function.h>
3400
#include <dolfin/function/FunctionSpace.h>
3402
class Poisson2D_3BilinearFormFunctionSpace0 : public dolfin::FunctionSpace
3406
Poisson2D_3BilinearFormFunctionSpace0(const dolfin::Mesh& mesh)
3407
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3408
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3409
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3416
class Poisson2D_3BilinearFormFunctionSpace1 : public dolfin::FunctionSpace
3420
Poisson2D_3BilinearFormFunctionSpace1(const dolfin::Mesh& mesh)
3421
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3422
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3423
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3430
class Poisson2D_3LinearFormFunctionSpace0 : public dolfin::FunctionSpace
3434
Poisson2D_3LinearFormFunctionSpace0(const dolfin::Mesh& mesh)
3435
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3436
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3437
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3444
class Poisson2D_3LinearFormCoefficientSpace0 : public dolfin::FunctionSpace
3448
Poisson2D_3LinearFormCoefficientSpace0(const dolfin::Mesh& mesh)
3449
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3450
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3451
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3458
class Poisson2D_3TestSpace : public dolfin::FunctionSpace
3462
Poisson2D_3TestSpace(const dolfin::Mesh& mesh)
3463
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3464
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3465
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3472
class Poisson2D_3TrialSpace : public dolfin::FunctionSpace
3476
Poisson2D_3TrialSpace(const dolfin::Mesh& mesh)
3477
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3478
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3479
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3486
class Poisson2D_3CoefficientSpace : public dolfin::FunctionSpace
3490
Poisson2D_3CoefficientSpace(const dolfin::Mesh& mesh)
3491
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3492
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3493
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3500
class Poisson2D_3FunctionSpace : public dolfin::FunctionSpace
3504
Poisson2D_3FunctionSpace(const dolfin::Mesh& mesh)
3505
: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
3506
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_3LinearForm_finite_element_1()))),
3507
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_3LinearForm_dof_map_1()), mesh)))
3514
class Poisson2D_3BilinearForm : public dolfin::Form
3518
// Create form on given function space(s)
3519
Poisson2D_3BilinearForm(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1) : dolfin::Form()
3521
boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());
3522
_function_spaces.push_back(_V0);
3523
boost::shared_ptr<const dolfin::FunctionSpace> _V1(&V1, dolfin::NoDeleter<const dolfin::FunctionSpace>());
3524
_function_spaces.push_back(_V1);
3526
_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_3BilinearForm());
3529
// Create form on given function space(s) (shared data)
3530
Poisson2D_3BilinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::FunctionSpace> V1) : dolfin::Form()
3532
_function_spaces.push_back(V0);
3533
_function_spaces.push_back(V1);
3535
_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_3BilinearForm());
3539
~Poisson2D_3BilinearForm() {}
3543
class Poisson2D_3LinearFormCoefficient0 : public dolfin::Coefficient
3548
Poisson2D_3LinearFormCoefficient0(dolfin::Form& form) : dolfin::Coefficient(form) {}
3551
~Poisson2D_3LinearFormCoefficient0() {}
3553
// Attach function to coefficient
3554
const Poisson2D_3LinearFormCoefficient0& operator= (dolfin::Function& v)
3560
/// Create function space for coefficient
3561
const dolfin::FunctionSpace* create_function_space() const
3563
return new Poisson2D_3LinearFormCoefficientSpace0(form.mesh());
3566
/// Return coefficient number
3567
dolfin::uint number() const
3572
/// Return coefficient name
3573
virtual std::string name() const
3579
class Poisson2D_3LinearForm : public dolfin::Form
3583
// Create form on given function space(s)
3584
Poisson2D_3LinearForm(const dolfin::FunctionSpace& V0) : dolfin::Form(), f(*this)
3586
boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());
3587
_function_spaces.push_back(_V0);
3589
_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
3591
_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_3LinearForm());
3594
// Create form on given function space(s) (shared data)
3595
Poisson2D_3LinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0) : dolfin::Form(), f(*this)
3597
_function_spaces.push_back(V0);
3599
_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
3601
_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_3LinearForm());
3604
// Create form on given function space(s) with given coefficient(s)
3605
Poisson2D_3LinearForm(const dolfin::FunctionSpace& V0, dolfin::Function& w0) : dolfin::Form(), f(*this)
3607
boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());
3608
_function_spaces.push_back(_V0);
3610
_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
3614
_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_3LinearForm());
3617
// Create form on given function space(s) with given coefficient(s) (shared data)
3618
Poisson2D_3LinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, dolfin::Function& w0) : dolfin::Form(), f(*this)
3620
_function_spaces.push_back(V0);
3622
_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
3626
_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_3LinearForm());
3630
~Poisson2D_3LinearForm() {}
3633
Poisson2D_3LinearFormCoefficient0 f;