11
11
#include <stdexcept>
15
/// This class defines the interface for a finite element.
17
class UFC_Poisson2D_4BilinearForm_finite_element_0: public ufc::finite_element
22
UFC_Poisson2D_4BilinearForm_finite_element_0() : ufc::finite_element()
28
virtual ~UFC_Poisson2D_4BilinearForm_finite_element_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "FiniteElement('Lagrange', 'triangle', 4)";
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);
109
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
111
// Compute psitilde_a
112
const double psitilde_a_0 = 1;
113
const double psitilde_a_1 = x;
114
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
115
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
116
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
118
// Compute psitilde_bs
119
const double psitilde_bs_0_0 = 1;
120
const double psitilde_bs_0_1 = 1.5*y + 0.5;
121
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;
122
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;
123
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
124
const double psitilde_bs_1_0 = 1;
125
const double psitilde_bs_1_1 = 2.5*y + 1.5;
126
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;
127
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
128
const double psitilde_bs_2_0 = 1;
129
const double psitilde_bs_2_1 = 3.5*y + 2.5;
130
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
131
const double psitilde_bs_3_0 = 1;
132
const double psitilde_bs_3_1 = 4.5*y + 3.5;
133
const double psitilde_bs_4_0 = 1;
135
// Compute basisvalues
136
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
137
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
138
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
139
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
140
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
141
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
142
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
143
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
144
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
145
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
146
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
147
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
148
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
149
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
150
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
152
// Table(s) of coefficients
153
const static double coefficients0[15][15] = \
154
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366312, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
155
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233825, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.023944356611608, -0.0185472188784818, 0.0107082418122104},
156
{0, 0, 0.0476190476190476, 0, 0, 0.0388807895679869, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
157
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543895, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406279, 0.0834624849531681, -0.0535412090610519},
158
{-0.0314269680527355, 0.0109971479845644, 0.00634920634920629, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
159
{0.125707872210942, 0.0439885919382572, 0.126984126984127, 0, 0.0359165349174119, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
160
{0.125707872210942, -0.131965775814772, -0.0253968253968255, 0.139104141588614, 0.0718330698348238, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.047888713223216, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
161
{-0.0314269680527357, -0.0109971479845642, 0.00634920634920626, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
162
{0.125707872210942, -0.0439885919382573, 0.126984126984127, 0, -0.035916534917412, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
163
{0.125707872210942, -0.0879771838765143, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
164
{-0.0314269680527355, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951818, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.0239443566116079, 0, 0.0107082418122103},
165
{0.125707872210942, 0.0879771838765144, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277045, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
166
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406278, -0.0278208283177228, -0.0535412090610519},
167
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.026770604530526, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406279, 0.0278208283177227, -0.0535412090610519},
168
{0.251415744421884, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
170
// Extract relevant coefficients
171
const double coeff0_0 = coefficients0[dof][0];
172
const double coeff0_1 = coefficients0[dof][1];
173
const double coeff0_2 = coefficients0[dof][2];
174
const double coeff0_3 = coefficients0[dof][3];
175
const double coeff0_4 = coefficients0[dof][4];
176
const double coeff0_5 = coefficients0[dof][5];
177
const double coeff0_6 = coefficients0[dof][6];
178
const double coeff0_7 = coefficients0[dof][7];
179
const double coeff0_8 = coefficients0[dof][8];
180
const double coeff0_9 = coefficients0[dof][9];
181
const double coeff0_10 = coefficients0[dof][10];
182
const double coeff0_11 = coefficients0[dof][11];
183
const double coeff0_12 = coefficients0[dof][12];
184
const double coeff0_13 = coefficients0[dof][13];
185
const double coeff0_14 = coefficients0[dof][14];
188
*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 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
191
/// Evaluate all basis functions at given point in cell
192
virtual void evaluate_basis_all(double* values,
193
const double* coordinates,
194
const ufc::cell& c) const
196
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
199
/// Evaluate order n derivatives of basis function i at given point in cell
200
virtual void evaluate_basis_derivatives(unsigned int i,
203
const double* coordinates,
204
const ufc::cell& c) const
206
// Extract vertex coordinates
207
const double * const * element_coordinates = c.coordinates;
209
// Compute Jacobian of affine map from reference cell
210
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
211
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
212
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
213
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
215
// Compute determinant of Jacobian
216
const double detJ = J_00*J_11 - J_01*J_10;
218
// Compute inverse of Jacobian
220
// Get coordinates and map to the reference (UFC) element
221
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
222
element_coordinates[0][0]*element_coordinates[2][1] +\
223
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
224
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
225
element_coordinates[1][0]*element_coordinates[0][1] -\
226
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
228
// Map coordinates to the reference square
229
if (std::abs(y - 1.0) < 1e-14)
232
x = 2.0 *x/(1.0 - y) - 1.0;
235
// Compute number of derivatives
236
unsigned int num_derivatives = 1;
238
for (unsigned int j = 0; j < n; j++)
239
num_derivatives *= 2;
242
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
243
unsigned int **combinations = new unsigned int *[num_derivatives];
245
for (unsigned int j = 0; j < num_derivatives; j++)
247
combinations[j] = new unsigned int [n];
248
for (unsigned int k = 0; k < n; k++)
249
combinations[j][k] = 0;
252
// Generate combinations of derivatives
253
for (unsigned int row = 1; row < num_derivatives; row++)
255
for (unsigned int num = 0; num < row; num++)
257
for (unsigned int col = n-1; col+1 > 0; col--)
259
if (combinations[row][col] + 1 > 1)
260
combinations[row][col] = 0;
263
combinations[row][col] += 1;
270
// Compute inverse of Jacobian
271
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
273
// Declare transformation matrix
274
// Declare pointer to two dimensional array and initialise
275
double **transform = new double *[num_derivatives];
277
for (unsigned int j = 0; j < num_derivatives; j++)
279
transform[j] = new double [num_derivatives];
280
for (unsigned int k = 0; k < num_derivatives; k++)
284
// Construct transformation matrix
285
for (unsigned int row = 0; row < num_derivatives; row++)
287
for (unsigned int col = 0; col < num_derivatives; col++)
289
for (unsigned int k = 0; k < n; k++)
290
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
295
for (unsigned int j = 0; j < 1*num_derivatives; j++)
298
// Map degree of freedom to element degree of freedom
299
const unsigned int dof = i;
302
const double scalings_y_0 = 1;
303
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
304
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
305
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
306
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
308
// Compute psitilde_a
309
const double psitilde_a_0 = 1;
310
const double psitilde_a_1 = x;
311
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
312
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
313
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
315
// Compute psitilde_bs
316
const double psitilde_bs_0_0 = 1;
317
const double psitilde_bs_0_1 = 1.5*y + 0.5;
318
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;
319
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;
320
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
321
const double psitilde_bs_1_0 = 1;
322
const double psitilde_bs_1_1 = 2.5*y + 1.5;
323
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;
324
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
325
const double psitilde_bs_2_0 = 1;
326
const double psitilde_bs_2_1 = 3.5*y + 2.5;
327
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
328
const double psitilde_bs_3_0 = 1;
329
const double psitilde_bs_3_1 = 4.5*y + 3.5;
330
const double psitilde_bs_4_0 = 1;
332
// Compute basisvalues
333
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
334
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
335
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
336
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
337
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
338
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
339
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
340
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
341
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
342
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
343
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
344
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
345
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
346
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
347
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
349
// Table(s) of coefficients
350
const static double coefficients0[15][15] = \
351
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366312, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
352
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233825, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.023944356611608, -0.0185472188784818, 0.0107082418122104},
353
{0, 0, 0.0476190476190476, 0, 0, 0.0388807895679869, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
354
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543895, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406279, 0.0834624849531681, -0.0535412090610519},
355
{-0.0314269680527355, 0.0109971479845644, 0.00634920634920629, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
356
{0.125707872210942, 0.0439885919382572, 0.126984126984127, 0, 0.0359165349174119, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
357
{0.125707872210942, -0.131965775814772, -0.0253968253968255, 0.139104141588614, 0.0718330698348238, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.047888713223216, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
358
{-0.0314269680527357, -0.0109971479845642, 0.00634920634920626, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
359
{0.125707872210942, -0.0439885919382573, 0.126984126984127, 0, -0.035916534917412, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
360
{0.125707872210942, -0.0879771838765143, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
361
{-0.0314269680527355, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951818, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.0239443566116079, 0, 0.0107082418122103},
362
{0.125707872210942, 0.0879771838765144, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277045, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
363
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406278, -0.0278208283177228, -0.0535412090610519},
364
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.026770604530526, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406279, 0.0278208283177227, -0.0535412090610519},
365
{0.251415744421884, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
367
// Interesting (new) part
368
// Tables of derivatives of the polynomial base (transpose)
369
const static double dmats0[15][15] = \
370
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
371
{4.89897948556635, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
372
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
373
{0, 9.48683298050513, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
374
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
375
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
376
{5.29150262212917, 0, -2.99332590941916, 13.6626010212795, 0, 0.611010092660778, 0, 0, 0, 0, 0, 0, 0, 0, 0},
377
{0, 4.38178046004132, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
378
{3.46410161513776, 0, 7.83836717690617, 0, 0, 8.40000000000001, 0, 0, 0, 0, 0, 0, 0, 0, 0},
379
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
380
{0, 10.9544511501033, 0, 0, -3.83325938999965, 0, 17.7482393492988, 0, 0.553283335172492, 0, 0, 0, 0, 0, 0},
381
{4.73286382647968, 0, 3.3466401061363, 4.36435780471985, 0, -5.07468037933239, 0, 17.0084012854152, 0, 1.52127765851133, 0, 0, 0, 0, 0},
382
{0, 2.44948974278317, 0, 0, 9.14285714285714, 0, 0, 0, 14.8461497791618, 0, 0, 0, 0, 0, 0},
383
{3.09838667696594, 0, 7.66811580507233, 0, 0, 10.733126291999, 0, 0, 0, 9.2951600308978, 0, 0, 0, 0, 0},
384
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
386
const static double dmats1[15][15] = \
387
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
388
{2.44948974278317, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
389
{4.24264068711928, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
390
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
391
{2, 6.12372435695794, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
392
{-2.30940107675849, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
393
{2.64575131106458, 5.18459255872628, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330385, 0, 0, 0, 0, 0, 0, 0, 0, 0},
394
{2.23606797749978, 2.19089023002066, 2.5298221281347, 8.08290376865476, 6.26099033699942, -1.80739222823014, 0, 0, 0, 0, 0, 0, 0, 0, 0},
395
{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238572, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
396
{5, 0, -2.8284271247462, 0, 0, 12.1243556529821, 0, 0, 0, 0, 0, 0, 0, 0, 0},
397
{2.68328157299974, 5.47722557505166, -1.89736659610103, 7.4230748895809, -1.91662969499982, 0.663940002206987, 8.87411967464942, -1.07142857142857, 0.276641667586245, -0.095831484749991, 0, 0, 0, 0, 0},
398
{2.36643191323984, 2.89827534923788, 1.67332005306815, 2.18217890235993, 5.74704893215391, -2.53734018966619, 10.0623058987491, 8.50420064270761, -2.1957751641342, 0.760638829255663, 0, 0, 0, 0, 0},
399
{2, 1.22474487139159, 3.53553390593274, -7.37711113563317, 4.57142857142857, 1.64957219768464, 0, 11.4997781699989, 7.4230748895809, -2.57142857142858, 0, 0, 0, 0, 0},
400
{1.54919333848296, 6.6407830863536, 3.83405790253616, 0, -6.19677335393188, 5.3665631459995, 0, 0, 13.4164078649987, 4.6475800154489, 0, 0, 0, 0, 0},
401
{-3.57770876399967, 0, 8.85437744847147, 0, 0, -3.09838667696593, 0, 0, 0, 16.0996894379985, 0, 0, 0, 0, 0}};
403
// Compute reference derivatives
404
// Declare pointer to array of derivatives on FIAT element
405
double *derivatives = new double [num_derivatives];
407
// Declare coefficients
418
double coeff0_10 = 0;
419
double coeff0_11 = 0;
420
double coeff0_12 = 0;
421
double coeff0_13 = 0;
422
double coeff0_14 = 0;
424
// Declare new coefficients
425
double new_coeff0_0 = 0;
426
double new_coeff0_1 = 0;
427
double new_coeff0_2 = 0;
428
double new_coeff0_3 = 0;
429
double new_coeff0_4 = 0;
430
double new_coeff0_5 = 0;
431
double new_coeff0_6 = 0;
432
double new_coeff0_7 = 0;
433
double new_coeff0_8 = 0;
434
double new_coeff0_9 = 0;
435
double new_coeff0_10 = 0;
436
double new_coeff0_11 = 0;
437
double new_coeff0_12 = 0;
438
double new_coeff0_13 = 0;
439
double new_coeff0_14 = 0;
441
// Loop possible derivatives
442
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
444
// Get values from coefficients array
445
new_coeff0_0 = coefficients0[dof][0];
446
new_coeff0_1 = coefficients0[dof][1];
447
new_coeff0_2 = coefficients0[dof][2];
448
new_coeff0_3 = coefficients0[dof][3];
449
new_coeff0_4 = coefficients0[dof][4];
450
new_coeff0_5 = coefficients0[dof][5];
451
new_coeff0_6 = coefficients0[dof][6];
452
new_coeff0_7 = coefficients0[dof][7];
453
new_coeff0_8 = coefficients0[dof][8];
454
new_coeff0_9 = coefficients0[dof][9];
455
new_coeff0_10 = coefficients0[dof][10];
456
new_coeff0_11 = coefficients0[dof][11];
457
new_coeff0_12 = coefficients0[dof][12];
458
new_coeff0_13 = coefficients0[dof][13];
459
new_coeff0_14 = coefficients0[dof][14];
461
// Loop derivative order
462
for (unsigned int j = 0; j < n; j++)
464
// Update old coefficients
465
coeff0_0 = new_coeff0_0;
466
coeff0_1 = new_coeff0_1;
467
coeff0_2 = new_coeff0_2;
468
coeff0_3 = new_coeff0_3;
469
coeff0_4 = new_coeff0_4;
470
coeff0_5 = new_coeff0_5;
471
coeff0_6 = new_coeff0_6;
472
coeff0_7 = new_coeff0_7;
473
coeff0_8 = new_coeff0_8;
474
coeff0_9 = new_coeff0_9;
475
coeff0_10 = new_coeff0_10;
476
coeff0_11 = new_coeff0_11;
477
coeff0_12 = new_coeff0_12;
478
coeff0_13 = new_coeff0_13;
479
coeff0_14 = new_coeff0_14;
481
if(combinations[deriv_num][j] == 0)
483
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] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
484
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] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
485
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] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
486
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] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
487
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] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
488
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] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
489
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] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
490
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] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
491
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] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
492
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] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
493
new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
494
new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
495
new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
496
new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
497
new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
499
if(combinations[deriv_num][j] == 1)
501
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] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
502
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] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
503
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] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
504
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] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
505
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] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
506
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] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
507
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] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
508
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] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
509
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] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
510
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] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
511
new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
512
new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
513
new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
514
new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
515
new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
519
// Compute derivatives on reference element as dot product of coefficients and basisvalues
520
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 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
523
// Transform derivatives back to physical element
524
for (unsigned int row = 0; row < num_derivatives; row++)
526
for (unsigned int col = 0; col < num_derivatives; col++)
528
values[row] += transform[row][col]*derivatives[col];
531
// Delete pointer to array of derivatives on FIAT element
532
delete [] derivatives;
534
// Delete pointer to array of combinations of derivatives and transform
535
for (unsigned int row = 0; row < num_derivatives; row++)
537
delete [] combinations[row];
538
delete [] transform[row];
541
delete [] combinations;
545
/// Evaluate order n derivatives of all basis functions at given point in cell
546
virtual void evaluate_basis_derivatives_all(unsigned int n,
548
const double* coordinates,
549
const ufc::cell& c) const
551
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
554
/// Evaluate linear functional for dof i on the function f
555
virtual double evaluate_dof(unsigned int i,
556
const ufc::function& f,
557
const ufc::cell& c) const
559
// The reference points, direction and weights:
560
const static double X[15][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}, {{0.25, 0.25}}, {{0.5, 0.25}}, {{0.25, 0.5}}};
561
const static double W[15][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
562
const static double D[15][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
564
const double * const * x = c.coordinates;
566
// Iterate over the points:
567
// Evaluate basis functions for affine mapping
568
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
569
const double w1 = X[i][0][0];
570
const double w2 = X[i][0][1];
572
// Compute affine mapping y = F(X)
574
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
575
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
577
// Evaluate function at physical points
579
f.evaluate(values, y, c);
581
// Map function values using appropriate mapping
582
// Affine map: Do nothing
584
// Note that we do not map the weights (yet).
586
// Take directional components
587
for(int k = 0; k < 1; k++)
588
result += values[k]*D[i][0][k];
589
// Multiply by weights
595
/// Evaluate linear functionals for all dofs on the function f
596
virtual void evaluate_dofs(double* values,
597
const ufc::function& f,
598
const ufc::cell& c) const
600
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
603
/// Interpolate vertex values from dof values
604
virtual void interpolate_vertex_values(double* vertex_values,
605
const double* dof_values,
606
const ufc::cell& c) const
608
// Evaluate at vertices and use affine mapping
609
vertex_values[0] = dof_values[0];
610
vertex_values[1] = dof_values[1];
611
vertex_values[2] = dof_values[2];
614
/// Return the number of sub elements (for a mixed element)
615
virtual unsigned int num_sub_elements() const
620
/// Create a new finite element for sub element i (for a mixed element)
621
virtual ufc::finite_element* create_sub_element(unsigned int i) const
623
return new UFC_Poisson2D_4BilinearForm_finite_element_0();
628
/// This class defines the interface for a finite element.
630
class UFC_Poisson2D_4BilinearForm_finite_element_1: public ufc::finite_element
635
UFC_Poisson2D_4BilinearForm_finite_element_1() : ufc::finite_element()
641
virtual ~UFC_Poisson2D_4BilinearForm_finite_element_1()
646
/// Return a string identifying the finite element
647
virtual const char* signature() const
649
return "FiniteElement('Lagrange', 'triangle', 4)";
652
/// Return the cell shape
653
virtual ufc::shape cell_shape() const
655
return ufc::triangle;
658
/// Return the dimension of the finite element function space
659
virtual unsigned int space_dimension() const
664
/// Return the rank of the value space
665
virtual unsigned int value_rank() const
670
/// Return the dimension of the value space for axis i
671
virtual unsigned int value_dimension(unsigned int i) const
676
/// Evaluate basis function i at given point in cell
677
virtual void evaluate_basis(unsigned int i,
679
const double* coordinates,
680
const ufc::cell& c) const
682
// Extract vertex coordinates
683
const double * const * element_coordinates = c.coordinates;
685
// Compute Jacobian of affine map from reference cell
686
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
687
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
688
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
689
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
691
// Compute determinant of Jacobian
692
const double detJ = J_00*J_11 - J_01*J_10;
694
// Compute inverse of Jacobian
696
// Get coordinates and map to the reference (UFC) element
697
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
698
element_coordinates[0][0]*element_coordinates[2][1] +\
699
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
700
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
701
element_coordinates[1][0]*element_coordinates[0][1] -\
702
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
704
// Map coordinates to the reference square
705
if (std::abs(y - 1.0) < 1e-14)
708
x = 2.0 *x/(1.0 - y) - 1.0;
714
// Map degree of freedom to element degree of freedom
715
const unsigned int dof = i;
718
const double scalings_y_0 = 1;
719
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
720
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
721
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
722
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
724
// Compute psitilde_a
725
const double psitilde_a_0 = 1;
726
const double psitilde_a_1 = x;
727
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
728
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
729
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
731
// Compute psitilde_bs
732
const double psitilde_bs_0_0 = 1;
733
const double psitilde_bs_0_1 = 1.5*y + 0.5;
734
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;
735
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;
736
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
737
const double psitilde_bs_1_0 = 1;
738
const double psitilde_bs_1_1 = 2.5*y + 1.5;
739
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;
740
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
741
const double psitilde_bs_2_0 = 1;
742
const double psitilde_bs_2_1 = 3.5*y + 2.5;
743
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
744
const double psitilde_bs_3_0 = 1;
745
const double psitilde_bs_3_1 = 4.5*y + 3.5;
746
const double psitilde_bs_4_0 = 1;
748
// Compute basisvalues
749
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
750
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
751
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
752
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
753
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
754
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
755
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
756
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
757
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
758
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
759
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
760
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
761
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
762
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
763
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
765
// Table(s) of coefficients
766
const static double coefficients0[15][15] = \
767
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366312, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
768
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233825, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.023944356611608, -0.0185472188784818, 0.0107082418122104},
769
{0, 0, 0.0476190476190476, 0, 0, 0.0388807895679869, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
770
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543895, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406279, 0.0834624849531681, -0.0535412090610519},
771
{-0.0314269680527355, 0.0109971479845644, 0.00634920634920629, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
772
{0.125707872210942, 0.0439885919382572, 0.126984126984127, 0, 0.0359165349174119, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
773
{0.125707872210942, -0.131965775814772, -0.0253968253968255, 0.139104141588614, 0.0718330698348238, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.047888713223216, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
774
{-0.0314269680527357, -0.0109971479845642, 0.00634920634920626, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
775
{0.125707872210942, -0.0439885919382573, 0.126984126984127, 0, -0.035916534917412, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
776
{0.125707872210942, -0.0879771838765143, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
777
{-0.0314269680527355, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951818, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.0239443566116079, 0, 0.0107082418122103},
778
{0.125707872210942, 0.0879771838765144, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277045, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
779
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406278, -0.0278208283177228, -0.0535412090610519},
780
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.026770604530526, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406279, 0.0278208283177227, -0.0535412090610519},
781
{0.251415744421884, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
783
// Extract relevant coefficients
784
const double coeff0_0 = coefficients0[dof][0];
785
const double coeff0_1 = coefficients0[dof][1];
786
const double coeff0_2 = coefficients0[dof][2];
787
const double coeff0_3 = coefficients0[dof][3];
788
const double coeff0_4 = coefficients0[dof][4];
789
const double coeff0_5 = coefficients0[dof][5];
790
const double coeff0_6 = coefficients0[dof][6];
791
const double coeff0_7 = coefficients0[dof][7];
792
const double coeff0_8 = coefficients0[dof][8];
793
const double coeff0_9 = coefficients0[dof][9];
794
const double coeff0_10 = coefficients0[dof][10];
795
const double coeff0_11 = coefficients0[dof][11];
796
const double coeff0_12 = coefficients0[dof][12];
797
const double coeff0_13 = coefficients0[dof][13];
798
const double coeff0_14 = coefficients0[dof][14];
801
*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 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
804
/// Evaluate all basis functions at given point in cell
805
virtual void evaluate_basis_all(double* values,
806
const double* coordinates,
807
const ufc::cell& c) const
809
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
812
/// Evaluate order n derivatives of basis function i at given point in cell
813
virtual void evaluate_basis_derivatives(unsigned int i,
816
const double* coordinates,
817
const ufc::cell& c) const
819
// Extract vertex coordinates
820
const double * const * element_coordinates = c.coordinates;
822
// Compute Jacobian of affine map from reference cell
823
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
824
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
825
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
826
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
828
// Compute determinant of Jacobian
829
const double detJ = J_00*J_11 - J_01*J_10;
831
// Compute inverse of Jacobian
833
// Get coordinates and map to the reference (UFC) element
834
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
835
element_coordinates[0][0]*element_coordinates[2][1] +\
836
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
837
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
838
element_coordinates[1][0]*element_coordinates[0][1] -\
839
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
841
// Map coordinates to the reference square
842
if (std::abs(y - 1.0) < 1e-14)
845
x = 2.0 *x/(1.0 - y) - 1.0;
848
// Compute number of derivatives
849
unsigned int num_derivatives = 1;
851
for (unsigned int j = 0; j < n; j++)
852
num_derivatives *= 2;
855
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
856
unsigned int **combinations = new unsigned int *[num_derivatives];
858
for (unsigned int j = 0; j < num_derivatives; j++)
860
combinations[j] = new unsigned int [n];
861
for (unsigned int k = 0; k < n; k++)
862
combinations[j][k] = 0;
865
// Generate combinations of derivatives
866
for (unsigned int row = 1; row < num_derivatives; row++)
868
for (unsigned int num = 0; num < row; num++)
870
for (unsigned int col = n-1; col+1 > 0; col--)
872
if (combinations[row][col] + 1 > 1)
873
combinations[row][col] = 0;
876
combinations[row][col] += 1;
883
// Compute inverse of Jacobian
884
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
886
// Declare transformation matrix
887
// Declare pointer to two dimensional array and initialise
888
double **transform = new double *[num_derivatives];
890
for (unsigned int j = 0; j < num_derivatives; j++)
892
transform[j] = new double [num_derivatives];
893
for (unsigned int k = 0; k < num_derivatives; k++)
897
// Construct transformation matrix
898
for (unsigned int row = 0; row < num_derivatives; row++)
900
for (unsigned int col = 0; col < num_derivatives; col++)
902
for (unsigned int k = 0; k < n; k++)
903
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
908
for (unsigned int j = 0; j < 1*num_derivatives; j++)
911
// Map degree of freedom to element degree of freedom
912
const unsigned int dof = i;
915
const double scalings_y_0 = 1;
916
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
917
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
918
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
919
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
921
// Compute psitilde_a
922
const double psitilde_a_0 = 1;
923
const double psitilde_a_1 = x;
924
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
925
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
926
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
928
// Compute psitilde_bs
929
const double psitilde_bs_0_0 = 1;
930
const double psitilde_bs_0_1 = 1.5*y + 0.5;
931
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;
932
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;
933
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
934
const double psitilde_bs_1_0 = 1;
935
const double psitilde_bs_1_1 = 2.5*y + 1.5;
936
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;
937
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
938
const double psitilde_bs_2_0 = 1;
939
const double psitilde_bs_2_1 = 3.5*y + 2.5;
940
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
941
const double psitilde_bs_3_0 = 1;
942
const double psitilde_bs_3_1 = 4.5*y + 3.5;
943
const double psitilde_bs_4_0 = 1;
945
// Compute basisvalues
946
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
947
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
948
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
949
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
950
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
951
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
952
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
953
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
954
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
955
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
956
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
957
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
958
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
959
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
960
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
962
// Table(s) of coefficients
963
const static double coefficients0[15][15] = \
964
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366312, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
965
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233825, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.023944356611608, -0.0185472188784818, 0.0107082418122104},
966
{0, 0, 0.0476190476190476, 0, 0, 0.0388807895679869, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
967
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543895, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406279, 0.0834624849531681, -0.0535412090610519},
968
{-0.0314269680527355, 0.0109971479845644, 0.00634920634920629, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
969
{0.125707872210942, 0.0439885919382572, 0.126984126984127, 0, 0.0359165349174119, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
970
{0.125707872210942, -0.131965775814772, -0.0253968253968255, 0.139104141588614, 0.0718330698348238, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.047888713223216, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
971
{-0.0314269680527357, -0.0109971479845642, 0.00634920634920626, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
972
{0.125707872210942, -0.0439885919382573, 0.126984126984127, 0, -0.035916534917412, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
973
{0.125707872210942, -0.0879771838765143, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
974
{-0.0314269680527355, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951818, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.0239443566116079, 0, 0.0107082418122103},
975
{0.125707872210942, 0.0879771838765144, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277045, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
976
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406278, -0.0278208283177228, -0.0535412090610519},
977
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.026770604530526, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406279, 0.0278208283177227, -0.0535412090610519},
978
{0.251415744421884, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
980
// Interesting (new) part
981
// Tables of derivatives of the polynomial base (transpose)
982
const static double dmats0[15][15] = \
983
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
984
{4.89897948556635, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
985
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
986
{0, 9.48683298050513, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
987
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
988
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
989
{5.29150262212917, 0, -2.99332590941916, 13.6626010212795, 0, 0.611010092660778, 0, 0, 0, 0, 0, 0, 0, 0, 0},
990
{0, 4.38178046004132, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
991
{3.46410161513776, 0, 7.83836717690617, 0, 0, 8.40000000000001, 0, 0, 0, 0, 0, 0, 0, 0, 0},
992
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
993
{0, 10.9544511501033, 0, 0, -3.83325938999965, 0, 17.7482393492988, 0, 0.553283335172492, 0, 0, 0, 0, 0, 0},
994
{4.73286382647968, 0, 3.3466401061363, 4.36435780471985, 0, -5.07468037933239, 0, 17.0084012854152, 0, 1.52127765851133, 0, 0, 0, 0, 0},
995
{0, 2.44948974278317, 0, 0, 9.14285714285714, 0, 0, 0, 14.8461497791618, 0, 0, 0, 0, 0, 0},
996
{3.09838667696594, 0, 7.66811580507233, 0, 0, 10.733126291999, 0, 0, 0, 9.2951600308978, 0, 0, 0, 0, 0},
997
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
999
const static double dmats1[15][15] = \
1000
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1001
{2.44948974278317, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1002
{4.24264068711928, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1003
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1004
{2, 6.12372435695794, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1005
{-2.30940107675849, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1006
{2.64575131106458, 5.18459255872628, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330385, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1007
{2.23606797749978, 2.19089023002066, 2.5298221281347, 8.08290376865476, 6.26099033699942, -1.80739222823014, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1008
{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238572, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1009
{5, 0, -2.8284271247462, 0, 0, 12.1243556529821, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1010
{2.68328157299974, 5.47722557505166, -1.89736659610103, 7.4230748895809, -1.91662969499982, 0.663940002206987, 8.87411967464942, -1.07142857142857, 0.276641667586245, -0.095831484749991, 0, 0, 0, 0, 0},
1011
{2.36643191323984, 2.89827534923788, 1.67332005306815, 2.18217890235993, 5.74704893215391, -2.53734018966619, 10.0623058987491, 8.50420064270761, -2.1957751641342, 0.760638829255663, 0, 0, 0, 0, 0},
1012
{2, 1.22474487139159, 3.53553390593274, -7.37711113563317, 4.57142857142857, 1.64957219768464, 0, 11.4997781699989, 7.4230748895809, -2.57142857142858, 0, 0, 0, 0, 0},
1013
{1.54919333848296, 6.6407830863536, 3.83405790253616, 0, -6.19677335393188, 5.3665631459995, 0, 0, 13.4164078649987, 4.6475800154489, 0, 0, 0, 0, 0},
1014
{-3.57770876399967, 0, 8.85437744847147, 0, 0, -3.09838667696593, 0, 0, 0, 16.0996894379985, 0, 0, 0, 0, 0}};
1016
// Compute reference derivatives
1017
// Declare pointer to array of derivatives on FIAT element
1018
double *derivatives = new double [num_derivatives];
1020
// Declare coefficients
1021
double coeff0_0 = 0;
1022
double coeff0_1 = 0;
1023
double coeff0_2 = 0;
1024
double coeff0_3 = 0;
1025
double coeff0_4 = 0;
1026
double coeff0_5 = 0;
1027
double coeff0_6 = 0;
1028
double coeff0_7 = 0;
1029
double coeff0_8 = 0;
1030
double coeff0_9 = 0;
1031
double coeff0_10 = 0;
1032
double coeff0_11 = 0;
1033
double coeff0_12 = 0;
1034
double coeff0_13 = 0;
1035
double coeff0_14 = 0;
1037
// Declare new coefficients
1038
double new_coeff0_0 = 0;
1039
double new_coeff0_1 = 0;
1040
double new_coeff0_2 = 0;
1041
double new_coeff0_3 = 0;
1042
double new_coeff0_4 = 0;
1043
double new_coeff0_5 = 0;
1044
double new_coeff0_6 = 0;
1045
double new_coeff0_7 = 0;
1046
double new_coeff0_8 = 0;
1047
double new_coeff0_9 = 0;
1048
double new_coeff0_10 = 0;
1049
double new_coeff0_11 = 0;
1050
double new_coeff0_12 = 0;
1051
double new_coeff0_13 = 0;
1052
double new_coeff0_14 = 0;
1054
// Loop possible derivatives
1055
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1057
// Get values from coefficients array
1058
new_coeff0_0 = coefficients0[dof][0];
1059
new_coeff0_1 = coefficients0[dof][1];
1060
new_coeff0_2 = coefficients0[dof][2];
1061
new_coeff0_3 = coefficients0[dof][3];
1062
new_coeff0_4 = coefficients0[dof][4];
1063
new_coeff0_5 = coefficients0[dof][5];
1064
new_coeff0_6 = coefficients0[dof][6];
1065
new_coeff0_7 = coefficients0[dof][7];
1066
new_coeff0_8 = coefficients0[dof][8];
1067
new_coeff0_9 = coefficients0[dof][9];
1068
new_coeff0_10 = coefficients0[dof][10];
1069
new_coeff0_11 = coefficients0[dof][11];
1070
new_coeff0_12 = coefficients0[dof][12];
1071
new_coeff0_13 = coefficients0[dof][13];
1072
new_coeff0_14 = coefficients0[dof][14];
1074
// Loop derivative order
1075
for (unsigned int j = 0; j < n; j++)
1077
// Update old coefficients
1078
coeff0_0 = new_coeff0_0;
1079
coeff0_1 = new_coeff0_1;
1080
coeff0_2 = new_coeff0_2;
1081
coeff0_3 = new_coeff0_3;
1082
coeff0_4 = new_coeff0_4;
1083
coeff0_5 = new_coeff0_5;
1084
coeff0_6 = new_coeff0_6;
1085
coeff0_7 = new_coeff0_7;
1086
coeff0_8 = new_coeff0_8;
1087
coeff0_9 = new_coeff0_9;
1088
coeff0_10 = new_coeff0_10;
1089
coeff0_11 = new_coeff0_11;
1090
coeff0_12 = new_coeff0_12;
1091
coeff0_13 = new_coeff0_13;
1092
coeff0_14 = new_coeff0_14;
1094
if(combinations[deriv_num][j] == 0)
1096
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] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
1097
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] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
1098
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] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
1099
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] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
1100
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] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
1101
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] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
1102
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] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
1103
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] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
1104
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] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
1105
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] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
1106
new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
1107
new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
1108
new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
1109
new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
1110
new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
1112
if(combinations[deriv_num][j] == 1)
1114
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] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
1115
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] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
1116
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] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
1117
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] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
1118
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] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
1119
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] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
1120
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] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
1121
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] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
1122
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] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
1123
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] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
1124
new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
1125
new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
1126
new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
1127
new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
1128
new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
1132
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1133
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 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
1136
// Transform derivatives back to physical element
1137
for (unsigned int row = 0; row < num_derivatives; row++)
1139
for (unsigned int col = 0; col < num_derivatives; col++)
1141
values[row] += transform[row][col]*derivatives[col];
1144
// Delete pointer to array of derivatives on FIAT element
1145
delete [] derivatives;
1147
// Delete pointer to array of combinations of derivatives and transform
1148
for (unsigned int row = 0; row < num_derivatives; row++)
1150
delete [] combinations[row];
1151
delete [] transform[row];
1154
delete [] combinations;
1155
delete [] transform;
1158
/// Evaluate order n derivatives of all basis functions at given point in cell
1159
virtual void evaluate_basis_derivatives_all(unsigned int n,
1161
const double* coordinates,
1162
const ufc::cell& c) const
1164
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1167
/// Evaluate linear functional for dof i on the function f
1168
virtual double evaluate_dof(unsigned int i,
1169
const ufc::function& f,
1170
const ufc::cell& c) const
1172
// The reference points, direction and weights:
1173
const static double X[15][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}, {{0.25, 0.25}}, {{0.5, 0.25}}, {{0.25, 0.5}}};
1174
const static double W[15][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
1175
const static double D[15][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
1177
const double * const * x = c.coordinates;
1178
double result = 0.0;
1179
// Iterate over the points:
1180
// Evaluate basis functions for affine mapping
1181
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1182
const double w1 = X[i][0][0];
1183
const double w2 = X[i][0][1];
1185
// Compute affine mapping y = F(X)
1187
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1188
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1190
// Evaluate function at physical points
1192
f.evaluate(values, y, c);
1194
// Map function values using appropriate mapping
1195
// Affine map: Do nothing
1197
// Note that we do not map the weights (yet).
1199
// Take directional components
1200
for(int k = 0; k < 1; k++)
1201
result += values[k]*D[i][0][k];
1202
// Multiply by weights
1208
/// Evaluate linear functionals for all dofs on the function f
1209
virtual void evaluate_dofs(double* values,
1210
const ufc::function& f,
1211
const ufc::cell& c) const
1213
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1216
/// Interpolate vertex values from dof values
1217
virtual void interpolate_vertex_values(double* vertex_values,
1218
const double* dof_values,
1219
const ufc::cell& c) const
1221
// Evaluate at vertices and use affine mapping
1222
vertex_values[0] = dof_values[0];
1223
vertex_values[1] = dof_values[1];
1224
vertex_values[2] = dof_values[2];
1227
/// Return the number of sub elements (for a mixed element)
1228
virtual unsigned int num_sub_elements() const
1233
/// Create a new finite element for sub element i (for a mixed element)
1234
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1236
return new UFC_Poisson2D_4BilinearForm_finite_element_1();
1241
/// This class defines the interface for a local-to-global mapping of
1242
/// degrees of freedom (dofs).
1244
class UFC_Poisson2D_4BilinearForm_dof_map_0: public ufc::dof_map
1248
unsigned int __global_dimension;
1253
UFC_Poisson2D_4BilinearForm_dof_map_0() : ufc::dof_map()
1255
__global_dimension = 0;
1259
virtual ~UFC_Poisson2D_4BilinearForm_dof_map_0()
1264
/// Return a string identifying the dof map
1265
virtual const char* signature() const
1267
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 4)";
1270
/// Return true iff mesh entities of topological dimension d are needed
1271
virtual bool needs_mesh_entities(unsigned int d) const
1288
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1289
virtual bool init_mesh(const ufc::mesh& m)
1291
__global_dimension = m.num_entities[0] + 3*m.num_entities[1] + 3*m.num_entities[2];
1295
/// Initialize dof map for given cell
1296
virtual void init_cell(const ufc::mesh& m,
1302
/// Finish initialization of dof map for cells
1303
virtual void init_cell_finalize()
1308
/// Return the dimension of the global finite element function space
1309
virtual unsigned int global_dimension() const
1311
return __global_dimension;
1314
/// Return the dimension of the local finite element function space
1315
virtual unsigned int local_dimension() const
1320
// Return the geometric dimension of the coordinates this dof map provides
1321
virtual unsigned int geometric_dimension() const
1326
/// Return the number of dofs on each cell facet
1327
virtual unsigned int num_facet_dofs() const
1332
/// Return the number of dofs associated with each cell entity of dimension d
1333
virtual unsigned int num_entity_dofs(unsigned int d) const
1335
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1338
/// Tabulate the local-to-global mapping of dofs on a cell
1339
virtual void tabulate_dofs(unsigned int* dofs,
1341
const ufc::cell& c) const
1343
dofs[0] = c.entity_indices[0][0];
1344
dofs[1] = c.entity_indices[0][1];
1345
dofs[2] = c.entity_indices[0][2];
1346
unsigned int offset = m.num_entities[0];
1347
dofs[3] = offset + 3*c.entity_indices[1][0];
1348
dofs[4] = offset + 3*c.entity_indices[1][0] + 1;
1349
dofs[5] = offset + 3*c.entity_indices[1][0] + 2;
1350
dofs[6] = offset + 3*c.entity_indices[1][1];
1351
dofs[7] = offset + 3*c.entity_indices[1][1] + 1;
1352
dofs[8] = offset + 3*c.entity_indices[1][1] + 2;
1353
dofs[9] = offset + 3*c.entity_indices[1][2];
1354
dofs[10] = offset + 3*c.entity_indices[1][2] + 1;
1355
dofs[11] = offset + 3*c.entity_indices[1][2] + 2;
1356
offset = offset + 3*m.num_entities[1];
1357
dofs[12] = offset + 3*c.entity_indices[2][0];
1358
dofs[13] = offset + 3*c.entity_indices[2][0] + 1;
1359
dofs[14] = offset + 3*c.entity_indices[2][0] + 2;
1362
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1363
virtual void tabulate_facet_dofs(unsigned int* dofs,
1364
unsigned int facet) const
1392
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1393
virtual void tabulate_entity_dofs(unsigned int* dofs,
1394
unsigned int d, unsigned int i) const
1396
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1399
/// Tabulate the coordinates of all dofs on a cell
1400
virtual void tabulate_coordinates(double** coordinates,
1401
const ufc::cell& c) const
1403
const double * const * x = c.coordinates;
1404
coordinates[0][0] = x[0][0];
1405
coordinates[0][1] = x[0][1];
1406
coordinates[1][0] = x[1][0];
1407
coordinates[1][1] = x[1][1];
1408
coordinates[2][0] = x[2][0];
1409
coordinates[2][1] = x[2][1];
1410
coordinates[3][0] = 0.75*x[1][0] + 0.25*x[2][0];
1411
coordinates[3][1] = 0.75*x[1][1] + 0.25*x[2][1];
1412
coordinates[4][0] = 0.5*x[1][0] + 0.5*x[2][0];
1413
coordinates[4][1] = 0.5*x[1][1] + 0.5*x[2][1];
1414
coordinates[5][0] = 0.25*x[1][0] + 0.75*x[2][0];
1415
coordinates[5][1] = 0.25*x[1][1] + 0.75*x[2][1];
1416
coordinates[6][0] = 0.75*x[0][0] + 0.25*x[2][0];
1417
coordinates[6][1] = 0.75*x[0][1] + 0.25*x[2][1];
1418
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[2][0];
1419
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[2][1];
1420
coordinates[8][0] = 0.25*x[0][0] + 0.75*x[2][0];
1421
coordinates[8][1] = 0.25*x[0][1] + 0.75*x[2][1];
1422
coordinates[9][0] = 0.75*x[0][0] + 0.25*x[1][0];
1423
coordinates[9][1] = 0.75*x[0][1] + 0.25*x[1][1];
1424
coordinates[10][0] = 0.5*x[0][0] + 0.5*x[1][0];
1425
coordinates[10][1] = 0.5*x[0][1] + 0.5*x[1][1];
1426
coordinates[11][0] = 0.25*x[0][0] + 0.75*x[1][0];
1427
coordinates[11][1] = 0.25*x[0][1] + 0.75*x[1][1];
1428
coordinates[12][0] = 0.5*x[0][0] + 0.25*x[1][0] + 0.25*x[2][0];
1429
coordinates[12][1] = 0.5*x[0][1] + 0.25*x[1][1] + 0.25*x[2][1];
1430
coordinates[13][0] = 0.25*x[0][0] + 0.5*x[1][0] + 0.25*x[2][0];
1431
coordinates[13][1] = 0.25*x[0][1] + 0.5*x[1][1] + 0.25*x[2][1];
1432
coordinates[14][0] = 0.25*x[0][0] + 0.25*x[1][0] + 0.5*x[2][0];
1433
coordinates[14][1] = 0.25*x[0][1] + 0.25*x[1][1] + 0.5*x[2][1];
1436
/// Return the number of sub dof maps (for a mixed element)
1437
virtual unsigned int num_sub_dof_maps() const
1442
/// Create a new dof_map for sub dof map i (for a mixed element)
1443
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1445
return new UFC_Poisson2D_4BilinearForm_dof_map_0();
1450
/// This class defines the interface for a local-to-global mapping of
1451
/// degrees of freedom (dofs).
1453
class UFC_Poisson2D_4BilinearForm_dof_map_1: public ufc::dof_map
1457
unsigned int __global_dimension;
1462
UFC_Poisson2D_4BilinearForm_dof_map_1() : ufc::dof_map()
1464
__global_dimension = 0;
1468
virtual ~UFC_Poisson2D_4BilinearForm_dof_map_1()
1473
/// Return a string identifying the dof map
1474
virtual const char* signature() const
1476
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 4)";
1479
/// Return true iff mesh entities of topological dimension d are needed
1480
virtual bool needs_mesh_entities(unsigned int d) const
1497
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1498
virtual bool init_mesh(const ufc::mesh& m)
1500
__global_dimension = m.num_entities[0] + 3*m.num_entities[1] + 3*m.num_entities[2];
1504
/// Initialize dof map for given cell
1505
virtual void init_cell(const ufc::mesh& m,
1511
/// Finish initialization of dof map for cells
1512
virtual void init_cell_finalize()
1517
/// Return the dimension of the global finite element function space
1518
virtual unsigned int global_dimension() const
1520
return __global_dimension;
1523
/// Return the dimension of the local finite element function space
1524
virtual unsigned int local_dimension() const
1529
// Return the geometric dimension of the coordinates this dof map provides
1530
virtual unsigned int geometric_dimension() const
1535
/// Return the number of dofs on each cell facet
1536
virtual unsigned int num_facet_dofs() const
1541
/// Return the number of dofs associated with each cell entity of dimension d
1542
virtual unsigned int num_entity_dofs(unsigned int d) const
1544
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1547
/// Tabulate the local-to-global mapping of dofs on a cell
1548
virtual void tabulate_dofs(unsigned int* dofs,
1550
const ufc::cell& c) const
1552
dofs[0] = c.entity_indices[0][0];
1553
dofs[1] = c.entity_indices[0][1];
1554
dofs[2] = c.entity_indices[0][2];
1555
unsigned int offset = m.num_entities[0];
1556
dofs[3] = offset + 3*c.entity_indices[1][0];
1557
dofs[4] = offset + 3*c.entity_indices[1][0] + 1;
1558
dofs[5] = offset + 3*c.entity_indices[1][0] + 2;
1559
dofs[6] = offset + 3*c.entity_indices[1][1];
1560
dofs[7] = offset + 3*c.entity_indices[1][1] + 1;
1561
dofs[8] = offset + 3*c.entity_indices[1][1] + 2;
1562
dofs[9] = offset + 3*c.entity_indices[1][2];
1563
dofs[10] = offset + 3*c.entity_indices[1][2] + 1;
1564
dofs[11] = offset + 3*c.entity_indices[1][2] + 2;
1565
offset = offset + 3*m.num_entities[1];
1566
dofs[12] = offset + 3*c.entity_indices[2][0];
1567
dofs[13] = offset + 3*c.entity_indices[2][0] + 1;
1568
dofs[14] = offset + 3*c.entity_indices[2][0] + 2;
1571
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1572
virtual void tabulate_facet_dofs(unsigned int* dofs,
1573
unsigned int facet) const
1601
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1602
virtual void tabulate_entity_dofs(unsigned int* dofs,
1603
unsigned int d, unsigned int i) const
1605
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1608
/// Tabulate the coordinates of all dofs on a cell
1609
virtual void tabulate_coordinates(double** coordinates,
1610
const ufc::cell& c) const
1612
const double * const * x = c.coordinates;
1613
coordinates[0][0] = x[0][0];
1614
coordinates[0][1] = x[0][1];
1615
coordinates[1][0] = x[1][0];
1616
coordinates[1][1] = x[1][1];
1617
coordinates[2][0] = x[2][0];
1618
coordinates[2][1] = x[2][1];
1619
coordinates[3][0] = 0.75*x[1][0] + 0.25*x[2][0];
1620
coordinates[3][1] = 0.75*x[1][1] + 0.25*x[2][1];
1621
coordinates[4][0] = 0.5*x[1][0] + 0.5*x[2][0];
1622
coordinates[4][1] = 0.5*x[1][1] + 0.5*x[2][1];
1623
coordinates[5][0] = 0.25*x[1][0] + 0.75*x[2][0];
1624
coordinates[5][1] = 0.25*x[1][1] + 0.75*x[2][1];
1625
coordinates[6][0] = 0.75*x[0][0] + 0.25*x[2][0];
1626
coordinates[6][1] = 0.75*x[0][1] + 0.25*x[2][1];
1627
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[2][0];
1628
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[2][1];
1629
coordinates[8][0] = 0.25*x[0][0] + 0.75*x[2][0];
1630
coordinates[8][1] = 0.25*x[0][1] + 0.75*x[2][1];
1631
coordinates[9][0] = 0.75*x[0][0] + 0.25*x[1][0];
1632
coordinates[9][1] = 0.75*x[0][1] + 0.25*x[1][1];
1633
coordinates[10][0] = 0.5*x[0][0] + 0.5*x[1][0];
1634
coordinates[10][1] = 0.5*x[0][1] + 0.5*x[1][1];
1635
coordinates[11][0] = 0.25*x[0][0] + 0.75*x[1][0];
1636
coordinates[11][1] = 0.25*x[0][1] + 0.75*x[1][1];
1637
coordinates[12][0] = 0.5*x[0][0] + 0.25*x[1][0] + 0.25*x[2][0];
1638
coordinates[12][1] = 0.5*x[0][1] + 0.25*x[1][1] + 0.25*x[2][1];
1639
coordinates[13][0] = 0.25*x[0][0] + 0.5*x[1][0] + 0.25*x[2][0];
1640
coordinates[13][1] = 0.25*x[0][1] + 0.5*x[1][1] + 0.25*x[2][1];
1641
coordinates[14][0] = 0.25*x[0][0] + 0.25*x[1][0] + 0.5*x[2][0];
1642
coordinates[14][1] = 0.25*x[0][1] + 0.25*x[1][1] + 0.5*x[2][1];
1645
/// Return the number of sub dof maps (for a mixed element)
1646
virtual unsigned int num_sub_dof_maps() const
1651
/// Create a new dof_map for sub dof map i (for a mixed element)
1652
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1654
return new UFC_Poisson2D_4BilinearForm_dof_map_1();
15
/// This class defines the interface for a finite element.
17
class poisson2d_4_0_finite_element_0: public ufc::finite_element
22
poisson2d_4_0_finite_element_0() : ufc::finite_element()
28
virtual ~poisson2d_4_0_finite_element_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "FiniteElement('Lagrange', 'triangle', 4)";
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);
109
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
111
// Compute psitilde_a
112
const double psitilde_a_0 = 1;
113
const double psitilde_a_1 = x;
114
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
115
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
116
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
118
// Compute psitilde_bs
119
const double psitilde_bs_0_0 = 1;
120
const double psitilde_bs_0_1 = 1.5*y + 0.5;
121
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;
122
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;
123
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
124
const double psitilde_bs_1_0 = 1;
125
const double psitilde_bs_1_1 = 2.5*y + 1.5;
126
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;
127
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
128
const double psitilde_bs_2_0 = 1;
129
const double psitilde_bs_2_1 = 3.5*y + 2.5;
130
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
131
const double psitilde_bs_3_0 = 1;
132
const double psitilde_bs_3_1 = 4.5*y + 3.5;
133
const double psitilde_bs_4_0 = 1;
135
// Compute basisvalues
136
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
137
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
138
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
139
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
140
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
141
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
142
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
143
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
144
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
145
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
146
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
147
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
148
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
149
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
150
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
152
// Table(s) of coefficients
153
static const double coefficients0[15][15] = \
154
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366311, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
155
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233824, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.0239443566116079, -0.0185472188784818, 0.0107082418122104},
156
{0, 0, 0.0476190476190476, 0, 0, 0.038880789567987, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
157
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543896, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406278, 0.0834624849531682, -0.0535412090610519},
158
{-0.0314269680527354, 0.0109971479845642, 0.00634920634920638, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
159
{0.125707872210942, 0.0439885919382573, 0.126984126984127, 0, 0.035916534917412, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
160
{0.125707872210942, -0.131965775814772, -0.0253968253968254, 0.139104141588614, 0.0718330698348239, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.0478887132232159, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
161
{-0.0314269680527354, -0.0109971479845643, 0.00634920634920633, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
162
{0.125707872210942, -0.0439885919382572, 0.126984126984127, 0, -0.0359165349174119, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
163
{0.125707872210942, -0.0879771838765144, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
164
{-0.0314269680527354, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951817, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.023944356611608, 0, 0.0107082418122104},
165
{0.125707872210942, 0.0879771838765145, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277046, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
166
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406279, -0.0278208283177228, -0.0535412090610519},
167
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.0267706045305259, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406278, 0.0278208283177228, -0.0535412090610519},
168
{0.251415744421883, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
170
// Extract relevant coefficients
171
const double coeff0_0 = coefficients0[dof][0];
172
const double coeff0_1 = coefficients0[dof][1];
173
const double coeff0_2 = coefficients0[dof][2];
174
const double coeff0_3 = coefficients0[dof][3];
175
const double coeff0_4 = coefficients0[dof][4];
176
const double coeff0_5 = coefficients0[dof][5];
177
const double coeff0_6 = coefficients0[dof][6];
178
const double coeff0_7 = coefficients0[dof][7];
179
const double coeff0_8 = coefficients0[dof][8];
180
const double coeff0_9 = coefficients0[dof][9];
181
const double coeff0_10 = coefficients0[dof][10];
182
const double coeff0_11 = coefficients0[dof][11];
183
const double coeff0_12 = coefficients0[dof][12];
184
const double coeff0_13 = coefficients0[dof][13];
185
const double coeff0_14 = coefficients0[dof][14];
188
*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 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
191
/// Evaluate all basis functions at given point in cell
192
virtual void evaluate_basis_all(double* values,
193
const double* coordinates,
194
const ufc::cell& c) const
196
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
199
/// Evaluate order n derivatives of basis function i at given point in cell
200
virtual void evaluate_basis_derivatives(unsigned int i,
203
const double* coordinates,
204
const ufc::cell& c) const
206
// Extract vertex coordinates
207
const double * const * element_coordinates = c.coordinates;
209
// Compute Jacobian of affine map from reference cell
210
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
211
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
212
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
213
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
215
// Compute determinant of Jacobian
216
const double detJ = J_00*J_11 - J_01*J_10;
218
// Compute inverse of Jacobian
220
// Get coordinates and map to the reference (UFC) element
221
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
222
element_coordinates[0][0]*element_coordinates[2][1] +\
223
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
224
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
225
element_coordinates[1][0]*element_coordinates[0][1] -\
226
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
228
// Map coordinates to the reference square
229
if (std::abs(y - 1.0) < 1e-14)
232
x = 2.0 *x/(1.0 - y) - 1.0;
235
// Compute number of derivatives
236
unsigned int num_derivatives = 1;
238
for (unsigned int j = 0; j < n; j++)
239
num_derivatives *= 2;
242
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
243
unsigned int **combinations = new unsigned int *[num_derivatives];
245
for (unsigned int j = 0; j < num_derivatives; j++)
247
combinations[j] = new unsigned int [n];
248
for (unsigned int k = 0; k < n; k++)
249
combinations[j][k] = 0;
252
// Generate combinations of derivatives
253
for (unsigned int row = 1; row < num_derivatives; row++)
255
for (unsigned int num = 0; num < row; num++)
257
for (unsigned int col = n-1; col+1 > 0; col--)
259
if (combinations[row][col] + 1 > 1)
260
combinations[row][col] = 0;
263
combinations[row][col] += 1;
270
// Compute inverse of Jacobian
271
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
273
// Declare transformation matrix
274
// Declare pointer to two dimensional array and initialise
275
double **transform = new double *[num_derivatives];
277
for (unsigned int j = 0; j < num_derivatives; j++)
279
transform[j] = new double [num_derivatives];
280
for (unsigned int k = 0; k < num_derivatives; k++)
284
// Construct transformation matrix
285
for (unsigned int row = 0; row < num_derivatives; row++)
287
for (unsigned int col = 0; col < num_derivatives; col++)
289
for (unsigned int k = 0; k < n; k++)
290
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
295
for (unsigned int j = 0; j < 1*num_derivatives; j++)
298
// Map degree of freedom to element degree of freedom
299
const unsigned int dof = i;
302
const double scalings_y_0 = 1;
303
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
304
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
305
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
306
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
308
// Compute psitilde_a
309
const double psitilde_a_0 = 1;
310
const double psitilde_a_1 = x;
311
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
312
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
313
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
315
// Compute psitilde_bs
316
const double psitilde_bs_0_0 = 1;
317
const double psitilde_bs_0_1 = 1.5*y + 0.5;
318
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;
319
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;
320
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
321
const double psitilde_bs_1_0 = 1;
322
const double psitilde_bs_1_1 = 2.5*y + 1.5;
323
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;
324
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
325
const double psitilde_bs_2_0 = 1;
326
const double psitilde_bs_2_1 = 3.5*y + 2.5;
327
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
328
const double psitilde_bs_3_0 = 1;
329
const double psitilde_bs_3_1 = 4.5*y + 3.5;
330
const double psitilde_bs_4_0 = 1;
332
// Compute basisvalues
333
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
334
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
335
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
336
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
337
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
338
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
339
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
340
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
341
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
342
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
343
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
344
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
345
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
346
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
347
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
349
// Table(s) of coefficients
350
static const double coefficients0[15][15] = \
351
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366311, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
352
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233824, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.0239443566116079, -0.0185472188784818, 0.0107082418122104},
353
{0, 0, 0.0476190476190476, 0, 0, 0.038880789567987, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
354
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543896, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406278, 0.0834624849531682, -0.0535412090610519},
355
{-0.0314269680527354, 0.0109971479845642, 0.00634920634920638, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
356
{0.125707872210942, 0.0439885919382573, 0.126984126984127, 0, 0.035916534917412, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
357
{0.125707872210942, -0.131965775814772, -0.0253968253968254, 0.139104141588614, 0.0718330698348239, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.0478887132232159, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
358
{-0.0314269680527354, -0.0109971479845643, 0.00634920634920633, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
359
{0.125707872210942, -0.0439885919382572, 0.126984126984127, 0, -0.0359165349174119, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
360
{0.125707872210942, -0.0879771838765144, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
361
{-0.0314269680527354, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951817, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.023944356611608, 0, 0.0107082418122104},
362
{0.125707872210942, 0.0879771838765145, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277046, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
363
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406279, -0.0278208283177228, -0.0535412090610519},
364
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.0267706045305259, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406278, 0.0278208283177228, -0.0535412090610519},
365
{0.251415744421883, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
367
// Interesting (new) part
368
// Tables of derivatives of the polynomial base (transpose)
369
static const double dmats0[15][15] = \
370
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
371
{4.89897948556635, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
372
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
373
{0, 9.48683298050514, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
374
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
375
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
376
{5.29150262212918, 0, -2.99332590941916, 13.6626010212795, 0, 0.611010092660776, 0, 0, 0, 0, 0, 0, 0, 0, 0},
377
{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
378
{3.46410161513776, 0, 7.83836717690617, 0, 0, 8.40000000000001, 0, 0, 0, 0, 0, 0, 0, 0, 0},
379
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
380
{0, 10.9544511501033, 0, 0, -3.83325938999965, 0, 17.7482393492988, 0, 0.55328333517249, 0, 0, 0, 0, 0, 0},
381
{4.73286382647969, 0, 3.3466401061363, 4.36435780471985, 0, -5.07468037933238, 0, 17.0084012854152, 0, 1.52127765851133, 0, 0, 0, 0, 0},
382
{0, 2.44948974278318, 0, 0, 9.14285714285715, 0, 0, 0, 14.8461497791618, 0, 0, 0, 0, 0, 0},
383
{3.09838667696593, 0, 7.66811580507233, 0, 0, 10.733126291999, 0, 0, 0, 9.2951600308978, 0, 0, 0, 0, 0},
384
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
386
static const double dmats1[15][15] = \
387
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
388
{2.44948974278318, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
389
{4.24264068711928, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
390
{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
391
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
392
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
393
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330388, 0, 0, 0, 0, 0, 0, 0, 0, 0},
394
{2.23606797749978, 2.19089023002066, 2.5298221281347, 8.08290376865476, 6.26099033699942, -1.80739222823014, 0, 0, 0, 0, 0, 0, 0, 0, 0},
395
{1.73205080756888, -5.09116882454315, 3.91918358845308, 0, 9.69948452238572, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
396
{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0, 0, 0, 0, 0, 0},
397
{2.68328157299975, 5.47722557505166, -1.89736659610103, 7.4230748895809, -1.91662969499982, 0.663940002206986, 8.87411967464942, -1.07142857142857, 0.276641667586244, -0.0958314847499919, 0, 0, 0, 0, 0},
398
{2.36643191323985, 2.89827534923789, 1.67332005306815, 2.18217890235992, 5.74704893215392, -2.53734018966619, 10.0623058987491, 8.50420064270761, -2.1957751641342, 0.760638829255664, 0, 0, 0, 0, 0},
399
{2, 1.22474487139159, 3.53553390593274, -7.37711113563317, 4.57142857142858, 1.64957219768464, 0, 11.4997781699989, 7.4230748895809, -2.57142857142857, 0, 0, 0, 0, 0},
400
{1.54919333848297, 6.6407830863536, 3.83405790253617, 0, -6.19677335393188, 5.3665631459995, 0, 0, 13.4164078649987, 4.6475800154489, 0, 0, 0, 0, 0},
401
{-3.57770876399967, 0, 8.85437744847147, 0, 0, -3.09838667696594, 0, 0, 0, 16.0996894379985, 0, 0, 0, 0, 0}};
403
// Compute reference derivatives
404
// Declare pointer to array of derivatives on FIAT element
405
double *derivatives = new double [num_derivatives];
407
// Declare coefficients
418
double coeff0_10 = 0;
419
double coeff0_11 = 0;
420
double coeff0_12 = 0;
421
double coeff0_13 = 0;
422
double coeff0_14 = 0;
424
// Declare new coefficients
425
double new_coeff0_0 = 0;
426
double new_coeff0_1 = 0;
427
double new_coeff0_2 = 0;
428
double new_coeff0_3 = 0;
429
double new_coeff0_4 = 0;
430
double new_coeff0_5 = 0;
431
double new_coeff0_6 = 0;
432
double new_coeff0_7 = 0;
433
double new_coeff0_8 = 0;
434
double new_coeff0_9 = 0;
435
double new_coeff0_10 = 0;
436
double new_coeff0_11 = 0;
437
double new_coeff0_12 = 0;
438
double new_coeff0_13 = 0;
439
double new_coeff0_14 = 0;
441
// Loop possible derivatives
442
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
444
// Get values from coefficients array
445
new_coeff0_0 = coefficients0[dof][0];
446
new_coeff0_1 = coefficients0[dof][1];
447
new_coeff0_2 = coefficients0[dof][2];
448
new_coeff0_3 = coefficients0[dof][3];
449
new_coeff0_4 = coefficients0[dof][4];
450
new_coeff0_5 = coefficients0[dof][5];
451
new_coeff0_6 = coefficients0[dof][6];
452
new_coeff0_7 = coefficients0[dof][7];
453
new_coeff0_8 = coefficients0[dof][8];
454
new_coeff0_9 = coefficients0[dof][9];
455
new_coeff0_10 = coefficients0[dof][10];
456
new_coeff0_11 = coefficients0[dof][11];
457
new_coeff0_12 = coefficients0[dof][12];
458
new_coeff0_13 = coefficients0[dof][13];
459
new_coeff0_14 = coefficients0[dof][14];
461
// Loop derivative order
462
for (unsigned int j = 0; j < n; j++)
464
// Update old coefficients
465
coeff0_0 = new_coeff0_0;
466
coeff0_1 = new_coeff0_1;
467
coeff0_2 = new_coeff0_2;
468
coeff0_3 = new_coeff0_3;
469
coeff0_4 = new_coeff0_4;
470
coeff0_5 = new_coeff0_5;
471
coeff0_6 = new_coeff0_6;
472
coeff0_7 = new_coeff0_7;
473
coeff0_8 = new_coeff0_8;
474
coeff0_9 = new_coeff0_9;
475
coeff0_10 = new_coeff0_10;
476
coeff0_11 = new_coeff0_11;
477
coeff0_12 = new_coeff0_12;
478
coeff0_13 = new_coeff0_13;
479
coeff0_14 = new_coeff0_14;
481
if(combinations[deriv_num][j] == 0)
483
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] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
484
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] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
485
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] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
486
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] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
487
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] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
488
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] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
489
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] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
490
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] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
491
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] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
492
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] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
493
new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
494
new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
495
new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
496
new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
497
new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
499
if(combinations[deriv_num][j] == 1)
501
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] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
502
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] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
503
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] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
504
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] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
505
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] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
506
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] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
507
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] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
508
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] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
509
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] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
510
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] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
511
new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
512
new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
513
new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
514
new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
515
new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
519
// Compute derivatives on reference element as dot product of coefficients and basisvalues
520
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 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
523
// Transform derivatives back to physical element
524
for (unsigned int row = 0; row < num_derivatives; row++)
526
for (unsigned int col = 0; col < num_derivatives; col++)
528
values[row] += transform[row][col]*derivatives[col];
531
// Delete pointer to array of derivatives on FIAT element
532
delete [] derivatives;
534
// Delete pointer to array of combinations of derivatives and transform
535
for (unsigned int row = 0; row < num_derivatives; row++)
537
delete [] combinations[row];
538
delete [] transform[row];
541
delete [] combinations;
545
/// Evaluate order n derivatives of all basis functions at given point in cell
546
virtual void evaluate_basis_derivatives_all(unsigned int n,
548
const double* coordinates,
549
const ufc::cell& c) const
551
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
554
/// Evaluate linear functional for dof i on the function f
555
virtual double evaluate_dof(unsigned int i,
556
const ufc::function& f,
557
const ufc::cell& c) const
559
// The reference points, direction and weights:
560
static const double X[15][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}, {{0.25, 0.25}}, {{0.5, 0.25}}, {{0.25, 0.5}}};
561
static const double W[15][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
562
static const double D[15][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
564
const double * const * x = c.coordinates;
566
// Iterate over the points:
567
// Evaluate basis functions for affine mapping
568
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
569
const double w1 = X[i][0][0];
570
const double w2 = X[i][0][1];
572
// Compute affine mapping y = F(X)
574
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
575
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
577
// Evaluate function at physical points
579
f.evaluate(values, y, c);
581
// Map function values using appropriate mapping
582
// Affine map: Do nothing
584
// Note that we do not map the weights (yet).
586
// Take directional components
587
for(int k = 0; k < 1; k++)
588
result += values[k]*D[i][0][k];
589
// Multiply by weights
595
/// Evaluate linear functionals for all dofs on the function f
596
virtual void evaluate_dofs(double* values,
597
const ufc::function& f,
598
const ufc::cell& c) const
600
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
603
/// Interpolate vertex values from dof values
604
virtual void interpolate_vertex_values(double* vertex_values,
605
const double* dof_values,
606
const ufc::cell& c) const
608
// Evaluate at vertices and use affine mapping
609
vertex_values[0] = dof_values[0];
610
vertex_values[1] = dof_values[1];
611
vertex_values[2] = dof_values[2];
614
/// Return the number of sub elements (for a mixed element)
615
virtual unsigned int num_sub_elements() const
620
/// Create a new finite element for sub element i (for a mixed element)
621
virtual ufc::finite_element* create_sub_element(unsigned int i) const
623
return new poisson2d_4_0_finite_element_0();
628
/// This class defines the interface for a finite element.
630
class poisson2d_4_0_finite_element_1: public ufc::finite_element
635
poisson2d_4_0_finite_element_1() : ufc::finite_element()
641
virtual ~poisson2d_4_0_finite_element_1()
646
/// Return a string identifying the finite element
647
virtual const char* signature() const
649
return "FiniteElement('Lagrange', 'triangle', 4)";
652
/// Return the cell shape
653
virtual ufc::shape cell_shape() const
655
return ufc::triangle;
658
/// Return the dimension of the finite element function space
659
virtual unsigned int space_dimension() const
664
/// Return the rank of the value space
665
virtual unsigned int value_rank() const
670
/// Return the dimension of the value space for axis i
671
virtual unsigned int value_dimension(unsigned int i) const
676
/// Evaluate basis function i at given point in cell
677
virtual void evaluate_basis(unsigned int i,
679
const double* coordinates,
680
const ufc::cell& c) const
682
// Extract vertex coordinates
683
const double * const * element_coordinates = c.coordinates;
685
// Compute Jacobian of affine map from reference cell
686
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
687
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
688
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
689
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
691
// Compute determinant of Jacobian
692
const double detJ = J_00*J_11 - J_01*J_10;
694
// Compute inverse of Jacobian
696
// Get coordinates and map to the reference (UFC) element
697
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
698
element_coordinates[0][0]*element_coordinates[2][1] +\
699
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
700
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
701
element_coordinates[1][0]*element_coordinates[0][1] -\
702
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
704
// Map coordinates to the reference square
705
if (std::abs(y - 1.0) < 1e-14)
708
x = 2.0 *x/(1.0 - y) - 1.0;
714
// Map degree of freedom to element degree of freedom
715
const unsigned int dof = i;
718
const double scalings_y_0 = 1;
719
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
720
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
721
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
722
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
724
// Compute psitilde_a
725
const double psitilde_a_0 = 1;
726
const double psitilde_a_1 = x;
727
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
728
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
729
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
731
// Compute psitilde_bs
732
const double psitilde_bs_0_0 = 1;
733
const double psitilde_bs_0_1 = 1.5*y + 0.5;
734
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;
735
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;
736
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
737
const double psitilde_bs_1_0 = 1;
738
const double psitilde_bs_1_1 = 2.5*y + 1.5;
739
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;
740
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
741
const double psitilde_bs_2_0 = 1;
742
const double psitilde_bs_2_1 = 3.5*y + 2.5;
743
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
744
const double psitilde_bs_3_0 = 1;
745
const double psitilde_bs_3_1 = 4.5*y + 3.5;
746
const double psitilde_bs_4_0 = 1;
748
// Compute basisvalues
749
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
750
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
751
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
752
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
753
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
754
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
755
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
756
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
757
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
758
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
759
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
760
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
761
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
762
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
763
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
765
// Table(s) of coefficients
766
static const double coefficients0[15][15] = \
767
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366311, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
768
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233824, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.0239443566116079, -0.0185472188784818, 0.0107082418122104},
769
{0, 0, 0.0476190476190476, 0, 0, 0.038880789567987, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
770
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543896, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406278, 0.0834624849531682, -0.0535412090610519},
771
{-0.0314269680527354, 0.0109971479845642, 0.00634920634920638, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
772
{0.125707872210942, 0.0439885919382573, 0.126984126984127, 0, 0.035916534917412, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
773
{0.125707872210942, -0.131965775814772, -0.0253968253968254, 0.139104141588614, 0.0718330698348239, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.0478887132232159, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
774
{-0.0314269680527354, -0.0109971479845643, 0.00634920634920633, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
775
{0.125707872210942, -0.0439885919382572, 0.126984126984127, 0, -0.0359165349174119, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
776
{0.125707872210942, -0.0879771838765144, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
777
{-0.0314269680527354, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951817, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.023944356611608, 0, 0.0107082418122104},
778
{0.125707872210942, 0.0879771838765145, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277046, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
779
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406279, -0.0278208283177228, -0.0535412090610519},
780
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.0267706045305259, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406278, 0.0278208283177228, -0.0535412090610519},
781
{0.251415744421883, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
783
// Extract relevant coefficients
784
const double coeff0_0 = coefficients0[dof][0];
785
const double coeff0_1 = coefficients0[dof][1];
786
const double coeff0_2 = coefficients0[dof][2];
787
const double coeff0_3 = coefficients0[dof][3];
788
const double coeff0_4 = coefficients0[dof][4];
789
const double coeff0_5 = coefficients0[dof][5];
790
const double coeff0_6 = coefficients0[dof][6];
791
const double coeff0_7 = coefficients0[dof][7];
792
const double coeff0_8 = coefficients0[dof][8];
793
const double coeff0_9 = coefficients0[dof][9];
794
const double coeff0_10 = coefficients0[dof][10];
795
const double coeff0_11 = coefficients0[dof][11];
796
const double coeff0_12 = coefficients0[dof][12];
797
const double coeff0_13 = coefficients0[dof][13];
798
const double coeff0_14 = coefficients0[dof][14];
801
*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 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
804
/// Evaluate all basis functions at given point in cell
805
virtual void evaluate_basis_all(double* values,
806
const double* coordinates,
807
const ufc::cell& c) const
809
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
812
/// Evaluate order n derivatives of basis function i at given point in cell
813
virtual void evaluate_basis_derivatives(unsigned int i,
816
const double* coordinates,
817
const ufc::cell& c) const
819
// Extract vertex coordinates
820
const double * const * element_coordinates = c.coordinates;
822
// Compute Jacobian of affine map from reference cell
823
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
824
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
825
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
826
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
828
// Compute determinant of Jacobian
829
const double detJ = J_00*J_11 - J_01*J_10;
831
// Compute inverse of Jacobian
833
// Get coordinates and map to the reference (UFC) element
834
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
835
element_coordinates[0][0]*element_coordinates[2][1] +\
836
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
837
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
838
element_coordinates[1][0]*element_coordinates[0][1] -\
839
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
841
// Map coordinates to the reference square
842
if (std::abs(y - 1.0) < 1e-14)
845
x = 2.0 *x/(1.0 - y) - 1.0;
848
// Compute number of derivatives
849
unsigned int num_derivatives = 1;
851
for (unsigned int j = 0; j < n; j++)
852
num_derivatives *= 2;
855
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
856
unsigned int **combinations = new unsigned int *[num_derivatives];
858
for (unsigned int j = 0; j < num_derivatives; j++)
860
combinations[j] = new unsigned int [n];
861
for (unsigned int k = 0; k < n; k++)
862
combinations[j][k] = 0;
865
// Generate combinations of derivatives
866
for (unsigned int row = 1; row < num_derivatives; row++)
868
for (unsigned int num = 0; num < row; num++)
870
for (unsigned int col = n-1; col+1 > 0; col--)
872
if (combinations[row][col] + 1 > 1)
873
combinations[row][col] = 0;
876
combinations[row][col] += 1;
883
// Compute inverse of Jacobian
884
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
886
// Declare transformation matrix
887
// Declare pointer to two dimensional array and initialise
888
double **transform = new double *[num_derivatives];
890
for (unsigned int j = 0; j < num_derivatives; j++)
892
transform[j] = new double [num_derivatives];
893
for (unsigned int k = 0; k < num_derivatives; k++)
897
// Construct transformation matrix
898
for (unsigned int row = 0; row < num_derivatives; row++)
900
for (unsigned int col = 0; col < num_derivatives; col++)
902
for (unsigned int k = 0; k < n; k++)
903
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
908
for (unsigned int j = 0; j < 1*num_derivatives; j++)
911
// Map degree of freedom to element degree of freedom
912
const unsigned int dof = i;
915
const double scalings_y_0 = 1;
916
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
917
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
918
const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
919
const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
921
// Compute psitilde_a
922
const double psitilde_a_0 = 1;
923
const double psitilde_a_1 = x;
924
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
925
const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
926
const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
928
// Compute psitilde_bs
929
const double psitilde_bs_0_0 = 1;
930
const double psitilde_bs_0_1 = 1.5*y + 0.5;
931
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;
932
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;
933
const double psitilde_bs_0_4 = 0.0285714285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571428571*psitilde_bs_0_2;
934
const double psitilde_bs_1_0 = 1;
935
const double psitilde_bs_1_1 = 2.5*y + 1.5;
936
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;
937
const double psitilde_bs_1_3 = 0.285714285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714285714*psitilde_bs_1_1;
938
const double psitilde_bs_2_0 = 1;
939
const double psitilde_bs_2_1 = 3.5*y + 2.5;
940
const double psitilde_bs_2_2 = 1.02040816326531*psitilde_bs_2_1 + 2.57142857142857*y*psitilde_bs_2_1 - 0.551020408163265*psitilde_bs_2_0;
941
const double psitilde_bs_3_0 = 1;
942
const double psitilde_bs_3_1 = 4.5*y + 3.5;
943
const double psitilde_bs_4_0 = 1;
945
// Compute basisvalues
946
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
947
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
948
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
949
const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
950
const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
951
const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
952
const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
953
const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
954
const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
955
const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
956
const double basisvalue10 = 4.74341649025257*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
957
const double basisvalue11 = 4.18330013267038*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
958
const double basisvalue12 = 3.53553390593274*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
959
const double basisvalue13 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
960
const double basisvalue14 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
962
// Table(s) of coefficients
963
static const double coefficients0[15][15] = \
964
{{0, -0.0412393049421161, -0.0238095238095238, 0.0289800294976279, 0.0224478343233825, 0.012960263189329, -0.0395942580610999, -0.0334632556631574, -0.025920526378658, -0.014965222882255, 0.0321247254366311, 0.0283313448138523, 0.0239443566116079, 0.0185472188784818, 0.0107082418122104},
965
{0, 0.0412393049421161, -0.0238095238095238, 0.0289800294976278, -0.0224478343233824, 0.012960263189329, 0.0395942580610999, -0.0334632556631574, 0.025920526378658, -0.014965222882255, 0.0321247254366312, -0.0283313448138523, 0.0239443566116079, -0.0185472188784818, 0.0107082418122104},
966
{0, 0, 0.0476190476190476, 0, 0, 0.038880789567987, 0, 0, 0, 0.0598608915290199, 0, 0, 0, 0, 0.0535412090610519},
967
{0.125707872210942, 0.131965775814772, -0.0253968253968254, 0.139104141588614, -0.0718330698348239, 0.0311046316543896, 0.0633508128977598, 0.0267706045305259, -0.0622092633087792, 0.0478887132232159, 0, 0.0566626896277046, -0.0838052481406278, 0.0834624849531682, -0.0535412090610519},
968
{-0.0314269680527354, 0.0109971479845642, 0.00634920634920638, 0, 0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, -0.139104141588614, 0.107082418122104},
969
{0.125707872210942, 0.0439885919382573, 0.126984126984127, 0, 0.035916534917412, 0.155523158271948, 0, 0, 0.103682105514632, -0.011972178305804, 0, 0, 0, 0.0927360943924091, -0.107082418122104},
970
{0.125707872210942, -0.131965775814772, -0.0253968253968254, 0.139104141588614, 0.0718330698348239, 0.0311046316543895, -0.0633508128977598, 0.0267706045305259, 0.0622092633087791, 0.0478887132232159, 0, -0.0566626896277046, -0.0838052481406278, -0.0834624849531682, -0.0535412090610519},
971
{-0.0314269680527354, -0.0109971479845643, 0.00634920634920633, 0, -0.188561808316413, -0.163299316185545, 0, 0.0936971158568409, 0, -0.0419026240703139, 0, 0, 0.0838052481406278, 0.139104141588614, 0.107082418122104},
972
{0.125707872210942, -0.0439885919382572, 0.126984126984127, 0, -0.0359165349174119, 0.155523158271948, 0, 0, -0.103682105514632, -0.011972178305804, 0, 0, 0, -0.0927360943924091, -0.107082418122104},
973
{0.125707872210942, -0.0879771838765144, -0.101587301587302, 0.0927360943924091, 0.107749604752236, 0.0725774738602423, 0.0791885161221998, -0.013385302265263, -0.0518410527573159, -0.0419026240703139, -0.128498901746525, -0.0566626896277046, -0.011972178305804, 0.00927360943924092, 0.0107082418122104},
974
{-0.0314269680527354, 0, -0.0126984126984127, -0.243432247780074, 0, 0.0544331053951817, 0, 0.0936971158568408, 0, -0.0419026240703139, 0.192748352619787, 0, -0.023944356611608, 0, 0.0107082418122104},
975
{0.125707872210942, 0.0879771838765145, -0.101587301587302, 0.0927360943924091, -0.107749604752236, 0.0725774738602423, -0.0791885161221998, -0.013385302265263, 0.0518410527573159, -0.0419026240703139, -0.128498901746525, 0.0566626896277046, -0.011972178305804, -0.0092736094392409, 0.0107082418122104},
976
{0.251415744421884, -0.351908735506058, -0.203174603174603, -0.139104141588614, -0.107749604752236, -0.0622092633087791, 0.19005243869328, -0.0267706045305259, 0.124418526617558, 0.155638317975452, 0, 0.169988068883114, 0.0838052481406279, -0.0278208283177228, -0.0535412090610519},
977
{0.251415744421884, 0.351908735506058, -0.203174603174603, -0.139104141588614, 0.107749604752236, -0.0622092633087791, -0.19005243869328, -0.0267706045305259, -0.124418526617558, 0.155638317975452, 0, -0.169988068883114, 0.0838052481406278, 0.0278208283177228, -0.0535412090610519},
978
{0.251415744421883, 0, 0.406349206349206, 0, 0, -0.186627789926337, 0, -0.187394231713682, 0, -0.203527031198668, 0, 0, -0.167610496281256, 0, 0.107082418122104}};
980
// Interesting (new) part
981
// Tables of derivatives of the polynomial base (transpose)
982
static const double dmats0[15][15] = \
983
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
984
{4.89897948556635, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
985
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
986
{0, 9.48683298050514, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
987
{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
988
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
989
{5.29150262212918, 0, -2.99332590941916, 13.6626010212795, 0, 0.611010092660776, 0, 0, 0, 0, 0, 0, 0, 0, 0},
990
{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
991
{3.46410161513776, 0, 7.83836717690617, 0, 0, 8.40000000000001, 0, 0, 0, 0, 0, 0, 0, 0, 0},
992
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
993
{0, 10.9544511501033, 0, 0, -3.83325938999965, 0, 17.7482393492988, 0, 0.55328333517249, 0, 0, 0, 0, 0, 0},
994
{4.73286382647969, 0, 3.3466401061363, 4.36435780471985, 0, -5.07468037933238, 0, 17.0084012854152, 0, 1.52127765851133, 0, 0, 0, 0, 0},
995
{0, 2.44948974278318, 0, 0, 9.14285714285715, 0, 0, 0, 14.8461497791618, 0, 0, 0, 0, 0, 0},
996
{3.09838667696593, 0, 7.66811580507233, 0, 0, 10.733126291999, 0, 0, 0, 9.2951600308978, 0, 0, 0, 0, 0},
997
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
999
static const double dmats1[15][15] = \
1000
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1001
{2.44948974278318, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1002
{4.24264068711928, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1003
{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1004
{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1005
{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1006
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330388, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1007
{2.23606797749978, 2.19089023002066, 2.5298221281347, 8.08290376865476, 6.26099033699942, -1.80739222823014, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1008
{1.73205080756888, -5.09116882454315, 3.91918358845308, 0, 9.69948452238572, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1009
{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1010
{2.68328157299975, 5.47722557505166, -1.89736659610103, 7.4230748895809, -1.91662969499982, 0.663940002206986, 8.87411967464942, -1.07142857142857, 0.276641667586244, -0.0958314847499919, 0, 0, 0, 0, 0},
1011
{2.36643191323985, 2.89827534923789, 1.67332005306815, 2.18217890235992, 5.74704893215392, -2.53734018966619, 10.0623058987491, 8.50420064270761, -2.1957751641342, 0.760638829255664, 0, 0, 0, 0, 0},
1012
{2, 1.22474487139159, 3.53553390593274, -7.37711113563317, 4.57142857142858, 1.64957219768464, 0, 11.4997781699989, 7.4230748895809, -2.57142857142857, 0, 0, 0, 0, 0},
1013
{1.54919333848297, 6.6407830863536, 3.83405790253617, 0, -6.19677335393188, 5.3665631459995, 0, 0, 13.4164078649987, 4.6475800154489, 0, 0, 0, 0, 0},
1014
{-3.57770876399967, 0, 8.85437744847147, 0, 0, -3.09838667696594, 0, 0, 0, 16.0996894379985, 0, 0, 0, 0, 0}};
1016
// Compute reference derivatives
1017
// Declare pointer to array of derivatives on FIAT element
1018
double *derivatives = new double [num_derivatives];
1020
// Declare coefficients
1021
double coeff0_0 = 0;
1022
double coeff0_1 = 0;
1023
double coeff0_2 = 0;
1024
double coeff0_3 = 0;
1025
double coeff0_4 = 0;
1026
double coeff0_5 = 0;
1027
double coeff0_6 = 0;
1028
double coeff0_7 = 0;
1029
double coeff0_8 = 0;
1030
double coeff0_9 = 0;
1031
double coeff0_10 = 0;
1032
double coeff0_11 = 0;
1033
double coeff0_12 = 0;
1034
double coeff0_13 = 0;
1035
double coeff0_14 = 0;
1037
// Declare new coefficients
1038
double new_coeff0_0 = 0;
1039
double new_coeff0_1 = 0;
1040
double new_coeff0_2 = 0;
1041
double new_coeff0_3 = 0;
1042
double new_coeff0_4 = 0;
1043
double new_coeff0_5 = 0;
1044
double new_coeff0_6 = 0;
1045
double new_coeff0_7 = 0;
1046
double new_coeff0_8 = 0;
1047
double new_coeff0_9 = 0;
1048
double new_coeff0_10 = 0;
1049
double new_coeff0_11 = 0;
1050
double new_coeff0_12 = 0;
1051
double new_coeff0_13 = 0;
1052
double new_coeff0_14 = 0;
1054
// Loop possible derivatives
1055
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1057
// Get values from coefficients array
1058
new_coeff0_0 = coefficients0[dof][0];
1059
new_coeff0_1 = coefficients0[dof][1];
1060
new_coeff0_2 = coefficients0[dof][2];
1061
new_coeff0_3 = coefficients0[dof][3];
1062
new_coeff0_4 = coefficients0[dof][4];
1063
new_coeff0_5 = coefficients0[dof][5];
1064
new_coeff0_6 = coefficients0[dof][6];
1065
new_coeff0_7 = coefficients0[dof][7];
1066
new_coeff0_8 = coefficients0[dof][8];
1067
new_coeff0_9 = coefficients0[dof][9];
1068
new_coeff0_10 = coefficients0[dof][10];
1069
new_coeff0_11 = coefficients0[dof][11];
1070
new_coeff0_12 = coefficients0[dof][12];
1071
new_coeff0_13 = coefficients0[dof][13];
1072
new_coeff0_14 = coefficients0[dof][14];
1074
// Loop derivative order
1075
for (unsigned int j = 0; j < n; j++)
1077
// Update old coefficients
1078
coeff0_0 = new_coeff0_0;
1079
coeff0_1 = new_coeff0_1;
1080
coeff0_2 = new_coeff0_2;
1081
coeff0_3 = new_coeff0_3;
1082
coeff0_4 = new_coeff0_4;
1083
coeff0_5 = new_coeff0_5;
1084
coeff0_6 = new_coeff0_6;
1085
coeff0_7 = new_coeff0_7;
1086
coeff0_8 = new_coeff0_8;
1087
coeff0_9 = new_coeff0_9;
1088
coeff0_10 = new_coeff0_10;
1089
coeff0_11 = new_coeff0_11;
1090
coeff0_12 = new_coeff0_12;
1091
coeff0_13 = new_coeff0_13;
1092
coeff0_14 = new_coeff0_14;
1094
if(combinations[deriv_num][j] == 0)
1096
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] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
1097
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] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
1098
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] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
1099
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] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
1100
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] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
1101
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] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
1102
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] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
1103
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] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
1104
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] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
1105
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] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
1106
new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
1107
new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
1108
new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
1109
new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
1110
new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
1112
if(combinations[deriv_num][j] == 1)
1114
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] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
1115
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] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
1116
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] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
1117
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] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
1118
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] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
1119
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] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
1120
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] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
1121
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] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
1122
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] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
1123
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] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
1124
new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
1125
new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
1126
new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
1127
new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
1128
new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
1132
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1133
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 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
1136
// Transform derivatives back to physical element
1137
for (unsigned int row = 0; row < num_derivatives; row++)
1139
for (unsigned int col = 0; col < num_derivatives; col++)
1141
values[row] += transform[row][col]*derivatives[col];
1144
// Delete pointer to array of derivatives on FIAT element
1145
delete [] derivatives;
1147
// Delete pointer to array of combinations of derivatives and transform
1148
for (unsigned int row = 0; row < num_derivatives; row++)
1150
delete [] combinations[row];
1151
delete [] transform[row];
1154
delete [] combinations;
1155
delete [] transform;
1158
/// Evaluate order n derivatives of all basis functions at given point in cell
1159
virtual void evaluate_basis_derivatives_all(unsigned int n,
1161
const double* coordinates,
1162
const ufc::cell& c) const
1164
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1167
/// Evaluate linear functional for dof i on the function f
1168
virtual double evaluate_dof(unsigned int i,
1169
const ufc::function& f,
1170
const ufc::cell& c) const
1172
// The reference points, direction and weights:
1173
static const double X[15][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}, {{0.25, 0.25}}, {{0.5, 0.25}}, {{0.25, 0.5}}};
1174
static const double W[15][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
1175
static const double D[15][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
1177
const double * const * x = c.coordinates;
1178
double result = 0.0;
1179
// Iterate over the points:
1180
// Evaluate basis functions for affine mapping
1181
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1182
const double w1 = X[i][0][0];
1183
const double w2 = X[i][0][1];
1185
// Compute affine mapping y = F(X)
1187
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1188
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1190
// Evaluate function at physical points
1192
f.evaluate(values, y, c);
1194
// Map function values using appropriate mapping
1195
// Affine map: Do nothing
1197
// Note that we do not map the weights (yet).
1199
// Take directional components
1200
for(int k = 0; k < 1; k++)
1201
result += values[k]*D[i][0][k];
1202
// Multiply by weights
1208
/// Evaluate linear functionals for all dofs on the function f
1209
virtual void evaluate_dofs(double* values,
1210
const ufc::function& f,
1211
const ufc::cell& c) const
1213
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1216
/// Interpolate vertex values from dof values
1217
virtual void interpolate_vertex_values(double* vertex_values,
1218
const double* dof_values,
1219
const ufc::cell& c) const
1221
// Evaluate at vertices and use affine mapping
1222
vertex_values[0] = dof_values[0];
1223
vertex_values[1] = dof_values[1];
1224
vertex_values[2] = dof_values[2];
1227
/// Return the number of sub elements (for a mixed element)
1228
virtual unsigned int num_sub_elements() const
1233
/// Create a new finite element for sub element i (for a mixed element)
1234
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1236
return new poisson2d_4_0_finite_element_1();
1241
/// This class defines the interface for a local-to-global mapping of
1242
/// degrees of freedom (dofs).
1244
class poisson2d_4_0_dof_map_0: public ufc::dof_map
1248
unsigned int __global_dimension;
1253
poisson2d_4_0_dof_map_0() : ufc::dof_map()
1255
__global_dimension = 0;
1259
virtual ~poisson2d_4_0_dof_map_0()
1264
/// Return a string identifying the dof map
1265
virtual const char* signature() const
1267
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 4)";
1270
/// Return true iff mesh entities of topological dimension d are needed
1271
virtual bool needs_mesh_entities(unsigned int d) const
1288
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1289
virtual bool init_mesh(const ufc::mesh& m)
1291
__global_dimension = m.num_entities[0] + 3*m.num_entities[1] + 3*m.num_entities[2];
1295
/// Initialize dof map for given cell
1296
virtual void init_cell(const ufc::mesh& m,
1302
/// Finish initialization of dof map for cells
1303
virtual void init_cell_finalize()
1308
/// Return the dimension of the global finite element function space
1309
virtual unsigned int global_dimension() const
1311
return __global_dimension;
1314
/// Return the dimension of the local finite element function space for a cell
1315
virtual unsigned int local_dimension(const ufc::cell& c) const
1320
/// Return the maximum dimension of the local finite element function space
1321
virtual unsigned int max_local_dimension() const
1326
// Return the geometric dimension of the coordinates this dof map provides
1327
virtual unsigned int geometric_dimension() const
1332
/// Return the number of dofs on each cell facet
1333
virtual unsigned int num_facet_dofs() const
1338
/// Return the number of dofs associated with each cell entity of dimension d
1339
virtual unsigned int num_entity_dofs(unsigned int d) const
1341
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1344
/// Tabulate the local-to-global mapping of dofs on a cell
1345
virtual void tabulate_dofs(unsigned int* dofs,
1347
const ufc::cell& c) const
1349
dofs[0] = c.entity_indices[0][0];
1350
dofs[1] = c.entity_indices[0][1];
1351
dofs[2] = c.entity_indices[0][2];
1352
unsigned int offset = m.num_entities[0];
1353
dofs[3] = offset + 3*c.entity_indices[1][0];
1354
dofs[4] = offset + 3*c.entity_indices[1][0] + 1;
1355
dofs[5] = offset + 3*c.entity_indices[1][0] + 2;
1356
dofs[6] = offset + 3*c.entity_indices[1][1];
1357
dofs[7] = offset + 3*c.entity_indices[1][1] + 1;
1358
dofs[8] = offset + 3*c.entity_indices[1][1] + 2;
1359
dofs[9] = offset + 3*c.entity_indices[1][2];
1360
dofs[10] = offset + 3*c.entity_indices[1][2] + 1;
1361
dofs[11] = offset + 3*c.entity_indices[1][2] + 2;
1362
offset = offset + 3*m.num_entities[1];
1363
dofs[12] = offset + 3*c.entity_indices[2][0];
1364
dofs[13] = offset + 3*c.entity_indices[2][0] + 1;
1365
dofs[14] = offset + 3*c.entity_indices[2][0] + 2;
1368
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1369
virtual void tabulate_facet_dofs(unsigned int* dofs,
1370
unsigned int facet) const
1398
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1399
virtual void tabulate_entity_dofs(unsigned int* dofs,
1400
unsigned int d, unsigned int i) const
1402
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1405
/// Tabulate the coordinates of all dofs on a cell
1406
virtual void tabulate_coordinates(double** coordinates,
1407
const ufc::cell& c) const
1409
const double * const * x = c.coordinates;
1410
coordinates[0][0] = x[0][0];
1411
coordinates[0][1] = x[0][1];
1412
coordinates[1][0] = x[1][0];
1413
coordinates[1][1] = x[1][1];
1414
coordinates[2][0] = x[2][0];
1415
coordinates[2][1] = x[2][1];
1416
coordinates[3][0] = 0.75*x[1][0] + 0.25*x[2][0];
1417
coordinates[3][1] = 0.75*x[1][1] + 0.25*x[2][1];
1418
coordinates[4][0] = 0.5*x[1][0] + 0.5*x[2][0];
1419
coordinates[4][1] = 0.5*x[1][1] + 0.5*x[2][1];
1420
coordinates[5][0] = 0.25*x[1][0] + 0.75*x[2][0];
1421
coordinates[5][1] = 0.25*x[1][1] + 0.75*x[2][1];
1422
coordinates[6][0] = 0.75*x[0][0] + 0.25*x[2][0];
1423
coordinates[6][1] = 0.75*x[0][1] + 0.25*x[2][1];
1424
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[2][0];
1425
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[2][1];
1426
coordinates[8][0] = 0.25*x[0][0] + 0.75*x[2][0];
1427
coordinates[8][1] = 0.25*x[0][1] + 0.75*x[2][1];
1428
coordinates[9][0] = 0.75*x[0][0] + 0.25*x[1][0];
1429
coordinates[9][1] = 0.75*x[0][1] + 0.25*x[1][1];
1430
coordinates[10][0] = 0.5*x[0][0] + 0.5*x[1][0];
1431
coordinates[10][1] = 0.5*x[0][1] + 0.5*x[1][1];
1432
coordinates[11][0] = 0.25*x[0][0] + 0.75*x[1][0];
1433
coordinates[11][1] = 0.25*x[0][1] + 0.75*x[1][1];
1434
coordinates[12][0] = 0.5*x[0][0] + 0.25*x[1][0] + 0.25*x[2][0];
1435
coordinates[12][1] = 0.5*x[0][1] + 0.25*x[1][1] + 0.25*x[2][1];
1436
coordinates[13][0] = 0.25*x[0][0] + 0.5*x[1][0] + 0.25*x[2][0];
1437
coordinates[13][1] = 0.25*x[0][1] + 0.5*x[1][1] + 0.25*x[2][1];
1438
coordinates[14][0] = 0.25*x[0][0] + 0.25*x[1][0] + 0.5*x[2][0];
1439
coordinates[14][1] = 0.25*x[0][1] + 0.25*x[1][1] + 0.5*x[2][1];
1442
/// Return the number of sub dof maps (for a mixed element)
1443
virtual unsigned int num_sub_dof_maps() const
1448
/// Create a new dof_map for sub dof map i (for a mixed element)
1449
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1451
return new poisson2d_4_0_dof_map_0();
1456
/// This class defines the interface for a local-to-global mapping of
1457
/// degrees of freedom (dofs).
1459
class poisson2d_4_0_dof_map_1: public ufc::dof_map
1463
unsigned int __global_dimension;
1468
poisson2d_4_0_dof_map_1() : ufc::dof_map()
1470
__global_dimension = 0;
1474
virtual ~poisson2d_4_0_dof_map_1()
1479
/// Return a string identifying the dof map
1480
virtual const char* signature() const
1482
return "FFC dof map for FiniteElement('Lagrange', 'triangle', 4)";
1485
/// Return true iff mesh entities of topological dimension d are needed
1486
virtual bool needs_mesh_entities(unsigned int d) const
1503
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1504
virtual bool init_mesh(const ufc::mesh& m)
1506
__global_dimension = m.num_entities[0] + 3*m.num_entities[1] + 3*m.num_entities[2];
1510
/// Initialize dof map for given cell
1511
virtual void init_cell(const ufc::mesh& m,
1517
/// Finish initialization of dof map for cells
1518
virtual void init_cell_finalize()
1523
/// Return the dimension of the global finite element function space
1524
virtual unsigned int global_dimension() const
1526
return __global_dimension;
1529
/// Return the dimension of the local finite element function space for a cell
1530
virtual unsigned int local_dimension(const ufc::cell& c) const
1535
/// Return the maximum dimension of the local finite element function space
1536
virtual unsigned int max_local_dimension() const
1541
// Return the geometric dimension of the coordinates this dof map provides
1542
virtual unsigned int geometric_dimension() const
1547
/// Return the number of dofs on each cell facet
1548
virtual unsigned int num_facet_dofs() const
1553
/// Return the number of dofs associated with each cell entity of dimension d
1554
virtual unsigned int num_entity_dofs(unsigned int d) const
1556
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1559
/// Tabulate the local-to-global mapping of dofs on a cell
1560
virtual void tabulate_dofs(unsigned int* dofs,
1562
const ufc::cell& c) const
1564
dofs[0] = c.entity_indices[0][0];
1565
dofs[1] = c.entity_indices[0][1];
1566
dofs[2] = c.entity_indices[0][2];
1567
unsigned int offset = m.num_entities[0];
1568
dofs[3] = offset + 3*c.entity_indices[1][0];
1569
dofs[4] = offset + 3*c.entity_indices[1][0] + 1;
1570
dofs[5] = offset + 3*c.entity_indices[1][0] + 2;
1571
dofs[6] = offset + 3*c.entity_indices[1][1];
1572
dofs[7] = offset + 3*c.entity_indices[1][1] + 1;
1573
dofs[8] = offset + 3*c.entity_indices[1][1] + 2;
1574
dofs[9] = offset + 3*c.entity_indices[1][2];
1575
dofs[10] = offset + 3*c.entity_indices[1][2] + 1;
1576
dofs[11] = offset + 3*c.entity_indices[1][2] + 2;
1577
offset = offset + 3*m.num_entities[1];
1578
dofs[12] = offset + 3*c.entity_indices[2][0];
1579
dofs[13] = offset + 3*c.entity_indices[2][0] + 1;
1580
dofs[14] = offset + 3*c.entity_indices[2][0] + 2;
1583
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1584
virtual void tabulate_facet_dofs(unsigned int* dofs,
1585
unsigned int facet) const
1613
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1614
virtual void tabulate_entity_dofs(unsigned int* dofs,
1615
unsigned int d, unsigned int i) const
1617
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1620
/// Tabulate the coordinates of all dofs on a cell
1621
virtual void tabulate_coordinates(double** coordinates,
1622
const ufc::cell& c) const
1624
const double * const * x = c.coordinates;
1625
coordinates[0][0] = x[0][0];
1626
coordinates[0][1] = x[0][1];
1627
coordinates[1][0] = x[1][0];
1628
coordinates[1][1] = x[1][1];
1629
coordinates[2][0] = x[2][0];
1630
coordinates[2][1] = x[2][1];
1631
coordinates[3][0] = 0.75*x[1][0] + 0.25*x[2][0];
1632
coordinates[3][1] = 0.75*x[1][1] + 0.25*x[2][1];
1633
coordinates[4][0] = 0.5*x[1][0] + 0.5*x[2][0];
1634
coordinates[4][1] = 0.5*x[1][1] + 0.5*x[2][1];
1635
coordinates[5][0] = 0.25*x[1][0] + 0.75*x[2][0];
1636
coordinates[5][1] = 0.25*x[1][1] + 0.75*x[2][1];
1637
coordinates[6][0] = 0.75*x[0][0] + 0.25*x[2][0];
1638
coordinates[6][1] = 0.75*x[0][1] + 0.25*x[2][1];
1639
coordinates[7][0] = 0.5*x[0][0] + 0.5*x[2][0];
1640
coordinates[7][1] = 0.5*x[0][1] + 0.5*x[2][1];
1641
coordinates[8][0] = 0.25*x[0][0] + 0.75*x[2][0];
1642
coordinates[8][1] = 0.25*x[0][1] + 0.75*x[2][1];
1643
coordinates[9][0] = 0.75*x[0][0] + 0.25*x[1][0];
1644
coordinates[9][1] = 0.75*x[0][1] + 0.25*x[1][1];
1645
coordinates[10][0] = 0.5*x[0][0] + 0.5*x[1][0];
1646
coordinates[10][1] = 0.5*x[0][1] + 0.5*x[1][1];
1647
coordinates[11][0] = 0.25*x[0][0] + 0.75*x[1][0];
1648
coordinates[11][1] = 0.25*x[0][1] + 0.75*x[1][1];
1649
coordinates[12][0] = 0.5*x[0][0] + 0.25*x[1][0] + 0.25*x[2][0];
1650
coordinates[12][1] = 0.5*x[0][1] + 0.25*x[1][1] + 0.25*x[2][1];
1651
coordinates[13][0] = 0.25*x[0][0] + 0.5*x[1][0] + 0.25*x[2][0];
1652
coordinates[13][1] = 0.25*x[0][1] + 0.5*x[1][1] + 0.25*x[2][1];
1653
coordinates[14][0] = 0.25*x[0][0] + 0.25*x[1][0] + 0.5*x[2][0];
1654
coordinates[14][1] = 0.25*x[0][1] + 0.25*x[1][1] + 0.5*x[2][1];
1657
/// Return the number of sub dof maps (for a mixed element)
1658
virtual unsigned int num_sub_dof_maps() const
1663
/// Create a new dof_map for sub dof map i (for a mixed element)
1664
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1666
return new poisson2d_4_0_dof_map_1();