1
// This code conforms with the UFC specification version 1.0
2
// and was automatically generated by FFC version 0.4.5.
4
// Warning: This code was generated with the option '-l dolfin'
5
// and contains DOLFIN-specific wrappers that depend on DOLFIN.
15
/// This class defines the interface for a finite element.
17
class UFC_Poisson2DBilinearForm_finite_element_0: public ufc::finite_element
22
UFC_Poisson2DBilinearForm_finite_element_0() : ufc::finite_element()
28
virtual ~UFC_Poisson2DBilinearForm_finite_element_0()
33
/// Return a string identifying the finite element
34
virtual const char* signature() const
36
return "Lagrange finite element of degree 1 on a triangle";
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);
108
// Compute psitilde_a
109
const double psitilde_a_0 = 1;
110
const double psitilde_a_1 = x;
112
// Compute psitilde_bs
113
const double psitilde_bs_0_0 = 1;
114
const double psitilde_bs_0_1 = 1.5*y + 0.5;
115
const double psitilde_bs_1_0 = 1;
117
// Compute basisvalues
118
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
119
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
120
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
122
// Table(s) of coefficients
123
const static double coefficients0[3][3] = \
124
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
125
{0.471404520791032, 0.288675134594813, -0.166666666666667},
126
{0.471404520791032, 0, 0.333333333333333}};
128
// Extract relevant coefficients
129
const double coeff0_0 = coefficients0[dof][0];
130
const double coeff0_1 = coefficients0[dof][1];
131
const double coeff0_2 = coefficients0[dof][2];
134
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
137
/// Evaluate all basis functions at given point in cell
138
virtual void evaluate_basis_all(double* values,
139
const double* coordinates,
140
const ufc::cell& c) const
142
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
145
/// Evaluate order n derivatives of basis function i at given point in cell
146
virtual void evaluate_basis_derivatives(unsigned int i,
149
const double* coordinates,
150
const ufc::cell& c) const
152
// Extract vertex coordinates
153
const double * const * element_coordinates = c.coordinates;
155
// Compute Jacobian of affine map from reference cell
156
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
157
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
158
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
159
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
161
// Compute determinant of Jacobian
162
const double detJ = J_00*J_11 - J_01*J_10;
164
// Compute inverse of Jacobian
166
// Get coordinates and map to the reference (UFC) element
167
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
168
element_coordinates[0][0]*element_coordinates[2][1] +\
169
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
170
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
171
element_coordinates[1][0]*element_coordinates[0][1] -\
172
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
174
// Map coordinates to the reference square
175
if (std::abs(y - 1.0) < 1e-14)
178
x = 2.0 *x/(1.0 - y) - 1.0;
181
// Compute number of derivatives
182
unsigned int num_derivatives = 1;
184
for (unsigned int j = 0; j < n; j++)
185
num_derivatives *= 2;
188
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
189
unsigned int **combinations = new unsigned int *[num_derivatives];
191
for (unsigned int j = 0; j < num_derivatives; j++)
193
combinations[j] = new unsigned int [n];
194
for (unsigned int k = 0; k < n; k++)
195
combinations[j][k] = 0;
198
// Generate combinations of derivatives
199
for (unsigned int row = 1; row < num_derivatives; row++)
201
for (unsigned int num = 0; num < row; num++)
203
for (unsigned int col = n-1; col+1 > 0; col--)
205
if (combinations[row][col] + 1 > 1)
206
combinations[row][col] = 0;
209
combinations[row][col] += 1;
216
// Compute inverse of Jacobian
217
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
219
// Declare transformation matrix
220
// Declare pointer to two dimensional array and initialise
221
double **transform = new double *[num_derivatives];
223
for (unsigned int j = 0; j < num_derivatives; j++)
225
transform[j] = new double [num_derivatives];
226
for (unsigned int k = 0; k < num_derivatives; k++)
230
// Construct transformation matrix
231
for (unsigned int row = 0; row < num_derivatives; row++)
233
for (unsigned int col = 0; col < num_derivatives; col++)
235
for (unsigned int k = 0; k < n; k++)
236
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
241
for (unsigned int j = 0; j < 1*num_derivatives; j++)
244
// Map degree of freedom to element degree of freedom
245
const unsigned int dof = i;
248
const double scalings_y_0 = 1;
249
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
251
// Compute psitilde_a
252
const double psitilde_a_0 = 1;
253
const double psitilde_a_1 = x;
255
// Compute psitilde_bs
256
const double psitilde_bs_0_0 = 1;
257
const double psitilde_bs_0_1 = 1.5*y + 0.5;
258
const double psitilde_bs_1_0 = 1;
260
// Compute basisvalues
261
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
262
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
263
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
265
// Table(s) of coefficients
266
const static double coefficients0[3][3] = \
267
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
268
{0.471404520791032, 0.288675134594813, -0.166666666666667},
269
{0.471404520791032, 0, 0.333333333333333}};
271
// Interesting (new) part
272
// Tables of derivatives of the polynomial base (transpose)
273
const static double dmats0[3][3] = \
275
{4.89897948556636, 0, 0},
278
const static double dmats1[3][3] = \
280
{2.44948974278318, 0, 0},
281
{4.24264068711928, 0, 0}};
283
// Compute reference derivatives
284
// Declare pointer to array of derivatives on FIAT element
285
double *derivatives = new double [num_derivatives];
287
// Declare coefficients
292
// Declare new coefficients
293
double new_coeff0_0 = 0;
294
double new_coeff0_1 = 0;
295
double new_coeff0_2 = 0;
297
// Loop possible derivatives
298
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
300
// Get values from coefficients array
301
new_coeff0_0 = coefficients0[dof][0];
302
new_coeff0_1 = coefficients0[dof][1];
303
new_coeff0_2 = coefficients0[dof][2];
305
// Loop derivative order
306
for (unsigned int j = 0; j < n; j++)
308
// Update old coefficients
309
coeff0_0 = new_coeff0_0;
310
coeff0_1 = new_coeff0_1;
311
coeff0_2 = new_coeff0_2;
313
if(combinations[deriv_num][j] == 0)
315
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
316
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
317
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
319
if(combinations[deriv_num][j] == 1)
321
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
322
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
323
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
327
// Compute derivatives on reference element as dot product of coefficients and basisvalues
328
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
331
// Transform derivatives back to physical element
332
for (unsigned int row = 0; row < num_derivatives; row++)
334
for (unsigned int col = 0; col < num_derivatives; col++)
336
values[row] += transform[row][col]*derivatives[col];
339
// Delete pointer to array of derivatives on FIAT element
340
delete [] derivatives;
342
// Delete pointer to array of combinations of derivatives and transform
343
for (unsigned int row = 0; row < num_derivatives; row++)
345
delete [] combinations[row];
346
delete [] transform[row];
349
delete [] combinations;
353
/// Evaluate order n derivatives of all basis functions at given point in cell
354
virtual void evaluate_basis_derivatives_all(unsigned int n,
356
const double* coordinates,
357
const ufc::cell& c) const
359
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
362
/// Evaluate linear functional for dof i on the function f
363
virtual double evaluate_dof(unsigned int i,
364
const ufc::function& f,
365
const ufc::cell& c) const
367
// The reference points, direction and weights:
368
const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
369
const static double W[3][1] = {{1}, {1}, {1}};
370
const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
372
const double * const * x = c.coordinates;
374
// Iterate over the points:
375
// Evaluate basis functions for affine mapping
376
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
377
const double w1 = X[i][0][0];
378
const double w2 = X[i][0][1];
380
// Compute affine mapping y = F(X)
382
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
383
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
385
// Evaluate function at physical points
387
f.evaluate(values, y, c);
389
// Map function values using appropriate mapping
390
// Affine map: Do nothing
392
// Note that we do not map the weights (yet).
394
// Take directional components
395
for(int k = 0; k < 1; k++)
396
result += values[k]*D[i][0][k];
397
// Multiply by weights
403
/// Evaluate linear functionals for all dofs on the function f
404
virtual void evaluate_dofs(double* values,
405
const ufc::function& f,
406
const ufc::cell& c) const
408
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
411
/// Interpolate vertex values from dof values
412
virtual void interpolate_vertex_values(double* vertex_values,
413
const double* dof_values,
414
const ufc::cell& c) const
416
// Evaluate at vertices and use affine mapping
417
vertex_values[0] = dof_values[0];
418
vertex_values[1] = dof_values[1];
419
vertex_values[2] = dof_values[2];
422
/// Return the number of sub elements (for a mixed element)
423
virtual unsigned int num_sub_elements() const
428
/// Create a new finite element for sub element i (for a mixed element)
429
virtual ufc::finite_element* create_sub_element(unsigned int i) const
431
return new UFC_Poisson2DBilinearForm_finite_element_0();
436
/// This class defines the interface for a finite element.
438
class UFC_Poisson2DBilinearForm_finite_element_1: public ufc::finite_element
443
UFC_Poisson2DBilinearForm_finite_element_1() : ufc::finite_element()
449
virtual ~UFC_Poisson2DBilinearForm_finite_element_1()
454
/// Return a string identifying the finite element
455
virtual const char* signature() const
457
return "Lagrange finite element of degree 1 on a triangle";
460
/// Return the cell shape
461
virtual ufc::shape cell_shape() const
463
return ufc::triangle;
466
/// Return the dimension of the finite element function space
467
virtual unsigned int space_dimension() const
472
/// Return the rank of the value space
473
virtual unsigned int value_rank() const
478
/// Return the dimension of the value space for axis i
479
virtual unsigned int value_dimension(unsigned int i) const
484
/// Evaluate basis function i at given point in cell
485
virtual void evaluate_basis(unsigned int i,
487
const double* coordinates,
488
const ufc::cell& c) const
490
// Extract vertex coordinates
491
const double * const * element_coordinates = c.coordinates;
493
// Compute Jacobian of affine map from reference cell
494
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
495
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
496
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
497
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
499
// Compute determinant of Jacobian
500
const double detJ = J_00*J_11 - J_01*J_10;
502
// Compute inverse of Jacobian
504
// Get coordinates and map to the reference (UFC) element
505
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
506
element_coordinates[0][0]*element_coordinates[2][1] +\
507
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
508
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
509
element_coordinates[1][0]*element_coordinates[0][1] -\
510
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
512
// Map coordinates to the reference square
513
if (std::abs(y - 1.0) < 1e-14)
516
x = 2.0 *x/(1.0 - y) - 1.0;
522
// Map degree of freedom to element degree of freedom
523
const unsigned int dof = i;
526
const double scalings_y_0 = 1;
527
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
529
// Compute psitilde_a
530
const double psitilde_a_0 = 1;
531
const double psitilde_a_1 = x;
533
// Compute psitilde_bs
534
const double psitilde_bs_0_0 = 1;
535
const double psitilde_bs_0_1 = 1.5*y + 0.5;
536
const double psitilde_bs_1_0 = 1;
538
// Compute basisvalues
539
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
540
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
541
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
543
// Table(s) of coefficients
544
const static double coefficients0[3][3] = \
545
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
546
{0.471404520791032, 0.288675134594813, -0.166666666666667},
547
{0.471404520791032, 0, 0.333333333333333}};
549
// Extract relevant coefficients
550
const double coeff0_0 = coefficients0[dof][0];
551
const double coeff0_1 = coefficients0[dof][1];
552
const double coeff0_2 = coefficients0[dof][2];
555
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
558
/// Evaluate all basis functions at given point in cell
559
virtual void evaluate_basis_all(double* values,
560
const double* coordinates,
561
const ufc::cell& c) const
563
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
566
/// Evaluate order n derivatives of basis function i at given point in cell
567
virtual void evaluate_basis_derivatives(unsigned int i,
570
const double* coordinates,
571
const ufc::cell& c) const
573
// Extract vertex coordinates
574
const double * const * element_coordinates = c.coordinates;
576
// Compute Jacobian of affine map from reference cell
577
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
578
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
579
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
580
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
582
// Compute determinant of Jacobian
583
const double detJ = J_00*J_11 - J_01*J_10;
585
// Compute inverse of Jacobian
587
// Get coordinates and map to the reference (UFC) element
588
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
589
element_coordinates[0][0]*element_coordinates[2][1] +\
590
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
591
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
592
element_coordinates[1][0]*element_coordinates[0][1] -\
593
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
595
// Map coordinates to the reference square
596
if (std::abs(y - 1.0) < 1e-14)
599
x = 2.0 *x/(1.0 - y) - 1.0;
602
// Compute number of derivatives
603
unsigned int num_derivatives = 1;
605
for (unsigned int j = 0; j < n; j++)
606
num_derivatives *= 2;
609
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
610
unsigned int **combinations = new unsigned int *[num_derivatives];
612
for (unsigned int j = 0; j < num_derivatives; j++)
614
combinations[j] = new unsigned int [n];
615
for (unsigned int k = 0; k < n; k++)
616
combinations[j][k] = 0;
619
// Generate combinations of derivatives
620
for (unsigned int row = 1; row < num_derivatives; row++)
622
for (unsigned int num = 0; num < row; num++)
624
for (unsigned int col = n-1; col+1 > 0; col--)
626
if (combinations[row][col] + 1 > 1)
627
combinations[row][col] = 0;
630
combinations[row][col] += 1;
637
// Compute inverse of Jacobian
638
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
640
// Declare transformation matrix
641
// Declare pointer to two dimensional array and initialise
642
double **transform = new double *[num_derivatives];
644
for (unsigned int j = 0; j < num_derivatives; j++)
646
transform[j] = new double [num_derivatives];
647
for (unsigned int k = 0; k < num_derivatives; k++)
651
// Construct transformation matrix
652
for (unsigned int row = 0; row < num_derivatives; row++)
654
for (unsigned int col = 0; col < num_derivatives; col++)
656
for (unsigned int k = 0; k < n; k++)
657
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
662
for (unsigned int j = 0; j < 1*num_derivatives; j++)
665
// Map degree of freedom to element degree of freedom
666
const unsigned int dof = i;
669
const double scalings_y_0 = 1;
670
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
672
// Compute psitilde_a
673
const double psitilde_a_0 = 1;
674
const double psitilde_a_1 = x;
676
// Compute psitilde_bs
677
const double psitilde_bs_0_0 = 1;
678
const double psitilde_bs_0_1 = 1.5*y + 0.5;
679
const double psitilde_bs_1_0 = 1;
681
// Compute basisvalues
682
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
683
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
684
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
686
// Table(s) of coefficients
687
const static double coefficients0[3][3] = \
688
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
689
{0.471404520791032, 0.288675134594813, -0.166666666666667},
690
{0.471404520791032, 0, 0.333333333333333}};
692
// Interesting (new) part
693
// Tables of derivatives of the polynomial base (transpose)
694
const static double dmats0[3][3] = \
696
{4.89897948556636, 0, 0},
699
const static double dmats1[3][3] = \
701
{2.44948974278318, 0, 0},
702
{4.24264068711928, 0, 0}};
704
// Compute reference derivatives
705
// Declare pointer to array of derivatives on FIAT element
706
double *derivatives = new double [num_derivatives];
708
// Declare coefficients
713
// Declare new coefficients
714
double new_coeff0_0 = 0;
715
double new_coeff0_1 = 0;
716
double new_coeff0_2 = 0;
718
// Loop possible derivatives
719
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
721
// Get values from coefficients array
722
new_coeff0_0 = coefficients0[dof][0];
723
new_coeff0_1 = coefficients0[dof][1];
724
new_coeff0_2 = coefficients0[dof][2];
726
// Loop derivative order
727
for (unsigned int j = 0; j < n; j++)
729
// Update old coefficients
730
coeff0_0 = new_coeff0_0;
731
coeff0_1 = new_coeff0_1;
732
coeff0_2 = new_coeff0_2;
734
if(combinations[deriv_num][j] == 0)
736
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
737
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
738
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
740
if(combinations[deriv_num][j] == 1)
742
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
743
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
744
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
748
// Compute derivatives on reference element as dot product of coefficients and basisvalues
749
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
752
// Transform derivatives back to physical element
753
for (unsigned int row = 0; row < num_derivatives; row++)
755
for (unsigned int col = 0; col < num_derivatives; col++)
757
values[row] += transform[row][col]*derivatives[col];
760
// Delete pointer to array of derivatives on FIAT element
761
delete [] derivatives;
763
// Delete pointer to array of combinations of derivatives and transform
764
for (unsigned int row = 0; row < num_derivatives; row++)
766
delete [] combinations[row];
767
delete [] transform[row];
770
delete [] combinations;
774
/// Evaluate order n derivatives of all basis functions at given point in cell
775
virtual void evaluate_basis_derivatives_all(unsigned int n,
777
const double* coordinates,
778
const ufc::cell& c) const
780
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
783
/// Evaluate linear functional for dof i on the function f
784
virtual double evaluate_dof(unsigned int i,
785
const ufc::function& f,
786
const ufc::cell& c) const
788
// The reference points, direction and weights:
789
const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
790
const static double W[3][1] = {{1}, {1}, {1}};
791
const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
793
const double * const * x = c.coordinates;
795
// Iterate over the points:
796
// Evaluate basis functions for affine mapping
797
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
798
const double w1 = X[i][0][0];
799
const double w2 = X[i][0][1];
801
// Compute affine mapping y = F(X)
803
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
804
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
806
// Evaluate function at physical points
808
f.evaluate(values, y, c);
810
// Map function values using appropriate mapping
811
// Affine map: Do nothing
813
// Note that we do not map the weights (yet).
815
// Take directional components
816
for(int k = 0; k < 1; k++)
817
result += values[k]*D[i][0][k];
818
// Multiply by weights
824
/// Evaluate linear functionals for all dofs on the function f
825
virtual void evaluate_dofs(double* values,
826
const ufc::function& f,
827
const ufc::cell& c) const
829
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
832
/// Interpolate vertex values from dof values
833
virtual void interpolate_vertex_values(double* vertex_values,
834
const double* dof_values,
835
const ufc::cell& c) const
837
// Evaluate at vertices and use affine mapping
838
vertex_values[0] = dof_values[0];
839
vertex_values[1] = dof_values[1];
840
vertex_values[2] = dof_values[2];
843
/// Return the number of sub elements (for a mixed element)
844
virtual unsigned int num_sub_elements() const
849
/// Create a new finite element for sub element i (for a mixed element)
850
virtual ufc::finite_element* create_sub_element(unsigned int i) const
852
return new UFC_Poisson2DBilinearForm_finite_element_1();
857
/// This class defines the interface for a local-to-global mapping of
858
/// degrees of freedom (dofs).
860
class UFC_Poisson2DBilinearForm_dof_map_0: public ufc::dof_map
864
unsigned int __global_dimension;
869
UFC_Poisson2DBilinearForm_dof_map_0() : ufc::dof_map()
871
__global_dimension = 0;
875
virtual ~UFC_Poisson2DBilinearForm_dof_map_0()
880
/// Return a string identifying the dof map
881
virtual const char* signature() const
883
return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
886
/// Return true iff mesh entities of topological dimension d are needed
887
virtual bool needs_mesh_entities(unsigned int d) const
904
/// Initialize dof map for mesh (return true iff init_cell() is needed)
905
virtual bool init_mesh(const ufc::mesh& m)
907
__global_dimension = m.num_entities[0];
911
/// Initialize dof map for given cell
912
virtual void init_cell(const ufc::mesh& m,
918
/// Finish initialization of dof map for cells
919
virtual void init_cell_finalize()
924
/// Return the dimension of the global finite element function space
925
virtual unsigned int global_dimension() const
927
return __global_dimension;
930
/// Return the dimension of the local finite element function space
931
virtual unsigned int local_dimension() const
936
// Return the geometric dimension of the coordinates this dof map provides
937
virtual unsigned int geometric_dimension() const
942
/// Return the number of dofs on each cell facet
943
virtual unsigned int num_facet_dofs() const
948
/// Return the number of dofs associated with each cell entity of dimension d
949
virtual unsigned int num_entity_dofs(unsigned int d) const
951
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
954
/// Tabulate the local-to-global mapping of dofs on a cell
955
virtual void tabulate_dofs(unsigned int* dofs,
957
const ufc::cell& c) const
959
dofs[0] = c.entity_indices[0][0];
960
dofs[1] = c.entity_indices[0][1];
961
dofs[2] = c.entity_indices[0][2];
964
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
965
virtual void tabulate_facet_dofs(unsigned int* dofs,
966
unsigned int facet) const
985
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
986
virtual void tabulate_entity_dofs(unsigned int* dofs,
987
unsigned int d, unsigned int i) const
989
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
992
/// Tabulate the coordinates of all dofs on a cell
993
virtual void tabulate_coordinates(double** coordinates,
994
const ufc::cell& c) const
996
const double * const * x = c.coordinates;
997
coordinates[0][0] = x[0][0];
998
coordinates[0][1] = x[0][1];
999
coordinates[1][0] = x[1][0];
1000
coordinates[1][1] = x[1][1];
1001
coordinates[2][0] = x[2][0];
1002
coordinates[2][1] = x[2][1];
1005
/// Return the number of sub dof maps (for a mixed element)
1006
virtual unsigned int num_sub_dof_maps() const
1011
/// Create a new dof_map for sub dof map i (for a mixed element)
1012
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1014
return new UFC_Poisson2DBilinearForm_dof_map_0();
1019
/// This class defines the interface for a local-to-global mapping of
1020
/// degrees of freedom (dofs).
1022
class UFC_Poisson2DBilinearForm_dof_map_1: public ufc::dof_map
1026
unsigned int __global_dimension;
1031
UFC_Poisson2DBilinearForm_dof_map_1() : ufc::dof_map()
1033
__global_dimension = 0;
1037
virtual ~UFC_Poisson2DBilinearForm_dof_map_1()
1042
/// Return a string identifying the dof map
1043
virtual const char* signature() const
1045
return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
1048
/// Return true iff mesh entities of topological dimension d are needed
1049
virtual bool needs_mesh_entities(unsigned int d) const
1066
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1067
virtual bool init_mesh(const ufc::mesh& m)
1069
__global_dimension = m.num_entities[0];
1073
/// Initialize dof map for given cell
1074
virtual void init_cell(const ufc::mesh& m,
1080
/// Finish initialization of dof map for cells
1081
virtual void init_cell_finalize()
1086
/// Return the dimension of the global finite element function space
1087
virtual unsigned int global_dimension() const
1089
return __global_dimension;
1092
/// Return the dimension of the local finite element function space
1093
virtual unsigned int local_dimension() const
1098
// Return the geometric dimension of the coordinates this dof map provides
1099
virtual unsigned int geometric_dimension() const
1104
/// Return the number of dofs on each cell facet
1105
virtual unsigned int num_facet_dofs() const
1110
/// Return the number of dofs associated with each cell entity of dimension d
1111
virtual unsigned int num_entity_dofs(unsigned int d) const
1113
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1116
/// Tabulate the local-to-global mapping of dofs on a cell
1117
virtual void tabulate_dofs(unsigned int* dofs,
1119
const ufc::cell& c) const
1121
dofs[0] = c.entity_indices[0][0];
1122
dofs[1] = c.entity_indices[0][1];
1123
dofs[2] = c.entity_indices[0][2];
1126
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1127
virtual void tabulate_facet_dofs(unsigned int* dofs,
1128
unsigned int facet) const
1147
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1148
virtual void tabulate_entity_dofs(unsigned int* dofs,
1149
unsigned int d, unsigned int i) const
1151
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1154
/// Tabulate the coordinates of all dofs on a cell
1155
virtual void tabulate_coordinates(double** coordinates,
1156
const ufc::cell& c) const
1158
const double * const * x = c.coordinates;
1159
coordinates[0][0] = x[0][0];
1160
coordinates[0][1] = x[0][1];
1161
coordinates[1][0] = x[1][0];
1162
coordinates[1][1] = x[1][1];
1163
coordinates[2][0] = x[2][0];
1164
coordinates[2][1] = x[2][1];
1167
/// Return the number of sub dof maps (for a mixed element)
1168
virtual unsigned int num_sub_dof_maps() const
1173
/// Create a new dof_map for sub dof map i (for a mixed element)
1174
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1176
return new UFC_Poisson2DBilinearForm_dof_map_1();
1181
/// This class defines the interface for the tabulation of the cell
1182
/// tensor corresponding to the local contribution to a form from
1183
/// the integral over a cell.
1185
class UFC_Poisson2DBilinearForm_cell_integral_0: public ufc::cell_integral
1190
UFC_Poisson2DBilinearForm_cell_integral_0() : ufc::cell_integral()
1196
virtual ~UFC_Poisson2DBilinearForm_cell_integral_0()
1201
/// Tabulate the tensor for the contribution from a local cell
1202
virtual void tabulate_tensor(double* A,
1203
const double * const * w,
1204
const ufc::cell& c) const
1206
// Extract vertex coordinates
1207
const double * const * x = c.coordinates;
1209
// Compute Jacobian of affine map from reference cell
1210
const double J_00 = x[1][0] - x[0][0];
1211
const double J_01 = x[2][0] - x[0][0];
1212
const double J_10 = x[1][1] - x[0][1];
1213
const double J_11 = x[2][1] - x[0][1];
1215
// Compute determinant of Jacobian
1216
double detJ = J_00*J_11 - J_01*J_10;
1218
// Compute inverse of Jacobian
1219
const double Jinv_00 = J_11 / detJ;
1220
const double Jinv_01 = -J_01 / detJ;
1221
const double Jinv_10 = -J_10 / detJ;
1222
const double Jinv_11 = J_00 / detJ;
1225
const double det = std::abs(detJ);
1227
// Compute geometry tensors
1228
const double G0_0_0 = det*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
1229
const double G0_0_1 = det*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
1230
const double G0_1_0 = det*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
1231
const double G0_1_1 = det*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
1233
// Compute element tensor
1234
A[0] = 0.5*G0_0_0 + 0.5*G0_0_1 + 0.5*G0_1_0 + 0.5*G0_1_1;
1235
A[1] = -0.5*G0_0_0 - 0.5*G0_1_0;
1236
A[2] = -0.5*G0_0_1 - 0.5*G0_1_1;
1237
A[3] = -0.5*G0_0_0 - 0.5*G0_0_1;
1240
A[6] = -0.5*G0_1_0 - 0.5*G0_1_1;
1247
/// This class defines the interface for the assembly of the global
1248
/// tensor corresponding to a form with r + n arguments, that is, a
1251
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
1253
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
1254
/// global tensor A is defined by
1256
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
1258
/// where each argument Vj represents the application to the
1259
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
1260
/// fixed functions (coefficients).
1262
class UFC_Poisson2DBilinearForm: public ufc::form
1267
UFC_Poisson2DBilinearForm() : ufc::form()
1273
virtual ~UFC_Poisson2DBilinearForm()
1278
/// Return a string identifying the form
1279
virtual const char* signature() const
1281
return "(dXa0[0, 1]/dxb0[0, 1])(dXa1[0, 1]/dxb0[0, 1]) | ((d/dXa0[0, 1])vi0[0, 1, 2])*((d/dXa1[0, 1])vi1[0, 1, 2])*dX(0)";
1284
/// Return the rank of the global tensor (r)
1285
virtual unsigned int rank() const
1290
/// Return the number of coefficients (n)
1291
virtual unsigned int num_coefficients() const
1296
/// Return the number of cell integrals
1297
virtual unsigned int num_cell_integrals() const
1302
/// Return the number of exterior facet integrals
1303
virtual unsigned int num_exterior_facet_integrals() const
1308
/// Return the number of interior facet integrals
1309
virtual unsigned int num_interior_facet_integrals() const
1314
/// Create a new finite element for argument function i
1315
virtual ufc::finite_element* create_finite_element(unsigned int i) const
1320
return new UFC_Poisson2DBilinearForm_finite_element_0();
1323
return new UFC_Poisson2DBilinearForm_finite_element_1();
1329
/// Create a new dof map for argument function i
1330
virtual ufc::dof_map* create_dof_map(unsigned int i) const
1335
return new UFC_Poisson2DBilinearForm_dof_map_0();
1338
return new UFC_Poisson2DBilinearForm_dof_map_1();
1344
/// Create a new cell integral on sub domain i
1345
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
1347
return new UFC_Poisson2DBilinearForm_cell_integral_0();
1350
/// Create a new exterior facet integral on sub domain i
1351
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
1356
/// Create a new interior facet integral on sub domain i
1357
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
1364
/// This class defines the interface for a finite element.
1366
class UFC_Poisson2DLinearForm_finite_element_0: public ufc::finite_element
1371
UFC_Poisson2DLinearForm_finite_element_0() : ufc::finite_element()
1377
virtual ~UFC_Poisson2DLinearForm_finite_element_0()
1382
/// Return a string identifying the finite element
1383
virtual const char* signature() const
1385
return "Lagrange finite element of degree 1 on a triangle";
1388
/// Return the cell shape
1389
virtual ufc::shape cell_shape() const
1391
return ufc::triangle;
1394
/// Return the dimension of the finite element function space
1395
virtual unsigned int space_dimension() const
1400
/// Return the rank of the value space
1401
virtual unsigned int value_rank() const
1406
/// Return the dimension of the value space for axis i
1407
virtual unsigned int value_dimension(unsigned int i) const
1412
/// Evaluate basis function i at given point in cell
1413
virtual void evaluate_basis(unsigned int i,
1415
const double* coordinates,
1416
const ufc::cell& c) const
1418
// Extract vertex coordinates
1419
const double * const * element_coordinates = c.coordinates;
1421
// Compute Jacobian of affine map from reference cell
1422
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1423
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1424
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1425
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1427
// Compute determinant of Jacobian
1428
const double detJ = J_00*J_11 - J_01*J_10;
1430
// Compute inverse of Jacobian
1432
// Get coordinates and map to the reference (UFC) element
1433
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1434
element_coordinates[0][0]*element_coordinates[2][1] +\
1435
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1436
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1437
element_coordinates[1][0]*element_coordinates[0][1] -\
1438
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1440
// Map coordinates to the reference square
1441
if (std::abs(y - 1.0) < 1e-14)
1444
x = 2.0 *x/(1.0 - y) - 1.0;
1450
// Map degree of freedom to element degree of freedom
1451
const unsigned int dof = i;
1453
// Generate scalings
1454
const double scalings_y_0 = 1;
1455
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1457
// Compute psitilde_a
1458
const double psitilde_a_0 = 1;
1459
const double psitilde_a_1 = x;
1461
// Compute psitilde_bs
1462
const double psitilde_bs_0_0 = 1;
1463
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1464
const double psitilde_bs_1_0 = 1;
1466
// Compute basisvalues
1467
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1468
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1469
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1471
// Table(s) of coefficients
1472
const static double coefficients0[3][3] = \
1473
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
1474
{0.471404520791032, 0.288675134594813, -0.166666666666667},
1475
{0.471404520791032, 0, 0.333333333333333}};
1477
// Extract relevant coefficients
1478
const double coeff0_0 = coefficients0[dof][0];
1479
const double coeff0_1 = coefficients0[dof][1];
1480
const double coeff0_2 = coefficients0[dof][2];
1483
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
1486
/// Evaluate all basis functions at given point in cell
1487
virtual void evaluate_basis_all(double* values,
1488
const double* coordinates,
1489
const ufc::cell& c) const
1491
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1494
/// Evaluate order n derivatives of basis function i at given point in cell
1495
virtual void evaluate_basis_derivatives(unsigned int i,
1498
const double* coordinates,
1499
const ufc::cell& c) const
1501
// Extract vertex coordinates
1502
const double * const * element_coordinates = c.coordinates;
1504
// Compute Jacobian of affine map from reference cell
1505
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1506
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1507
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1508
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1510
// Compute determinant of Jacobian
1511
const double detJ = J_00*J_11 - J_01*J_10;
1513
// Compute inverse of Jacobian
1515
// Get coordinates and map to the reference (UFC) element
1516
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1517
element_coordinates[0][0]*element_coordinates[2][1] +\
1518
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1519
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1520
element_coordinates[1][0]*element_coordinates[0][1] -\
1521
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1523
// Map coordinates to the reference square
1524
if (std::abs(y - 1.0) < 1e-14)
1527
x = 2.0 *x/(1.0 - y) - 1.0;
1530
// Compute number of derivatives
1531
unsigned int num_derivatives = 1;
1533
for (unsigned int j = 0; j < n; j++)
1534
num_derivatives *= 2;
1537
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1538
unsigned int **combinations = new unsigned int *[num_derivatives];
1540
for (unsigned int j = 0; j < num_derivatives; j++)
1542
combinations[j] = new unsigned int [n];
1543
for (unsigned int k = 0; k < n; k++)
1544
combinations[j][k] = 0;
1547
// Generate combinations of derivatives
1548
for (unsigned int row = 1; row < num_derivatives; row++)
1550
for (unsigned int num = 0; num < row; num++)
1552
for (unsigned int col = n-1; col+1 > 0; col--)
1554
if (combinations[row][col] + 1 > 1)
1555
combinations[row][col] = 0;
1558
combinations[row][col] += 1;
1565
// Compute inverse of Jacobian
1566
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1568
// Declare transformation matrix
1569
// Declare pointer to two dimensional array and initialise
1570
double **transform = new double *[num_derivatives];
1572
for (unsigned int j = 0; j < num_derivatives; j++)
1574
transform[j] = new double [num_derivatives];
1575
for (unsigned int k = 0; k < num_derivatives; k++)
1576
transform[j][k] = 1;
1579
// Construct transformation matrix
1580
for (unsigned int row = 0; row < num_derivatives; row++)
1582
for (unsigned int col = 0; col < num_derivatives; col++)
1584
for (unsigned int k = 0; k < n; k++)
1585
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1590
for (unsigned int j = 0; j < 1*num_derivatives; j++)
1593
// Map degree of freedom to element degree of freedom
1594
const unsigned int dof = i;
1596
// Generate scalings
1597
const double scalings_y_0 = 1;
1598
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1600
// Compute psitilde_a
1601
const double psitilde_a_0 = 1;
1602
const double psitilde_a_1 = x;
1604
// Compute psitilde_bs
1605
const double psitilde_bs_0_0 = 1;
1606
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1607
const double psitilde_bs_1_0 = 1;
1609
// Compute basisvalues
1610
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1611
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1612
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1614
// Table(s) of coefficients
1615
const static double coefficients0[3][3] = \
1616
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
1617
{0.471404520791032, 0.288675134594813, -0.166666666666667},
1618
{0.471404520791032, 0, 0.333333333333333}};
1620
// Interesting (new) part
1621
// Tables of derivatives of the polynomial base (transpose)
1622
const static double dmats0[3][3] = \
1624
{4.89897948556636, 0, 0},
1627
const static double dmats1[3][3] = \
1629
{2.44948974278318, 0, 0},
1630
{4.24264068711928, 0, 0}};
1632
// Compute reference derivatives
1633
// Declare pointer to array of derivatives on FIAT element
1634
double *derivatives = new double [num_derivatives];
1636
// Declare coefficients
1637
double coeff0_0 = 0;
1638
double coeff0_1 = 0;
1639
double coeff0_2 = 0;
1641
// Declare new coefficients
1642
double new_coeff0_0 = 0;
1643
double new_coeff0_1 = 0;
1644
double new_coeff0_2 = 0;
1646
// Loop possible derivatives
1647
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1649
// Get values from coefficients array
1650
new_coeff0_0 = coefficients0[dof][0];
1651
new_coeff0_1 = coefficients0[dof][1];
1652
new_coeff0_2 = coefficients0[dof][2];
1654
// Loop derivative order
1655
for (unsigned int j = 0; j < n; j++)
1657
// Update old coefficients
1658
coeff0_0 = new_coeff0_0;
1659
coeff0_1 = new_coeff0_1;
1660
coeff0_2 = new_coeff0_2;
1662
if(combinations[deriv_num][j] == 0)
1664
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
1665
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
1666
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
1668
if(combinations[deriv_num][j] == 1)
1670
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
1671
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
1672
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
1676
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1677
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
1680
// Transform derivatives back to physical element
1681
for (unsigned int row = 0; row < num_derivatives; row++)
1683
for (unsigned int col = 0; col < num_derivatives; col++)
1685
values[row] += transform[row][col]*derivatives[col];
1688
// Delete pointer to array of derivatives on FIAT element
1689
delete [] derivatives;
1691
// Delete pointer to array of combinations of derivatives and transform
1692
for (unsigned int row = 0; row < num_derivatives; row++)
1694
delete [] combinations[row];
1695
delete [] transform[row];
1698
delete [] combinations;
1699
delete [] transform;
1702
/// Evaluate order n derivatives of all basis functions at given point in cell
1703
virtual void evaluate_basis_derivatives_all(unsigned int n,
1705
const double* coordinates,
1706
const ufc::cell& c) const
1708
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1711
/// Evaluate linear functional for dof i on the function f
1712
virtual double evaluate_dof(unsigned int i,
1713
const ufc::function& f,
1714
const ufc::cell& c) const
1716
// The reference points, direction and weights:
1717
const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
1718
const static double W[3][1] = {{1}, {1}, {1}};
1719
const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
1721
const double * const * x = c.coordinates;
1722
double result = 0.0;
1723
// Iterate over the points:
1724
// Evaluate basis functions for affine mapping
1725
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1726
const double w1 = X[i][0][0];
1727
const double w2 = X[i][0][1];
1729
// Compute affine mapping y = F(X)
1731
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1732
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1734
// Evaluate function at physical points
1736
f.evaluate(values, y, c);
1738
// Map function values using appropriate mapping
1739
// Affine map: Do nothing
1741
// Note that we do not map the weights (yet).
1743
// Take directional components
1744
for(int k = 0; k < 1; k++)
1745
result += values[k]*D[i][0][k];
1746
// Multiply by weights
1752
/// Evaluate linear functionals for all dofs on the function f
1753
virtual void evaluate_dofs(double* values,
1754
const ufc::function& f,
1755
const ufc::cell& c) const
1757
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1760
/// Interpolate vertex values from dof values
1761
virtual void interpolate_vertex_values(double* vertex_values,
1762
const double* dof_values,
1763
const ufc::cell& c) const
1765
// Evaluate at vertices and use affine mapping
1766
vertex_values[0] = dof_values[0];
1767
vertex_values[1] = dof_values[1];
1768
vertex_values[2] = dof_values[2];
1771
/// Return the number of sub elements (for a mixed element)
1772
virtual unsigned int num_sub_elements() const
1777
/// Create a new finite element for sub element i (for a mixed element)
1778
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1780
return new UFC_Poisson2DLinearForm_finite_element_0();
1785
/// This class defines the interface for a finite element.
1787
class UFC_Poisson2DLinearForm_finite_element_1: public ufc::finite_element
1792
UFC_Poisson2DLinearForm_finite_element_1() : ufc::finite_element()
1798
virtual ~UFC_Poisson2DLinearForm_finite_element_1()
1803
/// Return a string identifying the finite element
1804
virtual const char* signature() const
1806
return "Lagrange finite element of degree 1 on a triangle";
1809
/// Return the cell shape
1810
virtual ufc::shape cell_shape() const
1812
return ufc::triangle;
1815
/// Return the dimension of the finite element function space
1816
virtual unsigned int space_dimension() const
1821
/// Return the rank of the value space
1822
virtual unsigned int value_rank() const
1827
/// Return the dimension of the value space for axis i
1828
virtual unsigned int value_dimension(unsigned int i) const
1833
/// Evaluate basis function i at given point in cell
1834
virtual void evaluate_basis(unsigned int i,
1836
const double* coordinates,
1837
const ufc::cell& c) const
1839
// Extract vertex coordinates
1840
const double * const * element_coordinates = c.coordinates;
1842
// Compute Jacobian of affine map from reference cell
1843
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1844
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1845
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1846
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1848
// Compute determinant of Jacobian
1849
const double detJ = J_00*J_11 - J_01*J_10;
1851
// Compute inverse of Jacobian
1853
// Get coordinates and map to the reference (UFC) element
1854
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1855
element_coordinates[0][0]*element_coordinates[2][1] +\
1856
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1857
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1858
element_coordinates[1][0]*element_coordinates[0][1] -\
1859
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1861
// Map coordinates to the reference square
1862
if (std::abs(y - 1.0) < 1e-14)
1865
x = 2.0 *x/(1.0 - y) - 1.0;
1871
// Map degree of freedom to element degree of freedom
1872
const unsigned int dof = i;
1874
// Generate scalings
1875
const double scalings_y_0 = 1;
1876
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1878
// Compute psitilde_a
1879
const double psitilde_a_0 = 1;
1880
const double psitilde_a_1 = x;
1882
// Compute psitilde_bs
1883
const double psitilde_bs_0_0 = 1;
1884
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1885
const double psitilde_bs_1_0 = 1;
1887
// Compute basisvalues
1888
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1889
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1890
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1892
// Table(s) of coefficients
1893
const static double coefficients0[3][3] = \
1894
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
1895
{0.471404520791032, 0.288675134594813, -0.166666666666667},
1896
{0.471404520791032, 0, 0.333333333333333}};
1898
// Extract relevant coefficients
1899
const double coeff0_0 = coefficients0[dof][0];
1900
const double coeff0_1 = coefficients0[dof][1];
1901
const double coeff0_2 = coefficients0[dof][2];
1904
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
1907
/// Evaluate all basis functions at given point in cell
1908
virtual void evaluate_basis_all(double* values,
1909
const double* coordinates,
1910
const ufc::cell& c) const
1912
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1915
/// Evaluate order n derivatives of basis function i at given point in cell
1916
virtual void evaluate_basis_derivatives(unsigned int i,
1919
const double* coordinates,
1920
const ufc::cell& c) const
1922
// Extract vertex coordinates
1923
const double * const * element_coordinates = c.coordinates;
1925
// Compute Jacobian of affine map from reference cell
1926
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1927
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1928
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1929
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1931
// Compute determinant of Jacobian
1932
const double detJ = J_00*J_11 - J_01*J_10;
1934
// Compute inverse of Jacobian
1936
// Get coordinates and map to the reference (UFC) element
1937
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1938
element_coordinates[0][0]*element_coordinates[2][1] +\
1939
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1940
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1941
element_coordinates[1][0]*element_coordinates[0][1] -\
1942
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1944
// Map coordinates to the reference square
1945
if (std::abs(y - 1.0) < 1e-14)
1948
x = 2.0 *x/(1.0 - y) - 1.0;
1951
// Compute number of derivatives
1952
unsigned int num_derivatives = 1;
1954
for (unsigned int j = 0; j < n; j++)
1955
num_derivatives *= 2;
1958
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1959
unsigned int **combinations = new unsigned int *[num_derivatives];
1961
for (unsigned int j = 0; j < num_derivatives; j++)
1963
combinations[j] = new unsigned int [n];
1964
for (unsigned int k = 0; k < n; k++)
1965
combinations[j][k] = 0;
1968
// Generate combinations of derivatives
1969
for (unsigned int row = 1; row < num_derivatives; row++)
1971
for (unsigned int num = 0; num < row; num++)
1973
for (unsigned int col = n-1; col+1 > 0; col--)
1975
if (combinations[row][col] + 1 > 1)
1976
combinations[row][col] = 0;
1979
combinations[row][col] += 1;
1986
// Compute inverse of Jacobian
1987
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1989
// Declare transformation matrix
1990
// Declare pointer to two dimensional array and initialise
1991
double **transform = new double *[num_derivatives];
1993
for (unsigned int j = 0; j < num_derivatives; j++)
1995
transform[j] = new double [num_derivatives];
1996
for (unsigned int k = 0; k < num_derivatives; k++)
1997
transform[j][k] = 1;
2000
// Construct transformation matrix
2001
for (unsigned int row = 0; row < num_derivatives; row++)
2003
for (unsigned int col = 0; col < num_derivatives; col++)
2005
for (unsigned int k = 0; k < n; k++)
2006
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2011
for (unsigned int j = 0; j < 1*num_derivatives; j++)
2014
// Map degree of freedom to element degree of freedom
2015
const unsigned int dof = i;
2017
// Generate scalings
2018
const double scalings_y_0 = 1;
2019
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
2021
// Compute psitilde_a
2022
const double psitilde_a_0 = 1;
2023
const double psitilde_a_1 = x;
2025
// Compute psitilde_bs
2026
const double psitilde_bs_0_0 = 1;
2027
const double psitilde_bs_0_1 = 1.5*y + 0.5;
2028
const double psitilde_bs_1_0 = 1;
2030
// Compute basisvalues
2031
const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2032
const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
2033
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
2035
// Table(s) of coefficients
2036
const static double coefficients0[3][3] = \
2037
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
2038
{0.471404520791032, 0.288675134594813, -0.166666666666667},
2039
{0.471404520791032, 0, 0.333333333333333}};
2041
// Interesting (new) part
2042
// Tables of derivatives of the polynomial base (transpose)
2043
const static double dmats0[3][3] = \
2045
{4.89897948556636, 0, 0},
2048
const static double dmats1[3][3] = \
2050
{2.44948974278318, 0, 0},
2051
{4.24264068711928, 0, 0}};
2053
// Compute reference derivatives
2054
// Declare pointer to array of derivatives on FIAT element
2055
double *derivatives = new double [num_derivatives];
2057
// Declare coefficients
2058
double coeff0_0 = 0;
2059
double coeff0_1 = 0;
2060
double coeff0_2 = 0;
2062
// Declare new coefficients
2063
double new_coeff0_0 = 0;
2064
double new_coeff0_1 = 0;
2065
double new_coeff0_2 = 0;
2067
// Loop possible derivatives
2068
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
2070
// Get values from coefficients array
2071
new_coeff0_0 = coefficients0[dof][0];
2072
new_coeff0_1 = coefficients0[dof][1];
2073
new_coeff0_2 = coefficients0[dof][2];
2075
// Loop derivative order
2076
for (unsigned int j = 0; j < n; j++)
2078
// Update old coefficients
2079
coeff0_0 = new_coeff0_0;
2080
coeff0_1 = new_coeff0_1;
2081
coeff0_2 = new_coeff0_2;
2083
if(combinations[deriv_num][j] == 0)
2085
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
2086
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
2087
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
2089
if(combinations[deriv_num][j] == 1)
2091
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
2092
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
2093
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
2097
// Compute derivatives on reference element as dot product of coefficients and basisvalues
2098
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
2101
// Transform derivatives back to physical element
2102
for (unsigned int row = 0; row < num_derivatives; row++)
2104
for (unsigned int col = 0; col < num_derivatives; col++)
2106
values[row] += transform[row][col]*derivatives[col];
2109
// Delete pointer to array of derivatives on FIAT element
2110
delete [] derivatives;
2112
// Delete pointer to array of combinations of derivatives and transform
2113
for (unsigned int row = 0; row < num_derivatives; row++)
2115
delete [] combinations[row];
2116
delete [] transform[row];
2119
delete [] combinations;
2120
delete [] transform;
2123
/// Evaluate order n derivatives of all basis functions at given point in cell
2124
virtual void evaluate_basis_derivatives_all(unsigned int n,
2126
const double* coordinates,
2127
const ufc::cell& c) const
2129
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
2132
/// Evaluate linear functional for dof i on the function f
2133
virtual double evaluate_dof(unsigned int i,
2134
const ufc::function& f,
2135
const ufc::cell& c) const
2137
// The reference points, direction and weights:
2138
const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
2139
const static double W[3][1] = {{1}, {1}, {1}};
2140
const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
2142
const double * const * x = c.coordinates;
2143
double result = 0.0;
2144
// Iterate over the points:
2145
// Evaluate basis functions for affine mapping
2146
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
2147
const double w1 = X[i][0][0];
2148
const double w2 = X[i][0][1];
2150
// Compute affine mapping y = F(X)
2152
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
2153
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
2155
// Evaluate function at physical points
2157
f.evaluate(values, y, c);
2159
// Map function values using appropriate mapping
2160
// Affine map: Do nothing
2162
// Note that we do not map the weights (yet).
2164
// Take directional components
2165
for(int k = 0; k < 1; k++)
2166
result += values[k]*D[i][0][k];
2167
// Multiply by weights
2173
/// Evaluate linear functionals for all dofs on the function f
2174
virtual void evaluate_dofs(double* values,
2175
const ufc::function& f,
2176
const ufc::cell& c) const
2178
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2181
/// Interpolate vertex values from dof values
2182
virtual void interpolate_vertex_values(double* vertex_values,
2183
const double* dof_values,
2184
const ufc::cell& c) const
2186
// Evaluate at vertices and use affine mapping
2187
vertex_values[0] = dof_values[0];
2188
vertex_values[1] = dof_values[1];
2189
vertex_values[2] = dof_values[2];
2192
/// Return the number of sub elements (for a mixed element)
2193
virtual unsigned int num_sub_elements() const
2198
/// Create a new finite element for sub element i (for a mixed element)
2199
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2201
return new UFC_Poisson2DLinearForm_finite_element_1();
2206
/// This class defines the interface for a local-to-global mapping of
2207
/// degrees of freedom (dofs).
2209
class UFC_Poisson2DLinearForm_dof_map_0: public ufc::dof_map
2213
unsigned int __global_dimension;
2218
UFC_Poisson2DLinearForm_dof_map_0() : ufc::dof_map()
2220
__global_dimension = 0;
2224
virtual ~UFC_Poisson2DLinearForm_dof_map_0()
2229
/// Return a string identifying the dof map
2230
virtual const char* signature() const
2232
return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
2235
/// Return true iff mesh entities of topological dimension d are needed
2236
virtual bool needs_mesh_entities(unsigned int d) const
2253
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2254
virtual bool init_mesh(const ufc::mesh& m)
2256
__global_dimension = m.num_entities[0];
2260
/// Initialize dof map for given cell
2261
virtual void init_cell(const ufc::mesh& m,
2267
/// Finish initialization of dof map for cells
2268
virtual void init_cell_finalize()
2273
/// Return the dimension of the global finite element function space
2274
virtual unsigned int global_dimension() const
2276
return __global_dimension;
2279
/// Return the dimension of the local finite element function space
2280
virtual unsigned int local_dimension() const
2285
// Return the geometric dimension of the coordinates this dof map provides
2286
virtual unsigned int geometric_dimension() const
2291
/// Return the number of dofs on each cell facet
2292
virtual unsigned int num_facet_dofs() const
2297
/// Return the number of dofs associated with each cell entity of dimension d
2298
virtual unsigned int num_entity_dofs(unsigned int d) const
2300
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2303
/// Tabulate the local-to-global mapping of dofs on a cell
2304
virtual void tabulate_dofs(unsigned int* dofs,
2306
const ufc::cell& c) const
2308
dofs[0] = c.entity_indices[0][0];
2309
dofs[1] = c.entity_indices[0][1];
2310
dofs[2] = c.entity_indices[0][2];
2313
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2314
virtual void tabulate_facet_dofs(unsigned int* dofs,
2315
unsigned int facet) const
2334
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2335
virtual void tabulate_entity_dofs(unsigned int* dofs,
2336
unsigned int d, unsigned int i) const
2338
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2341
/// Tabulate the coordinates of all dofs on a cell
2342
virtual void tabulate_coordinates(double** coordinates,
2343
const ufc::cell& c) const
2345
const double * const * x = c.coordinates;
2346
coordinates[0][0] = x[0][0];
2347
coordinates[0][1] = x[0][1];
2348
coordinates[1][0] = x[1][0];
2349
coordinates[1][1] = x[1][1];
2350
coordinates[2][0] = x[2][0];
2351
coordinates[2][1] = x[2][1];
2354
/// Return the number of sub dof maps (for a mixed element)
2355
virtual unsigned int num_sub_dof_maps() const
2360
/// Create a new dof_map for sub dof map i (for a mixed element)
2361
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2363
return new UFC_Poisson2DLinearForm_dof_map_0();
2368
/// This class defines the interface for a local-to-global mapping of
2369
/// degrees of freedom (dofs).
2371
class UFC_Poisson2DLinearForm_dof_map_1: public ufc::dof_map
2375
unsigned int __global_dimension;
2380
UFC_Poisson2DLinearForm_dof_map_1() : ufc::dof_map()
2382
__global_dimension = 0;
2386
virtual ~UFC_Poisson2DLinearForm_dof_map_1()
2391
/// Return a string identifying the dof map
2392
virtual const char* signature() const
2394
return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
2397
/// Return true iff mesh entities of topological dimension d are needed
2398
virtual bool needs_mesh_entities(unsigned int d) const
2415
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2416
virtual bool init_mesh(const ufc::mesh& m)
2418
__global_dimension = m.num_entities[0];
2422
/// Initialize dof map for given cell
2423
virtual void init_cell(const ufc::mesh& m,
2429
/// Finish initialization of dof map for cells
2430
virtual void init_cell_finalize()
2435
/// Return the dimension of the global finite element function space
2436
virtual unsigned int global_dimension() const
2438
return __global_dimension;
2441
/// Return the dimension of the local finite element function space
2442
virtual unsigned int local_dimension() const
2447
// Return the geometric dimension of the coordinates this dof map provides
2448
virtual unsigned int geometric_dimension() const
2453
/// Return the number of dofs on each cell facet
2454
virtual unsigned int num_facet_dofs() const
2459
/// Return the number of dofs associated with each cell entity of dimension d
2460
virtual unsigned int num_entity_dofs(unsigned int d) const
2462
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2465
/// Tabulate the local-to-global mapping of dofs on a cell
2466
virtual void tabulate_dofs(unsigned int* dofs,
2468
const ufc::cell& c) const
2470
dofs[0] = c.entity_indices[0][0];
2471
dofs[1] = c.entity_indices[0][1];
2472
dofs[2] = c.entity_indices[0][2];
2475
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2476
virtual void tabulate_facet_dofs(unsigned int* dofs,
2477
unsigned int facet) const
2496
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2497
virtual void tabulate_entity_dofs(unsigned int* dofs,
2498
unsigned int d, unsigned int i) const
2500
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2503
/// Tabulate the coordinates of all dofs on a cell
2504
virtual void tabulate_coordinates(double** coordinates,
2505
const ufc::cell& c) const
2507
const double * const * x = c.coordinates;
2508
coordinates[0][0] = x[0][0];
2509
coordinates[0][1] = x[0][1];
2510
coordinates[1][0] = x[1][0];
2511
coordinates[1][1] = x[1][1];
2512
coordinates[2][0] = x[2][0];
2513
coordinates[2][1] = x[2][1];
2516
/// Return the number of sub dof maps (for a mixed element)
2517
virtual unsigned int num_sub_dof_maps() const
2522
/// Create a new dof_map for sub dof map i (for a mixed element)
2523
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2525
return new UFC_Poisson2DLinearForm_dof_map_1();
2530
/// This class defines the interface for the tabulation of the cell
2531
/// tensor corresponding to the local contribution to a form from
2532
/// the integral over a cell.
2534
class UFC_Poisson2DLinearForm_cell_integral_0: public ufc::cell_integral
2539
UFC_Poisson2DLinearForm_cell_integral_0() : ufc::cell_integral()
2545
virtual ~UFC_Poisson2DLinearForm_cell_integral_0()
2550
/// Tabulate the tensor for the contribution from a local cell
2551
virtual void tabulate_tensor(double* A,
2552
const double * const * w,
2553
const ufc::cell& c) const
2555
// Extract vertex coordinates
2556
const double * const * x = c.coordinates;
2558
// Compute Jacobian of affine map from reference cell
2559
const double J_00 = x[1][0] - x[0][0];
2560
const double J_01 = x[2][0] - x[0][0];
2561
const double J_10 = x[1][1] - x[0][1];
2562
const double J_11 = x[2][1] - x[0][1];
2564
// Compute determinant of Jacobian
2565
double detJ = J_00*J_11 - J_01*J_10;
2567
// Compute inverse of Jacobian
2570
const double det = std::abs(detJ);
2572
// Compute coefficients
2573
const double c0_0_0_0 = w[0][0];
2574
const double c0_0_0_1 = w[0][1];
2575
const double c0_0_0_2 = w[0][2];
2577
// Compute geometry tensors
2578
const double G0_0 = det*c0_0_0_0;
2579
const double G0_1 = det*c0_0_0_1;
2580
const double G0_2 = det*c0_0_0_2;
2582
// Compute element tensor
2583
A[0] = 0.0833333333333332*G0_0 + 0.0416666666666666*G0_1 + 0.0416666666666666*G0_2;
2584
A[1] = 0.0416666666666666*G0_0 + 0.0833333333333332*G0_1 + 0.0416666666666666*G0_2;
2585
A[2] = 0.0416666666666666*G0_0 + 0.0416666666666666*G0_1 + 0.0833333333333331*G0_2;
2590
/// This class defines the interface for the assembly of the global
2591
/// tensor corresponding to a form with r + n arguments, that is, a
2594
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
2596
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
2597
/// global tensor A is defined by
2599
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
2601
/// where each argument Vj represents the application to the
2602
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
2603
/// fixed functions (coefficients).
2605
class UFC_Poisson2DLinearForm: public ufc::form
2610
UFC_Poisson2DLinearForm() : ufc::form()
2616
virtual ~UFC_Poisson2DLinearForm()
2621
/// Return a string identifying the form
2622
virtual const char* signature() const
2624
return "w0_a0[0, 1, 2] | vi0[0, 1, 2]*va0[0, 1, 2]*dX(0)";
2627
/// Return the rank of the global tensor (r)
2628
virtual unsigned int rank() const
2633
/// Return the number of coefficients (n)
2634
virtual unsigned int num_coefficients() const
2639
/// Return the number of cell integrals
2640
virtual unsigned int num_cell_integrals() const
2645
/// Return the number of exterior facet integrals
2646
virtual unsigned int num_exterior_facet_integrals() const
2651
/// Return the number of interior facet integrals
2652
virtual unsigned int num_interior_facet_integrals() const
2657
/// Create a new finite element for argument function i
2658
virtual ufc::finite_element* create_finite_element(unsigned int i) const
2663
return new UFC_Poisson2DLinearForm_finite_element_0();
2666
return new UFC_Poisson2DLinearForm_finite_element_1();
2672
/// Create a new dof map for argument function i
2673
virtual ufc::dof_map* create_dof_map(unsigned int i) const
2678
return new UFC_Poisson2DLinearForm_dof_map_0();
2681
return new UFC_Poisson2DLinearForm_dof_map_1();
2687
/// Create a new cell integral on sub domain i
2688
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
2690
return new UFC_Poisson2DLinearForm_cell_integral_0();
2693
/// Create a new exterior facet integral on sub domain i
2694
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
2699
/// Create a new interior facet integral on sub domain i
2700
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
2709
#include <dolfin/fem/Form.h>
2711
class Poisson2DBilinearForm : public dolfin::Form
2715
Poisson2DBilinearForm() : dolfin::Form()
2721
virtual const ufc::form& form() const
2726
/// Return array of coefficients
2727
virtual const dolfin::Array<dolfin::Function*>& coefficients() const
2729
return __coefficients;
2735
UFC_Poisson2DBilinearForm __form;
2737
/// Array of coefficients
2738
dolfin::Array<dolfin::Function*> __coefficients;
2742
class Poisson2DLinearForm : public dolfin::Form
2746
Poisson2DLinearForm(dolfin::Function& w0) : dolfin::Form()
2748
__coefficients.push_back(&w0);
2752
virtual const ufc::form& form() const
2757
/// Return array of coefficients
2758
virtual const dolfin::Array<dolfin::Function*>& coefficients() const
2760
return __coefficients;
2766
UFC_Poisson2DLinearForm __form;
2768
/// Array of coefficients
2769
dolfin::Array<dolfin::Function*> __coefficients;
added
removed
sandbox/passembly/Poisson2D.h
