1
// This code conforms with the UFC specification version 1.0
2
// and was automatically generated by FFC version 0.6.2.
4
#ifndef __MIXEDPOISSON_H
5
#define __MIXEDPOISSON_H
11
/// This class defines the interface for a finite element.
13
class mixedpoisson_0_finite_element_0_0: public ufc::finite_element
18
mixedpoisson_0_finite_element_0_0() : ufc::finite_element()
24
virtual ~mixedpoisson_0_finite_element_0_0()
29
/// Return a string identifying the finite element
30
virtual const char* signature() const
32
return "FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1)";
35
/// Return the cell shape
36
virtual ufc::shape cell_shape() const
41
/// Return the dimension of the finite element function space
42
virtual unsigned int space_dimension() const
47
/// Return the rank of the value space
48
virtual unsigned int value_rank() const
53
/// Return the dimension of the value space for axis i
54
virtual unsigned int value_dimension(unsigned int i) const
59
/// Evaluate basis function i at given point in cell
60
virtual void evaluate_basis(unsigned int i,
62
const double* coordinates,
63
const ufc::cell& c) const
65
// Extract vertex coordinates
66
const double * const * element_coordinates = c.coordinates;
68
// Compute Jacobian of affine map from reference cell
69
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
70
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
71
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
72
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
74
// Compute determinant of Jacobian
75
const double detJ = J_00*J_11 - J_01*J_10;
77
// Compute inverse of Jacobian
79
// Get coordinates and map to the reference (UFC) element
80
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
81
element_coordinates[0][0]*element_coordinates[2][1] +\
82
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
83
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
84
element_coordinates[1][0]*element_coordinates[0][1] -\
85
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
87
// Map coordinates to the reference square
88
if (std::abs(y - 1.0) < 1e-08)
91
x = 2.0 *x/(1.0 - y) - 1.0;
98
// Map degree of freedom to element degree of freedom
99
const unsigned int dof = i;
102
const double scalings_y_0 = 1;
103
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
105
// Compute psitilde_a
106
const double psitilde_a_0 = 1;
107
const double psitilde_a_1 = x;
109
// Compute psitilde_bs
110
const double psitilde_bs_0_0 = 1;
111
const double psitilde_bs_0_1 = 1.5*y + 0.5;
112
const double psitilde_bs_1_0 = 1;
114
// Compute basisvalues
115
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
116
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
117
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
119
// Table(s) of coefficients
120
static const double coefficients0[6][3] = \
121
{{0.942809042, 0.577350269, -0.333333333},
122
{-0.471404521, -0.288675135, 0.166666667},
123
{0.471404521, -0.577350269, -0.666666667},
124
{0.471404521, 0.288675135, 0.833333333},
125
{-0.471404521, -0.288675135, 0.166666667},
126
{0.942809042, 0.577350269, -0.333333333}};
128
static const double coefficients1[6][3] = \
129
{{-0.471404521, 0, -0.333333333},
130
{0.942809042, 0, 0.666666667},
131
{0.471404521, 0, 0.333333333},
132
{-0.942809042, 0, -0.666666667},
133
{-0.471404521, 0.866025404, 0.166666667},
134
{-0.471404521, -0.866025404, 0.166666667}};
136
// Extract relevant coefficients
137
const double coeff0_0 = coefficients0[dof][0];
138
const double coeff0_1 = coefficients0[dof][1];
139
const double coeff0_2 = coefficients0[dof][2];
140
const double coeff1_0 = coefficients1[dof][0];
141
const double coeff1_1 = coefficients1[dof][1];
142
const double coeff1_2 = coefficients1[dof][2];
145
const double tmp0_0 = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
146
const double tmp0_1 = coeff1_0*basisvalue0 + coeff1_1*basisvalue1 + coeff1_2*basisvalue2;
147
// Using contravariant Piola transform to map values back to the physical element
148
values[0] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
149
values[1] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
152
/// Evaluate all basis functions at given point in cell
153
virtual void evaluate_basis_all(double* values,
154
const double* coordinates,
155
const ufc::cell& c) const
157
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
160
/// Evaluate order n derivatives of basis function i at given point in cell
161
virtual void evaluate_basis_derivatives(unsigned int i,
164
const double* coordinates,
165
const ufc::cell& c) const
167
// Extract vertex coordinates
168
const double * const * element_coordinates = c.coordinates;
170
// Compute Jacobian of affine map from reference cell
171
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
172
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
173
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
174
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
176
// Compute determinant of Jacobian
177
const double detJ = J_00*J_11 - J_01*J_10;
179
// Compute inverse of Jacobian
181
// Get coordinates and map to the reference (UFC) element
182
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
183
element_coordinates[0][0]*element_coordinates[2][1] +\
184
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
185
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
186
element_coordinates[1][0]*element_coordinates[0][1] -\
187
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
189
// Map coordinates to the reference square
190
if (std::abs(y - 1.0) < 1e-08)
193
x = 2.0 *x/(1.0 - y) - 1.0;
196
// Compute number of derivatives
197
unsigned int num_derivatives = 1;
199
for (unsigned int j = 0; j < n; j++)
200
num_derivatives *= 2;
203
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
204
unsigned int **combinations = new unsigned int *[num_derivatives];
206
for (unsigned int j = 0; j < num_derivatives; j++)
208
combinations[j] = new unsigned int [n];
209
for (unsigned int k = 0; k < n; k++)
210
combinations[j][k] = 0;
213
// Generate combinations of derivatives
214
for (unsigned int row = 1; row < num_derivatives; row++)
216
for (unsigned int num = 0; num < row; num++)
218
for (unsigned int col = n-1; col+1 > 0; col--)
220
if (combinations[row][col] + 1 > 1)
221
combinations[row][col] = 0;
224
combinations[row][col] += 1;
231
// Compute inverse of Jacobian
232
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
234
// Declare transformation matrix
235
// Declare pointer to two dimensional array and initialise
236
double **transform = new double *[num_derivatives];
238
for (unsigned int j = 0; j < num_derivatives; j++)
240
transform[j] = new double [num_derivatives];
241
for (unsigned int k = 0; k < num_derivatives; k++)
245
// Construct transformation matrix
246
for (unsigned int row = 0; row < num_derivatives; row++)
248
for (unsigned int col = 0; col < num_derivatives; col++)
250
for (unsigned int k = 0; k < n; k++)
251
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
256
for (unsigned int j = 0; j < 2*num_derivatives; j++)
259
// Map degree of freedom to element degree of freedom
260
const unsigned int dof = i;
263
const double scalings_y_0 = 1;
264
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
266
// Compute psitilde_a
267
const double psitilde_a_0 = 1;
268
const double psitilde_a_1 = x;
270
// Compute psitilde_bs
271
const double psitilde_bs_0_0 = 1;
272
const double psitilde_bs_0_1 = 1.5*y + 0.5;
273
const double psitilde_bs_1_0 = 1;
275
// Compute basisvalues
276
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
277
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
278
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
280
// Table(s) of coefficients
281
static const double coefficients0[6][3] = \
282
{{0.942809042, 0.577350269, -0.333333333},
283
{-0.471404521, -0.288675135, 0.166666667},
284
{0.471404521, -0.577350269, -0.666666667},
285
{0.471404521, 0.288675135, 0.833333333},
286
{-0.471404521, -0.288675135, 0.166666667},
287
{0.942809042, 0.577350269, -0.333333333}};
289
static const double coefficients1[6][3] = \
290
{{-0.471404521, 0, -0.333333333},
291
{0.942809042, 0, 0.666666667},
292
{0.471404521, 0, 0.333333333},
293
{-0.942809042, 0, -0.666666667},
294
{-0.471404521, 0.866025404, 0.166666667},
295
{-0.471404521, -0.866025404, 0.166666667}};
297
// Interesting (new) part
298
// Tables of derivatives of the polynomial base (transpose)
299
static const double dmats0[3][3] = \
304
static const double dmats1[3][3] = \
309
// Compute reference derivatives
310
// Declare pointer to array of derivatives on FIAT element
311
double *derivatives = new double [2*num_derivatives];
313
// Declare coefficients
321
// Declare new coefficients
322
double new_coeff0_0 = 0;
323
double new_coeff0_1 = 0;
324
double new_coeff0_2 = 0;
325
double new_coeff1_0 = 0;
326
double new_coeff1_1 = 0;
327
double new_coeff1_2 = 0;
329
// Loop possible derivatives
330
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
332
// Get values from coefficients array
333
new_coeff0_0 = coefficients0[dof][0];
334
new_coeff0_1 = coefficients0[dof][1];
335
new_coeff0_2 = coefficients0[dof][2];
336
new_coeff1_0 = coefficients1[dof][0];
337
new_coeff1_1 = coefficients1[dof][1];
338
new_coeff1_2 = coefficients1[dof][2];
340
// Loop derivative order
341
for (unsigned int j = 0; j < n; j++)
343
// Update old coefficients
344
coeff0_0 = new_coeff0_0;
345
coeff0_1 = new_coeff0_1;
346
coeff0_2 = new_coeff0_2;
347
coeff1_0 = new_coeff1_0;
348
coeff1_1 = new_coeff1_1;
349
coeff1_2 = new_coeff1_2;
351
if(combinations[deriv_num][j] == 0)
353
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
354
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
355
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
356
new_coeff1_0 = coeff1_0*dmats0[0][0] + coeff1_1*dmats0[1][0] + coeff1_2*dmats0[2][0];
357
new_coeff1_1 = coeff1_0*dmats0[0][1] + coeff1_1*dmats0[1][1] + coeff1_2*dmats0[2][1];
358
new_coeff1_2 = coeff1_0*dmats0[0][2] + coeff1_1*dmats0[1][2] + coeff1_2*dmats0[2][2];
360
if(combinations[deriv_num][j] == 1)
362
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
363
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
364
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
365
new_coeff1_0 = coeff1_0*dmats1[0][0] + coeff1_1*dmats1[1][0] + coeff1_2*dmats1[2][0];
366
new_coeff1_1 = coeff1_0*dmats1[0][1] + coeff1_1*dmats1[1][1] + coeff1_2*dmats1[2][1];
367
new_coeff1_2 = coeff1_0*dmats1[0][2] + coeff1_1*dmats1[1][2] + coeff1_2*dmats1[2][2];
371
// Compute derivatives on reference element as dot product of coefficients and basisvalues
372
// Correct values by the contravariant Piola transform
373
const double tmp0_0 = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
374
const double tmp0_1 = new_coeff1_0*basisvalue0 + new_coeff1_1*basisvalue1 + new_coeff1_2*basisvalue2;
375
derivatives[deriv_num] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
376
derivatives[num_derivatives + deriv_num] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
379
// Transform derivatives back to physical element
380
for (unsigned int row = 0; row < num_derivatives; row++)
382
for (unsigned int col = 0; col < num_derivatives; col++)
384
values[row] += transform[row][col]*derivatives[col];
385
values[num_derivatives + row] += transform[row][col]*derivatives[num_derivatives + col];
388
// Delete pointer to array of derivatives on FIAT element
389
delete [] derivatives;
391
// Delete pointer to array of combinations of derivatives and transform
392
for (unsigned int row = 0; row < num_derivatives; row++)
394
delete [] combinations[row];
395
delete [] transform[row];
398
delete [] combinations;
402
/// Evaluate order n derivatives of all basis functions at given point in cell
403
virtual void evaluate_basis_derivatives_all(unsigned int n,
405
const double* coordinates,
406
const ufc::cell& c) const
408
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
411
/// Evaluate linear functional for dof i on the function f
412
virtual double evaluate_dof(unsigned int i,
413
const ufc::function& f,
414
const ufc::cell& c) const
416
// The reference points, direction and weights:
417
static const double X[6][1][2] = {{{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}};
418
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
419
static const double D[6][1][2] = {{{1, 1}}, {{1, 1}}, {{1, 0}}, {{1, 0}}, {{0, -1}}, {{0, -1}}};
421
// Extract vertex coordinates
422
const double * const * x = c.coordinates;
424
// Compute Jacobian of affine map from reference cell
425
const double J_00 = x[1][0] - x[0][0];
426
const double J_01 = x[2][0] - x[0][0];
427
const double J_10 = x[1][1] - x[0][1];
428
const double J_11 = x[2][1] - x[0][1];
430
// Compute determinant of Jacobian
431
double detJ = J_00*J_11 - J_01*J_10;
433
// Compute inverse of Jacobian
434
const double Jinv_00 = J_11 / detJ;
435
const double Jinv_01 = -J_01 / detJ;
436
const double Jinv_10 = -J_10 / detJ;
437
const double Jinv_11 = J_00 / detJ;
439
double copyofvalues[2];
441
// Iterate over the points:
442
// Evaluate basis functions for affine mapping
443
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
444
const double w1 = X[i][0][0];
445
const double w2 = X[i][0][1];
447
// Compute affine mapping y = F(X)
449
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
450
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
452
// Evaluate function at physical points
454
f.evaluate(values, y, c);
456
// Map function values using appropriate mapping
458
copyofvalues[0] = values[0];
459
copyofvalues[1] = values[1];
460
// Do the inverse of div piola
461
values[0] = detJ*(Jinv_00*copyofvalues[0]+Jinv_01*copyofvalues[1]);
462
values[1] = detJ*(Jinv_10*copyofvalues[0]+Jinv_11*copyofvalues[1]);
464
// Note that we do not map the weights (yet).
466
// Take directional components
467
for(int k = 0; k < 2; k++)
468
result += values[k]*D[i][0][k];
469
// Multiply by weights
475
/// Evaluate linear functionals for all dofs on the function f
476
virtual void evaluate_dofs(double* values,
477
const ufc::function& f,
478
const ufc::cell& c) const
480
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
483
/// Interpolate vertex values from dof values
484
virtual void interpolate_vertex_values(double* vertex_values,
485
const double* dof_values,
486
const ufc::cell& c) const
488
// Extract vertex coordinates
489
const double * const * x = c.coordinates;
491
// Compute Jacobian of affine map from reference cell
492
const double J_00 = x[1][0] - x[0][0];
493
const double J_01 = x[2][0] - x[0][0];
494
const double J_10 = x[1][1] - x[0][1];
495
const double J_11 = x[2][1] - x[0][1];
497
// Compute determinant of Jacobian
498
double detJ = J_00*J_11 - J_01*J_10;
500
// Compute inverse of Jacobian
501
// Evaluate at vertices and use Piola mapping
502
vertex_values[0] = (1.0/detJ)*(dof_values[2]*2*J_00 + dof_values[3]*J_00 + dof_values[4]*(-2*J_01) + dof_values[5]*J_01);
503
vertex_values[2] = (1.0/detJ)*(dof_values[0]*2*J_00 + dof_values[1]*J_00 + dof_values[4]*(J_00 + J_01) + dof_values[5]*(2*J_00 - 2*J_01));
504
vertex_values[4] = (1.0/detJ)*(dof_values[0]*J_01 + dof_values[1]*2*J_01 + dof_values[2]*(J_00 + J_01) + dof_values[3]*(2*J_00 - 2*J_01));
505
vertex_values[1] = (1.0/detJ)*(dof_values[2]*2*J_10 + dof_values[3]*J_10 + dof_values[4]*(-2*J_11) + dof_values[5]*J_11);
506
vertex_values[3] = (1.0/detJ)*(dof_values[0]*2*J_10 + dof_values[1]*J_10 + dof_values[4]*(J_10 + J_11) + dof_values[5]*(2*J_10 - 2*J_11));
507
vertex_values[5] = (1.0/detJ)*(dof_values[0]*J_11 + dof_values[1]*2*J_11 + dof_values[2]*(J_10 + J_11) + dof_values[3]*(2*J_10 - 2*J_11));
510
/// Return the number of sub elements (for a mixed element)
511
virtual unsigned int num_sub_elements() const
516
/// Create a new finite element for sub element i (for a mixed element)
517
virtual ufc::finite_element* create_sub_element(unsigned int i) const
519
return new mixedpoisson_0_finite_element_0_0();
524
/// This class defines the interface for a finite element.
526
class mixedpoisson_0_finite_element_0_1: public ufc::finite_element
531
mixedpoisson_0_finite_element_0_1() : ufc::finite_element()
537
virtual ~mixedpoisson_0_finite_element_0_1()
542
/// Return a string identifying the finite element
543
virtual const char* signature() const
545
return "FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
548
/// Return the cell shape
549
virtual ufc::shape cell_shape() const
551
return ufc::triangle;
554
/// Return the dimension of the finite element function space
555
virtual unsigned int space_dimension() const
560
/// Return the rank of the value space
561
virtual unsigned int value_rank() const
566
/// Return the dimension of the value space for axis i
567
virtual unsigned int value_dimension(unsigned int i) const
572
/// Evaluate basis function i at given point in cell
573
virtual void evaluate_basis(unsigned int i,
575
const double* coordinates,
576
const ufc::cell& c) const
578
// Extract vertex coordinates
579
const double * const * element_coordinates = c.coordinates;
581
// Compute Jacobian of affine map from reference cell
582
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
583
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
584
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
585
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
587
// Compute determinant of Jacobian
588
const double detJ = J_00*J_11 - J_01*J_10;
590
// Compute inverse of Jacobian
592
// Get coordinates and map to the reference (UFC) element
593
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
594
element_coordinates[0][0]*element_coordinates[2][1] +\
595
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
596
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
597
element_coordinates[1][0]*element_coordinates[0][1] -\
598
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
600
// Map coordinates to the reference square
601
if (std::abs(y - 1.0) < 1e-08)
604
x = 2.0 *x/(1.0 - y) - 1.0;
610
// Map degree of freedom to element degree of freedom
611
const unsigned int dof = i;
614
const double scalings_y_0 = 1;
616
// Compute psitilde_a
617
const double psitilde_a_0 = 1;
619
// Compute psitilde_bs
620
const double psitilde_bs_0_0 = 1;
622
// Compute basisvalues
623
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
625
// Table(s) of coefficients
626
static const double coefficients0[1][1] = \
629
// Extract relevant coefficients
630
const double coeff0_0 = coefficients0[dof][0];
633
*values = coeff0_0*basisvalue0;
636
/// Evaluate all basis functions at given point in cell
637
virtual void evaluate_basis_all(double* values,
638
const double* coordinates,
639
const ufc::cell& c) const
641
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
644
/// Evaluate order n derivatives of basis function i at given point in cell
645
virtual void evaluate_basis_derivatives(unsigned int i,
648
const double* coordinates,
649
const ufc::cell& c) const
651
// Extract vertex coordinates
652
const double * const * element_coordinates = c.coordinates;
654
// Compute Jacobian of affine map from reference cell
655
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
656
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
657
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
658
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
660
// Compute determinant of Jacobian
661
const double detJ = J_00*J_11 - J_01*J_10;
663
// Compute inverse of Jacobian
665
// Get coordinates and map to the reference (UFC) element
666
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
667
element_coordinates[0][0]*element_coordinates[2][1] +\
668
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
669
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
670
element_coordinates[1][0]*element_coordinates[0][1] -\
671
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
673
// Map coordinates to the reference square
674
if (std::abs(y - 1.0) < 1e-08)
677
x = 2.0 *x/(1.0 - y) - 1.0;
680
// Compute number of derivatives
681
unsigned int num_derivatives = 1;
683
for (unsigned int j = 0; j < n; j++)
684
num_derivatives *= 2;
687
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
688
unsigned int **combinations = new unsigned int *[num_derivatives];
690
for (unsigned int j = 0; j < num_derivatives; j++)
692
combinations[j] = new unsigned int [n];
693
for (unsigned int k = 0; k < n; k++)
694
combinations[j][k] = 0;
697
// Generate combinations of derivatives
698
for (unsigned int row = 1; row < num_derivatives; row++)
700
for (unsigned int num = 0; num < row; num++)
702
for (unsigned int col = n-1; col+1 > 0; col--)
704
if (combinations[row][col] + 1 > 1)
705
combinations[row][col] = 0;
708
combinations[row][col] += 1;
715
// Compute inverse of Jacobian
716
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
718
// Declare transformation matrix
719
// Declare pointer to two dimensional array and initialise
720
double **transform = new double *[num_derivatives];
722
for (unsigned int j = 0; j < num_derivatives; j++)
724
transform[j] = new double [num_derivatives];
725
for (unsigned int k = 0; k < num_derivatives; k++)
729
// Construct transformation matrix
730
for (unsigned int row = 0; row < num_derivatives; row++)
732
for (unsigned int col = 0; col < num_derivatives; col++)
734
for (unsigned int k = 0; k < n; k++)
735
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
740
for (unsigned int j = 0; j < 1*num_derivatives; j++)
743
// Map degree of freedom to element degree of freedom
744
const unsigned int dof = i;
747
const double scalings_y_0 = 1;
749
// Compute psitilde_a
750
const double psitilde_a_0 = 1;
752
// Compute psitilde_bs
753
const double psitilde_bs_0_0 = 1;
755
// Compute basisvalues
756
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
758
// Table(s) of coefficients
759
static const double coefficients0[1][1] = \
762
// Interesting (new) part
763
// Tables of derivatives of the polynomial base (transpose)
764
static const double dmats0[1][1] = \
767
static const double dmats1[1][1] = \
770
// Compute reference derivatives
771
// Declare pointer to array of derivatives on FIAT element
772
double *derivatives = new double [num_derivatives];
774
// Declare coefficients
777
// Declare new coefficients
778
double new_coeff0_0 = 0;
780
// Loop possible derivatives
781
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
783
// Get values from coefficients array
784
new_coeff0_0 = coefficients0[dof][0];
786
// Loop derivative order
787
for (unsigned int j = 0; j < n; j++)
789
// Update old coefficients
790
coeff0_0 = new_coeff0_0;
792
if(combinations[deriv_num][j] == 0)
794
new_coeff0_0 = coeff0_0*dmats0[0][0];
796
if(combinations[deriv_num][j] == 1)
798
new_coeff0_0 = coeff0_0*dmats1[0][0];
802
// Compute derivatives on reference element as dot product of coefficients and basisvalues
803
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
806
// Transform derivatives back to physical element
807
for (unsigned int row = 0; row < num_derivatives; row++)
809
for (unsigned int col = 0; col < num_derivatives; col++)
811
values[row] += transform[row][col]*derivatives[col];
814
// Delete pointer to array of derivatives on FIAT element
815
delete [] derivatives;
817
// Delete pointer to array of combinations of derivatives and transform
818
for (unsigned int row = 0; row < num_derivatives; row++)
820
delete [] combinations[row];
821
delete [] transform[row];
824
delete [] combinations;
828
/// Evaluate order n derivatives of all basis functions at given point in cell
829
virtual void evaluate_basis_derivatives_all(unsigned int n,
831
const double* coordinates,
832
const ufc::cell& c) const
834
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
837
/// Evaluate linear functional for dof i on the function f
838
virtual double evaluate_dof(unsigned int i,
839
const ufc::function& f,
840
const ufc::cell& c) const
842
// The reference points, direction and weights:
843
static const double X[1][1][2] = {{{0.333333333, 0.333333333}}};
844
static const double W[1][1] = {{1}};
845
static const double D[1][1][1] = {{{1}}};
847
const double * const * x = c.coordinates;
849
// Iterate over the points:
850
// Evaluate basis functions for affine mapping
851
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
852
const double w1 = X[i][0][0];
853
const double w2 = X[i][0][1];
855
// Compute affine mapping y = F(X)
857
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
858
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
860
// Evaluate function at physical points
862
f.evaluate(values, y, c);
864
// Map function values using appropriate mapping
865
// Affine map: Do nothing
867
// Note that we do not map the weights (yet).
869
// Take directional components
870
for(int k = 0; k < 1; k++)
871
result += values[k]*D[i][0][k];
872
// Multiply by weights
878
/// Evaluate linear functionals for all dofs on the function f
879
virtual void evaluate_dofs(double* values,
880
const ufc::function& f,
881
const ufc::cell& c) const
883
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
886
/// Interpolate vertex values from dof values
887
virtual void interpolate_vertex_values(double* vertex_values,
888
const double* dof_values,
889
const ufc::cell& c) const
891
// Evaluate at vertices and use affine mapping
892
vertex_values[0] = dof_values[0];
893
vertex_values[1] = dof_values[0];
894
vertex_values[2] = dof_values[0];
897
/// Return the number of sub elements (for a mixed element)
898
virtual unsigned int num_sub_elements() const
903
/// Create a new finite element for sub element i (for a mixed element)
904
virtual ufc::finite_element* create_sub_element(unsigned int i) const
906
return new mixedpoisson_0_finite_element_0_1();
911
/// This class defines the interface for a finite element.
913
class mixedpoisson_0_finite_element_0: public ufc::finite_element
918
mixedpoisson_0_finite_element_0() : ufc::finite_element()
924
virtual ~mixedpoisson_0_finite_element_0()
929
/// Return a string identifying the finite element
930
virtual const char* signature() const
932
return "MixedElement([FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1), FiniteElement('Discontinuous Lagrange', 'triangle', 0)])";
935
/// Return the cell shape
936
virtual ufc::shape cell_shape() const
938
return ufc::triangle;
941
/// Return the dimension of the finite element function space
942
virtual unsigned int space_dimension() const
947
/// Return the rank of the value space
948
virtual unsigned int value_rank() const
953
/// Return the dimension of the value space for axis i
954
virtual unsigned int value_dimension(unsigned int i) const
959
/// Evaluate basis function i at given point in cell
960
virtual void evaluate_basis(unsigned int i,
962
const double* coordinates,
963
const ufc::cell& c) const
965
// Extract vertex coordinates
966
const double * const * element_coordinates = c.coordinates;
968
// Compute Jacobian of affine map from reference cell
969
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
970
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
971
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
972
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
974
// Compute determinant of Jacobian
975
const double detJ = J_00*J_11 - J_01*J_10;
977
// Compute inverse of Jacobian
979
// Get coordinates and map to the reference (UFC) element
980
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
981
element_coordinates[0][0]*element_coordinates[2][1] +\
982
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
983
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
984
element_coordinates[1][0]*element_coordinates[0][1] -\
985
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
987
// Map coordinates to the reference square
988
if (std::abs(y - 1.0) < 1e-08)
991
x = 2.0 *x/(1.0 - y) - 1.0;
999
if (0 <= i && i <= 5)
1001
// Map degree of freedom to element degree of freedom
1002
const unsigned int dof = i;
1004
// Generate scalings
1005
const double scalings_y_0 = 1;
1006
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1008
// Compute psitilde_a
1009
const double psitilde_a_0 = 1;
1010
const double psitilde_a_1 = x;
1012
// Compute psitilde_bs
1013
const double psitilde_bs_0_0 = 1;
1014
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1015
const double psitilde_bs_1_0 = 1;
1017
// Compute basisvalues
1018
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1019
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1020
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1022
// Table(s) of coefficients
1023
static const double coefficients0[6][3] = \
1024
{{0.942809042, 0.577350269, -0.333333333},
1025
{-0.471404521, -0.288675135, 0.166666667},
1026
{0.471404521, -0.577350269, -0.666666667},
1027
{0.471404521, 0.288675135, 0.833333333},
1028
{-0.471404521, -0.288675135, 0.166666667},
1029
{0.942809042, 0.577350269, -0.333333333}};
1031
static const double coefficients1[6][3] = \
1032
{{-0.471404521, 0, -0.333333333},
1033
{0.942809042, 0, 0.666666667},
1034
{0.471404521, 0, 0.333333333},
1035
{-0.942809042, 0, -0.666666667},
1036
{-0.471404521, 0.866025404, 0.166666667},
1037
{-0.471404521, -0.866025404, 0.166666667}};
1039
// Extract relevant coefficients
1040
const double coeff0_0 = coefficients0[dof][0];
1041
const double coeff0_1 = coefficients0[dof][1];
1042
const double coeff0_2 = coefficients0[dof][2];
1043
const double coeff1_0 = coefficients1[dof][0];
1044
const double coeff1_1 = coefficients1[dof][1];
1045
const double coeff1_2 = coefficients1[dof][2];
1048
const double tmp0_0 = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
1049
const double tmp0_1 = coeff1_0*basisvalue0 + coeff1_1*basisvalue1 + coeff1_2*basisvalue2;
1050
// Using contravariant Piola transform to map values back to the physical element
1051
values[0] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
1052
values[1] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
1055
if (6 <= i && i <= 6)
1057
// Map degree of freedom to element degree of freedom
1058
const unsigned int dof = i - 6;
1060
// Generate scalings
1061
const double scalings_y_0 = 1;
1063
// Compute psitilde_a
1064
const double psitilde_a_0 = 1;
1066
// Compute psitilde_bs
1067
const double psitilde_bs_0_0 = 1;
1069
// Compute basisvalues
1070
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1072
// Table(s) of coefficients
1073
static const double coefficients0[1][1] = \
1076
// Extract relevant coefficients
1077
const double coeff0_0 = coefficients0[dof][0];
1080
values[2] = coeff0_0*basisvalue0;
1085
/// Evaluate all basis functions at given point in cell
1086
virtual void evaluate_basis_all(double* values,
1087
const double* coordinates,
1088
const ufc::cell& c) const
1090
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1093
/// Evaluate order n derivatives of basis function i at given point in cell
1094
virtual void evaluate_basis_derivatives(unsigned int i,
1097
const double* coordinates,
1098
const ufc::cell& c) const
1100
// Extract vertex coordinates
1101
const double * const * element_coordinates = c.coordinates;
1103
// Compute Jacobian of affine map from reference cell
1104
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1105
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1106
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1107
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1109
// Compute determinant of Jacobian
1110
const double detJ = J_00*J_11 - J_01*J_10;
1112
// Compute inverse of Jacobian
1114
// Get coordinates and map to the reference (UFC) element
1115
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1116
element_coordinates[0][0]*element_coordinates[2][1] +\
1117
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1118
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1119
element_coordinates[1][0]*element_coordinates[0][1] -\
1120
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1122
// Map coordinates to the reference square
1123
if (std::abs(y - 1.0) < 1e-08)
1126
x = 2.0 *x/(1.0 - y) - 1.0;
1129
// Compute number of derivatives
1130
unsigned int num_derivatives = 1;
1132
for (unsigned int j = 0; j < n; j++)
1133
num_derivatives *= 2;
1136
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1137
unsigned int **combinations = new unsigned int *[num_derivatives];
1139
for (unsigned int j = 0; j < num_derivatives; j++)
1141
combinations[j] = new unsigned int [n];
1142
for (unsigned int k = 0; k < n; k++)
1143
combinations[j][k] = 0;
1146
// Generate combinations of derivatives
1147
for (unsigned int row = 1; row < num_derivatives; row++)
1149
for (unsigned int num = 0; num < row; num++)
1151
for (unsigned int col = n-1; col+1 > 0; col--)
1153
if (combinations[row][col] + 1 > 1)
1154
combinations[row][col] = 0;
1157
combinations[row][col] += 1;
1164
// Compute inverse of Jacobian
1165
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1167
// Declare transformation matrix
1168
// Declare pointer to two dimensional array and initialise
1169
double **transform = new double *[num_derivatives];
1171
for (unsigned int j = 0; j < num_derivatives; j++)
1173
transform[j] = new double [num_derivatives];
1174
for (unsigned int k = 0; k < num_derivatives; k++)
1175
transform[j][k] = 1;
1178
// Construct transformation matrix
1179
for (unsigned int row = 0; row < num_derivatives; row++)
1181
for (unsigned int col = 0; col < num_derivatives; col++)
1183
for (unsigned int k = 0; k < n; k++)
1184
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1189
for (unsigned int j = 0; j < 3*num_derivatives; j++)
1192
if (0 <= i && i <= 5)
1194
// Map degree of freedom to element degree of freedom
1195
const unsigned int dof = i;
1197
// Generate scalings
1198
const double scalings_y_0 = 1;
1199
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1201
// Compute psitilde_a
1202
const double psitilde_a_0 = 1;
1203
const double psitilde_a_1 = x;
1205
// Compute psitilde_bs
1206
const double psitilde_bs_0_0 = 1;
1207
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1208
const double psitilde_bs_1_0 = 1;
1210
// Compute basisvalues
1211
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1212
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1213
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1215
// Table(s) of coefficients
1216
static const double coefficients0[6][3] = \
1217
{{0.942809042, 0.577350269, -0.333333333},
1218
{-0.471404521, -0.288675135, 0.166666667},
1219
{0.471404521, -0.577350269, -0.666666667},
1220
{0.471404521, 0.288675135, 0.833333333},
1221
{-0.471404521, -0.288675135, 0.166666667},
1222
{0.942809042, 0.577350269, -0.333333333}};
1224
static const double coefficients1[6][3] = \
1225
{{-0.471404521, 0, -0.333333333},
1226
{0.942809042, 0, 0.666666667},
1227
{0.471404521, 0, 0.333333333},
1228
{-0.942809042, 0, -0.666666667},
1229
{-0.471404521, 0.866025404, 0.166666667},
1230
{-0.471404521, -0.866025404, 0.166666667}};
1232
// Interesting (new) part
1233
// Tables of derivatives of the polynomial base (transpose)
1234
static const double dmats0[3][3] = \
1239
static const double dmats1[3][3] = \
1242
{4.24264069, 0, 0}};
1244
// Compute reference derivatives
1245
// Declare pointer to array of derivatives on FIAT element
1246
double *derivatives = new double [2*num_derivatives];
1248
// Declare coefficients
1249
double coeff0_0 = 0;
1250
double coeff0_1 = 0;
1251
double coeff0_2 = 0;
1252
double coeff1_0 = 0;
1253
double coeff1_1 = 0;
1254
double coeff1_2 = 0;
1256
// Declare new coefficients
1257
double new_coeff0_0 = 0;
1258
double new_coeff0_1 = 0;
1259
double new_coeff0_2 = 0;
1260
double new_coeff1_0 = 0;
1261
double new_coeff1_1 = 0;
1262
double new_coeff1_2 = 0;
1264
// Loop possible derivatives
1265
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1267
// Get values from coefficients array
1268
new_coeff0_0 = coefficients0[dof][0];
1269
new_coeff0_1 = coefficients0[dof][1];
1270
new_coeff0_2 = coefficients0[dof][2];
1271
new_coeff1_0 = coefficients1[dof][0];
1272
new_coeff1_1 = coefficients1[dof][1];
1273
new_coeff1_2 = coefficients1[dof][2];
1275
// Loop derivative order
1276
for (unsigned int j = 0; j < n; j++)
1278
// Update old coefficients
1279
coeff0_0 = new_coeff0_0;
1280
coeff0_1 = new_coeff0_1;
1281
coeff0_2 = new_coeff0_2;
1282
coeff1_0 = new_coeff1_0;
1283
coeff1_1 = new_coeff1_1;
1284
coeff1_2 = new_coeff1_2;
1286
if(combinations[deriv_num][j] == 0)
1288
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
1289
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
1290
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
1291
new_coeff1_0 = coeff1_0*dmats0[0][0] + coeff1_1*dmats0[1][0] + coeff1_2*dmats0[2][0];
1292
new_coeff1_1 = coeff1_0*dmats0[0][1] + coeff1_1*dmats0[1][1] + coeff1_2*dmats0[2][1];
1293
new_coeff1_2 = coeff1_0*dmats0[0][2] + coeff1_1*dmats0[1][2] + coeff1_2*dmats0[2][2];
1295
if(combinations[deriv_num][j] == 1)
1297
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
1298
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
1299
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
1300
new_coeff1_0 = coeff1_0*dmats1[0][0] + coeff1_1*dmats1[1][0] + coeff1_2*dmats1[2][0];
1301
new_coeff1_1 = coeff1_0*dmats1[0][1] + coeff1_1*dmats1[1][1] + coeff1_2*dmats1[2][1];
1302
new_coeff1_2 = coeff1_0*dmats1[0][2] + coeff1_1*dmats1[1][2] + coeff1_2*dmats1[2][2];
1306
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1307
// Correct values by the contravariant Piola transform
1308
const double tmp0_0 = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
1309
const double tmp0_1 = new_coeff1_0*basisvalue0 + new_coeff1_1*basisvalue1 + new_coeff1_2*basisvalue2;
1310
derivatives[deriv_num] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
1311
derivatives[num_derivatives + deriv_num] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
1314
// Transform derivatives back to physical element
1315
for (unsigned int row = 0; row < num_derivatives; row++)
1317
for (unsigned int col = 0; col < num_derivatives; col++)
1319
values[row] += transform[row][col]*derivatives[col];
1320
values[num_derivatives + row] += transform[row][col]*derivatives[num_derivatives + col];
1323
// Delete pointer to array of derivatives on FIAT element
1324
delete [] derivatives;
1326
// Delete pointer to array of combinations of derivatives and transform
1327
for (unsigned int row = 0; row < num_derivatives; row++)
1329
delete [] combinations[row];
1330
delete [] transform[row];
1333
delete [] combinations;
1334
delete [] transform;
1337
if (6 <= i && i <= 6)
1339
// Map degree of freedom to element degree of freedom
1340
const unsigned int dof = i - 6;
1342
// Generate scalings
1343
const double scalings_y_0 = 1;
1345
// Compute psitilde_a
1346
const double psitilde_a_0 = 1;
1348
// Compute psitilde_bs
1349
const double psitilde_bs_0_0 = 1;
1351
// Compute basisvalues
1352
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1354
// Table(s) of coefficients
1355
static const double coefficients0[1][1] = \
1358
// Interesting (new) part
1359
// Tables of derivatives of the polynomial base (transpose)
1360
static const double dmats0[1][1] = \
1363
static const double dmats1[1][1] = \
1366
// Compute reference derivatives
1367
// Declare pointer to array of derivatives on FIAT element
1368
double *derivatives = new double [num_derivatives];
1370
// Declare coefficients
1371
double coeff0_0 = 0;
1373
// Declare new coefficients
1374
double new_coeff0_0 = 0;
1376
// Loop possible derivatives
1377
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1379
// Get values from coefficients array
1380
new_coeff0_0 = coefficients0[dof][0];
1382
// Loop derivative order
1383
for (unsigned int j = 0; j < n; j++)
1385
// Update old coefficients
1386
coeff0_0 = new_coeff0_0;
1388
if(combinations[deriv_num][j] == 0)
1390
new_coeff0_0 = coeff0_0*dmats0[0][0];
1392
if(combinations[deriv_num][j] == 1)
1394
new_coeff0_0 = coeff0_0*dmats1[0][0];
1398
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1399
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
1402
// Transform derivatives back to physical element
1403
for (unsigned int row = 0; row < num_derivatives; row++)
1405
for (unsigned int col = 0; col < num_derivatives; col++)
1407
values[2*num_derivatives + row] += transform[row][col]*derivatives[col];
1410
// Delete pointer to array of derivatives on FIAT element
1411
delete [] derivatives;
1413
// Delete pointer to array of combinations of derivatives and transform
1414
for (unsigned int row = 0; row < num_derivatives; row++)
1416
delete [] combinations[row];
1417
delete [] transform[row];
1420
delete [] combinations;
1421
delete [] transform;
1426
/// Evaluate order n derivatives of all basis functions at given point in cell
1427
virtual void evaluate_basis_derivatives_all(unsigned int n,
1429
const double* coordinates,
1430
const ufc::cell& c) const
1432
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1435
/// Evaluate linear functional for dof i on the function f
1436
virtual double evaluate_dof(unsigned int i,
1437
const ufc::function& f,
1438
const ufc::cell& c) const
1440
// The reference points, direction and weights:
1441
static const double X[7][1][2] = {{{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}, {{0.333333333, 0.333333333}}};
1442
static const double W[7][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}};
1443
static const double D[7][1][3] = {{{1, 1, 0}}, {{1, 1, 0}}, {{1, 0, 0}}, {{1, 0, 0}}, {{0, -1, 0}}, {{0, -1, 0}}, {{0, 0, 1}}};
1445
static const unsigned int mappings[7] = {1, 1, 1, 1, 1, 1, 0};
1446
// Extract vertex coordinates
1447
const double * const * x = c.coordinates;
1449
// Compute Jacobian of affine map from reference cell
1450
const double J_00 = x[1][0] - x[0][0];
1451
const double J_01 = x[2][0] - x[0][0];
1452
const double J_10 = x[1][1] - x[0][1];
1453
const double J_11 = x[2][1] - x[0][1];
1455
// Compute determinant of Jacobian
1456
double detJ = J_00*J_11 - J_01*J_10;
1458
// Compute inverse of Jacobian
1459
const double Jinv_00 = J_11 / detJ;
1460
const double Jinv_01 = -J_01 / detJ;
1461
const double Jinv_10 = -J_10 / detJ;
1462
const double Jinv_11 = J_00 / detJ;
1464
double copyofvalues[3];
1465
double result = 0.0;
1466
// Iterate over the points:
1467
// Evaluate basis functions for affine mapping
1468
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1469
const double w1 = X[i][0][0];
1470
const double w2 = X[i][0][1];
1472
// Compute affine mapping y = F(X)
1474
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1475
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1477
// Evaluate function at physical points
1479
f.evaluate(values, y, c);
1481
// Map function values using appropriate mapping
1482
if (mappings[i] == 0) {
1483
// Affine map: Do nothing
1484
} else if (mappings[i] == 1) {
1486
copyofvalues[0] = values[0];
1487
copyofvalues[1] = values[1];
1488
// Do the inverse of div piola
1489
values[0] = detJ*(Jinv_00*copyofvalues[0]+Jinv_01*copyofvalues[1]);
1490
values[1] = detJ*(Jinv_10*copyofvalues[0]+Jinv_11*copyofvalues[1]);
1492
// Other mappings not applicable.
1495
// Note that we do not map the weights (yet).
1497
// Take directional components
1498
for(int k = 0; k < 3; k++)
1499
result += values[k]*D[i][0][k];
1500
// Multiply by weights
1506
/// Evaluate linear functionals for all dofs on the function f
1507
virtual void evaluate_dofs(double* values,
1508
const ufc::function& f,
1509
const ufc::cell& c) const
1511
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1514
/// Interpolate vertex values from dof values
1515
virtual void interpolate_vertex_values(double* vertex_values,
1516
const double* dof_values,
1517
const ufc::cell& c) const
1519
// Extract vertex coordinates
1520
const double * const * x = c.coordinates;
1522
// Compute Jacobian of affine map from reference cell
1523
const double J_00 = x[1][0] - x[0][0];
1524
const double J_01 = x[2][0] - x[0][0];
1525
const double J_10 = x[1][1] - x[0][1];
1526
const double J_11 = x[2][1] - x[0][1];
1528
// Compute determinant of Jacobian
1529
double detJ = J_00*J_11 - J_01*J_10;
1531
// Compute inverse of Jacobian
1532
// Evaluate at vertices and use Piola mapping
1533
vertex_values[0] = (1.0/detJ)*(dof_values[2]*2*J_00 + dof_values[3]*J_00 + dof_values[4]*(-2*J_01) + dof_values[5]*J_01);
1534
vertex_values[3] = (1.0/detJ)*(dof_values[0]*2*J_00 + dof_values[1]*J_00 + dof_values[4]*(J_00 + J_01) + dof_values[5]*(2*J_00 - 2*J_01));
1535
vertex_values[6] = (1.0/detJ)*(dof_values[0]*J_01 + dof_values[1]*2*J_01 + dof_values[2]*(J_00 + J_01) + dof_values[3]*(2*J_00 - 2*J_01));
1536
vertex_values[1] = (1.0/detJ)*(dof_values[2]*2*J_10 + dof_values[3]*J_10 + dof_values[4]*(-2*J_11) + dof_values[5]*J_11);
1537
vertex_values[4] = (1.0/detJ)*(dof_values[0]*2*J_10 + dof_values[1]*J_10 + dof_values[4]*(J_10 + J_11) + dof_values[5]*(2*J_10 - 2*J_11));
1538
vertex_values[7] = (1.0/detJ)*(dof_values[0]*J_11 + dof_values[1]*2*J_11 + dof_values[2]*(J_10 + J_11) + dof_values[3]*(2*J_10 - 2*J_11));
1539
// Evaluate at vertices and use affine mapping
1540
vertex_values[2] = dof_values[6];
1541
vertex_values[5] = dof_values[6];
1542
vertex_values[8] = dof_values[6];
1545
/// Return the number of sub elements (for a mixed element)
1546
virtual unsigned int num_sub_elements() const
1551
/// Create a new finite element for sub element i (for a mixed element)
1552
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1557
return new mixedpoisson_0_finite_element_0_0();
1560
return new mixedpoisson_0_finite_element_0_1();
1568
/// This class defines the interface for a finite element.
1570
class mixedpoisson_0_finite_element_1_0: public ufc::finite_element
1575
mixedpoisson_0_finite_element_1_0() : ufc::finite_element()
1581
virtual ~mixedpoisson_0_finite_element_1_0()
1586
/// Return a string identifying the finite element
1587
virtual const char* signature() const
1589
return "FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1)";
1592
/// Return the cell shape
1593
virtual ufc::shape cell_shape() const
1595
return ufc::triangle;
1598
/// Return the dimension of the finite element function space
1599
virtual unsigned int space_dimension() const
1604
/// Return the rank of the value space
1605
virtual unsigned int value_rank() const
1610
/// Return the dimension of the value space for axis i
1611
virtual unsigned int value_dimension(unsigned int i) const
1616
/// Evaluate basis function i at given point in cell
1617
virtual void evaluate_basis(unsigned int i,
1619
const double* coordinates,
1620
const ufc::cell& c) const
1622
// Extract vertex coordinates
1623
const double * const * element_coordinates = c.coordinates;
1625
// Compute Jacobian of affine map from reference cell
1626
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1627
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1628
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1629
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1631
// Compute determinant of Jacobian
1632
const double detJ = J_00*J_11 - J_01*J_10;
1634
// Compute inverse of Jacobian
1636
// Get coordinates and map to the reference (UFC) element
1637
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1638
element_coordinates[0][0]*element_coordinates[2][1] +\
1639
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1640
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1641
element_coordinates[1][0]*element_coordinates[0][1] -\
1642
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1644
// Map coordinates to the reference square
1645
if (std::abs(y - 1.0) < 1e-08)
1648
x = 2.0 *x/(1.0 - y) - 1.0;
1655
// Map degree of freedom to element degree of freedom
1656
const unsigned int dof = i;
1658
// Generate scalings
1659
const double scalings_y_0 = 1;
1660
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1662
// Compute psitilde_a
1663
const double psitilde_a_0 = 1;
1664
const double psitilde_a_1 = x;
1666
// Compute psitilde_bs
1667
const double psitilde_bs_0_0 = 1;
1668
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1669
const double psitilde_bs_1_0 = 1;
1671
// Compute basisvalues
1672
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1673
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1674
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1676
// Table(s) of coefficients
1677
static const double coefficients0[6][3] = \
1678
{{0.942809042, 0.577350269, -0.333333333},
1679
{-0.471404521, -0.288675135, 0.166666667},
1680
{0.471404521, -0.577350269, -0.666666667},
1681
{0.471404521, 0.288675135, 0.833333333},
1682
{-0.471404521, -0.288675135, 0.166666667},
1683
{0.942809042, 0.577350269, -0.333333333}};
1685
static const double coefficients1[6][3] = \
1686
{{-0.471404521, 0, -0.333333333},
1687
{0.942809042, 0, 0.666666667},
1688
{0.471404521, 0, 0.333333333},
1689
{-0.942809042, 0, -0.666666667},
1690
{-0.471404521, 0.866025404, 0.166666667},
1691
{-0.471404521, -0.866025404, 0.166666667}};
1693
// Extract relevant coefficients
1694
const double coeff0_0 = coefficients0[dof][0];
1695
const double coeff0_1 = coefficients0[dof][1];
1696
const double coeff0_2 = coefficients0[dof][2];
1697
const double coeff1_0 = coefficients1[dof][0];
1698
const double coeff1_1 = coefficients1[dof][1];
1699
const double coeff1_2 = coefficients1[dof][2];
1702
const double tmp0_0 = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
1703
const double tmp0_1 = coeff1_0*basisvalue0 + coeff1_1*basisvalue1 + coeff1_2*basisvalue2;
1704
// Using contravariant Piola transform to map values back to the physical element
1705
values[0] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
1706
values[1] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
1709
/// Evaluate all basis functions at given point in cell
1710
virtual void evaluate_basis_all(double* values,
1711
const double* coordinates,
1712
const ufc::cell& c) const
1714
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1717
/// Evaluate order n derivatives of basis function i at given point in cell
1718
virtual void evaluate_basis_derivatives(unsigned int i,
1721
const double* coordinates,
1722
const ufc::cell& c) const
1724
// Extract vertex coordinates
1725
const double * const * element_coordinates = c.coordinates;
1727
// Compute Jacobian of affine map from reference cell
1728
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1729
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1730
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1731
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1733
// Compute determinant of Jacobian
1734
const double detJ = J_00*J_11 - J_01*J_10;
1736
// Compute inverse of Jacobian
1738
// Get coordinates and map to the reference (UFC) element
1739
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1740
element_coordinates[0][0]*element_coordinates[2][1] +\
1741
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1742
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1743
element_coordinates[1][0]*element_coordinates[0][1] -\
1744
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1746
// Map coordinates to the reference square
1747
if (std::abs(y - 1.0) < 1e-08)
1750
x = 2.0 *x/(1.0 - y) - 1.0;
1753
// Compute number of derivatives
1754
unsigned int num_derivatives = 1;
1756
for (unsigned int j = 0; j < n; j++)
1757
num_derivatives *= 2;
1760
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1761
unsigned int **combinations = new unsigned int *[num_derivatives];
1763
for (unsigned int j = 0; j < num_derivatives; j++)
1765
combinations[j] = new unsigned int [n];
1766
for (unsigned int k = 0; k < n; k++)
1767
combinations[j][k] = 0;
1770
// Generate combinations of derivatives
1771
for (unsigned int row = 1; row < num_derivatives; row++)
1773
for (unsigned int num = 0; num < row; num++)
1775
for (unsigned int col = n-1; col+1 > 0; col--)
1777
if (combinations[row][col] + 1 > 1)
1778
combinations[row][col] = 0;
1781
combinations[row][col] += 1;
1788
// Compute inverse of Jacobian
1789
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1791
// Declare transformation matrix
1792
// Declare pointer to two dimensional array and initialise
1793
double **transform = new double *[num_derivatives];
1795
for (unsigned int j = 0; j < num_derivatives; j++)
1797
transform[j] = new double [num_derivatives];
1798
for (unsigned int k = 0; k < num_derivatives; k++)
1799
transform[j][k] = 1;
1802
// Construct transformation matrix
1803
for (unsigned int row = 0; row < num_derivatives; row++)
1805
for (unsigned int col = 0; col < num_derivatives; col++)
1807
for (unsigned int k = 0; k < n; k++)
1808
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1813
for (unsigned int j = 0; j < 2*num_derivatives; j++)
1816
// Map degree of freedom to element degree of freedom
1817
const unsigned int dof = i;
1819
// Generate scalings
1820
const double scalings_y_0 = 1;
1821
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1823
// Compute psitilde_a
1824
const double psitilde_a_0 = 1;
1825
const double psitilde_a_1 = x;
1827
// Compute psitilde_bs
1828
const double psitilde_bs_0_0 = 1;
1829
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1830
const double psitilde_bs_1_0 = 1;
1832
// Compute basisvalues
1833
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1834
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1835
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1837
// Table(s) of coefficients
1838
static const double coefficients0[6][3] = \
1839
{{0.942809042, 0.577350269, -0.333333333},
1840
{-0.471404521, -0.288675135, 0.166666667},
1841
{0.471404521, -0.577350269, -0.666666667},
1842
{0.471404521, 0.288675135, 0.833333333},
1843
{-0.471404521, -0.288675135, 0.166666667},
1844
{0.942809042, 0.577350269, -0.333333333}};
1846
static const double coefficients1[6][3] = \
1847
{{-0.471404521, 0, -0.333333333},
1848
{0.942809042, 0, 0.666666667},
1849
{0.471404521, 0, 0.333333333},
1850
{-0.942809042, 0, -0.666666667},
1851
{-0.471404521, 0.866025404, 0.166666667},
1852
{-0.471404521, -0.866025404, 0.166666667}};
1854
// Interesting (new) part
1855
// Tables of derivatives of the polynomial base (transpose)
1856
static const double dmats0[3][3] = \
1861
static const double dmats1[3][3] = \
1864
{4.24264069, 0, 0}};
1866
// Compute reference derivatives
1867
// Declare pointer to array of derivatives on FIAT element
1868
double *derivatives = new double [2*num_derivatives];
1870
// Declare coefficients
1871
double coeff0_0 = 0;
1872
double coeff0_1 = 0;
1873
double coeff0_2 = 0;
1874
double coeff1_0 = 0;
1875
double coeff1_1 = 0;
1876
double coeff1_2 = 0;
1878
// Declare new coefficients
1879
double new_coeff0_0 = 0;
1880
double new_coeff0_1 = 0;
1881
double new_coeff0_2 = 0;
1882
double new_coeff1_0 = 0;
1883
double new_coeff1_1 = 0;
1884
double new_coeff1_2 = 0;
1886
// Loop possible derivatives
1887
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1889
// Get values from coefficients array
1890
new_coeff0_0 = coefficients0[dof][0];
1891
new_coeff0_1 = coefficients0[dof][1];
1892
new_coeff0_2 = coefficients0[dof][2];
1893
new_coeff1_0 = coefficients1[dof][0];
1894
new_coeff1_1 = coefficients1[dof][1];
1895
new_coeff1_2 = coefficients1[dof][2];
1897
// Loop derivative order
1898
for (unsigned int j = 0; j < n; j++)
1900
// Update old coefficients
1901
coeff0_0 = new_coeff0_0;
1902
coeff0_1 = new_coeff0_1;
1903
coeff0_2 = new_coeff0_2;
1904
coeff1_0 = new_coeff1_0;
1905
coeff1_1 = new_coeff1_1;
1906
coeff1_2 = new_coeff1_2;
1908
if(combinations[deriv_num][j] == 0)
1910
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
1911
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
1912
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
1913
new_coeff1_0 = coeff1_0*dmats0[0][0] + coeff1_1*dmats0[1][0] + coeff1_2*dmats0[2][0];
1914
new_coeff1_1 = coeff1_0*dmats0[0][1] + coeff1_1*dmats0[1][1] + coeff1_2*dmats0[2][1];
1915
new_coeff1_2 = coeff1_0*dmats0[0][2] + coeff1_1*dmats0[1][2] + coeff1_2*dmats0[2][2];
1917
if(combinations[deriv_num][j] == 1)
1919
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
1920
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
1921
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
1922
new_coeff1_0 = coeff1_0*dmats1[0][0] + coeff1_1*dmats1[1][0] + coeff1_2*dmats1[2][0];
1923
new_coeff1_1 = coeff1_0*dmats1[0][1] + coeff1_1*dmats1[1][1] + coeff1_2*dmats1[2][1];
1924
new_coeff1_2 = coeff1_0*dmats1[0][2] + coeff1_1*dmats1[1][2] + coeff1_2*dmats1[2][2];
1928
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1929
// Correct values by the contravariant Piola transform
1930
const double tmp0_0 = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
1931
const double tmp0_1 = new_coeff1_0*basisvalue0 + new_coeff1_1*basisvalue1 + new_coeff1_2*basisvalue2;
1932
derivatives[deriv_num] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
1933
derivatives[num_derivatives + deriv_num] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
1936
// Transform derivatives back to physical element
1937
for (unsigned int row = 0; row < num_derivatives; row++)
1939
for (unsigned int col = 0; col < num_derivatives; col++)
1941
values[row] += transform[row][col]*derivatives[col];
1942
values[num_derivatives + row] += transform[row][col]*derivatives[num_derivatives + col];
1945
// Delete pointer to array of derivatives on FIAT element
1946
delete [] derivatives;
1948
// Delete pointer to array of combinations of derivatives and transform
1949
for (unsigned int row = 0; row < num_derivatives; row++)
1951
delete [] combinations[row];
1952
delete [] transform[row];
1955
delete [] combinations;
1956
delete [] transform;
1959
/// Evaluate order n derivatives of all basis functions at given point in cell
1960
virtual void evaluate_basis_derivatives_all(unsigned int n,
1962
const double* coordinates,
1963
const ufc::cell& c) const
1965
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1968
/// Evaluate linear functional for dof i on the function f
1969
virtual double evaluate_dof(unsigned int i,
1970
const ufc::function& f,
1971
const ufc::cell& c) const
1973
// The reference points, direction and weights:
1974
static const double X[6][1][2] = {{{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}};
1975
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
1976
static const double D[6][1][2] = {{{1, 1}}, {{1, 1}}, {{1, 0}}, {{1, 0}}, {{0, -1}}, {{0, -1}}};
1978
// Extract vertex coordinates
1979
const double * const * x = c.coordinates;
1981
// Compute Jacobian of affine map from reference cell
1982
const double J_00 = x[1][0] - x[0][0];
1983
const double J_01 = x[2][0] - x[0][0];
1984
const double J_10 = x[1][1] - x[0][1];
1985
const double J_11 = x[2][1] - x[0][1];
1987
// Compute determinant of Jacobian
1988
double detJ = J_00*J_11 - J_01*J_10;
1990
// Compute inverse of Jacobian
1991
const double Jinv_00 = J_11 / detJ;
1992
const double Jinv_01 = -J_01 / detJ;
1993
const double Jinv_10 = -J_10 / detJ;
1994
const double Jinv_11 = J_00 / detJ;
1996
double copyofvalues[2];
1997
double result = 0.0;
1998
// Iterate over the points:
1999
// Evaluate basis functions for affine mapping
2000
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
2001
const double w1 = X[i][0][0];
2002
const double w2 = X[i][0][1];
2004
// Compute affine mapping y = F(X)
2006
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
2007
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
2009
// Evaluate function at physical points
2011
f.evaluate(values, y, c);
2013
// Map function values using appropriate mapping
2015
copyofvalues[0] = values[0];
2016
copyofvalues[1] = values[1];
2017
// Do the inverse of div piola
2018
values[0] = detJ*(Jinv_00*copyofvalues[0]+Jinv_01*copyofvalues[1]);
2019
values[1] = detJ*(Jinv_10*copyofvalues[0]+Jinv_11*copyofvalues[1]);
2021
// Note that we do not map the weights (yet).
2023
// Take directional components
2024
for(int k = 0; k < 2; k++)
2025
result += values[k]*D[i][0][k];
2026
// Multiply by weights
2032
/// Evaluate linear functionals for all dofs on the function f
2033
virtual void evaluate_dofs(double* values,
2034
const ufc::function& f,
2035
const ufc::cell& c) const
2037
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2040
/// Interpolate vertex values from dof values
2041
virtual void interpolate_vertex_values(double* vertex_values,
2042
const double* dof_values,
2043
const ufc::cell& c) const
2045
// Extract vertex coordinates
2046
const double * const * x = c.coordinates;
2048
// Compute Jacobian of affine map from reference cell
2049
const double J_00 = x[1][0] - x[0][0];
2050
const double J_01 = x[2][0] - x[0][0];
2051
const double J_10 = x[1][1] - x[0][1];
2052
const double J_11 = x[2][1] - x[0][1];
2054
// Compute determinant of Jacobian
2055
double detJ = J_00*J_11 - J_01*J_10;
2057
// Compute inverse of Jacobian
2058
// Evaluate at vertices and use Piola mapping
2059
vertex_values[0] = (1.0/detJ)*(dof_values[2]*2*J_00 + dof_values[3]*J_00 + dof_values[4]*(-2*J_01) + dof_values[5]*J_01);
2060
vertex_values[2] = (1.0/detJ)*(dof_values[0]*2*J_00 + dof_values[1]*J_00 + dof_values[4]*(J_00 + J_01) + dof_values[5]*(2*J_00 - 2*J_01));
2061
vertex_values[4] = (1.0/detJ)*(dof_values[0]*J_01 + dof_values[1]*2*J_01 + dof_values[2]*(J_00 + J_01) + dof_values[3]*(2*J_00 - 2*J_01));
2062
vertex_values[1] = (1.0/detJ)*(dof_values[2]*2*J_10 + dof_values[3]*J_10 + dof_values[4]*(-2*J_11) + dof_values[5]*J_11);
2063
vertex_values[3] = (1.0/detJ)*(dof_values[0]*2*J_10 + dof_values[1]*J_10 + dof_values[4]*(J_10 + J_11) + dof_values[5]*(2*J_10 - 2*J_11));
2064
vertex_values[5] = (1.0/detJ)*(dof_values[0]*J_11 + dof_values[1]*2*J_11 + dof_values[2]*(J_10 + J_11) + dof_values[3]*(2*J_10 - 2*J_11));
2067
/// Return the number of sub elements (for a mixed element)
2068
virtual unsigned int num_sub_elements() const
2073
/// Create a new finite element for sub element i (for a mixed element)
2074
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2076
return new mixedpoisson_0_finite_element_1_0();
2081
/// This class defines the interface for a finite element.
2083
class mixedpoisson_0_finite_element_1_1: public ufc::finite_element
2088
mixedpoisson_0_finite_element_1_1() : ufc::finite_element()
2094
virtual ~mixedpoisson_0_finite_element_1_1()
2099
/// Return a string identifying the finite element
2100
virtual const char* signature() const
2102
return "FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
2105
/// Return the cell shape
2106
virtual ufc::shape cell_shape() const
2108
return ufc::triangle;
2111
/// Return the dimension of the finite element function space
2112
virtual unsigned int space_dimension() const
2117
/// Return the rank of the value space
2118
virtual unsigned int value_rank() const
2123
/// Return the dimension of the value space for axis i
2124
virtual unsigned int value_dimension(unsigned int i) const
2129
/// Evaluate basis function i at given point in cell
2130
virtual void evaluate_basis(unsigned int i,
2132
const double* coordinates,
2133
const ufc::cell& c) const
2135
// Extract vertex coordinates
2136
const double * const * element_coordinates = c.coordinates;
2138
// Compute Jacobian of affine map from reference cell
2139
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
2140
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
2141
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
2142
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
2144
// Compute determinant of Jacobian
2145
const double detJ = J_00*J_11 - J_01*J_10;
2147
// Compute inverse of Jacobian
2149
// Get coordinates and map to the reference (UFC) element
2150
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
2151
element_coordinates[0][0]*element_coordinates[2][1] +\
2152
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
2153
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
2154
element_coordinates[1][0]*element_coordinates[0][1] -\
2155
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
2157
// Map coordinates to the reference square
2158
if (std::abs(y - 1.0) < 1e-08)
2161
x = 2.0 *x/(1.0 - y) - 1.0;
2167
// Map degree of freedom to element degree of freedom
2168
const unsigned int dof = i;
2170
// Generate scalings
2171
const double scalings_y_0 = 1;
2173
// Compute psitilde_a
2174
const double psitilde_a_0 = 1;
2176
// Compute psitilde_bs
2177
const double psitilde_bs_0_0 = 1;
2179
// Compute basisvalues
2180
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2182
// Table(s) of coefficients
2183
static const double coefficients0[1][1] = \
2186
// Extract relevant coefficients
2187
const double coeff0_0 = coefficients0[dof][0];
2190
*values = coeff0_0*basisvalue0;
2193
/// Evaluate all basis functions at given point in cell
2194
virtual void evaluate_basis_all(double* values,
2195
const double* coordinates,
2196
const ufc::cell& c) const
2198
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
2201
/// Evaluate order n derivatives of basis function i at given point in cell
2202
virtual void evaluate_basis_derivatives(unsigned int i,
2205
const double* coordinates,
2206
const ufc::cell& c) const
2208
// Extract vertex coordinates
2209
const double * const * element_coordinates = c.coordinates;
2211
// Compute Jacobian of affine map from reference cell
2212
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
2213
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
2214
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
2215
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
2217
// Compute determinant of Jacobian
2218
const double detJ = J_00*J_11 - J_01*J_10;
2220
// Compute inverse of Jacobian
2222
// Get coordinates and map to the reference (UFC) element
2223
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
2224
element_coordinates[0][0]*element_coordinates[2][1] +\
2225
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
2226
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
2227
element_coordinates[1][0]*element_coordinates[0][1] -\
2228
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
2230
// Map coordinates to the reference square
2231
if (std::abs(y - 1.0) < 1e-08)
2234
x = 2.0 *x/(1.0 - y) - 1.0;
2237
// Compute number of derivatives
2238
unsigned int num_derivatives = 1;
2240
for (unsigned int j = 0; j < n; j++)
2241
num_derivatives *= 2;
2244
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
2245
unsigned int **combinations = new unsigned int *[num_derivatives];
2247
for (unsigned int j = 0; j < num_derivatives; j++)
2249
combinations[j] = new unsigned int [n];
2250
for (unsigned int k = 0; k < n; k++)
2251
combinations[j][k] = 0;
2254
// Generate combinations of derivatives
2255
for (unsigned int row = 1; row < num_derivatives; row++)
2257
for (unsigned int num = 0; num < row; num++)
2259
for (unsigned int col = n-1; col+1 > 0; col--)
2261
if (combinations[row][col] + 1 > 1)
2262
combinations[row][col] = 0;
2265
combinations[row][col] += 1;
2272
// Compute inverse of Jacobian
2273
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
2275
// Declare transformation matrix
2276
// Declare pointer to two dimensional array and initialise
2277
double **transform = new double *[num_derivatives];
2279
for (unsigned int j = 0; j < num_derivatives; j++)
2281
transform[j] = new double [num_derivatives];
2282
for (unsigned int k = 0; k < num_derivatives; k++)
2283
transform[j][k] = 1;
2286
// Construct transformation matrix
2287
for (unsigned int row = 0; row < num_derivatives; row++)
2289
for (unsigned int col = 0; col < num_derivatives; col++)
2291
for (unsigned int k = 0; k < n; k++)
2292
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2297
for (unsigned int j = 0; j < 1*num_derivatives; j++)
2300
// Map degree of freedom to element degree of freedom
2301
const unsigned int dof = i;
2303
// Generate scalings
2304
const double scalings_y_0 = 1;
2306
// Compute psitilde_a
2307
const double psitilde_a_0 = 1;
2309
// Compute psitilde_bs
2310
const double psitilde_bs_0_0 = 1;
2312
// Compute basisvalues
2313
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2315
// Table(s) of coefficients
2316
static const double coefficients0[1][1] = \
2319
// Interesting (new) part
2320
// Tables of derivatives of the polynomial base (transpose)
2321
static const double dmats0[1][1] = \
2324
static const double dmats1[1][1] = \
2327
// Compute reference derivatives
2328
// Declare pointer to array of derivatives on FIAT element
2329
double *derivatives = new double [num_derivatives];
2331
// Declare coefficients
2332
double coeff0_0 = 0;
2334
// Declare new coefficients
2335
double new_coeff0_0 = 0;
2337
// Loop possible derivatives
2338
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
2340
// Get values from coefficients array
2341
new_coeff0_0 = coefficients0[dof][0];
2343
// Loop derivative order
2344
for (unsigned int j = 0; j < n; j++)
2346
// Update old coefficients
2347
coeff0_0 = new_coeff0_0;
2349
if(combinations[deriv_num][j] == 0)
2351
new_coeff0_0 = coeff0_0*dmats0[0][0];
2353
if(combinations[deriv_num][j] == 1)
2355
new_coeff0_0 = coeff0_0*dmats1[0][0];
2359
// Compute derivatives on reference element as dot product of coefficients and basisvalues
2360
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
2363
// Transform derivatives back to physical element
2364
for (unsigned int row = 0; row < num_derivatives; row++)
2366
for (unsigned int col = 0; col < num_derivatives; col++)
2368
values[row] += transform[row][col]*derivatives[col];
2371
// Delete pointer to array of derivatives on FIAT element
2372
delete [] derivatives;
2374
// Delete pointer to array of combinations of derivatives and transform
2375
for (unsigned int row = 0; row < num_derivatives; row++)
2377
delete [] combinations[row];
2378
delete [] transform[row];
2381
delete [] combinations;
2382
delete [] transform;
2385
/// Evaluate order n derivatives of all basis functions at given point in cell
2386
virtual void evaluate_basis_derivatives_all(unsigned int n,
2388
const double* coordinates,
2389
const ufc::cell& c) const
2391
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
2394
/// Evaluate linear functional for dof i on the function f
2395
virtual double evaluate_dof(unsigned int i,
2396
const ufc::function& f,
2397
const ufc::cell& c) const
2399
// The reference points, direction and weights:
2400
static const double X[1][1][2] = {{{0.333333333, 0.333333333}}};
2401
static const double W[1][1] = {{1}};
2402
static const double D[1][1][1] = {{{1}}};
2404
const double * const * x = c.coordinates;
2405
double result = 0.0;
2406
// Iterate over the points:
2407
// Evaluate basis functions for affine mapping
2408
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
2409
const double w1 = X[i][0][0];
2410
const double w2 = X[i][0][1];
2412
// Compute affine mapping y = F(X)
2414
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
2415
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
2417
// Evaluate function at physical points
2419
f.evaluate(values, y, c);
2421
// Map function values using appropriate mapping
2422
// Affine map: Do nothing
2424
// Note that we do not map the weights (yet).
2426
// Take directional components
2427
for(int k = 0; k < 1; k++)
2428
result += values[k]*D[i][0][k];
2429
// Multiply by weights
2435
/// Evaluate linear functionals for all dofs on the function f
2436
virtual void evaluate_dofs(double* values,
2437
const ufc::function& f,
2438
const ufc::cell& c) const
2440
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2443
/// Interpolate vertex values from dof values
2444
virtual void interpolate_vertex_values(double* vertex_values,
2445
const double* dof_values,
2446
const ufc::cell& c) const
2448
// Evaluate at vertices and use affine mapping
2449
vertex_values[0] = dof_values[0];
2450
vertex_values[1] = dof_values[0];
2451
vertex_values[2] = dof_values[0];
2454
/// Return the number of sub elements (for a mixed element)
2455
virtual unsigned int num_sub_elements() const
2460
/// Create a new finite element for sub element i (for a mixed element)
2461
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2463
return new mixedpoisson_0_finite_element_1_1();
2468
/// This class defines the interface for a finite element.
2470
class mixedpoisson_0_finite_element_1: public ufc::finite_element
2475
mixedpoisson_0_finite_element_1() : ufc::finite_element()
2481
virtual ~mixedpoisson_0_finite_element_1()
2486
/// Return a string identifying the finite element
2487
virtual const char* signature() const
2489
return "MixedElement([FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1), FiniteElement('Discontinuous Lagrange', 'triangle', 0)])";
2492
/// Return the cell shape
2493
virtual ufc::shape cell_shape() const
2495
return ufc::triangle;
2498
/// Return the dimension of the finite element function space
2499
virtual unsigned int space_dimension() const
2504
/// Return the rank of the value space
2505
virtual unsigned int value_rank() const
2510
/// Return the dimension of the value space for axis i
2511
virtual unsigned int value_dimension(unsigned int i) const
2516
/// Evaluate basis function i at given point in cell
2517
virtual void evaluate_basis(unsigned int i,
2519
const double* coordinates,
2520
const ufc::cell& c) const
2522
// Extract vertex coordinates
2523
const double * const * element_coordinates = c.coordinates;
2525
// Compute Jacobian of affine map from reference cell
2526
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
2527
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
2528
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
2529
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
2531
// Compute determinant of Jacobian
2532
const double detJ = J_00*J_11 - J_01*J_10;
2534
// Compute inverse of Jacobian
2536
// Get coordinates and map to the reference (UFC) element
2537
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
2538
element_coordinates[0][0]*element_coordinates[2][1] +\
2539
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
2540
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
2541
element_coordinates[1][0]*element_coordinates[0][1] -\
2542
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
2544
// Map coordinates to the reference square
2545
if (std::abs(y - 1.0) < 1e-08)
2548
x = 2.0 *x/(1.0 - y) - 1.0;
2556
if (0 <= i && i <= 5)
2558
// Map degree of freedom to element degree of freedom
2559
const unsigned int dof = i;
2561
// Generate scalings
2562
const double scalings_y_0 = 1;
2563
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
2565
// Compute psitilde_a
2566
const double psitilde_a_0 = 1;
2567
const double psitilde_a_1 = x;
2569
// Compute psitilde_bs
2570
const double psitilde_bs_0_0 = 1;
2571
const double psitilde_bs_0_1 = 1.5*y + 0.5;
2572
const double psitilde_bs_1_0 = 1;
2574
// Compute basisvalues
2575
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2576
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
2577
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
2579
// Table(s) of coefficients
2580
static const double coefficients0[6][3] = \
2581
{{0.942809042, 0.577350269, -0.333333333},
2582
{-0.471404521, -0.288675135, 0.166666667},
2583
{0.471404521, -0.577350269, -0.666666667},
2584
{0.471404521, 0.288675135, 0.833333333},
2585
{-0.471404521, -0.288675135, 0.166666667},
2586
{0.942809042, 0.577350269, -0.333333333}};
2588
static const double coefficients1[6][3] = \
2589
{{-0.471404521, 0, -0.333333333},
2590
{0.942809042, 0, 0.666666667},
2591
{0.471404521, 0, 0.333333333},
2592
{-0.942809042, 0, -0.666666667},
2593
{-0.471404521, 0.866025404, 0.166666667},
2594
{-0.471404521, -0.866025404, 0.166666667}};
2596
// Extract relevant coefficients
2597
const double coeff0_0 = coefficients0[dof][0];
2598
const double coeff0_1 = coefficients0[dof][1];
2599
const double coeff0_2 = coefficients0[dof][2];
2600
const double coeff1_0 = coefficients1[dof][0];
2601
const double coeff1_1 = coefficients1[dof][1];
2602
const double coeff1_2 = coefficients1[dof][2];
2605
const double tmp0_0 = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
2606
const double tmp0_1 = coeff1_0*basisvalue0 + coeff1_1*basisvalue1 + coeff1_2*basisvalue2;
2607
// Using contravariant Piola transform to map values back to the physical element
2608
values[0] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
2609
values[1] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
2612
if (6 <= i && i <= 6)
2614
// Map degree of freedom to element degree of freedom
2615
const unsigned int dof = i - 6;
2617
// Generate scalings
2618
const double scalings_y_0 = 1;
2620
// Compute psitilde_a
2621
const double psitilde_a_0 = 1;
2623
// Compute psitilde_bs
2624
const double psitilde_bs_0_0 = 1;
2626
// Compute basisvalues
2627
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2629
// Table(s) of coefficients
2630
static const double coefficients0[1][1] = \
2633
// Extract relevant coefficients
2634
const double coeff0_0 = coefficients0[dof][0];
2637
values[2] = coeff0_0*basisvalue0;
2642
/// Evaluate all basis functions at given point in cell
2643
virtual void evaluate_basis_all(double* values,
2644
const double* coordinates,
2645
const ufc::cell& c) const
2647
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
2650
/// Evaluate order n derivatives of basis function i at given point in cell
2651
virtual void evaluate_basis_derivatives(unsigned int i,
2654
const double* coordinates,
2655
const ufc::cell& c) const
2657
// Extract vertex coordinates
2658
const double * const * element_coordinates = c.coordinates;
2660
// Compute Jacobian of affine map from reference cell
2661
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
2662
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
2663
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
2664
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
2666
// Compute determinant of Jacobian
2667
const double detJ = J_00*J_11 - J_01*J_10;
2669
// Compute inverse of Jacobian
2671
// Get coordinates and map to the reference (UFC) element
2672
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
2673
element_coordinates[0][0]*element_coordinates[2][1] +\
2674
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
2675
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
2676
element_coordinates[1][0]*element_coordinates[0][1] -\
2677
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
2679
// Map coordinates to the reference square
2680
if (std::abs(y - 1.0) < 1e-08)
2683
x = 2.0 *x/(1.0 - y) - 1.0;
2686
// Compute number of derivatives
2687
unsigned int num_derivatives = 1;
2689
for (unsigned int j = 0; j < n; j++)
2690
num_derivatives *= 2;
2693
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
2694
unsigned int **combinations = new unsigned int *[num_derivatives];
2696
for (unsigned int j = 0; j < num_derivatives; j++)
2698
combinations[j] = new unsigned int [n];
2699
for (unsigned int k = 0; k < n; k++)
2700
combinations[j][k] = 0;
2703
// Generate combinations of derivatives
2704
for (unsigned int row = 1; row < num_derivatives; row++)
2706
for (unsigned int num = 0; num < row; num++)
2708
for (unsigned int col = n-1; col+1 > 0; col--)
2710
if (combinations[row][col] + 1 > 1)
2711
combinations[row][col] = 0;
2714
combinations[row][col] += 1;
2721
// Compute inverse of Jacobian
2722
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
2724
// Declare transformation matrix
2725
// Declare pointer to two dimensional array and initialise
2726
double **transform = new double *[num_derivatives];
2728
for (unsigned int j = 0; j < num_derivatives; j++)
2730
transform[j] = new double [num_derivatives];
2731
for (unsigned int k = 0; k < num_derivatives; k++)
2732
transform[j][k] = 1;
2735
// Construct transformation matrix
2736
for (unsigned int row = 0; row < num_derivatives; row++)
2738
for (unsigned int col = 0; col < num_derivatives; col++)
2740
for (unsigned int k = 0; k < n; k++)
2741
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2746
for (unsigned int j = 0; j < 3*num_derivatives; j++)
2749
if (0 <= i && i <= 5)
2751
// Map degree of freedom to element degree of freedom
2752
const unsigned int dof = i;
2754
// Generate scalings
2755
const double scalings_y_0 = 1;
2756
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
2758
// Compute psitilde_a
2759
const double psitilde_a_0 = 1;
2760
const double psitilde_a_1 = x;
2762
// Compute psitilde_bs
2763
const double psitilde_bs_0_0 = 1;
2764
const double psitilde_bs_0_1 = 1.5*y + 0.5;
2765
const double psitilde_bs_1_0 = 1;
2767
// Compute basisvalues
2768
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2769
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
2770
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
2772
// Table(s) of coefficients
2773
static const double coefficients0[6][3] = \
2774
{{0.942809042, 0.577350269, -0.333333333},
2775
{-0.471404521, -0.288675135, 0.166666667},
2776
{0.471404521, -0.577350269, -0.666666667},
2777
{0.471404521, 0.288675135, 0.833333333},
2778
{-0.471404521, -0.288675135, 0.166666667},
2779
{0.942809042, 0.577350269, -0.333333333}};
2781
static const double coefficients1[6][3] = \
2782
{{-0.471404521, 0, -0.333333333},
2783
{0.942809042, 0, 0.666666667},
2784
{0.471404521, 0, 0.333333333},
2785
{-0.942809042, 0, -0.666666667},
2786
{-0.471404521, 0.866025404, 0.166666667},
2787
{-0.471404521, -0.866025404, 0.166666667}};
2789
// Interesting (new) part
2790
// Tables of derivatives of the polynomial base (transpose)
2791
static const double dmats0[3][3] = \
2796
static const double dmats1[3][3] = \
2799
{4.24264069, 0, 0}};
2801
// Compute reference derivatives
2802
// Declare pointer to array of derivatives on FIAT element
2803
double *derivatives = new double [2*num_derivatives];
2805
// Declare coefficients
2806
double coeff0_0 = 0;
2807
double coeff0_1 = 0;
2808
double coeff0_2 = 0;
2809
double coeff1_0 = 0;
2810
double coeff1_1 = 0;
2811
double coeff1_2 = 0;
2813
// Declare new coefficients
2814
double new_coeff0_0 = 0;
2815
double new_coeff0_1 = 0;
2816
double new_coeff0_2 = 0;
2817
double new_coeff1_0 = 0;
2818
double new_coeff1_1 = 0;
2819
double new_coeff1_2 = 0;
2821
// Loop possible derivatives
2822
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
2824
// Get values from coefficients array
2825
new_coeff0_0 = coefficients0[dof][0];
2826
new_coeff0_1 = coefficients0[dof][1];
2827
new_coeff0_2 = coefficients0[dof][2];
2828
new_coeff1_0 = coefficients1[dof][0];
2829
new_coeff1_1 = coefficients1[dof][1];
2830
new_coeff1_2 = coefficients1[dof][2];
2832
// Loop derivative order
2833
for (unsigned int j = 0; j < n; j++)
2835
// Update old coefficients
2836
coeff0_0 = new_coeff0_0;
2837
coeff0_1 = new_coeff0_1;
2838
coeff0_2 = new_coeff0_2;
2839
coeff1_0 = new_coeff1_0;
2840
coeff1_1 = new_coeff1_1;
2841
coeff1_2 = new_coeff1_2;
2843
if(combinations[deriv_num][j] == 0)
2845
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
2846
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
2847
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
2848
new_coeff1_0 = coeff1_0*dmats0[0][0] + coeff1_1*dmats0[1][0] + coeff1_2*dmats0[2][0];
2849
new_coeff1_1 = coeff1_0*dmats0[0][1] + coeff1_1*dmats0[1][1] + coeff1_2*dmats0[2][1];
2850
new_coeff1_2 = coeff1_0*dmats0[0][2] + coeff1_1*dmats0[1][2] + coeff1_2*dmats0[2][2];
2852
if(combinations[deriv_num][j] == 1)
2854
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
2855
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
2856
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
2857
new_coeff1_0 = coeff1_0*dmats1[0][0] + coeff1_1*dmats1[1][0] + coeff1_2*dmats1[2][0];
2858
new_coeff1_1 = coeff1_0*dmats1[0][1] + coeff1_1*dmats1[1][1] + coeff1_2*dmats1[2][1];
2859
new_coeff1_2 = coeff1_0*dmats1[0][2] + coeff1_1*dmats1[1][2] + coeff1_2*dmats1[2][2];
2863
// Compute derivatives on reference element as dot product of coefficients and basisvalues
2864
// Correct values by the contravariant Piola transform
2865
const double tmp0_0 = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
2866
const double tmp0_1 = new_coeff1_0*basisvalue0 + new_coeff1_1*basisvalue1 + new_coeff1_2*basisvalue2;
2867
derivatives[deriv_num] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
2868
derivatives[num_derivatives + deriv_num] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
2871
// Transform derivatives back to physical element
2872
for (unsigned int row = 0; row < num_derivatives; row++)
2874
for (unsigned int col = 0; col < num_derivatives; col++)
2876
values[row] += transform[row][col]*derivatives[col];
2877
values[num_derivatives + row] += transform[row][col]*derivatives[num_derivatives + col];
2880
// Delete pointer to array of derivatives on FIAT element
2881
delete [] derivatives;
2883
// Delete pointer to array of combinations of derivatives and transform
2884
for (unsigned int row = 0; row < num_derivatives; row++)
2886
delete [] combinations[row];
2887
delete [] transform[row];
2890
delete [] combinations;
2891
delete [] transform;
2894
if (6 <= i && i <= 6)
2896
// Map degree of freedom to element degree of freedom
2897
const unsigned int dof = i - 6;
2899
// Generate scalings
2900
const double scalings_y_0 = 1;
2902
// Compute psitilde_a
2903
const double psitilde_a_0 = 1;
2905
// Compute psitilde_bs
2906
const double psitilde_bs_0_0 = 1;
2908
// Compute basisvalues
2909
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2911
// Table(s) of coefficients
2912
static const double coefficients0[1][1] = \
2915
// Interesting (new) part
2916
// Tables of derivatives of the polynomial base (transpose)
2917
static const double dmats0[1][1] = \
2920
static const double dmats1[1][1] = \
2923
// Compute reference derivatives
2924
// Declare pointer to array of derivatives on FIAT element
2925
double *derivatives = new double [num_derivatives];
2927
// Declare coefficients
2928
double coeff0_0 = 0;
2930
// Declare new coefficients
2931
double new_coeff0_0 = 0;
2933
// Loop possible derivatives
2934
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
2936
// Get values from coefficients array
2937
new_coeff0_0 = coefficients0[dof][0];
2939
// Loop derivative order
2940
for (unsigned int j = 0; j < n; j++)
2942
// Update old coefficients
2943
coeff0_0 = new_coeff0_0;
2945
if(combinations[deriv_num][j] == 0)
2947
new_coeff0_0 = coeff0_0*dmats0[0][0];
2949
if(combinations[deriv_num][j] == 1)
2951
new_coeff0_0 = coeff0_0*dmats1[0][0];
2955
// Compute derivatives on reference element as dot product of coefficients and basisvalues
2956
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
2959
// Transform derivatives back to physical element
2960
for (unsigned int row = 0; row < num_derivatives; row++)
2962
for (unsigned int col = 0; col < num_derivatives; col++)
2964
values[2*num_derivatives + row] += transform[row][col]*derivatives[col];
2967
// Delete pointer to array of derivatives on FIAT element
2968
delete [] derivatives;
2970
// Delete pointer to array of combinations of derivatives and transform
2971
for (unsigned int row = 0; row < num_derivatives; row++)
2973
delete [] combinations[row];
2974
delete [] transform[row];
2977
delete [] combinations;
2978
delete [] transform;
2983
/// Evaluate order n derivatives of all basis functions at given point in cell
2984
virtual void evaluate_basis_derivatives_all(unsigned int n,
2986
const double* coordinates,
2987
const ufc::cell& c) const
2989
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
2992
/// Evaluate linear functional for dof i on the function f
2993
virtual double evaluate_dof(unsigned int i,
2994
const ufc::function& f,
2995
const ufc::cell& c) const
2997
// The reference points, direction and weights:
2998
static const double X[7][1][2] = {{{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}, {{0.333333333, 0.333333333}}};
2999
static const double W[7][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}};
3000
static const double D[7][1][3] = {{{1, 1, 0}}, {{1, 1, 0}}, {{1, 0, 0}}, {{1, 0, 0}}, {{0, -1, 0}}, {{0, -1, 0}}, {{0, 0, 1}}};
3002
static const unsigned int mappings[7] = {1, 1, 1, 1, 1, 1, 0};
3003
// Extract vertex coordinates
3004
const double * const * x = c.coordinates;
3006
// Compute Jacobian of affine map from reference cell
3007
const double J_00 = x[1][0] - x[0][0];
3008
const double J_01 = x[2][0] - x[0][0];
3009
const double J_10 = x[1][1] - x[0][1];
3010
const double J_11 = x[2][1] - x[0][1];
3012
// Compute determinant of Jacobian
3013
double detJ = J_00*J_11 - J_01*J_10;
3015
// Compute inverse of Jacobian
3016
const double Jinv_00 = J_11 / detJ;
3017
const double Jinv_01 = -J_01 / detJ;
3018
const double Jinv_10 = -J_10 / detJ;
3019
const double Jinv_11 = J_00 / detJ;
3021
double copyofvalues[3];
3022
double result = 0.0;
3023
// Iterate over the points:
3024
// Evaluate basis functions for affine mapping
3025
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
3026
const double w1 = X[i][0][0];
3027
const double w2 = X[i][0][1];
3029
// Compute affine mapping y = F(X)
3031
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
3032
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
3034
// Evaluate function at physical points
3036
f.evaluate(values, y, c);
3038
// Map function values using appropriate mapping
3039
if (mappings[i] == 0) {
3040
// Affine map: Do nothing
3041
} else if (mappings[i] == 1) {
3043
copyofvalues[0] = values[0];
3044
copyofvalues[1] = values[1];
3045
// Do the inverse of div piola
3046
values[0] = detJ*(Jinv_00*copyofvalues[0]+Jinv_01*copyofvalues[1]);
3047
values[1] = detJ*(Jinv_10*copyofvalues[0]+Jinv_11*copyofvalues[1]);
3049
// Other mappings not applicable.
3052
// Note that we do not map the weights (yet).
3054
// Take directional components
3055
for(int k = 0; k < 3; k++)
3056
result += values[k]*D[i][0][k];
3057
// Multiply by weights
3063
/// Evaluate linear functionals for all dofs on the function f
3064
virtual void evaluate_dofs(double* values,
3065
const ufc::function& f,
3066
const ufc::cell& c) const
3068
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3071
/// Interpolate vertex values from dof values
3072
virtual void interpolate_vertex_values(double* vertex_values,
3073
const double* dof_values,
3074
const ufc::cell& c) const
3076
// Extract vertex coordinates
3077
const double * const * x = c.coordinates;
3079
// Compute Jacobian of affine map from reference cell
3080
const double J_00 = x[1][0] - x[0][0];
3081
const double J_01 = x[2][0] - x[0][0];
3082
const double J_10 = x[1][1] - x[0][1];
3083
const double J_11 = x[2][1] - x[0][1];
3085
// Compute determinant of Jacobian
3086
double detJ = J_00*J_11 - J_01*J_10;
3088
// Compute inverse of Jacobian
3089
// Evaluate at vertices and use Piola mapping
3090
vertex_values[0] = (1.0/detJ)*(dof_values[2]*2*J_00 + dof_values[3]*J_00 + dof_values[4]*(-2*J_01) + dof_values[5]*J_01);
3091
vertex_values[3] = (1.0/detJ)*(dof_values[0]*2*J_00 + dof_values[1]*J_00 + dof_values[4]*(J_00 + J_01) + dof_values[5]*(2*J_00 - 2*J_01));
3092
vertex_values[6] = (1.0/detJ)*(dof_values[0]*J_01 + dof_values[1]*2*J_01 + dof_values[2]*(J_00 + J_01) + dof_values[3]*(2*J_00 - 2*J_01));
3093
vertex_values[1] = (1.0/detJ)*(dof_values[2]*2*J_10 + dof_values[3]*J_10 + dof_values[4]*(-2*J_11) + dof_values[5]*J_11);
3094
vertex_values[4] = (1.0/detJ)*(dof_values[0]*2*J_10 + dof_values[1]*J_10 + dof_values[4]*(J_10 + J_11) + dof_values[5]*(2*J_10 - 2*J_11));
3095
vertex_values[7] = (1.0/detJ)*(dof_values[0]*J_11 + dof_values[1]*2*J_11 + dof_values[2]*(J_10 + J_11) + dof_values[3]*(2*J_10 - 2*J_11));
3096
// Evaluate at vertices and use affine mapping
3097
vertex_values[2] = dof_values[6];
3098
vertex_values[5] = dof_values[6];
3099
vertex_values[8] = dof_values[6];
3102
/// Return the number of sub elements (for a mixed element)
3103
virtual unsigned int num_sub_elements() const
3108
/// Create a new finite element for sub element i (for a mixed element)
3109
virtual ufc::finite_element* create_sub_element(unsigned int i) const
3114
return new mixedpoisson_0_finite_element_1_0();
3117
return new mixedpoisson_0_finite_element_1_1();
3125
/// This class defines the interface for a local-to-global mapping of
3126
/// degrees of freedom (dofs).
3128
class mixedpoisson_0_dof_map_0_0: public ufc::dof_map
3132
unsigned int __global_dimension;
3137
mixedpoisson_0_dof_map_0_0() : ufc::dof_map()
3139
__global_dimension = 0;
3143
virtual ~mixedpoisson_0_dof_map_0_0()
3148
/// Return a string identifying the dof map
3149
virtual const char* signature() const
3151
return "FFC dof map for FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1)";
3154
/// Return true iff mesh entities of topological dimension d are needed
3155
virtual bool needs_mesh_entities(unsigned int d) const
3172
/// Initialize dof map for mesh (return true iff init_cell() is needed)
3173
virtual bool init_mesh(const ufc::mesh& m)
3175
__global_dimension = 2*m.num_entities[1];
3179
/// Initialize dof map for given cell
3180
virtual void init_cell(const ufc::mesh& m,
3186
/// Finish initialization of dof map for cells
3187
virtual void init_cell_finalize()
3192
/// Return the dimension of the global finite element function space
3193
virtual unsigned int global_dimension() const
3195
return __global_dimension;
3198
/// Return the dimension of the local finite element function space for a cell
3199
virtual unsigned int local_dimension(const ufc::cell& c) const
3204
/// Return the maximum dimension of the local finite element function space
3205
virtual unsigned int max_local_dimension() const
3210
// Return the geometric dimension of the coordinates this dof map provides
3211
virtual unsigned int geometric_dimension() const
3216
/// Return the number of dofs on each cell facet
3217
virtual unsigned int num_facet_dofs() const
3222
/// Return the number of dofs associated with each cell entity of dimension d
3223
virtual unsigned int num_entity_dofs(unsigned int d) const
3225
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3228
/// Tabulate the local-to-global mapping of dofs on a cell
3229
virtual void tabulate_dofs(unsigned int* dofs,
3231
const ufc::cell& c) const
3233
dofs[0] = 2*c.entity_indices[1][0];
3234
dofs[1] = 2*c.entity_indices[1][0] + 1;
3235
dofs[2] = 2*c.entity_indices[1][1];
3236
dofs[3] = 2*c.entity_indices[1][1] + 1;
3237
dofs[4] = 2*c.entity_indices[1][2];
3238
dofs[5] = 2*c.entity_indices[1][2] + 1;
3241
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3242
virtual void tabulate_facet_dofs(unsigned int* dofs,
3243
unsigned int facet) const
3262
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3263
virtual void tabulate_entity_dofs(unsigned int* dofs,
3264
unsigned int d, unsigned int i) const
3266
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3269
/// Tabulate the coordinates of all dofs on a cell
3270
virtual void tabulate_coordinates(double** coordinates,
3271
const ufc::cell& c) const
3273
const double * const * x = c.coordinates;
3274
coordinates[0][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];
3275
coordinates[0][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];
3276
coordinates[1][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];
3277
coordinates[1][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];
3278
coordinates[2][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];
3279
coordinates[2][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];
3280
coordinates[3][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];
3281
coordinates[3][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];
3282
coordinates[4][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];
3283
coordinates[4][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];
3284
coordinates[5][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];
3285
coordinates[5][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];
3288
/// Return the number of sub dof maps (for a mixed element)
3289
virtual unsigned int num_sub_dof_maps() const
3294
/// Create a new dof_map for sub dof map i (for a mixed element)
3295
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
3297
return new mixedpoisson_0_dof_map_0_0();
3302
/// This class defines the interface for a local-to-global mapping of
3303
/// degrees of freedom (dofs).
3305
class mixedpoisson_0_dof_map_0_1: public ufc::dof_map
3309
unsigned int __global_dimension;
3314
mixedpoisson_0_dof_map_0_1() : ufc::dof_map()
3316
__global_dimension = 0;
3320
virtual ~mixedpoisson_0_dof_map_0_1()
3325
/// Return a string identifying the dof map
3326
virtual const char* signature() const
3328
return "FFC dof map for FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
3331
/// Return true iff mesh entities of topological dimension d are needed
3332
virtual bool needs_mesh_entities(unsigned int d) const
3349
/// Initialize dof map for mesh (return true iff init_cell() is needed)
3350
virtual bool init_mesh(const ufc::mesh& m)
3352
__global_dimension = m.num_entities[2];
3356
/// Initialize dof map for given cell
3357
virtual void init_cell(const ufc::mesh& m,
3363
/// Finish initialization of dof map for cells
3364
virtual void init_cell_finalize()
3369
/// Return the dimension of the global finite element function space
3370
virtual unsigned int global_dimension() const
3372
return __global_dimension;
3375
/// Return the dimension of the local finite element function space for a cell
3376
virtual unsigned int local_dimension(const ufc::cell& c) const
3381
/// Return the maximum dimension of the local finite element function space
3382
virtual unsigned int max_local_dimension() const
3387
// Return the geometric dimension of the coordinates this dof map provides
3388
virtual unsigned int geometric_dimension() const
3393
/// Return the number of dofs on each cell facet
3394
virtual unsigned int num_facet_dofs() const
3399
/// Return the number of dofs associated with each cell entity of dimension d
3400
virtual unsigned int num_entity_dofs(unsigned int d) const
3402
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3405
/// Tabulate the local-to-global mapping of dofs on a cell
3406
virtual void tabulate_dofs(unsigned int* dofs,
3408
const ufc::cell& c) const
3410
dofs[0] = c.entity_indices[2][0];
3413
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3414
virtual void tabulate_facet_dofs(unsigned int* dofs,
3415
unsigned int facet) const
3431
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3432
virtual void tabulate_entity_dofs(unsigned int* dofs,
3433
unsigned int d, unsigned int i) const
3435
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3438
/// Tabulate the coordinates of all dofs on a cell
3439
virtual void tabulate_coordinates(double** coordinates,
3440
const ufc::cell& c) const
3442
const double * const * x = c.coordinates;
3443
coordinates[0][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
3444
coordinates[0][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
3447
/// Return the number of sub dof maps (for a mixed element)
3448
virtual unsigned int num_sub_dof_maps() const
3453
/// Create a new dof_map for sub dof map i (for a mixed element)
3454
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
3456
return new mixedpoisson_0_dof_map_0_1();
3461
/// This class defines the interface for a local-to-global mapping of
3462
/// degrees of freedom (dofs).
3464
class mixedpoisson_0_dof_map_0: public ufc::dof_map
3468
unsigned int __global_dimension;
3473
mixedpoisson_0_dof_map_0() : ufc::dof_map()
3475
__global_dimension = 0;
3479
virtual ~mixedpoisson_0_dof_map_0()
3484
/// Return a string identifying the dof map
3485
virtual const char* signature() const
3487
return "FFC dof map for MixedElement([FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1), FiniteElement('Discontinuous Lagrange', 'triangle', 0)])";
3490
/// Return true iff mesh entities of topological dimension d are needed
3491
virtual bool needs_mesh_entities(unsigned int d) const
3508
/// Initialize dof map for mesh (return true iff init_cell() is needed)
3509
virtual bool init_mesh(const ufc::mesh& m)
3511
__global_dimension = 2*m.num_entities[1] + m.num_entities[2];
3515
/// Initialize dof map for given cell
3516
virtual void init_cell(const ufc::mesh& m,
3522
/// Finish initialization of dof map for cells
3523
virtual void init_cell_finalize()
3528
/// Return the dimension of the global finite element function space
3529
virtual unsigned int global_dimension() const
3531
return __global_dimension;
3534
/// Return the dimension of the local finite element function space for a cell
3535
virtual unsigned int local_dimension(const ufc::cell& c) const
3540
/// Return the maximum dimension of the local finite element function space
3541
virtual unsigned int max_local_dimension() const
3546
// Return the geometric dimension of the coordinates this dof map provides
3547
virtual unsigned int geometric_dimension() const
3552
/// Return the number of dofs on each cell facet
3553
virtual unsigned int num_facet_dofs() const
3558
/// Return the number of dofs associated with each cell entity of dimension d
3559
virtual unsigned int num_entity_dofs(unsigned int d) const
3561
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3564
/// Tabulate the local-to-global mapping of dofs on a cell
3565
virtual void tabulate_dofs(unsigned int* dofs,
3567
const ufc::cell& c) const
3569
dofs[0] = 2*c.entity_indices[1][0];
3570
dofs[1] = 2*c.entity_indices[1][0] + 1;
3571
dofs[2] = 2*c.entity_indices[1][1];
3572
dofs[3] = 2*c.entity_indices[1][1] + 1;
3573
dofs[4] = 2*c.entity_indices[1][2];
3574
dofs[5] = 2*c.entity_indices[1][2] + 1;
3575
unsigned int offset = 2*m.num_entities[1];
3576
dofs[6] = offset + c.entity_indices[2][0];
3579
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3580
virtual void tabulate_facet_dofs(unsigned int* dofs,
3581
unsigned int facet) const
3600
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3601
virtual void tabulate_entity_dofs(unsigned int* dofs,
3602
unsigned int d, unsigned int i) const
3604
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3607
/// Tabulate the coordinates of all dofs on a cell
3608
virtual void tabulate_coordinates(double** coordinates,
3609
const ufc::cell& c) const
3611
const double * const * x = c.coordinates;
3612
coordinates[0][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];
3613
coordinates[0][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];
3614
coordinates[1][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];
3615
coordinates[1][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];
3616
coordinates[2][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];
3617
coordinates[2][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];
3618
coordinates[3][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];
3619
coordinates[3][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];
3620
coordinates[4][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];
3621
coordinates[4][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];
3622
coordinates[5][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];
3623
coordinates[5][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];
3624
coordinates[6][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
3625
coordinates[6][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
3628
/// Return the number of sub dof maps (for a mixed element)
3629
virtual unsigned int num_sub_dof_maps() const
3634
/// Create a new dof_map for sub dof map i (for a mixed element)
3635
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
3640
return new mixedpoisson_0_dof_map_0_0();
3643
return new mixedpoisson_0_dof_map_0_1();
3651
/// This class defines the interface for a local-to-global mapping of
3652
/// degrees of freedom (dofs).
3654
class mixedpoisson_0_dof_map_1_0: public ufc::dof_map
3658
unsigned int __global_dimension;
3663
mixedpoisson_0_dof_map_1_0() : ufc::dof_map()
3665
__global_dimension = 0;
3669
virtual ~mixedpoisson_0_dof_map_1_0()
3674
/// Return a string identifying the dof map
3675
virtual const char* signature() const
3677
return "FFC dof map for FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1)";
3680
/// Return true iff mesh entities of topological dimension d are needed
3681
virtual bool needs_mesh_entities(unsigned int d) const
3698
/// Initialize dof map for mesh (return true iff init_cell() is needed)
3699
virtual bool init_mesh(const ufc::mesh& m)
3701
__global_dimension = 2*m.num_entities[1];
3705
/// Initialize dof map for given cell
3706
virtual void init_cell(const ufc::mesh& m,
3712
/// Finish initialization of dof map for cells
3713
virtual void init_cell_finalize()
3718
/// Return the dimension of the global finite element function space
3719
virtual unsigned int global_dimension() const
3721
return __global_dimension;
3724
/// Return the dimension of the local finite element function space for a cell
3725
virtual unsigned int local_dimension(const ufc::cell& c) const
3730
/// Return the maximum dimension of the local finite element function space
3731
virtual unsigned int max_local_dimension() const
3736
// Return the geometric dimension of the coordinates this dof map provides
3737
virtual unsigned int geometric_dimension() const
3742
/// Return the number of dofs on each cell facet
3743
virtual unsigned int num_facet_dofs() const
3748
/// Return the number of dofs associated with each cell entity of dimension d
3749
virtual unsigned int num_entity_dofs(unsigned int d) const
3751
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3754
/// Tabulate the local-to-global mapping of dofs on a cell
3755
virtual void tabulate_dofs(unsigned int* dofs,
3757
const ufc::cell& c) const
3759
dofs[0] = 2*c.entity_indices[1][0];
3760
dofs[1] = 2*c.entity_indices[1][0] + 1;
3761
dofs[2] = 2*c.entity_indices[1][1];
3762
dofs[3] = 2*c.entity_indices[1][1] + 1;
3763
dofs[4] = 2*c.entity_indices[1][2];
3764
dofs[5] = 2*c.entity_indices[1][2] + 1;
3767
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3768
virtual void tabulate_facet_dofs(unsigned int* dofs,
3769
unsigned int facet) const
3788
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3789
virtual void tabulate_entity_dofs(unsigned int* dofs,
3790
unsigned int d, unsigned int i) const
3792
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3795
/// Tabulate the coordinates of all dofs on a cell
3796
virtual void tabulate_coordinates(double** coordinates,
3797
const ufc::cell& c) const
3799
const double * const * x = c.coordinates;
3800
coordinates[0][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];
3801
coordinates[0][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];
3802
coordinates[1][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];
3803
coordinates[1][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];
3804
coordinates[2][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];
3805
coordinates[2][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];
3806
coordinates[3][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];
3807
coordinates[3][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];
3808
coordinates[4][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];
3809
coordinates[4][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];
3810
coordinates[5][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];
3811
coordinates[5][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];
3814
/// Return the number of sub dof maps (for a mixed element)
3815
virtual unsigned int num_sub_dof_maps() const
3820
/// Create a new dof_map for sub dof map i (for a mixed element)
3821
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
3823
return new mixedpoisson_0_dof_map_1_0();
3828
/// This class defines the interface for a local-to-global mapping of
3829
/// degrees of freedom (dofs).
3831
class mixedpoisson_0_dof_map_1_1: public ufc::dof_map
3835
unsigned int __global_dimension;
3840
mixedpoisson_0_dof_map_1_1() : ufc::dof_map()
3842
__global_dimension = 0;
3846
virtual ~mixedpoisson_0_dof_map_1_1()
3851
/// Return a string identifying the dof map
3852
virtual const char* signature() const
3854
return "FFC dof map for FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
3857
/// Return true iff mesh entities of topological dimension d are needed
3858
virtual bool needs_mesh_entities(unsigned int d) const
3875
/// Initialize dof map for mesh (return true iff init_cell() is needed)
3876
virtual bool init_mesh(const ufc::mesh& m)
3878
__global_dimension = m.num_entities[2];
3882
/// Initialize dof map for given cell
3883
virtual void init_cell(const ufc::mesh& m,
3889
/// Finish initialization of dof map for cells
3890
virtual void init_cell_finalize()
3895
/// Return the dimension of the global finite element function space
3896
virtual unsigned int global_dimension() const
3898
return __global_dimension;
3901
/// Return the dimension of the local finite element function space for a cell
3902
virtual unsigned int local_dimension(const ufc::cell& c) const
3907
/// Return the maximum dimension of the local finite element function space
3908
virtual unsigned int max_local_dimension() const
3913
// Return the geometric dimension of the coordinates this dof map provides
3914
virtual unsigned int geometric_dimension() const
3919
/// Return the number of dofs on each cell facet
3920
virtual unsigned int num_facet_dofs() const
3925
/// Return the number of dofs associated with each cell entity of dimension d
3926
virtual unsigned int num_entity_dofs(unsigned int d) const
3928
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3931
/// Tabulate the local-to-global mapping of dofs on a cell
3932
virtual void tabulate_dofs(unsigned int* dofs,
3934
const ufc::cell& c) const
3936
dofs[0] = c.entity_indices[2][0];
3939
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3940
virtual void tabulate_facet_dofs(unsigned int* dofs,
3941
unsigned int facet) const
3957
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3958
virtual void tabulate_entity_dofs(unsigned int* dofs,
3959
unsigned int d, unsigned int i) const
3961
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3964
/// Tabulate the coordinates of all dofs on a cell
3965
virtual void tabulate_coordinates(double** coordinates,
3966
const ufc::cell& c) const
3968
const double * const * x = c.coordinates;
3969
coordinates[0][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
3970
coordinates[0][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
3973
/// Return the number of sub dof maps (for a mixed element)
3974
virtual unsigned int num_sub_dof_maps() const
3979
/// Create a new dof_map for sub dof map i (for a mixed element)
3980
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
3982
return new mixedpoisson_0_dof_map_1_1();
3987
/// This class defines the interface for a local-to-global mapping of
3988
/// degrees of freedom (dofs).
3990
class mixedpoisson_0_dof_map_1: public ufc::dof_map
3994
unsigned int __global_dimension;
3999
mixedpoisson_0_dof_map_1() : ufc::dof_map()
4001
__global_dimension = 0;
4005
virtual ~mixedpoisson_0_dof_map_1()
4010
/// Return a string identifying the dof map
4011
virtual const char* signature() const
4013
return "FFC dof map for MixedElement([FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1), FiniteElement('Discontinuous Lagrange', 'triangle', 0)])";
4016
/// Return true iff mesh entities of topological dimension d are needed
4017
virtual bool needs_mesh_entities(unsigned int d) const
4034
/// Initialize dof map for mesh (return true iff init_cell() is needed)
4035
virtual bool init_mesh(const ufc::mesh& m)
4037
__global_dimension = 2*m.num_entities[1] + m.num_entities[2];
4041
/// Initialize dof map for given cell
4042
virtual void init_cell(const ufc::mesh& m,
4048
/// Finish initialization of dof map for cells
4049
virtual void init_cell_finalize()
4054
/// Return the dimension of the global finite element function space
4055
virtual unsigned int global_dimension() const
4057
return __global_dimension;
4060
/// Return the dimension of the local finite element function space for a cell
4061
virtual unsigned int local_dimension(const ufc::cell& c) const
4066
/// Return the maximum dimension of the local finite element function space
4067
virtual unsigned int max_local_dimension() const
4072
// Return the geometric dimension of the coordinates this dof map provides
4073
virtual unsigned int geometric_dimension() const
4078
/// Return the number of dofs on each cell facet
4079
virtual unsigned int num_facet_dofs() const
4084
/// Return the number of dofs associated with each cell entity of dimension d
4085
virtual unsigned int num_entity_dofs(unsigned int d) const
4087
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4090
/// Tabulate the local-to-global mapping of dofs on a cell
4091
virtual void tabulate_dofs(unsigned int* dofs,
4093
const ufc::cell& c) const
4095
dofs[0] = 2*c.entity_indices[1][0];
4096
dofs[1] = 2*c.entity_indices[1][0] + 1;
4097
dofs[2] = 2*c.entity_indices[1][1];
4098
dofs[3] = 2*c.entity_indices[1][1] + 1;
4099
dofs[4] = 2*c.entity_indices[1][2];
4100
dofs[5] = 2*c.entity_indices[1][2] + 1;
4101
unsigned int offset = 2*m.num_entities[1];
4102
dofs[6] = offset + c.entity_indices[2][0];
4105
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4106
virtual void tabulate_facet_dofs(unsigned int* dofs,
4107
unsigned int facet) const
4126
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4127
virtual void tabulate_entity_dofs(unsigned int* dofs,
4128
unsigned int d, unsigned int i) const
4130
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4133
/// Tabulate the coordinates of all dofs on a cell
4134
virtual void tabulate_coordinates(double** coordinates,
4135
const ufc::cell& c) const
4137
const double * const * x = c.coordinates;
4138
coordinates[0][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];
4139
coordinates[0][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];
4140
coordinates[1][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];
4141
coordinates[1][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];
4142
coordinates[2][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];
4143
coordinates[2][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];
4144
coordinates[3][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];
4145
coordinates[3][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];
4146
coordinates[4][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];
4147
coordinates[4][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];
4148
coordinates[5][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];
4149
coordinates[5][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];
4150
coordinates[6][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
4151
coordinates[6][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
4154
/// Return the number of sub dof maps (for a mixed element)
4155
virtual unsigned int num_sub_dof_maps() const
4160
/// Create a new dof_map for sub dof map i (for a mixed element)
4161
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
4166
return new mixedpoisson_0_dof_map_1_0();
4169
return new mixedpoisson_0_dof_map_1_1();
4177
/// This class defines the interface for the tabulation of the cell
4178
/// tensor corresponding to the local contribution to a form from
4179
/// the integral over a cell.
4181
class mixedpoisson_0_cell_integral_0_tensor: public ufc::cell_integral
4186
mixedpoisson_0_cell_integral_0_tensor() : ufc::cell_integral()
4192
virtual ~mixedpoisson_0_cell_integral_0_tensor()
4197
/// Tabulate the tensor for the contribution from a local cell
4198
virtual void tabulate_tensor(double* A,
4199
const double * const * w,
4200
const ufc::cell& c) const
4202
// Number of operations to compute geometry tensor: 48
4203
// Number of operations to compute tensor contraction: 592
4204
// Total number of operations to compute cell tensor: 640
4206
// Extract vertex coordinates
4207
const double * const * x = c.coordinates;
4209
// Compute Jacobian of affine map from reference cell
4210
const double J_00 = x[1][0] - x[0][0];
4211
const double J_01 = x[2][0] - x[0][0];
4212
const double J_10 = x[1][1] - x[0][1];
4213
const double J_11 = x[2][1] - x[0][1];
4215
// Compute determinant of Jacobian
4216
double detJ = J_00*J_11 - J_01*J_10;
4218
// Compute inverse of Jacobian
4219
const double Jinv_00 = J_11 / detJ;
4220
const double Jinv_01 = -J_01 / detJ;
4221
const double Jinv_10 = -J_10 / detJ;
4222
const double Jinv_11 = J_00 / detJ;
4225
const double det = std::abs(detJ);
4227
// Compute geometry tensor
4228
const double G0_0_0 = 1.0/(detJ)*det*J_00*Jinv_00;
4229
const double G0_0_1 = 1.0/(detJ)*det*J_00*Jinv_10;
4230
const double G0_1_0 = 1.0/(detJ)*det*J_01*Jinv_00;
4231
const double G0_1_1 = 1.0/(detJ)*det*J_01*Jinv_10;
4232
const double G1_0_0 = 1.0/(detJ)*det*J_10*Jinv_01;
4233
const double G1_0_1 = 1.0/(detJ)*det*J_10*Jinv_11;
4234
const double G1_1_0 = 1.0/(detJ)*det*J_11*Jinv_01;
4235
const double G1_1_1 = 1.0/(detJ)*det*J_11*Jinv_11;
4236
const double G2_0_0 = 1.0/(detJ)*det*J_00*Jinv_00;
4237
const double G2_0_1 = 1.0/(detJ)*det*J_00*Jinv_10;
4238
const double G2_1_0 = 1.0/(detJ)*det*J_01*Jinv_00;
4239
const double G2_1_1 = 1.0/(detJ)*det*J_01*Jinv_10;
4240
const double G3_0_0 = 1.0/(detJ)*det*J_10*Jinv_01;
4241
const double G3_0_1 = 1.0/(detJ)*det*J_10*Jinv_11;
4242
const double G3_1_0 = 1.0/(detJ)*det*J_11*Jinv_01;
4243
const double G3_1_1 = 1.0/(detJ)*det*J_11*Jinv_11;
4244
const double G4_0_0 = 1.0/(detJ*detJ)*det*J_00*J_00;
4245
const double G4_0_1 = 1.0/(detJ*detJ)*det*J_00*J_01;
4246
const double G4_1_0 = 1.0/(detJ*detJ)*det*J_01*J_00;
4247
const double G4_1_1 = 1.0/(detJ*detJ)*det*J_01*J_01;
4248
const double G5_0_0 = 1.0/(detJ*detJ)*det*J_10*J_10;
4249
const double G5_0_1 = 1.0/(detJ*detJ)*det*J_10*J_11;
4250
const double G5_1_0 = 1.0/(detJ*detJ)*det*J_11*J_10;
4251
const double G5_1_1 = 1.0/(detJ*detJ)*det*J_11*J_11;
4253
// Compute element tensor
4254
A[0] += 0.333333333*G4_0_0 - 0.0833333333*G4_0_1 - 0.0833333333*G4_1_0 + 0.0833333333*G4_1_1 + 0.333333333*G5_0_0 - 0.0833333333*G5_0_1 - 0.0833333333*G5_1_0 + 0.0833333333*G5_1_1;
4255
A[1] += -0.166666667*G4_0_0 + 0.166666667*G4_0_1 + 0.0416666667*G4_1_0 - 0.166666667*G4_1_1 - 0.166666667*G5_0_0 + 0.166666667*G5_0_1 + 0.0416666667*G5_1_0 - 0.166666667*G5_1_1;
4256
A[2] += 0.0833333333*G4_0_0 + 0.0833333333*G4_0_1 - 0.0833333333*G4_1_1 + 0.0833333333*G5_0_0 + 0.0833333333*G5_0_1 - 0.0833333333*G5_1_1;
4257
A[3] += 0.0833333333*G4_0_0 - 0.166666667*G4_0_1 - 0.125*G4_1_0 + 0.166666667*G4_1_1 + 0.0833333333*G5_0_0 - 0.166666667*G5_0_1 - 0.125*G5_1_0 + 0.166666667*G5_1_1;
4258
A[4] += -0.166666667*G4_0_0 + 0.0416666667*G4_1_0 + 0.0416666667*G4_1_1 - 0.166666667*G5_0_0 + 0.0416666667*G5_1_0 + 0.0416666667*G5_1_1;
4259
A[5] += 0.333333333*G4_0_0 - 0.25*G4_0_1 - 0.0833333333*G4_1_0 + 0.0416666667*G4_1_1 + 0.333333333*G5_0_0 - 0.25*G5_0_1 - 0.0833333333*G5_1_0 + 0.0416666667*G5_1_1;
4260
A[6] += -1*G2_0_0 + 0.5*G2_1_1 - 1*G3_0_0 + 0.5*G3_1_1;
4261
A[7] += -0.166666667*G4_0_0 + 0.0416666667*G4_0_1 + 0.166666667*G4_1_0 - 0.166666667*G4_1_1 - 0.166666667*G5_0_0 + 0.0416666667*G5_0_1 + 0.166666667*G5_1_0 - 0.166666667*G5_1_1;
4262
A[8] += 0.0833333333*G4_0_0 - 0.0833333333*G4_0_1 - 0.0833333333*G4_1_0 + 0.333333333*G4_1_1 + 0.0833333333*G5_0_0 - 0.0833333333*G5_0_1 - 0.0833333333*G5_1_0 + 0.333333333*G5_1_1;
4263
A[9] += -0.0416666667*G4_0_0 - 0.0416666667*G4_0_1 + 0.166666667*G4_1_1 - 0.0416666667*G5_0_0 - 0.0416666667*G5_0_1 + 0.166666667*G5_1_1;
4264
A[10] += -0.0416666667*G4_0_0 + 0.0833333333*G4_0_1 + 0.25*G4_1_0 - 0.333333333*G4_1_1 - 0.0416666667*G5_0_0 + 0.0833333333*G5_0_1 + 0.25*G5_1_0 - 0.333333333*G5_1_1;
4265
A[11] += 0.0833333333*G4_0_0 - 0.0833333333*G4_1_0 - 0.0833333333*G4_1_1 + 0.0833333333*G5_0_0 - 0.0833333333*G5_1_0 - 0.0833333333*G5_1_1;
4266
A[12] += -0.166666667*G4_0_0 + 0.125*G4_0_1 + 0.166666667*G4_1_0 - 0.0833333333*G4_1_1 - 0.166666667*G5_0_0 + 0.125*G5_0_1 + 0.166666667*G5_1_0 - 0.0833333333*G5_1_1;
4267
A[13] += 0.5*G2_0_0 - 1*G2_1_1 + 0.5*G3_0_0 - 1*G3_1_1;
4268
A[14] += 0.0833333333*G4_0_0 + 0.0833333333*G4_1_0 - 0.0833333333*G4_1_1 + 0.0833333333*G5_0_0 + 0.0833333333*G5_1_0 - 0.0833333333*G5_1_1;
4269
A[15] += -0.0416666667*G4_0_0 - 0.0416666667*G4_1_0 + 0.166666667*G4_1_1 - 0.0416666667*G5_0_0 - 0.0416666667*G5_1_0 + 0.166666667*G5_1_1;
4270
A[16] += 0.25*G4_0_0 + 0.0833333333*G4_1_1 + 0.25*G5_0_0 + 0.0833333333*G5_1_1;
4271
A[17] += -0.125*G4_0_0 + 0.125*G4_1_0 - 0.166666667*G4_1_1 - 0.125*G5_0_0 + 0.125*G5_1_0 - 0.166666667*G5_1_1;
4272
A[18] += -0.0416666667*G4_0_0 - 0.208333333*G4_0_1 - 0.0416666667*G4_1_0 - 0.0416666667*G4_1_1 - 0.0416666667*G5_0_0 - 0.208333333*G5_0_1 - 0.0416666667*G5_1_0 - 0.0416666667*G5_1_1;
4273
A[19] += 0.0833333333*G4_0_0 + 0.0416666667*G4_0_1 + 0.0833333333*G4_1_0 - 0.0416666667*G4_1_1 + 0.0833333333*G5_0_0 + 0.0416666667*G5_0_1 + 0.0833333333*G5_1_0 - 0.0416666667*G5_1_1;
4274
A[20] += 1*G2_0_0 + 1.5*G2_0_1 - 0.5*G2_1_1 + 1*G3_0_0 + 1.5*G3_0_1 - 0.5*G3_1_1;
4275
A[21] += 0.0833333333*G4_0_0 - 0.125*G4_0_1 - 0.166666667*G4_1_0 + 0.166666667*G4_1_1 + 0.0833333333*G5_0_0 - 0.125*G5_0_1 - 0.166666667*G5_1_0 + 0.166666667*G5_1_1;
4276
A[22] += -0.0416666667*G4_0_0 + 0.25*G4_0_1 + 0.0833333333*G4_1_0 - 0.333333333*G4_1_1 - 0.0416666667*G5_0_0 + 0.25*G5_0_1 + 0.0833333333*G5_1_0 - 0.333333333*G5_1_1;
4277
A[23] += -0.125*G4_0_0 + 0.125*G4_0_1 - 0.166666667*G4_1_1 - 0.125*G5_0_0 + 0.125*G5_0_1 - 0.166666667*G5_1_1;
4278
A[24] += 0.25*G4_0_0 - 0.25*G4_0_1 - 0.25*G4_1_0 + 0.333333333*G4_1_1 + 0.25*G5_0_0 - 0.25*G5_0_1 - 0.25*G5_1_0 + 0.333333333*G5_1_1;
4279
A[25] += -0.0416666667*G4_0_0 + 0.0416666667*G4_0_1 + 0.0833333333*G4_1_0 + 0.0833333333*G4_1_1 - 0.0416666667*G5_0_0 + 0.0416666667*G5_0_1 + 0.0833333333*G5_1_0 + 0.0833333333*G5_1_1;
4280
A[26] += 0.0833333333*G4_0_0 - 0.0833333333*G4_0_1 - 0.166666667*G4_1_0 + 0.0833333333*G4_1_1 + 0.0833333333*G5_0_0 - 0.0833333333*G5_0_1 - 0.166666667*G5_1_0 + 0.0833333333*G5_1_1;
4281
A[27] += -0.5*G2_0_0 - 1.5*G2_0_1 + 1*G2_1_1 - 0.5*G3_0_0 - 1.5*G3_0_1 + 1*G3_1_1;
4282
A[28] += -0.166666667*G4_0_0 + 0.0416666667*G4_0_1 + 0.0416666667*G4_1_1 - 0.166666667*G5_0_0 + 0.0416666667*G5_0_1 + 0.0416666667*G5_1_1;
4283
A[29] += 0.0833333333*G4_0_0 - 0.0833333333*G4_0_1 - 0.0833333333*G4_1_1 + 0.0833333333*G5_0_0 - 0.0833333333*G5_0_1 - 0.0833333333*G5_1_1;
4284
A[30] += -0.0416666667*G4_0_0 - 0.0416666667*G4_0_1 - 0.208333333*G4_1_0 - 0.0416666667*G4_1_1 - 0.0416666667*G5_0_0 - 0.0416666667*G5_0_1 - 0.208333333*G5_1_0 - 0.0416666667*G5_1_1;
4285
A[31] += -0.0416666667*G4_0_0 + 0.0833333333*G4_0_1 + 0.0416666667*G4_1_0 + 0.0833333333*G4_1_1 - 0.0416666667*G5_0_0 + 0.0833333333*G5_0_1 + 0.0416666667*G5_1_0 + 0.0833333333*G5_1_1;
4286
A[32] += 0.0833333333*G4_0_0 + 0.25*G4_1_1 + 0.0833333333*G5_0_0 + 0.25*G5_1_1;
4287
A[33] += -0.166666667*G4_0_0 + 0.125*G4_0_1 - 0.125*G4_1_1 - 0.166666667*G5_0_0 + 0.125*G5_0_1 - 0.125*G5_1_1;
4288
A[34] += 0.5*G2_0_0 - 1.5*G2_1_0 - 1*G2_1_1 + 0.5*G3_0_0 - 1.5*G3_1_0 - 1*G3_1_1;
4289
A[35] += 0.333333333*G4_0_0 - 0.0833333333*G4_0_1 - 0.25*G4_1_0 + 0.0416666667*G4_1_1 + 0.333333333*G5_0_0 - 0.0833333333*G5_0_1 - 0.25*G5_1_0 + 0.0416666667*G5_1_1;
4290
A[36] += -0.166666667*G4_0_0 + 0.166666667*G4_0_1 + 0.125*G4_1_0 - 0.0833333333*G4_1_1 - 0.166666667*G5_0_0 + 0.166666667*G5_0_1 + 0.125*G5_1_0 - 0.0833333333*G5_1_1;
4291
A[37] += 0.0833333333*G4_0_0 + 0.0833333333*G4_0_1 + 0.0416666667*G4_1_0 - 0.0416666667*G4_1_1 + 0.0833333333*G5_0_0 + 0.0833333333*G5_0_1 + 0.0416666667*G5_1_0 - 0.0416666667*G5_1_1;
4292
A[38] += 0.0833333333*G4_0_0 - 0.166666667*G4_0_1 - 0.0833333333*G4_1_0 + 0.0833333333*G4_1_1 + 0.0833333333*G5_0_0 - 0.166666667*G5_0_1 - 0.0833333333*G5_1_0 + 0.0833333333*G5_1_1;
4293
A[39] += -0.166666667*G4_0_0 + 0.125*G4_1_0 - 0.125*G4_1_1 - 0.166666667*G5_0_0 + 0.125*G5_1_0 - 0.125*G5_1_1;
4294
A[40] += 0.333333333*G4_0_0 - 0.25*G4_0_1 - 0.25*G4_1_0 + 0.25*G4_1_1 + 0.333333333*G5_0_0 - 0.25*G5_0_1 - 0.25*G5_1_0 + 0.25*G5_1_1;
4295
A[41] += -1*G2_0_0 + 1.5*G2_1_0 + 0.5*G2_1_1 - 1*G3_0_0 + 1.5*G3_1_0 + 0.5*G3_1_1;
4296
A[42] += 1*G0_0_0 - 0.5*G0_1_1 + 1*G1_0_0 - 0.5*G1_1_1;
4297
A[43] += -0.5*G0_0_0 + 1*G0_1_1 - 0.5*G1_0_0 + 1*G1_1_1;
4298
A[44] += -1*G0_0_0 - 1.5*G0_0_1 + 0.5*G0_1_1 - 1*G1_0_0 - 1.5*G1_0_1 + 0.5*G1_1_1;
4299
A[45] += 0.5*G0_0_0 + 1.5*G0_0_1 - 1*G0_1_1 + 0.5*G1_0_0 + 1.5*G1_0_1 - 1*G1_1_1;
4300
A[46] += -0.5*G0_0_0 + 1.5*G0_1_0 + 1*G0_1_1 - 0.5*G1_0_0 + 1.5*G1_1_0 + 1*G1_1_1;
4301
A[47] += 1*G0_0_0 - 1.5*G0_1_0 - 0.5*G0_1_1 + 1*G1_0_0 - 1.5*G1_1_0 - 0.5*G1_1_1;
4307
/// This class defines the interface for the tabulation of the cell
4308
/// tensor corresponding to the local contribution to a form from
4309
/// the integral over a cell.
4311
class mixedpoisson_0_cell_integral_0: public ufc::cell_integral
4315
mixedpoisson_0_cell_integral_0_tensor integral_0_tensor;
4320
mixedpoisson_0_cell_integral_0() : ufc::cell_integral()
4326
virtual ~mixedpoisson_0_cell_integral_0()
4331
/// Tabulate the tensor for the contribution from a local cell
4332
virtual void tabulate_tensor(double* A,
4333
const double * const * w,
4334
const ufc::cell& c) const
4336
// Reset values of the element tensor block
4337
for (unsigned int j = 0; j < 49; j++)
4340
// Add all contributions to element tensor
4341
integral_0_tensor.tabulate_tensor(A, w, c);
4346
/// This class defines the interface for the assembly of the global
4347
/// tensor corresponding to a form with r + n arguments, that is, a
4350
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
4352
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
4353
/// global tensor A is defined by
4355
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
4357
/// where each argument Vj represents the application to the
4358
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
4359
/// fixed functions (coefficients).
4361
class mixedpoisson_form_0: public ufc::form
4366
mixedpoisson_form_0() : ufc::form()
4372
virtual ~mixedpoisson_form_0()
4377
/// Return a string identifying the form
4378
virtual const char* signature() const
4380
return "Form([Integral(Sum(Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 1), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 1), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), Indexed(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(2),), {FixedIndex(2): 3}))), Sum(IndexSum(Product(Indexed(ListTensor(Indexed(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(0),), {FixedIndex(0): 3})), Indexed(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(1),), {FixedIndex(1): 3}))), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ListTensor(Indexed(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 1), MultiIndex((FixedIndex(0),), {FixedIndex(0): 3})), Indexed(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 1), MultiIndex((FixedIndex(1),), {FixedIndex(1): 3}))), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Product(IntValue(-1, (), (), {}), Product(IndexSum(Indexed(ListTensor(Indexed(SpatialDerivative(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 0), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((FixedIndex(0),), {})), Indexed(SpatialDerivative(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 0), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((FixedIndex(1),), {}))), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), Indexed(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 1), MultiIndex((FixedIndex(2),), {FixedIndex(2): 3})))))), Measure('cell', 0, None))])";
4383
/// Return the rank of the global tensor (r)
4384
virtual unsigned int rank() const
4389
/// Return the number of coefficients (n)
4390
virtual unsigned int num_coefficients() const
4395
/// Return the number of cell integrals
4396
virtual unsigned int num_cell_integrals() const
4401
/// Return the number of exterior facet integrals
4402
virtual unsigned int num_exterior_facet_integrals() const
4407
/// Return the number of interior facet integrals
4408
virtual unsigned int num_interior_facet_integrals() const
4413
/// Create a new finite element for argument function i
4414
virtual ufc::finite_element* create_finite_element(unsigned int i) const
4419
return new mixedpoisson_0_finite_element_0();
4422
return new mixedpoisson_0_finite_element_1();
4428
/// Create a new dof map for argument function i
4429
virtual ufc::dof_map* create_dof_map(unsigned int i) const
4434
return new mixedpoisson_0_dof_map_0();
4437
return new mixedpoisson_0_dof_map_1();
4443
/// Create a new cell integral on sub domain i
4444
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
4446
return new mixedpoisson_0_cell_integral_0();
4449
/// Create a new exterior facet integral on sub domain i
4450
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
4455
/// Create a new interior facet integral on sub domain i
4456
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
4463
/// This class defines the interface for a finite element.
4465
class mixedpoisson_1_finite_element_0_0: public ufc::finite_element
4470
mixedpoisson_1_finite_element_0_0() : ufc::finite_element()
4476
virtual ~mixedpoisson_1_finite_element_0_0()
4481
/// Return a string identifying the finite element
4482
virtual const char* signature() const
4484
return "FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1)";
4487
/// Return the cell shape
4488
virtual ufc::shape cell_shape() const
4490
return ufc::triangle;
4493
/// Return the dimension of the finite element function space
4494
virtual unsigned int space_dimension() const
4499
/// Return the rank of the value space
4500
virtual unsigned int value_rank() const
4505
/// Return the dimension of the value space for axis i
4506
virtual unsigned int value_dimension(unsigned int i) const
4511
/// Evaluate basis function i at given point in cell
4512
virtual void evaluate_basis(unsigned int i,
4514
const double* coordinates,
4515
const ufc::cell& c) const
4517
// Extract vertex coordinates
4518
const double * const * element_coordinates = c.coordinates;
4520
// Compute Jacobian of affine map from reference cell
4521
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
4522
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
4523
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
4524
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
4526
// Compute determinant of Jacobian
4527
const double detJ = J_00*J_11 - J_01*J_10;
4529
// Compute inverse of Jacobian
4531
// Get coordinates and map to the reference (UFC) element
4532
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
4533
element_coordinates[0][0]*element_coordinates[2][1] +\
4534
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
4535
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
4536
element_coordinates[1][0]*element_coordinates[0][1] -\
4537
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
4539
// Map coordinates to the reference square
4540
if (std::abs(y - 1.0) < 1e-08)
4543
x = 2.0 *x/(1.0 - y) - 1.0;
4550
// Map degree of freedom to element degree of freedom
4551
const unsigned int dof = i;
4553
// Generate scalings
4554
const double scalings_y_0 = 1;
4555
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
4557
// Compute psitilde_a
4558
const double psitilde_a_0 = 1;
4559
const double psitilde_a_1 = x;
4561
// Compute psitilde_bs
4562
const double psitilde_bs_0_0 = 1;
4563
const double psitilde_bs_0_1 = 1.5*y + 0.5;
4564
const double psitilde_bs_1_0 = 1;
4566
// Compute basisvalues
4567
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
4568
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
4569
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
4571
// Table(s) of coefficients
4572
static const double coefficients0[6][3] = \
4573
{{0.942809042, 0.577350269, -0.333333333},
4574
{-0.471404521, -0.288675135, 0.166666667},
4575
{0.471404521, -0.577350269, -0.666666667},
4576
{0.471404521, 0.288675135, 0.833333333},
4577
{-0.471404521, -0.288675135, 0.166666667},
4578
{0.942809042, 0.577350269, -0.333333333}};
4580
static const double coefficients1[6][3] = \
4581
{{-0.471404521, 0, -0.333333333},
4582
{0.942809042, 0, 0.666666667},
4583
{0.471404521, 0, 0.333333333},
4584
{-0.942809042, 0, -0.666666667},
4585
{-0.471404521, 0.866025404, 0.166666667},
4586
{-0.471404521, -0.866025404, 0.166666667}};
4588
// Extract relevant coefficients
4589
const double coeff0_0 = coefficients0[dof][0];
4590
const double coeff0_1 = coefficients0[dof][1];
4591
const double coeff0_2 = coefficients0[dof][2];
4592
const double coeff1_0 = coefficients1[dof][0];
4593
const double coeff1_1 = coefficients1[dof][1];
4594
const double coeff1_2 = coefficients1[dof][2];
4597
const double tmp0_0 = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
4598
const double tmp0_1 = coeff1_0*basisvalue0 + coeff1_1*basisvalue1 + coeff1_2*basisvalue2;
4599
// Using contravariant Piola transform to map values back to the physical element
4600
values[0] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
4601
values[1] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
4604
/// Evaluate all basis functions at given point in cell
4605
virtual void evaluate_basis_all(double* values,
4606
const double* coordinates,
4607
const ufc::cell& c) const
4609
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
4612
/// Evaluate order n derivatives of basis function i at given point in cell
4613
virtual void evaluate_basis_derivatives(unsigned int i,
4616
const double* coordinates,
4617
const ufc::cell& c) const
4619
// Extract vertex coordinates
4620
const double * const * element_coordinates = c.coordinates;
4622
// Compute Jacobian of affine map from reference cell
4623
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
4624
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
4625
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
4626
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
4628
// Compute determinant of Jacobian
4629
const double detJ = J_00*J_11 - J_01*J_10;
4631
// Compute inverse of Jacobian
4633
// Get coordinates and map to the reference (UFC) element
4634
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
4635
element_coordinates[0][0]*element_coordinates[2][1] +\
4636
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
4637
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
4638
element_coordinates[1][0]*element_coordinates[0][1] -\
4639
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
4641
// Map coordinates to the reference square
4642
if (std::abs(y - 1.0) < 1e-08)
4645
x = 2.0 *x/(1.0 - y) - 1.0;
4648
// Compute number of derivatives
4649
unsigned int num_derivatives = 1;
4651
for (unsigned int j = 0; j < n; j++)
4652
num_derivatives *= 2;
4655
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
4656
unsigned int **combinations = new unsigned int *[num_derivatives];
4658
for (unsigned int j = 0; j < num_derivatives; j++)
4660
combinations[j] = new unsigned int [n];
4661
for (unsigned int k = 0; k < n; k++)
4662
combinations[j][k] = 0;
4665
// Generate combinations of derivatives
4666
for (unsigned int row = 1; row < num_derivatives; row++)
4668
for (unsigned int num = 0; num < row; num++)
4670
for (unsigned int col = n-1; col+1 > 0; col--)
4672
if (combinations[row][col] + 1 > 1)
4673
combinations[row][col] = 0;
4676
combinations[row][col] += 1;
4683
// Compute inverse of Jacobian
4684
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
4686
// Declare transformation matrix
4687
// Declare pointer to two dimensional array and initialise
4688
double **transform = new double *[num_derivatives];
4690
for (unsigned int j = 0; j < num_derivatives; j++)
4692
transform[j] = new double [num_derivatives];
4693
for (unsigned int k = 0; k < num_derivatives; k++)
4694
transform[j][k] = 1;
4697
// Construct transformation matrix
4698
for (unsigned int row = 0; row < num_derivatives; row++)
4700
for (unsigned int col = 0; col < num_derivatives; col++)
4702
for (unsigned int k = 0; k < n; k++)
4703
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
4708
for (unsigned int j = 0; j < 2*num_derivatives; j++)
4711
// Map degree of freedom to element degree of freedom
4712
const unsigned int dof = i;
4714
// Generate scalings
4715
const double scalings_y_0 = 1;
4716
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
4718
// Compute psitilde_a
4719
const double psitilde_a_0 = 1;
4720
const double psitilde_a_1 = x;
4722
// Compute psitilde_bs
4723
const double psitilde_bs_0_0 = 1;
4724
const double psitilde_bs_0_1 = 1.5*y + 0.5;
4725
const double psitilde_bs_1_0 = 1;
4727
// Compute basisvalues
4728
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
4729
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
4730
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
4732
// Table(s) of coefficients
4733
static const double coefficients0[6][3] = \
4734
{{0.942809042, 0.577350269, -0.333333333},
4735
{-0.471404521, -0.288675135, 0.166666667},
4736
{0.471404521, -0.577350269, -0.666666667},
4737
{0.471404521, 0.288675135, 0.833333333},
4738
{-0.471404521, -0.288675135, 0.166666667},
4739
{0.942809042, 0.577350269, -0.333333333}};
4741
static const double coefficients1[6][3] = \
4742
{{-0.471404521, 0, -0.333333333},
4743
{0.942809042, 0, 0.666666667},
4744
{0.471404521, 0, 0.333333333},
4745
{-0.942809042, 0, -0.666666667},
4746
{-0.471404521, 0.866025404, 0.166666667},
4747
{-0.471404521, -0.866025404, 0.166666667}};
4749
// Interesting (new) part
4750
// Tables of derivatives of the polynomial base (transpose)
4751
static const double dmats0[3][3] = \
4756
static const double dmats1[3][3] = \
4759
{4.24264069, 0, 0}};
4761
// Compute reference derivatives
4762
// Declare pointer to array of derivatives on FIAT element
4763
double *derivatives = new double [2*num_derivatives];
4765
// Declare coefficients
4766
double coeff0_0 = 0;
4767
double coeff0_1 = 0;
4768
double coeff0_2 = 0;
4769
double coeff1_0 = 0;
4770
double coeff1_1 = 0;
4771
double coeff1_2 = 0;
4773
// Declare new coefficients
4774
double new_coeff0_0 = 0;
4775
double new_coeff0_1 = 0;
4776
double new_coeff0_2 = 0;
4777
double new_coeff1_0 = 0;
4778
double new_coeff1_1 = 0;
4779
double new_coeff1_2 = 0;
4781
// Loop possible derivatives
4782
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
4784
// Get values from coefficients array
4785
new_coeff0_0 = coefficients0[dof][0];
4786
new_coeff0_1 = coefficients0[dof][1];
4787
new_coeff0_2 = coefficients0[dof][2];
4788
new_coeff1_0 = coefficients1[dof][0];
4789
new_coeff1_1 = coefficients1[dof][1];
4790
new_coeff1_2 = coefficients1[dof][2];
4792
// Loop derivative order
4793
for (unsigned int j = 0; j < n; j++)
4795
// Update old coefficients
4796
coeff0_0 = new_coeff0_0;
4797
coeff0_1 = new_coeff0_1;
4798
coeff0_2 = new_coeff0_2;
4799
coeff1_0 = new_coeff1_0;
4800
coeff1_1 = new_coeff1_1;
4801
coeff1_2 = new_coeff1_2;
4803
if(combinations[deriv_num][j] == 0)
4805
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
4806
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
4807
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
4808
new_coeff1_0 = coeff1_0*dmats0[0][0] + coeff1_1*dmats0[1][0] + coeff1_2*dmats0[2][0];
4809
new_coeff1_1 = coeff1_0*dmats0[0][1] + coeff1_1*dmats0[1][1] + coeff1_2*dmats0[2][1];
4810
new_coeff1_2 = coeff1_0*dmats0[0][2] + coeff1_1*dmats0[1][2] + coeff1_2*dmats0[2][2];
4812
if(combinations[deriv_num][j] == 1)
4814
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
4815
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
4816
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
4817
new_coeff1_0 = coeff1_0*dmats1[0][0] + coeff1_1*dmats1[1][0] + coeff1_2*dmats1[2][0];
4818
new_coeff1_1 = coeff1_0*dmats1[0][1] + coeff1_1*dmats1[1][1] + coeff1_2*dmats1[2][1];
4819
new_coeff1_2 = coeff1_0*dmats1[0][2] + coeff1_1*dmats1[1][2] + coeff1_2*dmats1[2][2];
4823
// Compute derivatives on reference element as dot product of coefficients and basisvalues
4824
// Correct values by the contravariant Piola transform
4825
const double tmp0_0 = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
4826
const double tmp0_1 = new_coeff1_0*basisvalue0 + new_coeff1_1*basisvalue1 + new_coeff1_2*basisvalue2;
4827
derivatives[deriv_num] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
4828
derivatives[num_derivatives + deriv_num] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
4831
// Transform derivatives back to physical element
4832
for (unsigned int row = 0; row < num_derivatives; row++)
4834
for (unsigned int col = 0; col < num_derivatives; col++)
4836
values[row] += transform[row][col]*derivatives[col];
4837
values[num_derivatives + row] += transform[row][col]*derivatives[num_derivatives + col];
4840
// Delete pointer to array of derivatives on FIAT element
4841
delete [] derivatives;
4843
// Delete pointer to array of combinations of derivatives and transform
4844
for (unsigned int row = 0; row < num_derivatives; row++)
4846
delete [] combinations[row];
4847
delete [] transform[row];
4850
delete [] combinations;
4851
delete [] transform;
4854
/// Evaluate order n derivatives of all basis functions at given point in cell
4855
virtual void evaluate_basis_derivatives_all(unsigned int n,
4857
const double* coordinates,
4858
const ufc::cell& c) const
4860
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
4863
/// Evaluate linear functional for dof i on the function f
4864
virtual double evaluate_dof(unsigned int i,
4865
const ufc::function& f,
4866
const ufc::cell& c) const
4868
// The reference points, direction and weights:
4869
static const double X[6][1][2] = {{{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}};
4870
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
4871
static const double D[6][1][2] = {{{1, 1}}, {{1, 1}}, {{1, 0}}, {{1, 0}}, {{0, -1}}, {{0, -1}}};
4873
// Extract vertex coordinates
4874
const double * const * x = c.coordinates;
4876
// Compute Jacobian of affine map from reference cell
4877
const double J_00 = x[1][0] - x[0][0];
4878
const double J_01 = x[2][0] - x[0][0];
4879
const double J_10 = x[1][1] - x[0][1];
4880
const double J_11 = x[2][1] - x[0][1];
4882
// Compute determinant of Jacobian
4883
double detJ = J_00*J_11 - J_01*J_10;
4885
// Compute inverse of Jacobian
4886
const double Jinv_00 = J_11 / detJ;
4887
const double Jinv_01 = -J_01 / detJ;
4888
const double Jinv_10 = -J_10 / detJ;
4889
const double Jinv_11 = J_00 / detJ;
4891
double copyofvalues[2];
4892
double result = 0.0;
4893
// Iterate over the points:
4894
// Evaluate basis functions for affine mapping
4895
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
4896
const double w1 = X[i][0][0];
4897
const double w2 = X[i][0][1];
4899
// Compute affine mapping y = F(X)
4901
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
4902
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
4904
// Evaluate function at physical points
4906
f.evaluate(values, y, c);
4908
// Map function values using appropriate mapping
4910
copyofvalues[0] = values[0];
4911
copyofvalues[1] = values[1];
4912
// Do the inverse of div piola
4913
values[0] = detJ*(Jinv_00*copyofvalues[0]+Jinv_01*copyofvalues[1]);
4914
values[1] = detJ*(Jinv_10*copyofvalues[0]+Jinv_11*copyofvalues[1]);
4916
// Note that we do not map the weights (yet).
4918
// Take directional components
4919
for(int k = 0; k < 2; k++)
4920
result += values[k]*D[i][0][k];
4921
// Multiply by weights
4927
/// Evaluate linear functionals for all dofs on the function f
4928
virtual void evaluate_dofs(double* values,
4929
const ufc::function& f,
4930
const ufc::cell& c) const
4932
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4935
/// Interpolate vertex values from dof values
4936
virtual void interpolate_vertex_values(double* vertex_values,
4937
const double* dof_values,
4938
const ufc::cell& c) const
4940
// Extract vertex coordinates
4941
const double * const * x = c.coordinates;
4943
// Compute Jacobian of affine map from reference cell
4944
const double J_00 = x[1][0] - x[0][0];
4945
const double J_01 = x[2][0] - x[0][0];
4946
const double J_10 = x[1][1] - x[0][1];
4947
const double J_11 = x[2][1] - x[0][1];
4949
// Compute determinant of Jacobian
4950
double detJ = J_00*J_11 - J_01*J_10;
4952
// Compute inverse of Jacobian
4953
// Evaluate at vertices and use Piola mapping
4954
vertex_values[0] = (1.0/detJ)*(dof_values[2]*2*J_00 + dof_values[3]*J_00 + dof_values[4]*(-2*J_01) + dof_values[5]*J_01);
4955
vertex_values[2] = (1.0/detJ)*(dof_values[0]*2*J_00 + dof_values[1]*J_00 + dof_values[4]*(J_00 + J_01) + dof_values[5]*(2*J_00 - 2*J_01));
4956
vertex_values[4] = (1.0/detJ)*(dof_values[0]*J_01 + dof_values[1]*2*J_01 + dof_values[2]*(J_00 + J_01) + dof_values[3]*(2*J_00 - 2*J_01));
4957
vertex_values[1] = (1.0/detJ)*(dof_values[2]*2*J_10 + dof_values[3]*J_10 + dof_values[4]*(-2*J_11) + dof_values[5]*J_11);
4958
vertex_values[3] = (1.0/detJ)*(dof_values[0]*2*J_10 + dof_values[1]*J_10 + dof_values[4]*(J_10 + J_11) + dof_values[5]*(2*J_10 - 2*J_11));
4959
vertex_values[5] = (1.0/detJ)*(dof_values[0]*J_11 + dof_values[1]*2*J_11 + dof_values[2]*(J_10 + J_11) + dof_values[3]*(2*J_10 - 2*J_11));
4962
/// Return the number of sub elements (for a mixed element)
4963
virtual unsigned int num_sub_elements() const
4968
/// Create a new finite element for sub element i (for a mixed element)
4969
virtual ufc::finite_element* create_sub_element(unsigned int i) const
4971
return new mixedpoisson_1_finite_element_0_0();
4976
/// This class defines the interface for a finite element.
4978
class mixedpoisson_1_finite_element_0_1: public ufc::finite_element
4983
mixedpoisson_1_finite_element_0_1() : ufc::finite_element()
4989
virtual ~mixedpoisson_1_finite_element_0_1()
4994
/// Return a string identifying the finite element
4995
virtual const char* signature() const
4997
return "FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
5000
/// Return the cell shape
5001
virtual ufc::shape cell_shape() const
5003
return ufc::triangle;
5006
/// Return the dimension of the finite element function space
5007
virtual unsigned int space_dimension() const
5012
/// Return the rank of the value space
5013
virtual unsigned int value_rank() const
5018
/// Return the dimension of the value space for axis i
5019
virtual unsigned int value_dimension(unsigned int i) const
5024
/// Evaluate basis function i at given point in cell
5025
virtual void evaluate_basis(unsigned int i,
5027
const double* coordinates,
5028
const ufc::cell& c) const
5030
// Extract vertex coordinates
5031
const double * const * element_coordinates = c.coordinates;
5033
// Compute Jacobian of affine map from reference cell
5034
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
5035
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
5036
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
5037
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
5039
// Compute determinant of Jacobian
5040
const double detJ = J_00*J_11 - J_01*J_10;
5042
// Compute inverse of Jacobian
5044
// Get coordinates and map to the reference (UFC) element
5045
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
5046
element_coordinates[0][0]*element_coordinates[2][1] +\
5047
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
5048
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
5049
element_coordinates[1][0]*element_coordinates[0][1] -\
5050
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
5052
// Map coordinates to the reference square
5053
if (std::abs(y - 1.0) < 1e-08)
5056
x = 2.0 *x/(1.0 - y) - 1.0;
5062
// Map degree of freedom to element degree of freedom
5063
const unsigned int dof = i;
5065
// Generate scalings
5066
const double scalings_y_0 = 1;
5068
// Compute psitilde_a
5069
const double psitilde_a_0 = 1;
5071
// Compute psitilde_bs
5072
const double psitilde_bs_0_0 = 1;
5074
// Compute basisvalues
5075
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
5077
// Table(s) of coefficients
5078
static const double coefficients0[1][1] = \
5081
// Extract relevant coefficients
5082
const double coeff0_0 = coefficients0[dof][0];
5085
*values = coeff0_0*basisvalue0;
5088
/// Evaluate all basis functions at given point in cell
5089
virtual void evaluate_basis_all(double* values,
5090
const double* coordinates,
5091
const ufc::cell& c) const
5093
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
5096
/// Evaluate order n derivatives of basis function i at given point in cell
5097
virtual void evaluate_basis_derivatives(unsigned int i,
5100
const double* coordinates,
5101
const ufc::cell& c) const
5103
// Extract vertex coordinates
5104
const double * const * element_coordinates = c.coordinates;
5106
// Compute Jacobian of affine map from reference cell
5107
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
5108
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
5109
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
5110
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
5112
// Compute determinant of Jacobian
5113
const double detJ = J_00*J_11 - J_01*J_10;
5115
// Compute inverse of Jacobian
5117
// Get coordinates and map to the reference (UFC) element
5118
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
5119
element_coordinates[0][0]*element_coordinates[2][1] +\
5120
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
5121
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
5122
element_coordinates[1][0]*element_coordinates[0][1] -\
5123
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
5125
// Map coordinates to the reference square
5126
if (std::abs(y - 1.0) < 1e-08)
5129
x = 2.0 *x/(1.0 - y) - 1.0;
5132
// Compute number of derivatives
5133
unsigned int num_derivatives = 1;
5135
for (unsigned int j = 0; j < n; j++)
5136
num_derivatives *= 2;
5139
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
5140
unsigned int **combinations = new unsigned int *[num_derivatives];
5142
for (unsigned int j = 0; j < num_derivatives; j++)
5144
combinations[j] = new unsigned int [n];
5145
for (unsigned int k = 0; k < n; k++)
5146
combinations[j][k] = 0;
5149
// Generate combinations of derivatives
5150
for (unsigned int row = 1; row < num_derivatives; row++)
5152
for (unsigned int num = 0; num < row; num++)
5154
for (unsigned int col = n-1; col+1 > 0; col--)
5156
if (combinations[row][col] + 1 > 1)
5157
combinations[row][col] = 0;
5160
combinations[row][col] += 1;
5167
// Compute inverse of Jacobian
5168
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
5170
// Declare transformation matrix
5171
// Declare pointer to two dimensional array and initialise
5172
double **transform = new double *[num_derivatives];
5174
for (unsigned int j = 0; j < num_derivatives; j++)
5176
transform[j] = new double [num_derivatives];
5177
for (unsigned int k = 0; k < num_derivatives; k++)
5178
transform[j][k] = 1;
5181
// Construct transformation matrix
5182
for (unsigned int row = 0; row < num_derivatives; row++)
5184
for (unsigned int col = 0; col < num_derivatives; col++)
5186
for (unsigned int k = 0; k < n; k++)
5187
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
5192
for (unsigned int j = 0; j < 1*num_derivatives; j++)
5195
// Map degree of freedom to element degree of freedom
5196
const unsigned int dof = i;
5198
// Generate scalings
5199
const double scalings_y_0 = 1;
5201
// Compute psitilde_a
5202
const double psitilde_a_0 = 1;
5204
// Compute psitilde_bs
5205
const double psitilde_bs_0_0 = 1;
5207
// Compute basisvalues
5208
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
5210
// Table(s) of coefficients
5211
static const double coefficients0[1][1] = \
5214
// Interesting (new) part
5215
// Tables of derivatives of the polynomial base (transpose)
5216
static const double dmats0[1][1] = \
5219
static const double dmats1[1][1] = \
5222
// Compute reference derivatives
5223
// Declare pointer to array of derivatives on FIAT element
5224
double *derivatives = new double [num_derivatives];
5226
// Declare coefficients
5227
double coeff0_0 = 0;
5229
// Declare new coefficients
5230
double new_coeff0_0 = 0;
5232
// Loop possible derivatives
5233
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
5235
// Get values from coefficients array
5236
new_coeff0_0 = coefficients0[dof][0];
5238
// Loop derivative order
5239
for (unsigned int j = 0; j < n; j++)
5241
// Update old coefficients
5242
coeff0_0 = new_coeff0_0;
5244
if(combinations[deriv_num][j] == 0)
5246
new_coeff0_0 = coeff0_0*dmats0[0][0];
5248
if(combinations[deriv_num][j] == 1)
5250
new_coeff0_0 = coeff0_0*dmats1[0][0];
5254
// Compute derivatives on reference element as dot product of coefficients and basisvalues
5255
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
5258
// Transform derivatives back to physical element
5259
for (unsigned int row = 0; row < num_derivatives; row++)
5261
for (unsigned int col = 0; col < num_derivatives; col++)
5263
values[row] += transform[row][col]*derivatives[col];
5266
// Delete pointer to array of derivatives on FIAT element
5267
delete [] derivatives;
5269
// Delete pointer to array of combinations of derivatives and transform
5270
for (unsigned int row = 0; row < num_derivatives; row++)
5272
delete [] combinations[row];
5273
delete [] transform[row];
5276
delete [] combinations;
5277
delete [] transform;
5280
/// Evaluate order n derivatives of all basis functions at given point in cell
5281
virtual void evaluate_basis_derivatives_all(unsigned int n,
5283
const double* coordinates,
5284
const ufc::cell& c) const
5286
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
5289
/// Evaluate linear functional for dof i on the function f
5290
virtual double evaluate_dof(unsigned int i,
5291
const ufc::function& f,
5292
const ufc::cell& c) const
5294
// The reference points, direction and weights:
5295
static const double X[1][1][2] = {{{0.333333333, 0.333333333}}};
5296
static const double W[1][1] = {{1}};
5297
static const double D[1][1][1] = {{{1}}};
5299
const double * const * x = c.coordinates;
5300
double result = 0.0;
5301
// Iterate over the points:
5302
// Evaluate basis functions for affine mapping
5303
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
5304
const double w1 = X[i][0][0];
5305
const double w2 = X[i][0][1];
5307
// Compute affine mapping y = F(X)
5309
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
5310
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
5312
// Evaluate function at physical points
5314
f.evaluate(values, y, c);
5316
// Map function values using appropriate mapping
5317
// Affine map: Do nothing
5319
// Note that we do not map the weights (yet).
5321
// Take directional components
5322
for(int k = 0; k < 1; k++)
5323
result += values[k]*D[i][0][k];
5324
// Multiply by weights
5330
/// Evaluate linear functionals for all dofs on the function f
5331
virtual void evaluate_dofs(double* values,
5332
const ufc::function& f,
5333
const ufc::cell& c) const
5335
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
5338
/// Interpolate vertex values from dof values
5339
virtual void interpolate_vertex_values(double* vertex_values,
5340
const double* dof_values,
5341
const ufc::cell& c) const
5343
// Evaluate at vertices and use affine mapping
5344
vertex_values[0] = dof_values[0];
5345
vertex_values[1] = dof_values[0];
5346
vertex_values[2] = dof_values[0];
5349
/// Return the number of sub elements (for a mixed element)
5350
virtual unsigned int num_sub_elements() const
5355
/// Create a new finite element for sub element i (for a mixed element)
5356
virtual ufc::finite_element* create_sub_element(unsigned int i) const
5358
return new mixedpoisson_1_finite_element_0_1();
5363
/// This class defines the interface for a finite element.
5365
class mixedpoisson_1_finite_element_0: public ufc::finite_element
5370
mixedpoisson_1_finite_element_0() : ufc::finite_element()
5376
virtual ~mixedpoisson_1_finite_element_0()
5381
/// Return a string identifying the finite element
5382
virtual const char* signature() const
5384
return "MixedElement([FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1), FiniteElement('Discontinuous Lagrange', 'triangle', 0)])";
5387
/// Return the cell shape
5388
virtual ufc::shape cell_shape() const
5390
return ufc::triangle;
5393
/// Return the dimension of the finite element function space
5394
virtual unsigned int space_dimension() const
5399
/// Return the rank of the value space
5400
virtual unsigned int value_rank() const
5405
/// Return the dimension of the value space for axis i
5406
virtual unsigned int value_dimension(unsigned int i) const
5411
/// Evaluate basis function i at given point in cell
5412
virtual void evaluate_basis(unsigned int i,
5414
const double* coordinates,
5415
const ufc::cell& c) const
5417
// Extract vertex coordinates
5418
const double * const * element_coordinates = c.coordinates;
5420
// Compute Jacobian of affine map from reference cell
5421
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
5422
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
5423
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
5424
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
5426
// Compute determinant of Jacobian
5427
const double detJ = J_00*J_11 - J_01*J_10;
5429
// Compute inverse of Jacobian
5431
// Get coordinates and map to the reference (UFC) element
5432
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
5433
element_coordinates[0][0]*element_coordinates[2][1] +\
5434
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
5435
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
5436
element_coordinates[1][0]*element_coordinates[0][1] -\
5437
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
5439
// Map coordinates to the reference square
5440
if (std::abs(y - 1.0) < 1e-08)
5443
x = 2.0 *x/(1.0 - y) - 1.0;
5451
if (0 <= i && i <= 5)
5453
// Map degree of freedom to element degree of freedom
5454
const unsigned int dof = i;
5456
// Generate scalings
5457
const double scalings_y_0 = 1;
5458
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
5460
// Compute psitilde_a
5461
const double psitilde_a_0 = 1;
5462
const double psitilde_a_1 = x;
5464
// Compute psitilde_bs
5465
const double psitilde_bs_0_0 = 1;
5466
const double psitilde_bs_0_1 = 1.5*y + 0.5;
5467
const double psitilde_bs_1_0 = 1;
5469
// Compute basisvalues
5470
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
5471
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
5472
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
5474
// Table(s) of coefficients
5475
static const double coefficients0[6][3] = \
5476
{{0.942809042, 0.577350269, -0.333333333},
5477
{-0.471404521, -0.288675135, 0.166666667},
5478
{0.471404521, -0.577350269, -0.666666667},
5479
{0.471404521, 0.288675135, 0.833333333},
5480
{-0.471404521, -0.288675135, 0.166666667},
5481
{0.942809042, 0.577350269, -0.333333333}};
5483
static const double coefficients1[6][3] = \
5484
{{-0.471404521, 0, -0.333333333},
5485
{0.942809042, 0, 0.666666667},
5486
{0.471404521, 0, 0.333333333},
5487
{-0.942809042, 0, -0.666666667},
5488
{-0.471404521, 0.866025404, 0.166666667},
5489
{-0.471404521, -0.866025404, 0.166666667}};
5491
// Extract relevant coefficients
5492
const double coeff0_0 = coefficients0[dof][0];
5493
const double coeff0_1 = coefficients0[dof][1];
5494
const double coeff0_2 = coefficients0[dof][2];
5495
const double coeff1_0 = coefficients1[dof][0];
5496
const double coeff1_1 = coefficients1[dof][1];
5497
const double coeff1_2 = coefficients1[dof][2];
5500
const double tmp0_0 = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
5501
const double tmp0_1 = coeff1_0*basisvalue0 + coeff1_1*basisvalue1 + coeff1_2*basisvalue2;
5502
// Using contravariant Piola transform to map values back to the physical element
5503
values[0] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
5504
values[1] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
5507
if (6 <= i && i <= 6)
5509
// Map degree of freedom to element degree of freedom
5510
const unsigned int dof = i - 6;
5512
// Generate scalings
5513
const double scalings_y_0 = 1;
5515
// Compute psitilde_a
5516
const double psitilde_a_0 = 1;
5518
// Compute psitilde_bs
5519
const double psitilde_bs_0_0 = 1;
5521
// Compute basisvalues
5522
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
5524
// Table(s) of coefficients
5525
static const double coefficients0[1][1] = \
5528
// Extract relevant coefficients
5529
const double coeff0_0 = coefficients0[dof][0];
5532
values[2] = coeff0_0*basisvalue0;
5537
/// Evaluate all basis functions at given point in cell
5538
virtual void evaluate_basis_all(double* values,
5539
const double* coordinates,
5540
const ufc::cell& c) const
5542
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
5545
/// Evaluate order n derivatives of basis function i at given point in cell
5546
virtual void evaluate_basis_derivatives(unsigned int i,
5549
const double* coordinates,
5550
const ufc::cell& c) const
5552
// Extract vertex coordinates
5553
const double * const * element_coordinates = c.coordinates;
5555
// Compute Jacobian of affine map from reference cell
5556
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
5557
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
5558
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
5559
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
5561
// Compute determinant of Jacobian
5562
const double detJ = J_00*J_11 - J_01*J_10;
5564
// Compute inverse of Jacobian
5566
// Get coordinates and map to the reference (UFC) element
5567
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
5568
element_coordinates[0][0]*element_coordinates[2][1] +\
5569
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
5570
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
5571
element_coordinates[1][0]*element_coordinates[0][1] -\
5572
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
5574
// Map coordinates to the reference square
5575
if (std::abs(y - 1.0) < 1e-08)
5578
x = 2.0 *x/(1.0 - y) - 1.0;
5581
// Compute number of derivatives
5582
unsigned int num_derivatives = 1;
5584
for (unsigned int j = 0; j < n; j++)
5585
num_derivatives *= 2;
5588
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
5589
unsigned int **combinations = new unsigned int *[num_derivatives];
5591
for (unsigned int j = 0; j < num_derivatives; j++)
5593
combinations[j] = new unsigned int [n];
5594
for (unsigned int k = 0; k < n; k++)
5595
combinations[j][k] = 0;
5598
// Generate combinations of derivatives
5599
for (unsigned int row = 1; row < num_derivatives; row++)
5601
for (unsigned int num = 0; num < row; num++)
5603
for (unsigned int col = n-1; col+1 > 0; col--)
5605
if (combinations[row][col] + 1 > 1)
5606
combinations[row][col] = 0;
5609
combinations[row][col] += 1;
5616
// Compute inverse of Jacobian
5617
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
5619
// Declare transformation matrix
5620
// Declare pointer to two dimensional array and initialise
5621
double **transform = new double *[num_derivatives];
5623
for (unsigned int j = 0; j < num_derivatives; j++)
5625
transform[j] = new double [num_derivatives];
5626
for (unsigned int k = 0; k < num_derivatives; k++)
5627
transform[j][k] = 1;
5630
// Construct transformation matrix
5631
for (unsigned int row = 0; row < num_derivatives; row++)
5633
for (unsigned int col = 0; col < num_derivatives; col++)
5635
for (unsigned int k = 0; k < n; k++)
5636
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
5641
for (unsigned int j = 0; j < 3*num_derivatives; j++)
5644
if (0 <= i && i <= 5)
5646
// Map degree of freedom to element degree of freedom
5647
const unsigned int dof = i;
5649
// Generate scalings
5650
const double scalings_y_0 = 1;
5651
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
5653
// Compute psitilde_a
5654
const double psitilde_a_0 = 1;
5655
const double psitilde_a_1 = x;
5657
// Compute psitilde_bs
5658
const double psitilde_bs_0_0 = 1;
5659
const double psitilde_bs_0_1 = 1.5*y + 0.5;
5660
const double psitilde_bs_1_0 = 1;
5662
// Compute basisvalues
5663
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
5664
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
5665
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
5667
// Table(s) of coefficients
5668
static const double coefficients0[6][3] = \
5669
{{0.942809042, 0.577350269, -0.333333333},
5670
{-0.471404521, -0.288675135, 0.166666667},
5671
{0.471404521, -0.577350269, -0.666666667},
5672
{0.471404521, 0.288675135, 0.833333333},
5673
{-0.471404521, -0.288675135, 0.166666667},
5674
{0.942809042, 0.577350269, -0.333333333}};
5676
static const double coefficients1[6][3] = \
5677
{{-0.471404521, 0, -0.333333333},
5678
{0.942809042, 0, 0.666666667},
5679
{0.471404521, 0, 0.333333333},
5680
{-0.942809042, 0, -0.666666667},
5681
{-0.471404521, 0.866025404, 0.166666667},
5682
{-0.471404521, -0.866025404, 0.166666667}};
5684
// Interesting (new) part
5685
// Tables of derivatives of the polynomial base (transpose)
5686
static const double dmats0[3][3] = \
5691
static const double dmats1[3][3] = \
5694
{4.24264069, 0, 0}};
5696
// Compute reference derivatives
5697
// Declare pointer to array of derivatives on FIAT element
5698
double *derivatives = new double [2*num_derivatives];
5700
// Declare coefficients
5701
double coeff0_0 = 0;
5702
double coeff0_1 = 0;
5703
double coeff0_2 = 0;
5704
double coeff1_0 = 0;
5705
double coeff1_1 = 0;
5706
double coeff1_2 = 0;
5708
// Declare new coefficients
5709
double new_coeff0_0 = 0;
5710
double new_coeff0_1 = 0;
5711
double new_coeff0_2 = 0;
5712
double new_coeff1_0 = 0;
5713
double new_coeff1_1 = 0;
5714
double new_coeff1_2 = 0;
5716
// Loop possible derivatives
5717
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
5719
// Get values from coefficients array
5720
new_coeff0_0 = coefficients0[dof][0];
5721
new_coeff0_1 = coefficients0[dof][1];
5722
new_coeff0_2 = coefficients0[dof][2];
5723
new_coeff1_0 = coefficients1[dof][0];
5724
new_coeff1_1 = coefficients1[dof][1];
5725
new_coeff1_2 = coefficients1[dof][2];
5727
// Loop derivative order
5728
for (unsigned int j = 0; j < n; j++)
5730
// Update old coefficients
5731
coeff0_0 = new_coeff0_0;
5732
coeff0_1 = new_coeff0_1;
5733
coeff0_2 = new_coeff0_2;
5734
coeff1_0 = new_coeff1_0;
5735
coeff1_1 = new_coeff1_1;
5736
coeff1_2 = new_coeff1_2;
5738
if(combinations[deriv_num][j] == 0)
5740
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
5741
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
5742
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
5743
new_coeff1_0 = coeff1_0*dmats0[0][0] + coeff1_1*dmats0[1][0] + coeff1_2*dmats0[2][0];
5744
new_coeff1_1 = coeff1_0*dmats0[0][1] + coeff1_1*dmats0[1][1] + coeff1_2*dmats0[2][1];
5745
new_coeff1_2 = coeff1_0*dmats0[0][2] + coeff1_1*dmats0[1][2] + coeff1_2*dmats0[2][2];
5747
if(combinations[deriv_num][j] == 1)
5749
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
5750
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
5751
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
5752
new_coeff1_0 = coeff1_0*dmats1[0][0] + coeff1_1*dmats1[1][0] + coeff1_2*dmats1[2][0];
5753
new_coeff1_1 = coeff1_0*dmats1[0][1] + coeff1_1*dmats1[1][1] + coeff1_2*dmats1[2][1];
5754
new_coeff1_2 = coeff1_0*dmats1[0][2] + coeff1_1*dmats1[1][2] + coeff1_2*dmats1[2][2];
5758
// Compute derivatives on reference element as dot product of coefficients and basisvalues
5759
// Correct values by the contravariant Piola transform
5760
const double tmp0_0 = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
5761
const double tmp0_1 = new_coeff1_0*basisvalue0 + new_coeff1_1*basisvalue1 + new_coeff1_2*basisvalue2;
5762
derivatives[deriv_num] = (1.0/detJ)*(J_00*tmp0_0 + J_01*tmp0_1);
5763
derivatives[num_derivatives + deriv_num] = (1.0/detJ)*(J_10*tmp0_0 + J_11*tmp0_1);
5766
// Transform derivatives back to physical element
5767
for (unsigned int row = 0; row < num_derivatives; row++)
5769
for (unsigned int col = 0; col < num_derivatives; col++)
5771
values[row] += transform[row][col]*derivatives[col];
5772
values[num_derivatives + row] += transform[row][col]*derivatives[num_derivatives + col];
5775
// Delete pointer to array of derivatives on FIAT element
5776
delete [] derivatives;
5778
// Delete pointer to array of combinations of derivatives and transform
5779
for (unsigned int row = 0; row < num_derivatives; row++)
5781
delete [] combinations[row];
5782
delete [] transform[row];
5785
delete [] combinations;
5786
delete [] transform;
5789
if (6 <= i && i <= 6)
5791
// Map degree of freedom to element degree of freedom
5792
const unsigned int dof = i - 6;
5794
// Generate scalings
5795
const double scalings_y_0 = 1;
5797
// Compute psitilde_a
5798
const double psitilde_a_0 = 1;
5800
// Compute psitilde_bs
5801
const double psitilde_bs_0_0 = 1;
5803
// Compute basisvalues
5804
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
5806
// Table(s) of coefficients
5807
static const double coefficients0[1][1] = \
5810
// Interesting (new) part
5811
// Tables of derivatives of the polynomial base (transpose)
5812
static const double dmats0[1][1] = \
5815
static const double dmats1[1][1] = \
5818
// Compute reference derivatives
5819
// Declare pointer to array of derivatives on FIAT element
5820
double *derivatives = new double [num_derivatives];
5822
// Declare coefficients
5823
double coeff0_0 = 0;
5825
// Declare new coefficients
5826
double new_coeff0_0 = 0;
5828
// Loop possible derivatives
5829
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
5831
// Get values from coefficients array
5832
new_coeff0_0 = coefficients0[dof][0];
5834
// Loop derivative order
5835
for (unsigned int j = 0; j < n; j++)
5837
// Update old coefficients
5838
coeff0_0 = new_coeff0_0;
5840
if(combinations[deriv_num][j] == 0)
5842
new_coeff0_0 = coeff0_0*dmats0[0][0];
5844
if(combinations[deriv_num][j] == 1)
5846
new_coeff0_0 = coeff0_0*dmats1[0][0];
5850
// Compute derivatives on reference element as dot product of coefficients and basisvalues
5851
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
5854
// Transform derivatives back to physical element
5855
for (unsigned int row = 0; row < num_derivatives; row++)
5857
for (unsigned int col = 0; col < num_derivatives; col++)
5859
values[2*num_derivatives + row] += transform[row][col]*derivatives[col];
5862
// Delete pointer to array of derivatives on FIAT element
5863
delete [] derivatives;
5865
// Delete pointer to array of combinations of derivatives and transform
5866
for (unsigned int row = 0; row < num_derivatives; row++)
5868
delete [] combinations[row];
5869
delete [] transform[row];
5872
delete [] combinations;
5873
delete [] transform;
5878
/// Evaluate order n derivatives of all basis functions at given point in cell
5879
virtual void evaluate_basis_derivatives_all(unsigned int n,
5881
const double* coordinates,
5882
const ufc::cell& c) const
5884
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
5887
/// Evaluate linear functional for dof i on the function f
5888
virtual double evaluate_dof(unsigned int i,
5889
const ufc::function& f,
5890
const ufc::cell& c) const
5892
// The reference points, direction and weights:
5893
static const double X[7][1][2] = {{{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}, {{0.333333333, 0.333333333}}};
5894
static const double W[7][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}};
5895
static const double D[7][1][3] = {{{1, 1, 0}}, {{1, 1, 0}}, {{1, 0, 0}}, {{1, 0, 0}}, {{0, -1, 0}}, {{0, -1, 0}}, {{0, 0, 1}}};
5897
static const unsigned int mappings[7] = {1, 1, 1, 1, 1, 1, 0};
5898
// Extract vertex coordinates
5899
const double * const * x = c.coordinates;
5901
// Compute Jacobian of affine map from reference cell
5902
const double J_00 = x[1][0] - x[0][0];
5903
const double J_01 = x[2][0] - x[0][0];
5904
const double J_10 = x[1][1] - x[0][1];
5905
const double J_11 = x[2][1] - x[0][1];
5907
// Compute determinant of Jacobian
5908
double detJ = J_00*J_11 - J_01*J_10;
5910
// Compute inverse of Jacobian
5911
const double Jinv_00 = J_11 / detJ;
5912
const double Jinv_01 = -J_01 / detJ;
5913
const double Jinv_10 = -J_10 / detJ;
5914
const double Jinv_11 = J_00 / detJ;
5916
double copyofvalues[3];
5917
double result = 0.0;
5918
// Iterate over the points:
5919
// Evaluate basis functions for affine mapping
5920
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
5921
const double w1 = X[i][0][0];
5922
const double w2 = X[i][0][1];
5924
// Compute affine mapping y = F(X)
5926
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
5927
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
5929
// Evaluate function at physical points
5931
f.evaluate(values, y, c);
5933
// Map function values using appropriate mapping
5934
if (mappings[i] == 0) {
5935
// Affine map: Do nothing
5936
} else if (mappings[i] == 1) {
5938
copyofvalues[0] = values[0];
5939
copyofvalues[1] = values[1];
5940
// Do the inverse of div piola
5941
values[0] = detJ*(Jinv_00*copyofvalues[0]+Jinv_01*copyofvalues[1]);
5942
values[1] = detJ*(Jinv_10*copyofvalues[0]+Jinv_11*copyofvalues[1]);
5944
// Other mappings not applicable.
5947
// Note that we do not map the weights (yet).
5949
// Take directional components
5950
for(int k = 0; k < 3; k++)
5951
result += values[k]*D[i][0][k];
5952
// Multiply by weights
5958
/// Evaluate linear functionals for all dofs on the function f
5959
virtual void evaluate_dofs(double* values,
5960
const ufc::function& f,
5961
const ufc::cell& c) const
5963
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
5966
/// Interpolate vertex values from dof values
5967
virtual void interpolate_vertex_values(double* vertex_values,
5968
const double* dof_values,
5969
const ufc::cell& c) const
5971
// Extract vertex coordinates
5972
const double * const * x = c.coordinates;
5974
// Compute Jacobian of affine map from reference cell
5975
const double J_00 = x[1][0] - x[0][0];
5976
const double J_01 = x[2][0] - x[0][0];
5977
const double J_10 = x[1][1] - x[0][1];
5978
const double J_11 = x[2][1] - x[0][1];
5980
// Compute determinant of Jacobian
5981
double detJ = J_00*J_11 - J_01*J_10;
5983
// Compute inverse of Jacobian
5984
// Evaluate at vertices and use Piola mapping
5985
vertex_values[0] = (1.0/detJ)*(dof_values[2]*2*J_00 + dof_values[3]*J_00 + dof_values[4]*(-2*J_01) + dof_values[5]*J_01);
5986
vertex_values[3] = (1.0/detJ)*(dof_values[0]*2*J_00 + dof_values[1]*J_00 + dof_values[4]*(J_00 + J_01) + dof_values[5]*(2*J_00 - 2*J_01));
5987
vertex_values[6] = (1.0/detJ)*(dof_values[0]*J_01 + dof_values[1]*2*J_01 + dof_values[2]*(J_00 + J_01) + dof_values[3]*(2*J_00 - 2*J_01));
5988
vertex_values[1] = (1.0/detJ)*(dof_values[2]*2*J_10 + dof_values[3]*J_10 + dof_values[4]*(-2*J_11) + dof_values[5]*J_11);
5989
vertex_values[4] = (1.0/detJ)*(dof_values[0]*2*J_10 + dof_values[1]*J_10 + dof_values[4]*(J_10 + J_11) + dof_values[5]*(2*J_10 - 2*J_11));
5990
vertex_values[7] = (1.0/detJ)*(dof_values[0]*J_11 + dof_values[1]*2*J_11 + dof_values[2]*(J_10 + J_11) + dof_values[3]*(2*J_10 - 2*J_11));
5991
// Evaluate at vertices and use affine mapping
5992
vertex_values[2] = dof_values[6];
5993
vertex_values[5] = dof_values[6];
5994
vertex_values[8] = dof_values[6];
5997
/// Return the number of sub elements (for a mixed element)
5998
virtual unsigned int num_sub_elements() const
6003
/// Create a new finite element for sub element i (for a mixed element)
6004
virtual ufc::finite_element* create_sub_element(unsigned int i) const
6009
return new mixedpoisson_1_finite_element_0_0();
6012
return new mixedpoisson_1_finite_element_0_1();
6020
/// This class defines the interface for a finite element.
6022
class mixedpoisson_1_finite_element_1: public ufc::finite_element
6027
mixedpoisson_1_finite_element_1() : ufc::finite_element()
6033
virtual ~mixedpoisson_1_finite_element_1()
6038
/// Return a string identifying the finite element
6039
virtual const char* signature() const
6041
return "FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
6044
/// Return the cell shape
6045
virtual ufc::shape cell_shape() const
6047
return ufc::triangle;
6050
/// Return the dimension of the finite element function space
6051
virtual unsigned int space_dimension() const
6056
/// Return the rank of the value space
6057
virtual unsigned int value_rank() const
6062
/// Return the dimension of the value space for axis i
6063
virtual unsigned int value_dimension(unsigned int i) const
6068
/// Evaluate basis function i at given point in cell
6069
virtual void evaluate_basis(unsigned int i,
6071
const double* coordinates,
6072
const ufc::cell& c) const
6074
// Extract vertex coordinates
6075
const double * const * element_coordinates = c.coordinates;
6077
// Compute Jacobian of affine map from reference cell
6078
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
6079
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
6080
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
6081
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
6083
// Compute determinant of Jacobian
6084
const double detJ = J_00*J_11 - J_01*J_10;
6086
// Compute inverse of Jacobian
6088
// Get coordinates and map to the reference (UFC) element
6089
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
6090
element_coordinates[0][0]*element_coordinates[2][1] +\
6091
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
6092
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
6093
element_coordinates[1][0]*element_coordinates[0][1] -\
6094
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
6096
// Map coordinates to the reference square
6097
if (std::abs(y - 1.0) < 1e-08)
6100
x = 2.0 *x/(1.0 - y) - 1.0;
6106
// Map degree of freedom to element degree of freedom
6107
const unsigned int dof = i;
6109
// Generate scalings
6110
const double scalings_y_0 = 1;
6112
// Compute psitilde_a
6113
const double psitilde_a_0 = 1;
6115
// Compute psitilde_bs
6116
const double psitilde_bs_0_0 = 1;
6118
// Compute basisvalues
6119
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
6121
// Table(s) of coefficients
6122
static const double coefficients0[1][1] = \
6125
// Extract relevant coefficients
6126
const double coeff0_0 = coefficients0[dof][0];
6129
*values = coeff0_0*basisvalue0;
6132
/// Evaluate all basis functions at given point in cell
6133
virtual void evaluate_basis_all(double* values,
6134
const double* coordinates,
6135
const ufc::cell& c) const
6137
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
6140
/// Evaluate order n derivatives of basis function i at given point in cell
6141
virtual void evaluate_basis_derivatives(unsigned int i,
6144
const double* coordinates,
6145
const ufc::cell& c) const
6147
// Extract vertex coordinates
6148
const double * const * element_coordinates = c.coordinates;
6150
// Compute Jacobian of affine map from reference cell
6151
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
6152
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
6153
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
6154
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
6156
// Compute determinant of Jacobian
6157
const double detJ = J_00*J_11 - J_01*J_10;
6159
// Compute inverse of Jacobian
6161
// Get coordinates and map to the reference (UFC) element
6162
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
6163
element_coordinates[0][0]*element_coordinates[2][1] +\
6164
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
6165
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
6166
element_coordinates[1][0]*element_coordinates[0][1] -\
6167
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
6169
// Map coordinates to the reference square
6170
if (std::abs(y - 1.0) < 1e-08)
6173
x = 2.0 *x/(1.0 - y) - 1.0;
6176
// Compute number of derivatives
6177
unsigned int num_derivatives = 1;
6179
for (unsigned int j = 0; j < n; j++)
6180
num_derivatives *= 2;
6183
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
6184
unsigned int **combinations = new unsigned int *[num_derivatives];
6186
for (unsigned int j = 0; j < num_derivatives; j++)
6188
combinations[j] = new unsigned int [n];
6189
for (unsigned int k = 0; k < n; k++)
6190
combinations[j][k] = 0;
6193
// Generate combinations of derivatives
6194
for (unsigned int row = 1; row < num_derivatives; row++)
6196
for (unsigned int num = 0; num < row; num++)
6198
for (unsigned int col = n-1; col+1 > 0; col--)
6200
if (combinations[row][col] + 1 > 1)
6201
combinations[row][col] = 0;
6204
combinations[row][col] += 1;
6211
// Compute inverse of Jacobian
6212
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
6214
// Declare transformation matrix
6215
// Declare pointer to two dimensional array and initialise
6216
double **transform = new double *[num_derivatives];
6218
for (unsigned int j = 0; j < num_derivatives; j++)
6220
transform[j] = new double [num_derivatives];
6221
for (unsigned int k = 0; k < num_derivatives; k++)
6222
transform[j][k] = 1;
6225
// Construct transformation matrix
6226
for (unsigned int row = 0; row < num_derivatives; row++)
6228
for (unsigned int col = 0; col < num_derivatives; col++)
6230
for (unsigned int k = 0; k < n; k++)
6231
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
6236
for (unsigned int j = 0; j < 1*num_derivatives; j++)
6239
// Map degree of freedom to element degree of freedom
6240
const unsigned int dof = i;
6242
// Generate scalings
6243
const double scalings_y_0 = 1;
6245
// Compute psitilde_a
6246
const double psitilde_a_0 = 1;
6248
// Compute psitilde_bs
6249
const double psitilde_bs_0_0 = 1;
6251
// Compute basisvalues
6252
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
6254
// Table(s) of coefficients
6255
static const double coefficients0[1][1] = \
6258
// Interesting (new) part
6259
// Tables of derivatives of the polynomial base (transpose)
6260
static const double dmats0[1][1] = \
6263
static const double dmats1[1][1] = \
6266
// Compute reference derivatives
6267
// Declare pointer to array of derivatives on FIAT element
6268
double *derivatives = new double [num_derivatives];
6270
// Declare coefficients
6271
double coeff0_0 = 0;
6273
// Declare new coefficients
6274
double new_coeff0_0 = 0;
6276
// Loop possible derivatives
6277
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
6279
// Get values from coefficients array
6280
new_coeff0_0 = coefficients0[dof][0];
6282
// Loop derivative order
6283
for (unsigned int j = 0; j < n; j++)
6285
// Update old coefficients
6286
coeff0_0 = new_coeff0_0;
6288
if(combinations[deriv_num][j] == 0)
6290
new_coeff0_0 = coeff0_0*dmats0[0][0];
6292
if(combinations[deriv_num][j] == 1)
6294
new_coeff0_0 = coeff0_0*dmats1[0][0];
6298
// Compute derivatives on reference element as dot product of coefficients and basisvalues
6299
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
6302
// Transform derivatives back to physical element
6303
for (unsigned int row = 0; row < num_derivatives; row++)
6305
for (unsigned int col = 0; col < num_derivatives; col++)
6307
values[row] += transform[row][col]*derivatives[col];
6310
// Delete pointer to array of derivatives on FIAT element
6311
delete [] derivatives;
6313
// Delete pointer to array of combinations of derivatives and transform
6314
for (unsigned int row = 0; row < num_derivatives; row++)
6316
delete [] combinations[row];
6317
delete [] transform[row];
6320
delete [] combinations;
6321
delete [] transform;
6324
/// Evaluate order n derivatives of all basis functions at given point in cell
6325
virtual void evaluate_basis_derivatives_all(unsigned int n,
6327
const double* coordinates,
6328
const ufc::cell& c) const
6330
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
6333
/// Evaluate linear functional for dof i on the function f
6334
virtual double evaluate_dof(unsigned int i,
6335
const ufc::function& f,
6336
const ufc::cell& c) const
6338
// The reference points, direction and weights:
6339
static const double X[1][1][2] = {{{0.333333333, 0.333333333}}};
6340
static const double W[1][1] = {{1}};
6341
static const double D[1][1][1] = {{{1}}};
6343
const double * const * x = c.coordinates;
6344
double result = 0.0;
6345
// Iterate over the points:
6346
// Evaluate basis functions for affine mapping
6347
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
6348
const double w1 = X[i][0][0];
6349
const double w2 = X[i][0][1];
6351
// Compute affine mapping y = F(X)
6353
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
6354
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
6356
// Evaluate function at physical points
6358
f.evaluate(values, y, c);
6360
// Map function values using appropriate mapping
6361
// Affine map: Do nothing
6363
// Note that we do not map the weights (yet).
6365
// Take directional components
6366
for(int k = 0; k < 1; k++)
6367
result += values[k]*D[i][0][k];
6368
// Multiply by weights
6374
/// Evaluate linear functionals for all dofs on the function f
6375
virtual void evaluate_dofs(double* values,
6376
const ufc::function& f,
6377
const ufc::cell& c) const
6379
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
6382
/// Interpolate vertex values from dof values
6383
virtual void interpolate_vertex_values(double* vertex_values,
6384
const double* dof_values,
6385
const ufc::cell& c) const
6387
// Evaluate at vertices and use affine mapping
6388
vertex_values[0] = dof_values[0];
6389
vertex_values[1] = dof_values[0];
6390
vertex_values[2] = dof_values[0];
6393
/// Return the number of sub elements (for a mixed element)
6394
virtual unsigned int num_sub_elements() const
6399
/// Create a new finite element for sub element i (for a mixed element)
6400
virtual ufc::finite_element* create_sub_element(unsigned int i) const
6402
return new mixedpoisson_1_finite_element_1();
6407
/// This class defines the interface for a local-to-global mapping of
6408
/// degrees of freedom (dofs).
6410
class mixedpoisson_1_dof_map_0_0: public ufc::dof_map
6414
unsigned int __global_dimension;
6419
mixedpoisson_1_dof_map_0_0() : ufc::dof_map()
6421
__global_dimension = 0;
6425
virtual ~mixedpoisson_1_dof_map_0_0()
6430
/// Return a string identifying the dof map
6431
virtual const char* signature() const
6433
return "FFC dof map for FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1)";
6436
/// Return true iff mesh entities of topological dimension d are needed
6437
virtual bool needs_mesh_entities(unsigned int d) const
6454
/// Initialize dof map for mesh (return true iff init_cell() is needed)
6455
virtual bool init_mesh(const ufc::mesh& m)
6457
__global_dimension = 2*m.num_entities[1];
6461
/// Initialize dof map for given cell
6462
virtual void init_cell(const ufc::mesh& m,
6468
/// Finish initialization of dof map for cells
6469
virtual void init_cell_finalize()
6474
/// Return the dimension of the global finite element function space
6475
virtual unsigned int global_dimension() const
6477
return __global_dimension;
6480
/// Return the dimension of the local finite element function space for a cell
6481
virtual unsigned int local_dimension(const ufc::cell& c) const
6486
/// Return the maximum dimension of the local finite element function space
6487
virtual unsigned int max_local_dimension() const
6492
// Return the geometric dimension of the coordinates this dof map provides
6493
virtual unsigned int geometric_dimension() const
6498
/// Return the number of dofs on each cell facet
6499
virtual unsigned int num_facet_dofs() const
6504
/// Return the number of dofs associated with each cell entity of dimension d
6505
virtual unsigned int num_entity_dofs(unsigned int d) const
6507
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
6510
/// Tabulate the local-to-global mapping of dofs on a cell
6511
virtual void tabulate_dofs(unsigned int* dofs,
6513
const ufc::cell& c) const
6515
dofs[0] = 2*c.entity_indices[1][0];
6516
dofs[1] = 2*c.entity_indices[1][0] + 1;
6517
dofs[2] = 2*c.entity_indices[1][1];
6518
dofs[3] = 2*c.entity_indices[1][1] + 1;
6519
dofs[4] = 2*c.entity_indices[1][2];
6520
dofs[5] = 2*c.entity_indices[1][2] + 1;
6523
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
6524
virtual void tabulate_facet_dofs(unsigned int* dofs,
6525
unsigned int facet) const
6544
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
6545
virtual void tabulate_entity_dofs(unsigned int* dofs,
6546
unsigned int d, unsigned int i) const
6548
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
6551
/// Tabulate the coordinates of all dofs on a cell
6552
virtual void tabulate_coordinates(double** coordinates,
6553
const ufc::cell& c) const
6555
const double * const * x = c.coordinates;
6556
coordinates[0][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];
6557
coordinates[0][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];
6558
coordinates[1][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];
6559
coordinates[1][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];
6560
coordinates[2][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];
6561
coordinates[2][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];
6562
coordinates[3][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];
6563
coordinates[3][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];
6564
coordinates[4][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];
6565
coordinates[4][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];
6566
coordinates[5][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];
6567
coordinates[5][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];
6570
/// Return the number of sub dof maps (for a mixed element)
6571
virtual unsigned int num_sub_dof_maps() const
6576
/// Create a new dof_map for sub dof map i (for a mixed element)
6577
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
6579
return new mixedpoisson_1_dof_map_0_0();
6584
/// This class defines the interface for a local-to-global mapping of
6585
/// degrees of freedom (dofs).
6587
class mixedpoisson_1_dof_map_0_1: public ufc::dof_map
6591
unsigned int __global_dimension;
6596
mixedpoisson_1_dof_map_0_1() : ufc::dof_map()
6598
__global_dimension = 0;
6602
virtual ~mixedpoisson_1_dof_map_0_1()
6607
/// Return a string identifying the dof map
6608
virtual const char* signature() const
6610
return "FFC dof map for FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
6613
/// Return true iff mesh entities of topological dimension d are needed
6614
virtual bool needs_mesh_entities(unsigned int d) const
6631
/// Initialize dof map for mesh (return true iff init_cell() is needed)
6632
virtual bool init_mesh(const ufc::mesh& m)
6634
__global_dimension = m.num_entities[2];
6638
/// Initialize dof map for given cell
6639
virtual void init_cell(const ufc::mesh& m,
6645
/// Finish initialization of dof map for cells
6646
virtual void init_cell_finalize()
6651
/// Return the dimension of the global finite element function space
6652
virtual unsigned int global_dimension() const
6654
return __global_dimension;
6657
/// Return the dimension of the local finite element function space for a cell
6658
virtual unsigned int local_dimension(const ufc::cell& c) const
6663
/// Return the maximum dimension of the local finite element function space
6664
virtual unsigned int max_local_dimension() const
6669
// Return the geometric dimension of the coordinates this dof map provides
6670
virtual unsigned int geometric_dimension() const
6675
/// Return the number of dofs on each cell facet
6676
virtual unsigned int num_facet_dofs() const
6681
/// Return the number of dofs associated with each cell entity of dimension d
6682
virtual unsigned int num_entity_dofs(unsigned int d) const
6684
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
6687
/// Tabulate the local-to-global mapping of dofs on a cell
6688
virtual void tabulate_dofs(unsigned int* dofs,
6690
const ufc::cell& c) const
6692
dofs[0] = c.entity_indices[2][0];
6695
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
6696
virtual void tabulate_facet_dofs(unsigned int* dofs,
6697
unsigned int facet) const
6713
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
6714
virtual void tabulate_entity_dofs(unsigned int* dofs,
6715
unsigned int d, unsigned int i) const
6717
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
6720
/// Tabulate the coordinates of all dofs on a cell
6721
virtual void tabulate_coordinates(double** coordinates,
6722
const ufc::cell& c) const
6724
const double * const * x = c.coordinates;
6725
coordinates[0][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
6726
coordinates[0][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
6729
/// Return the number of sub dof maps (for a mixed element)
6730
virtual unsigned int num_sub_dof_maps() const
6735
/// Create a new dof_map for sub dof map i (for a mixed element)
6736
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
6738
return new mixedpoisson_1_dof_map_0_1();
6743
/// This class defines the interface for a local-to-global mapping of
6744
/// degrees of freedom (dofs).
6746
class mixedpoisson_1_dof_map_0: public ufc::dof_map
6750
unsigned int __global_dimension;
6755
mixedpoisson_1_dof_map_0() : ufc::dof_map()
6757
__global_dimension = 0;
6761
virtual ~mixedpoisson_1_dof_map_0()
6766
/// Return a string identifying the dof map
6767
virtual const char* signature() const
6769
return "FFC dof map for MixedElement([FiniteElement('Brezzi-Douglas-Marini', 'triangle', 1), FiniteElement('Discontinuous Lagrange', 'triangle', 0)])";
6772
/// Return true iff mesh entities of topological dimension d are needed
6773
virtual bool needs_mesh_entities(unsigned int d) const
6790
/// Initialize dof map for mesh (return true iff init_cell() is needed)
6791
virtual bool init_mesh(const ufc::mesh& m)
6793
__global_dimension = 2*m.num_entities[1] + m.num_entities[2];
6797
/// Initialize dof map for given cell
6798
virtual void init_cell(const ufc::mesh& m,
6804
/// Finish initialization of dof map for cells
6805
virtual void init_cell_finalize()
6810
/// Return the dimension of the global finite element function space
6811
virtual unsigned int global_dimension() const
6813
return __global_dimension;
6816
/// Return the dimension of the local finite element function space for a cell
6817
virtual unsigned int local_dimension(const ufc::cell& c) const
6822
/// Return the maximum dimension of the local finite element function space
6823
virtual unsigned int max_local_dimension() const
6828
// Return the geometric dimension of the coordinates this dof map provides
6829
virtual unsigned int geometric_dimension() const
6834
/// Return the number of dofs on each cell facet
6835
virtual unsigned int num_facet_dofs() const
6840
/// Return the number of dofs associated with each cell entity of dimension d
6841
virtual unsigned int num_entity_dofs(unsigned int d) const
6843
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
6846
/// Tabulate the local-to-global mapping of dofs on a cell
6847
virtual void tabulate_dofs(unsigned int* dofs,
6849
const ufc::cell& c) const
6851
dofs[0] = 2*c.entity_indices[1][0];
6852
dofs[1] = 2*c.entity_indices[1][0] + 1;
6853
dofs[2] = 2*c.entity_indices[1][1];
6854
dofs[3] = 2*c.entity_indices[1][1] + 1;
6855
dofs[4] = 2*c.entity_indices[1][2];
6856
dofs[5] = 2*c.entity_indices[1][2] + 1;
6857
unsigned int offset = 2*m.num_entities[1];
6858
dofs[6] = offset + c.entity_indices[2][0];
6861
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
6862
virtual void tabulate_facet_dofs(unsigned int* dofs,
6863
unsigned int facet) const
6882
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
6883
virtual void tabulate_entity_dofs(unsigned int* dofs,
6884
unsigned int d, unsigned int i) const
6886
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
6889
/// Tabulate the coordinates of all dofs on a cell
6890
virtual void tabulate_coordinates(double** coordinates,
6891
const ufc::cell& c) const
6893
const double * const * x = c.coordinates;
6894
coordinates[0][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];
6895
coordinates[0][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];
6896
coordinates[1][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];
6897
coordinates[1][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];
6898
coordinates[2][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];
6899
coordinates[2][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];
6900
coordinates[3][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];
6901
coordinates[3][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];
6902
coordinates[4][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];
6903
coordinates[4][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];
6904
coordinates[5][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];
6905
coordinates[5][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];
6906
coordinates[6][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
6907
coordinates[6][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
6910
/// Return the number of sub dof maps (for a mixed element)
6911
virtual unsigned int num_sub_dof_maps() const
6916
/// Create a new dof_map for sub dof map i (for a mixed element)
6917
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
6922
return new mixedpoisson_1_dof_map_0_0();
6925
return new mixedpoisson_1_dof_map_0_1();
6933
/// This class defines the interface for a local-to-global mapping of
6934
/// degrees of freedom (dofs).
6936
class mixedpoisson_1_dof_map_1: public ufc::dof_map
6940
unsigned int __global_dimension;
6945
mixedpoisson_1_dof_map_1() : ufc::dof_map()
6947
__global_dimension = 0;
6951
virtual ~mixedpoisson_1_dof_map_1()
6956
/// Return a string identifying the dof map
6957
virtual const char* signature() const
6959
return "FFC dof map for FiniteElement('Discontinuous Lagrange', 'triangle', 0)";
6962
/// Return true iff mesh entities of topological dimension d are needed
6963
virtual bool needs_mesh_entities(unsigned int d) const
6980
/// Initialize dof map for mesh (return true iff init_cell() is needed)
6981
virtual bool init_mesh(const ufc::mesh& m)
6983
__global_dimension = m.num_entities[2];
6987
/// Initialize dof map for given cell
6988
virtual void init_cell(const ufc::mesh& m,
6994
/// Finish initialization of dof map for cells
6995
virtual void init_cell_finalize()
7000
/// Return the dimension of the global finite element function space
7001
virtual unsigned int global_dimension() const
7003
return __global_dimension;
7006
/// Return the dimension of the local finite element function space for a cell
7007
virtual unsigned int local_dimension(const ufc::cell& c) const
7012
/// Return the maximum dimension of the local finite element function space
7013
virtual unsigned int max_local_dimension() const
7018
// Return the geometric dimension of the coordinates this dof map provides
7019
virtual unsigned int geometric_dimension() const
7024
/// Return the number of dofs on each cell facet
7025
virtual unsigned int num_facet_dofs() const
7030
/// Return the number of dofs associated with each cell entity of dimension d
7031
virtual unsigned int num_entity_dofs(unsigned int d) const
7033
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
7036
/// Tabulate the local-to-global mapping of dofs on a cell
7037
virtual void tabulate_dofs(unsigned int* dofs,
7039
const ufc::cell& c) const
7041
dofs[0] = c.entity_indices[2][0];
7044
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
7045
virtual void tabulate_facet_dofs(unsigned int* dofs,
7046
unsigned int facet) const
7062
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
7063
virtual void tabulate_entity_dofs(unsigned int* dofs,
7064
unsigned int d, unsigned int i) const
7066
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
7069
/// Tabulate the coordinates of all dofs on a cell
7070
virtual void tabulate_coordinates(double** coordinates,
7071
const ufc::cell& c) const
7073
const double * const * x = c.coordinates;
7074
coordinates[0][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
7075
coordinates[0][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
7078
/// Return the number of sub dof maps (for a mixed element)
7079
virtual unsigned int num_sub_dof_maps() const
7084
/// Create a new dof_map for sub dof map i (for a mixed element)
7085
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
7087
return new mixedpoisson_1_dof_map_1();
7092
/// This class defines the interface for the tabulation of the cell
7093
/// tensor corresponding to the local contribution to a form from
7094
/// the integral over a cell.
7096
class mixedpoisson_1_cell_integral_0_tensor: public ufc::cell_integral
7101
mixedpoisson_1_cell_integral_0_tensor() : ufc::cell_integral()
7107
virtual ~mixedpoisson_1_cell_integral_0_tensor()
7112
/// Tabulate the tensor for the contribution from a local cell
7113
virtual void tabulate_tensor(double* A,
7114
const double * const * w,
7115
const ufc::cell& c) const
7117
// Number of operations to compute geometry tensor: 1
7118
// Number of operations to compute tensor contraction: 1
7119
// Total number of operations to compute cell tensor: 2
7121
// Extract vertex coordinates
7122
const double * const * x = c.coordinates;
7124
// Compute Jacobian of affine map from reference cell
7125
const double J_00 = x[1][0] - x[0][0];
7126
const double J_01 = x[2][0] - x[0][0];
7127
const double J_10 = x[1][1] - x[0][1];
7128
const double J_11 = x[2][1] - x[0][1];
7130
// Compute determinant of Jacobian
7131
double detJ = J_00*J_11 - J_01*J_10;
7133
// Compute inverse of Jacobian
7136
const double det = std::abs(detJ);
7138
// Compute geometry tensor
7139
const double G0_0 = det*w[0][0];
7141
// Compute element tensor
7153
/// This class defines the interface for the tabulation of the cell
7154
/// tensor corresponding to the local contribution to a form from
7155
/// the integral over a cell.
7157
class mixedpoisson_1_cell_integral_0: public ufc::cell_integral
7161
mixedpoisson_1_cell_integral_0_tensor integral_0_tensor;
7166
mixedpoisson_1_cell_integral_0() : ufc::cell_integral()
7172
virtual ~mixedpoisson_1_cell_integral_0()
7177
/// Tabulate the tensor for the contribution from a local cell
7178
virtual void tabulate_tensor(double* A,
7179
const double * const * w,
7180
const ufc::cell& c) const
7182
// Reset values of the element tensor block
7183
for (unsigned int j = 0; j < 7; j++)
7186
// Add all contributions to element tensor
7187
integral_0_tensor.tabulate_tensor(A, w, c);
7192
/// This class defines the interface for the assembly of the global
7193
/// tensor corresponding to a form with r + n arguments, that is, a
7196
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
7198
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
7199
/// global tensor A is defined by
7201
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
7203
/// where each argument Vj represents the application to the
7204
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
7205
/// fixed functions (coefficients).
7207
class mixedpoisson_form_1: public ufc::form
7212
mixedpoisson_form_1() : ufc::form()
7218
virtual ~mixedpoisson_form_1()
7223
/// Return a string identifying the form
7224
virtual const char* signature() const
7226
return "Form([Integral(Product(Function(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0), 0), Indexed(BasisFunction(MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Cell('triangle', 1, Space(2)), 1), FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)], **{'value_shape': (3,) }), 0), MultiIndex((FixedIndex(2),), {FixedIndex(2): 3}))), Measure('cell', 0, None))])";
7229
/// Return the rank of the global tensor (r)
7230
virtual unsigned int rank() const
7235
/// Return the number of coefficients (n)
7236
virtual unsigned int num_coefficients() const
7241
/// Return the number of cell integrals
7242
virtual unsigned int num_cell_integrals() const
7247
/// Return the number of exterior facet integrals
7248
virtual unsigned int num_exterior_facet_integrals() const
7253
/// Return the number of interior facet integrals
7254
virtual unsigned int num_interior_facet_integrals() const
7259
/// Create a new finite element for argument function i
7260
virtual ufc::finite_element* create_finite_element(unsigned int i) const
7265
return new mixedpoisson_1_finite_element_0();
7268
return new mixedpoisson_1_finite_element_1();
7274
/// Create a new dof map for argument function i
7275
virtual ufc::dof_map* create_dof_map(unsigned int i) const
7280
return new mixedpoisson_1_dof_map_0();
7283
return new mixedpoisson_1_dof_map_1();
7289
/// Create a new cell integral on sub domain i
7290
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
7292
return new mixedpoisson_1_cell_integral_0();
7295
/// Create a new exterior facet integral on sub domain i
7296
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
7301
/// Create a new interior facet integral on sub domain i
7302
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const