← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to bench/fem/assembly/cpp/forms/Poisson2DP2.h

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2009-10-12 14:13:18 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20091012141318-hkbxl0sq555vqv7d
Tags: 0.9.4-1
* New upstream release. This version cleans up the design of the
  function class by adding a new abstraction for user-defined
  functions called Expression. A number of minor bugfixes and
  improvements have also been made.
* debian/watch: Update for new flat directory structure.
* Update debian/copyright.
* debian/rules: Use explicit paths to PETSc 3.0.0 and SLEPc 3.0.0.

Show diffs side-by-side

added added

removed removed

Lines of Context:
bench/fem/assembly/cpp/forms/Poisson2DP2.h
1
1
// This code conforms with the UFC specification version 1.0
2
 
// and was automatically generated by FFC version 0.6.0.
 
2
// and was automatically generated by FFC version 0.7.0.
3
3
//
4
4
// Warning: This code was generated with the option '-l dolfin'
5
5
// and contains DOLFIN-specific wrappers that depend on DOLFIN.
11
11
#include <stdexcept>
12
12
#include <fstream>
13
13
#include <ufc.h>
14
 
 
15
 
/// This class defines the interface for a finite element.
16
 
 
17
 
class UFC_Poisson2DP2BilinearForm_finite_element_0: public ufc::finite_element
18
 
{
19
 
public:
20
 
 
21
 
  /// Constructor
22
 
  UFC_Poisson2DP2BilinearForm_finite_element_0() : ufc::finite_element()
23
 
  {
24
 
    // Do nothing
25
 
  }
26
 
 
27
 
  /// Destructor
28
 
  virtual ~UFC_Poisson2DP2BilinearForm_finite_element_0()
29
 
  {
30
 
    // Do nothing
31
 
  }
32
 
 
33
 
  /// Return a string identifying the finite element
34
 
  virtual const char* signature() const
35
 
  {
36
 
    return "FiniteElement('Lagrange', 'triangle', 2)";
37
 
  }
38
 
 
39
 
  /// Return the cell shape
40
 
  virtual ufc::shape cell_shape() const
41
 
  {
42
 
    return ufc::triangle;
43
 
  }
44
 
 
45
 
  /// Return the dimension of the finite element function space
46
 
  virtual unsigned int space_dimension() const
47
 
  {
48
 
    return 6;
49
 
  }
50
 
 
51
 
  /// Return the rank of the value space
52
 
  virtual unsigned int value_rank() const
53
 
  {
54
 
    return 0;
55
 
  }
56
 
 
57
 
  /// Return the dimension of the value space for axis i
58
 
  virtual unsigned int value_dimension(unsigned int i) const
59
 
  {
60
 
    return 1;
61
 
  }
62
 
 
63
 
  /// Evaluate basis function i at given point in cell
64
 
  virtual void evaluate_basis(unsigned int i,
65
 
                              double* values,
66
 
                              const double* coordinates,
67
 
                              const ufc::cell& c) const
68
 
  {
69
 
    // Extract vertex coordinates
70
 
    const double * const * element_coordinates = c.coordinates;
71
 
    
72
 
    // Compute Jacobian of affine map from reference cell
73
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
74
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
75
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
76
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
77
 
      
78
 
    // Compute determinant of Jacobian
79
 
    const double detJ = J_00*J_11 - J_01*J_10;
80
 
    
81
 
    // Compute inverse of Jacobian
82
 
    
83
 
    // Get coordinates and map to the reference (UFC) element
84
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
85
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
86
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
87
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
88
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
89
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
90
 
    
91
 
    // Map coordinates to the reference square
92
 
    if (std::abs(y - 1.0) < 1e-14)
93
 
      x = -1.0;
94
 
    else
95
 
      x = 2.0 *x/(1.0 - y) - 1.0;
96
 
    y = 2.0*y - 1.0;
97
 
    
98
 
    // Reset values
99
 
    *values = 0;
100
 
    
101
 
    // Map degree of freedom to element degree of freedom
102
 
    const unsigned int dof = i;
103
 
    
104
 
    // Generate scalings
105
 
    const double scalings_y_0 = 1;
106
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
107
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
108
 
    
109
 
    // Compute psitilde_a
110
 
    const double psitilde_a_0 = 1;
111
 
    const double psitilde_a_1 = x;
112
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
113
 
    
114
 
    // Compute psitilde_bs
115
 
    const double psitilde_bs_0_0 = 1;
116
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
117
 
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
118
 
    const double psitilde_bs_1_0 = 1;
119
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
120
 
    const double psitilde_bs_2_0 = 1;
121
 
    
122
 
    // Compute basisvalues
123
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
124
 
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
125
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
126
 
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
127
 
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
128
 
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
129
 
    
130
 
    // Table(s) of coefficients
131
 
    const static double coefficients0[6][6] = \
132
 
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},
133
 
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},
134
 
    {0, 0, 0.2, 0, 0, 0.163299316185545},
135
 
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
136
 
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
137
 
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
138
 
    
139
 
    // Extract relevant coefficients
140
 
    const double coeff0_0 = coefficients0[dof][0];
141
 
    const double coeff0_1 = coefficients0[dof][1];
142
 
    const double coeff0_2 = coefficients0[dof][2];
143
 
    const double coeff0_3 = coefficients0[dof][3];
144
 
    const double coeff0_4 = coefficients0[dof][4];
145
 
    const double coeff0_5 = coefficients0[dof][5];
146
 
    
147
 
    // Compute value(s)
148
 
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
149
 
  }
150
 
 
151
 
  /// Evaluate all basis functions at given point in cell
152
 
  virtual void evaluate_basis_all(double* values,
153
 
                                  const double* coordinates,
154
 
                                  const ufc::cell& c) const
155
 
  {
156
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
157
 
  }
158
 
 
159
 
  /// Evaluate order n derivatives of basis function i at given point in cell
160
 
  virtual void evaluate_basis_derivatives(unsigned int i,
161
 
                                          unsigned int n,
162
 
                                          double* values,
163
 
                                          const double* coordinates,
164
 
                                          const ufc::cell& c) const
165
 
  {
166
 
    // Extract vertex coordinates
167
 
    const double * const * element_coordinates = c.coordinates;
168
 
    
169
 
    // Compute Jacobian of affine map from reference cell
170
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
171
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
172
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
173
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
174
 
      
175
 
    // Compute determinant of Jacobian
176
 
    const double detJ = J_00*J_11 - J_01*J_10;
177
 
    
178
 
    // Compute inverse of Jacobian
179
 
    
180
 
    // Get coordinates and map to the reference (UFC) element
181
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
182
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
183
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
184
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
185
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
186
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
187
 
    
188
 
    // Map coordinates to the reference square
189
 
    if (std::abs(y - 1.0) < 1e-14)
190
 
      x = -1.0;
191
 
    else
192
 
      x = 2.0 *x/(1.0 - y) - 1.0;
193
 
    y = 2.0*y - 1.0;
194
 
    
195
 
    // Compute number of derivatives
196
 
    unsigned int num_derivatives = 1;
197
 
    
198
 
    for (unsigned int j = 0; j < n; j++)
199
 
      num_derivatives *= 2;
200
 
    
201
 
    
202
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
203
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
204
 
        
205
 
    for (unsigned int j = 0; j < num_derivatives; j++)
206
 
    {
207
 
      combinations[j] = new unsigned int [n];
208
 
      for (unsigned int k = 0; k < n; k++)
209
 
        combinations[j][k] = 0;
210
 
    }
211
 
        
212
 
    // Generate combinations of derivatives
213
 
    for (unsigned int row = 1; row < num_derivatives; row++)
214
 
    {
215
 
      for (unsigned int num = 0; num < row; num++)
216
 
      {
217
 
        for (unsigned int col = n-1; col+1 > 0; col--)
218
 
        {
219
 
          if (combinations[row][col] + 1 > 1)
220
 
            combinations[row][col] = 0;
221
 
          else
222
 
          {
223
 
            combinations[row][col] += 1;
224
 
            break;
225
 
          }
226
 
        }
227
 
      }
228
 
    }
229
 
    
230
 
    // Compute inverse of Jacobian
231
 
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
232
 
    
233
 
    // Declare transformation matrix
234
 
    // Declare pointer to two dimensional array and initialise
235
 
    double **transform = new double *[num_derivatives];
236
 
        
237
 
    for (unsigned int j = 0; j < num_derivatives; j++)
238
 
    {
239
 
      transform[j] = new double [num_derivatives];
240
 
      for (unsigned int k = 0; k < num_derivatives; k++)
241
 
        transform[j][k] = 1;
242
 
    }
243
 
    
244
 
    // Construct transformation matrix
245
 
    for (unsigned int row = 0; row < num_derivatives; row++)
246
 
    {
247
 
      for (unsigned int col = 0; col < num_derivatives; col++)
248
 
      {
249
 
        for (unsigned int k = 0; k < n; k++)
250
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
251
 
      }
252
 
    }
253
 
    
254
 
    // Reset values
255
 
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
256
 
      values[j] = 0;
257
 
    
258
 
    // Map degree of freedom to element degree of freedom
259
 
    const unsigned int dof = i;
260
 
    
261
 
    // Generate scalings
262
 
    const double scalings_y_0 = 1;
263
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
264
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
265
 
    
266
 
    // Compute psitilde_a
267
 
    const double psitilde_a_0 = 1;
268
 
    const double psitilde_a_1 = x;
269
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
270
 
    
271
 
    // Compute psitilde_bs
272
 
    const double psitilde_bs_0_0 = 1;
273
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
274
 
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
275
 
    const double psitilde_bs_1_0 = 1;
276
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
277
 
    const double psitilde_bs_2_0 = 1;
278
 
    
279
 
    // Compute basisvalues
280
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
281
 
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
282
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
283
 
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
284
 
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
285
 
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
286
 
    
287
 
    // Table(s) of coefficients
288
 
    const static double coefficients0[6][6] = \
289
 
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},
290
 
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},
291
 
    {0, 0, 0.2, 0, 0, 0.163299316185545},
292
 
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
293
 
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
294
 
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
295
 
    
296
 
    // Interesting (new) part
297
 
    // Tables of derivatives of the polynomial base (transpose)
298
 
    const static double dmats0[6][6] = \
299
 
    {{0, 0, 0, 0, 0, 0},
300
 
    {4.89897948556636, 0, 0, 0, 0, 0},
301
 
    {0, 0, 0, 0, 0, 0},
302
 
    {0, 9.48683298050514, 0, 0, 0, 0},
303
 
    {4, 0, 7.07106781186548, 0, 0, 0},
304
 
    {0, 0, 0, 0, 0, 0}};
305
 
    
306
 
    const static double dmats1[6][6] = \
307
 
    {{0, 0, 0, 0, 0, 0},
308
 
    {2.44948974278318, 0, 0, 0, 0, 0},
309
 
    {4.24264068711928, 0, 0, 0, 0, 0},
310
 
    {2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
311
 
    {2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
312
 
    {-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
313
 
    
314
 
    // Compute reference derivatives
315
 
    // Declare pointer to array of derivatives on FIAT element
316
 
    double *derivatives = new double [num_derivatives];
317
 
    
318
 
    // Declare coefficients
319
 
    double coeff0_0 = 0;
320
 
    double coeff0_1 = 0;
321
 
    double coeff0_2 = 0;
322
 
    double coeff0_3 = 0;
323
 
    double coeff0_4 = 0;
324
 
    double coeff0_5 = 0;
325
 
    
326
 
    // Declare new coefficients
327
 
    double new_coeff0_0 = 0;
328
 
    double new_coeff0_1 = 0;
329
 
    double new_coeff0_2 = 0;
330
 
    double new_coeff0_3 = 0;
331
 
    double new_coeff0_4 = 0;
332
 
    double new_coeff0_5 = 0;
333
 
    
334
 
    // Loop possible derivatives
335
 
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
336
 
    {
337
 
      // Get values from coefficients array
338
 
      new_coeff0_0 = coefficients0[dof][0];
339
 
      new_coeff0_1 = coefficients0[dof][1];
340
 
      new_coeff0_2 = coefficients0[dof][2];
341
 
      new_coeff0_3 = coefficients0[dof][3];
342
 
      new_coeff0_4 = coefficients0[dof][4];
343
 
      new_coeff0_5 = coefficients0[dof][5];
344
 
    
345
 
      // Loop derivative order
346
 
      for (unsigned int j = 0; j < n; j++)
347
 
      {
348
 
        // Update old coefficients
349
 
        coeff0_0 = new_coeff0_0;
350
 
        coeff0_1 = new_coeff0_1;
351
 
        coeff0_2 = new_coeff0_2;
352
 
        coeff0_3 = new_coeff0_3;
353
 
        coeff0_4 = new_coeff0_4;
354
 
        coeff0_5 = new_coeff0_5;
355
 
    
356
 
        if(combinations[deriv_num][j] == 0)
357
 
        {
358
 
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
359
 
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
360
 
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
361
 
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
362
 
          new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
363
 
          new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
364
 
        }
365
 
        if(combinations[deriv_num][j] == 1)
366
 
        {
367
 
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
368
 
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
369
 
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
370
 
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
371
 
          new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
372
 
          new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
373
 
        }
374
 
    
375
 
      }
376
 
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
377
 
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
378
 
    }
379
 
    
380
 
    // Transform derivatives back to physical element
381
 
    for (unsigned int row = 0; row < num_derivatives; row++)
382
 
    {
383
 
      for (unsigned int col = 0; col < num_derivatives; col++)
384
 
      {
385
 
        values[row] += transform[row][col]*derivatives[col];
386
 
      }
387
 
    }
388
 
    // Delete pointer to array of derivatives on FIAT element
389
 
    delete [] derivatives;
390
 
    
391
 
    // Delete pointer to array of combinations of derivatives and transform
392
 
    for (unsigned int row = 0; row < num_derivatives; row++)
393
 
    {
394
 
      delete [] combinations[row];
395
 
      delete [] transform[row];
396
 
    }
397
 
    
398
 
    delete [] combinations;
399
 
    delete [] transform;
400
 
  }
401
 
 
402
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
403
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
404
 
                                              double* values,
405
 
                                              const double* coordinates,
406
 
                                              const ufc::cell& c) const
407
 
  {
408
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
409
 
  }
410
 
 
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
415
 
  {
416
 
    // The reference points, direction and weights:
417
 
    const static double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
418
 
    const static double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
419
 
    const static double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
420
 
    
421
 
    const double * const * x = c.coordinates;
422
 
    double result = 0.0;
423
 
    // Iterate over the points:
424
 
    // Evaluate basis functions for affine mapping
425
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
426
 
    const double w1 = X[i][0][0];
427
 
    const double w2 = X[i][0][1];
428
 
    
429
 
    // Compute affine mapping y = F(X)
430
 
    double y[2];
431
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
432
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
433
 
    
434
 
    // Evaluate function at physical points
435
 
    double values[1];
436
 
    f.evaluate(values, y, c);
437
 
    
438
 
    // Map function values using appropriate mapping
439
 
    // Affine map: Do nothing
440
 
    
441
 
    // Note that we do not map the weights (yet).
442
 
    
443
 
    // Take directional components
444
 
    for(int k = 0; k < 1; k++)
445
 
      result += values[k]*D[i][0][k];
446
 
    // Multiply by weights 
447
 
    result *= W[i][0];
448
 
    
449
 
    return result;
450
 
  }
451
 
 
452
 
  /// Evaluate linear functionals for all dofs on the function f
453
 
  virtual void evaluate_dofs(double* values,
454
 
                             const ufc::function& f,
455
 
                             const ufc::cell& c) const
456
 
  {
457
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
458
 
  }
459
 
 
460
 
  /// Interpolate vertex values from dof values
461
 
  virtual void interpolate_vertex_values(double* vertex_values,
462
 
                                         const double* dof_values,
463
 
                                         const ufc::cell& c) const
464
 
  {
465
 
    // Evaluate at vertices and use affine mapping
466
 
    vertex_values[0] = dof_values[0];
467
 
    vertex_values[1] = dof_values[1];
468
 
    vertex_values[2] = dof_values[2];
469
 
  }
470
 
 
471
 
  /// Return the number of sub elements (for a mixed element)
472
 
  virtual unsigned int num_sub_elements() const
473
 
  {
474
 
    return 1;
475
 
  }
476
 
 
477
 
  /// Create a new finite element for sub element i (for a mixed element)
478
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
479
 
  {
480
 
    return new UFC_Poisson2DP2BilinearForm_finite_element_0();
481
 
  }
482
 
 
483
 
};
484
 
 
485
 
/// This class defines the interface for a finite element.
486
 
 
487
 
class UFC_Poisson2DP2BilinearForm_finite_element_1: public ufc::finite_element
488
 
{
489
 
public:
490
 
 
491
 
  /// Constructor
492
 
  UFC_Poisson2DP2BilinearForm_finite_element_1() : ufc::finite_element()
493
 
  {
494
 
    // Do nothing
495
 
  }
496
 
 
497
 
  /// Destructor
498
 
  virtual ~UFC_Poisson2DP2BilinearForm_finite_element_1()
499
 
  {
500
 
    // Do nothing
501
 
  }
502
 
 
503
 
  /// Return a string identifying the finite element
504
 
  virtual const char* signature() const
505
 
  {
506
 
    return "FiniteElement('Lagrange', 'triangle', 2)";
507
 
  }
508
 
 
509
 
  /// Return the cell shape
510
 
  virtual ufc::shape cell_shape() const
511
 
  {
512
 
    return ufc::triangle;
513
 
  }
514
 
 
515
 
  /// Return the dimension of the finite element function space
516
 
  virtual unsigned int space_dimension() const
517
 
  {
518
 
    return 6;
519
 
  }
520
 
 
521
 
  /// Return the rank of the value space
522
 
  virtual unsigned int value_rank() const
523
 
  {
524
 
    return 0;
525
 
  }
526
 
 
527
 
  /// Return the dimension of the value space for axis i
528
 
  virtual unsigned int value_dimension(unsigned int i) const
529
 
  {
530
 
    return 1;
531
 
  }
532
 
 
533
 
  /// Evaluate basis function i at given point in cell
534
 
  virtual void evaluate_basis(unsigned int i,
535
 
                              double* values,
536
 
                              const double* coordinates,
537
 
                              const ufc::cell& c) const
538
 
  {
539
 
    // Extract vertex coordinates
540
 
    const double * const * element_coordinates = c.coordinates;
541
 
    
542
 
    // Compute Jacobian of affine map from reference cell
543
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
544
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
545
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
546
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
547
 
      
548
 
    // Compute determinant of Jacobian
549
 
    const double detJ = J_00*J_11 - J_01*J_10;
550
 
    
551
 
    // Compute inverse of Jacobian
552
 
    
553
 
    // Get coordinates and map to the reference (UFC) element
554
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
555
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
556
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
557
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
558
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
559
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
560
 
    
561
 
    // Map coordinates to the reference square
562
 
    if (std::abs(y - 1.0) < 1e-14)
563
 
      x = -1.0;
564
 
    else
565
 
      x = 2.0 *x/(1.0 - y) - 1.0;
566
 
    y = 2.0*y - 1.0;
567
 
    
568
 
    // Reset values
569
 
    *values = 0;
570
 
    
571
 
    // Map degree of freedom to element degree of freedom
572
 
    const unsigned int dof = i;
573
 
    
574
 
    // Generate scalings
575
 
    const double scalings_y_0 = 1;
576
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
577
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
578
 
    
579
 
    // Compute psitilde_a
580
 
    const double psitilde_a_0 = 1;
581
 
    const double psitilde_a_1 = x;
582
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
583
 
    
584
 
    // Compute psitilde_bs
585
 
    const double psitilde_bs_0_0 = 1;
586
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
587
 
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
588
 
    const double psitilde_bs_1_0 = 1;
589
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
590
 
    const double psitilde_bs_2_0 = 1;
591
 
    
592
 
    // Compute basisvalues
593
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
594
 
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
595
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
596
 
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
597
 
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
598
 
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
599
 
    
600
 
    // Table(s) of coefficients
601
 
    const static double coefficients0[6][6] = \
602
 
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},
603
 
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},
604
 
    {0, 0, 0.2, 0, 0, 0.163299316185545},
605
 
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
606
 
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
607
 
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
608
 
    
609
 
    // Extract relevant coefficients
610
 
    const double coeff0_0 = coefficients0[dof][0];
611
 
    const double coeff0_1 = coefficients0[dof][1];
612
 
    const double coeff0_2 = coefficients0[dof][2];
613
 
    const double coeff0_3 = coefficients0[dof][3];
614
 
    const double coeff0_4 = coefficients0[dof][4];
615
 
    const double coeff0_5 = coefficients0[dof][5];
616
 
    
617
 
    // Compute value(s)
618
 
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
619
 
  }
620
 
 
621
 
  /// Evaluate all basis functions at given point in cell
622
 
  virtual void evaluate_basis_all(double* values,
623
 
                                  const double* coordinates,
624
 
                                  const ufc::cell& c) const
625
 
  {
626
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
627
 
  }
628
 
 
629
 
  /// Evaluate order n derivatives of basis function i at given point in cell
630
 
  virtual void evaluate_basis_derivatives(unsigned int i,
631
 
                                          unsigned int n,
632
 
                                          double* values,
633
 
                                          const double* coordinates,
634
 
                                          const ufc::cell& c) const
635
 
  {
636
 
    // Extract vertex coordinates
637
 
    const double * const * element_coordinates = c.coordinates;
638
 
    
639
 
    // Compute Jacobian of affine map from reference cell
640
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
641
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
642
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
643
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
644
 
      
645
 
    // Compute determinant of Jacobian
646
 
    const double detJ = J_00*J_11 - J_01*J_10;
647
 
    
648
 
    // Compute inverse of Jacobian
649
 
    
650
 
    // Get coordinates and map to the reference (UFC) element
651
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
652
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
653
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
654
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
655
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
656
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
657
 
    
658
 
    // Map coordinates to the reference square
659
 
    if (std::abs(y - 1.0) < 1e-14)
660
 
      x = -1.0;
661
 
    else
662
 
      x = 2.0 *x/(1.0 - y) - 1.0;
663
 
    y = 2.0*y - 1.0;
664
 
    
665
 
    // Compute number of derivatives
666
 
    unsigned int num_derivatives = 1;
667
 
    
668
 
    for (unsigned int j = 0; j < n; j++)
669
 
      num_derivatives *= 2;
670
 
    
671
 
    
672
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
673
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
674
 
        
675
 
    for (unsigned int j = 0; j < num_derivatives; j++)
676
 
    {
677
 
      combinations[j] = new unsigned int [n];
678
 
      for (unsigned int k = 0; k < n; k++)
679
 
        combinations[j][k] = 0;
680
 
    }
681
 
        
682
 
    // Generate combinations of derivatives
683
 
    for (unsigned int row = 1; row < num_derivatives; row++)
684
 
    {
685
 
      for (unsigned int num = 0; num < row; num++)
686
 
      {
687
 
        for (unsigned int col = n-1; col+1 > 0; col--)
688
 
        {
689
 
          if (combinations[row][col] + 1 > 1)
690
 
            combinations[row][col] = 0;
691
 
          else
692
 
          {
693
 
            combinations[row][col] += 1;
694
 
            break;
695
 
          }
696
 
        }
697
 
      }
698
 
    }
699
 
    
700
 
    // Compute inverse of Jacobian
701
 
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
702
 
    
703
 
    // Declare transformation matrix
704
 
    // Declare pointer to two dimensional array and initialise
705
 
    double **transform = new double *[num_derivatives];
706
 
        
707
 
    for (unsigned int j = 0; j < num_derivatives; j++)
708
 
    {
709
 
      transform[j] = new double [num_derivatives];
710
 
      for (unsigned int k = 0; k < num_derivatives; k++)
711
 
        transform[j][k] = 1;
712
 
    }
713
 
    
714
 
    // Construct transformation matrix
715
 
    for (unsigned int row = 0; row < num_derivatives; row++)
716
 
    {
717
 
      for (unsigned int col = 0; col < num_derivatives; col++)
718
 
      {
719
 
        for (unsigned int k = 0; k < n; k++)
720
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
721
 
      }
722
 
    }
723
 
    
724
 
    // Reset values
725
 
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
726
 
      values[j] = 0;
727
 
    
728
 
    // Map degree of freedom to element degree of freedom
729
 
    const unsigned int dof = i;
730
 
    
731
 
    // Generate scalings
732
 
    const double scalings_y_0 = 1;
733
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
734
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
735
 
    
736
 
    // Compute psitilde_a
737
 
    const double psitilde_a_0 = 1;
738
 
    const double psitilde_a_1 = x;
739
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
740
 
    
741
 
    // Compute psitilde_bs
742
 
    const double psitilde_bs_0_0 = 1;
743
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
744
 
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
745
 
    const double psitilde_bs_1_0 = 1;
746
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
747
 
    const double psitilde_bs_2_0 = 1;
748
 
    
749
 
    // Compute basisvalues
750
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
751
 
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
752
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
753
 
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
754
 
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
755
 
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
756
 
    
757
 
    // Table(s) of coefficients
758
 
    const static double coefficients0[6][6] = \
759
 
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},
760
 
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},
761
 
    {0, 0, 0.2, 0, 0, 0.163299316185545},
762
 
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
763
 
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
764
 
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
765
 
    
766
 
    // Interesting (new) part
767
 
    // Tables of derivatives of the polynomial base (transpose)
768
 
    const static double dmats0[6][6] = \
769
 
    {{0, 0, 0, 0, 0, 0},
770
 
    {4.89897948556636, 0, 0, 0, 0, 0},
771
 
    {0, 0, 0, 0, 0, 0},
772
 
    {0, 9.48683298050514, 0, 0, 0, 0},
773
 
    {4, 0, 7.07106781186548, 0, 0, 0},
774
 
    {0, 0, 0, 0, 0, 0}};
775
 
    
776
 
    const static double dmats1[6][6] = \
777
 
    {{0, 0, 0, 0, 0, 0},
778
 
    {2.44948974278318, 0, 0, 0, 0, 0},
779
 
    {4.24264068711928, 0, 0, 0, 0, 0},
780
 
    {2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
781
 
    {2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
782
 
    {-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
783
 
    
784
 
    // Compute reference derivatives
785
 
    // Declare pointer to array of derivatives on FIAT element
786
 
    double *derivatives = new double [num_derivatives];
787
 
    
788
 
    // Declare coefficients
789
 
    double coeff0_0 = 0;
790
 
    double coeff0_1 = 0;
791
 
    double coeff0_2 = 0;
792
 
    double coeff0_3 = 0;
793
 
    double coeff0_4 = 0;
794
 
    double coeff0_5 = 0;
795
 
    
796
 
    // Declare new coefficients
797
 
    double new_coeff0_0 = 0;
798
 
    double new_coeff0_1 = 0;
799
 
    double new_coeff0_2 = 0;
800
 
    double new_coeff0_3 = 0;
801
 
    double new_coeff0_4 = 0;
802
 
    double new_coeff0_5 = 0;
803
 
    
804
 
    // Loop possible derivatives
805
 
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
806
 
    {
807
 
      // Get values from coefficients array
808
 
      new_coeff0_0 = coefficients0[dof][0];
809
 
      new_coeff0_1 = coefficients0[dof][1];
810
 
      new_coeff0_2 = coefficients0[dof][2];
811
 
      new_coeff0_3 = coefficients0[dof][3];
812
 
      new_coeff0_4 = coefficients0[dof][4];
813
 
      new_coeff0_5 = coefficients0[dof][5];
814
 
    
815
 
      // Loop derivative order
816
 
      for (unsigned int j = 0; j < n; j++)
817
 
      {
818
 
        // Update old coefficients
819
 
        coeff0_0 = new_coeff0_0;
820
 
        coeff0_1 = new_coeff0_1;
821
 
        coeff0_2 = new_coeff0_2;
822
 
        coeff0_3 = new_coeff0_3;
823
 
        coeff0_4 = new_coeff0_4;
824
 
        coeff0_5 = new_coeff0_5;
825
 
    
826
 
        if(combinations[deriv_num][j] == 0)
827
 
        {
828
 
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
829
 
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
830
 
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
831
 
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
832
 
          new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
833
 
          new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
834
 
        }
835
 
        if(combinations[deriv_num][j] == 1)
836
 
        {
837
 
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
838
 
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
839
 
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
840
 
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
841
 
          new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
842
 
          new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
843
 
        }
844
 
    
845
 
      }
846
 
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
847
 
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
848
 
    }
849
 
    
850
 
    // Transform derivatives back to physical element
851
 
    for (unsigned int row = 0; row < num_derivatives; row++)
852
 
    {
853
 
      for (unsigned int col = 0; col < num_derivatives; col++)
854
 
      {
855
 
        values[row] += transform[row][col]*derivatives[col];
856
 
      }
857
 
    }
858
 
    // Delete pointer to array of derivatives on FIAT element
859
 
    delete [] derivatives;
860
 
    
861
 
    // Delete pointer to array of combinations of derivatives and transform
862
 
    for (unsigned int row = 0; row < num_derivatives; row++)
863
 
    {
864
 
      delete [] combinations[row];
865
 
      delete [] transform[row];
866
 
    }
867
 
    
868
 
    delete [] combinations;
869
 
    delete [] transform;
870
 
  }
871
 
 
872
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
873
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
874
 
                                              double* values,
875
 
                                              const double* coordinates,
876
 
                                              const ufc::cell& c) const
877
 
  {
878
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
879
 
  }
880
 
 
881
 
  /// Evaluate linear functional for dof i on the function f
882
 
  virtual double evaluate_dof(unsigned int i,
883
 
                              const ufc::function& f,
884
 
                              const ufc::cell& c) const
885
 
  {
886
 
    // The reference points, direction and weights:
887
 
    const static double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
888
 
    const static double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
889
 
    const static double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
890
 
    
891
 
    const double * const * x = c.coordinates;
892
 
    double result = 0.0;
893
 
    // Iterate over the points:
894
 
    // Evaluate basis functions for affine mapping
895
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
896
 
    const double w1 = X[i][0][0];
897
 
    const double w2 = X[i][0][1];
898
 
    
899
 
    // Compute affine mapping y = F(X)
900
 
    double y[2];
901
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
902
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
903
 
    
904
 
    // Evaluate function at physical points
905
 
    double values[1];
906
 
    f.evaluate(values, y, c);
907
 
    
908
 
    // Map function values using appropriate mapping
909
 
    // Affine map: Do nothing
910
 
    
911
 
    // Note that we do not map the weights (yet).
912
 
    
913
 
    // Take directional components
914
 
    for(int k = 0; k < 1; k++)
915
 
      result += values[k]*D[i][0][k];
916
 
    // Multiply by weights 
917
 
    result *= W[i][0];
918
 
    
919
 
    return result;
920
 
  }
921
 
 
922
 
  /// Evaluate linear functionals for all dofs on the function f
923
 
  virtual void evaluate_dofs(double* values,
924
 
                             const ufc::function& f,
925
 
                             const ufc::cell& c) const
926
 
  {
927
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
928
 
  }
929
 
 
930
 
  /// Interpolate vertex values from dof values
931
 
  virtual void interpolate_vertex_values(double* vertex_values,
932
 
                                         const double* dof_values,
933
 
                                         const ufc::cell& c) const
934
 
  {
935
 
    // Evaluate at vertices and use affine mapping
936
 
    vertex_values[0] = dof_values[0];
937
 
    vertex_values[1] = dof_values[1];
938
 
    vertex_values[2] = dof_values[2];
939
 
  }
940
 
 
941
 
  /// Return the number of sub elements (for a mixed element)
942
 
  virtual unsigned int num_sub_elements() const
943
 
  {
944
 
    return 1;
945
 
  }
946
 
 
947
 
  /// Create a new finite element for sub element i (for a mixed element)
948
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
949
 
  {
950
 
    return new UFC_Poisson2DP2BilinearForm_finite_element_1();
951
 
  }
952
 
 
953
 
};
954
 
 
955
 
/// This class defines the interface for a local-to-global mapping of
956
 
/// degrees of freedom (dofs).
957
 
 
958
 
class UFC_Poisson2DP2BilinearForm_dof_map_0: public ufc::dof_map
959
 
{
960
 
private:
961
 
 
962
 
  unsigned int __global_dimension;
963
 
 
964
 
public:
965
 
 
966
 
  /// Constructor
967
 
  UFC_Poisson2DP2BilinearForm_dof_map_0() : ufc::dof_map()
968
 
  {
969
 
    __global_dimension = 0;
970
 
  }
971
 
 
972
 
  /// Destructor
973
 
  virtual ~UFC_Poisson2DP2BilinearForm_dof_map_0()
974
 
  {
975
 
    // Do nothing
976
 
  }
977
 
 
978
 
  /// Return a string identifying the dof map
979
 
  virtual const char* signature() const
980
 
  {
981
 
    return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
982
 
  }
983
 
 
984
 
  /// Return true iff mesh entities of topological dimension d are needed
985
 
  virtual bool needs_mesh_entities(unsigned int d) const
986
 
  {
987
 
    switch ( d )
988
 
    {
989
 
    case 0:
990
 
      return true;
991
 
      break;
992
 
    case 1:
993
 
      return true;
994
 
      break;
995
 
    case 2:
996
 
      return false;
997
 
      break;
998
 
    }
999
 
    return false;
1000
 
  }
1001
 
 
1002
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
1003
 
  virtual bool init_mesh(const ufc::mesh& m)
1004
 
  {
1005
 
    __global_dimension = m.num_entities[0] + m.num_entities[1];
1006
 
    return false;
1007
 
  }
1008
 
 
1009
 
  /// Initialize dof map for given cell
1010
 
  virtual void init_cell(const ufc::mesh& m,
1011
 
                         const ufc::cell& c)
1012
 
  {
1013
 
    // Do nothing
1014
 
  }
1015
 
 
1016
 
  /// Finish initialization of dof map for cells
1017
 
  virtual void init_cell_finalize()
1018
 
  {
1019
 
    // Do nothing
1020
 
  }
1021
 
 
1022
 
  /// Return the dimension of the global finite element function space
1023
 
  virtual unsigned int global_dimension() const
1024
 
  {
1025
 
    return __global_dimension;
1026
 
  }
1027
 
 
1028
 
  /// Return the dimension of the local finite element function space
1029
 
  virtual unsigned int local_dimension() const
1030
 
  {
1031
 
    return 6;
1032
 
  }
1033
 
 
1034
 
  // Return the geometric dimension of the coordinates this dof map provides
1035
 
  virtual unsigned int geometric_dimension() const
1036
 
  {
1037
 
    return 2;
1038
 
  }
1039
 
 
1040
 
  /// Return the number of dofs on each cell facet
1041
 
  virtual unsigned int num_facet_dofs() const
1042
 
  {
1043
 
    return 3;
1044
 
  }
1045
 
 
1046
 
  /// Return the number of dofs associated with each cell entity of dimension d
1047
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
1048
 
  {
1049
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1050
 
  }
1051
 
 
1052
 
  /// Tabulate the local-to-global mapping of dofs on a cell
1053
 
  virtual void tabulate_dofs(unsigned int* dofs,
1054
 
                             const ufc::mesh& m,
1055
 
                             const ufc::cell& c) const
1056
 
  {
1057
 
    dofs[0] = c.entity_indices[0][0];
1058
 
    dofs[1] = c.entity_indices[0][1];
1059
 
    dofs[2] = c.entity_indices[0][2];
1060
 
    unsigned int offset = m.num_entities[0];
1061
 
    dofs[3] = offset + c.entity_indices[1][0];
1062
 
    dofs[4] = offset + c.entity_indices[1][1];
1063
 
    dofs[5] = offset + c.entity_indices[1][2];
1064
 
  }
1065
 
 
1066
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
1067
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
1068
 
                                   unsigned int facet) const
1069
 
  {
1070
 
    switch ( facet )
1071
 
    {
1072
 
    case 0:
1073
 
      dofs[0] = 1;
1074
 
      dofs[1] = 2;
1075
 
      dofs[2] = 3;
1076
 
      break;
1077
 
    case 1:
1078
 
      dofs[0] = 0;
1079
 
      dofs[1] = 2;
1080
 
      dofs[2] = 4;
1081
 
      break;
1082
 
    case 2:
1083
 
      dofs[0] = 0;
1084
 
      dofs[1] = 1;
1085
 
      dofs[2] = 5;
1086
 
      break;
1087
 
    }
1088
 
  }
1089
 
 
1090
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
1091
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
1092
 
                                    unsigned int d, unsigned int i) const
1093
 
  {
1094
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1095
 
  }
1096
 
 
1097
 
  /// Tabulate the coordinates of all dofs on a cell
1098
 
  virtual void tabulate_coordinates(double** coordinates,
1099
 
                                    const ufc::cell& c) const
1100
 
  {
1101
 
    const double * const * x = c.coordinates;
1102
 
    coordinates[0][0] = x[0][0];
1103
 
    coordinates[0][1] = x[0][1];
1104
 
    coordinates[1][0] = x[1][0];
1105
 
    coordinates[1][1] = x[1][1];
1106
 
    coordinates[2][0] = x[2][0];
1107
 
    coordinates[2][1] = x[2][1];
1108
 
    coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1109
 
    coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1110
 
    coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1111
 
    coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1112
 
    coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1113
 
    coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
1114
 
  }
1115
 
 
1116
 
  /// Return the number of sub dof maps (for a mixed element)
1117
 
  virtual unsigned int num_sub_dof_maps() const
1118
 
  {
1119
 
    return 1;
1120
 
  }
1121
 
 
1122
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
1123
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1124
 
  {
1125
 
    return new UFC_Poisson2DP2BilinearForm_dof_map_0();
1126
 
  }
1127
 
 
1128
 
};
1129
 
 
1130
 
/// This class defines the interface for a local-to-global mapping of
1131
 
/// degrees of freedom (dofs).
1132
 
 
1133
 
class UFC_Poisson2DP2BilinearForm_dof_map_1: public ufc::dof_map
1134
 
{
1135
 
private:
1136
 
 
1137
 
  unsigned int __global_dimension;
1138
 
 
1139
 
public:
1140
 
 
1141
 
  /// Constructor
1142
 
  UFC_Poisson2DP2BilinearForm_dof_map_1() : ufc::dof_map()
1143
 
  {
1144
 
    __global_dimension = 0;
1145
 
  }
1146
 
 
1147
 
  /// Destructor
1148
 
  virtual ~UFC_Poisson2DP2BilinearForm_dof_map_1()
1149
 
  {
1150
 
    // Do nothing
1151
 
  }
1152
 
 
1153
 
  /// Return a string identifying the dof map
1154
 
  virtual const char* signature() const
1155
 
  {
1156
 
    return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
1157
 
  }
1158
 
 
1159
 
  /// Return true iff mesh entities of topological dimension d are needed
1160
 
  virtual bool needs_mesh_entities(unsigned int d) const
1161
 
  {
1162
 
    switch ( d )
1163
 
    {
1164
 
    case 0:
1165
 
      return true;
1166
 
      break;
1167
 
    case 1:
1168
 
      return true;
1169
 
      break;
1170
 
    case 2:
1171
 
      return false;
1172
 
      break;
1173
 
    }
1174
 
    return false;
1175
 
  }
1176
 
 
1177
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
1178
 
  virtual bool init_mesh(const ufc::mesh& m)
1179
 
  {
1180
 
    __global_dimension = m.num_entities[0] + m.num_entities[1];
1181
 
    return false;
1182
 
  }
1183
 
 
1184
 
  /// Initialize dof map for given cell
1185
 
  virtual void init_cell(const ufc::mesh& m,
1186
 
                         const ufc::cell& c)
1187
 
  {
1188
 
    // Do nothing
1189
 
  }
1190
 
 
1191
 
  /// Finish initialization of dof map for cells
1192
 
  virtual void init_cell_finalize()
1193
 
  {
1194
 
    // Do nothing
1195
 
  }
1196
 
 
1197
 
  /// Return the dimension of the global finite element function space
1198
 
  virtual unsigned int global_dimension() const
1199
 
  {
1200
 
    return __global_dimension;
1201
 
  }
1202
 
 
1203
 
  /// Return the dimension of the local finite element function space
1204
 
  virtual unsigned int local_dimension() const
1205
 
  {
1206
 
    return 6;
1207
 
  }
1208
 
 
1209
 
  // Return the geometric dimension of the coordinates this dof map provides
1210
 
  virtual unsigned int geometric_dimension() const
1211
 
  {
1212
 
    return 2;
1213
 
  }
1214
 
 
1215
 
  /// Return the number of dofs on each cell facet
1216
 
  virtual unsigned int num_facet_dofs() const
1217
 
  {
1218
 
    return 3;
1219
 
  }
1220
 
 
1221
 
  /// Return the number of dofs associated with each cell entity of dimension d
1222
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
1223
 
  {
1224
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1225
 
  }
1226
 
 
1227
 
  /// Tabulate the local-to-global mapping of dofs on a cell
1228
 
  virtual void tabulate_dofs(unsigned int* dofs,
1229
 
                             const ufc::mesh& m,
1230
 
                             const ufc::cell& c) const
1231
 
  {
1232
 
    dofs[0] = c.entity_indices[0][0];
1233
 
    dofs[1] = c.entity_indices[0][1];
1234
 
    dofs[2] = c.entity_indices[0][2];
1235
 
    unsigned int offset = m.num_entities[0];
1236
 
    dofs[3] = offset + c.entity_indices[1][0];
1237
 
    dofs[4] = offset + c.entity_indices[1][1];
1238
 
    dofs[5] = offset + c.entity_indices[1][2];
1239
 
  }
1240
 
 
1241
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
1242
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
1243
 
                                   unsigned int facet) const
1244
 
  {
1245
 
    switch ( facet )
1246
 
    {
1247
 
    case 0:
1248
 
      dofs[0] = 1;
1249
 
      dofs[1] = 2;
1250
 
      dofs[2] = 3;
1251
 
      break;
1252
 
    case 1:
1253
 
      dofs[0] = 0;
1254
 
      dofs[1] = 2;
1255
 
      dofs[2] = 4;
1256
 
      break;
1257
 
    case 2:
1258
 
      dofs[0] = 0;
1259
 
      dofs[1] = 1;
1260
 
      dofs[2] = 5;
1261
 
      break;
1262
 
    }
1263
 
  }
1264
 
 
1265
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
1266
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
1267
 
                                    unsigned int d, unsigned int i) const
1268
 
  {
1269
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1270
 
  }
1271
 
 
1272
 
  /// Tabulate the coordinates of all dofs on a cell
1273
 
  virtual void tabulate_coordinates(double** coordinates,
1274
 
                                    const ufc::cell& c) const
1275
 
  {
1276
 
    const double * const * x = c.coordinates;
1277
 
    coordinates[0][0] = x[0][0];
1278
 
    coordinates[0][1] = x[0][1];
1279
 
    coordinates[1][0] = x[1][0];
1280
 
    coordinates[1][1] = x[1][1];
1281
 
    coordinates[2][0] = x[2][0];
1282
 
    coordinates[2][1] = x[2][1];
1283
 
    coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1284
 
    coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1285
 
    coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1286
 
    coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1287
 
    coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1288
 
    coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
1289
 
  }
1290
 
 
1291
 
  /// Return the number of sub dof maps (for a mixed element)
1292
 
  virtual unsigned int num_sub_dof_maps() const
1293
 
  {
1294
 
    return 1;
1295
 
  }
1296
 
 
1297
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
1298
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1299
 
  {
1300
 
    return new UFC_Poisson2DP2BilinearForm_dof_map_1();
 
14
    
 
15
/// This class defines the interface for a finite element.
 
16
 
 
17
class poisson2dp2_0_finite_element_0: public ufc::finite_element
 
18
{
 
19
public:
 
20
 
 
21
  /// Constructor
 
22
  poisson2dp2_0_finite_element_0() : ufc::finite_element()
 
23
  {
 
24
    // Do nothing
 
25
  }
 
26
 
 
27
  /// Destructor
 
28
  virtual ~poisson2dp2_0_finite_element_0()
 
29
  {
 
30
    // Do nothing
 
31
  }
 
32
 
 
33
  /// Return a string identifying the finite element
 
34
  virtual const char* signature() const
 
35
  {
 
36
    return "FiniteElement('Lagrange', 'triangle', 2)";
 
37
  }
 
38
 
 
39
  /// Return the cell shape
 
40
  virtual ufc::shape cell_shape() const
 
41
  {
 
42
    return ufc::triangle;
 
43
  }
 
44
 
 
45
  /// Return the dimension of the finite element function space
 
46
  virtual unsigned int space_dimension() const
 
47
  {
 
48
    return 6;
 
49
  }
 
50
 
 
51
  /// Return the rank of the value space
 
52
  virtual unsigned int value_rank() const
 
53
  {
 
54
    return 0;
 
55
  }
 
56
 
 
57
  /// Return the dimension of the value space for axis i
 
58
  virtual unsigned int value_dimension(unsigned int i) const
 
59
  {
 
60
    return 1;
 
61
  }
 
62
 
 
63
  /// Evaluate basis function i at given point in cell
 
64
  virtual void evaluate_basis(unsigned int i,
 
65
                              double* values,
 
66
                              const double* coordinates,
 
67
                              const ufc::cell& c) const
 
68
  {
 
69
    // Extract vertex coordinates
 
70
    const double * const * element_coordinates = c.coordinates;
 
71
    
 
72
    // Compute Jacobian of affine map from reference cell
 
73
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
74
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
75
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
76
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
77
    
 
78
    // Compute determinant of Jacobian
 
79
    const double detJ = J_00*J_11 - J_01*J_10;
 
80
    
 
81
    // Compute inverse of Jacobian
 
82
    
 
83
    // Get coordinates and map to the reference (UFC) element
 
84
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
85
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
86
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
87
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
88
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
89
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
90
    
 
91
    // Map coordinates to the reference square
 
92
    if (std::abs(y - 1.0) < 1e-14)
 
93
      x = -1.0;
 
94
    else
 
95
      x = 2.0 *x/(1.0 - y) - 1.0;
 
96
    y = 2.0*y - 1.0;
 
97
    
 
98
    // Reset values
 
99
    *values = 0;
 
100
    
 
101
    // Map degree of freedom to element degree of freedom
 
102
    const unsigned int dof = i;
 
103
    
 
104
    // Generate scalings
 
105
    const double scalings_y_0 = 1;
 
106
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
107
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
 
108
    
 
109
    // Compute psitilde_a
 
110
    const double psitilde_a_0 = 1;
 
111
    const double psitilde_a_1 = x;
 
112
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
 
113
    
 
114
    // Compute psitilde_bs
 
115
    const double psitilde_bs_0_0 = 1;
 
116
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
117
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
 
118
    const double psitilde_bs_1_0 = 1;
 
119
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
 
120
    const double psitilde_bs_2_0 = 1;
 
121
    
 
122
    // Compute basisvalues
 
123
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
124
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
125
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
126
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
 
127
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
 
128
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
 
129
    
 
130
    // Table(s) of coefficients
 
131
    static const double coefficients0[6][6] = \
 
132
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
 
133
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
 
134
    {0, 0, 0.2, 0, 0, 0.163299316185545},
 
135
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
 
136
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
 
137
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
 
138
    
 
139
    // Extract relevant coefficients
 
140
    const double coeff0_0 = coefficients0[dof][0];
 
141
    const double coeff0_1 = coefficients0[dof][1];
 
142
    const double coeff0_2 = coefficients0[dof][2];
 
143
    const double coeff0_3 = coefficients0[dof][3];
 
144
    const double coeff0_4 = coefficients0[dof][4];
 
145
    const double coeff0_5 = coefficients0[dof][5];
 
146
    
 
147
    // Compute value(s)
 
148
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
 
149
  }
 
150
 
 
151
  /// Evaluate all basis functions at given point in cell
 
152
  virtual void evaluate_basis_all(double* values,
 
153
                                  const double* coordinates,
 
154
                                  const ufc::cell& c) const
 
155
  {
 
156
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
157
  }
 
158
 
 
159
  /// Evaluate order n derivatives of basis function i at given point in cell
 
160
  virtual void evaluate_basis_derivatives(unsigned int i,
 
161
                                          unsigned int n,
 
162
                                          double* values,
 
163
                                          const double* coordinates,
 
164
                                          const ufc::cell& c) const
 
165
  {
 
166
    // Extract vertex coordinates
 
167
    const double * const * element_coordinates = c.coordinates;
 
168
    
 
169
    // Compute Jacobian of affine map from reference cell
 
170
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
171
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
172
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
173
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
174
    
 
175
    // Compute determinant of Jacobian
 
176
    const double detJ = J_00*J_11 - J_01*J_10;
 
177
    
 
178
    // Compute inverse of Jacobian
 
179
    
 
180
    // Get coordinates and map to the reference (UFC) element
 
181
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
182
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
183
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
184
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
185
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
186
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
187
    
 
188
    // Map coordinates to the reference square
 
189
    if (std::abs(y - 1.0) < 1e-14)
 
190
      x = -1.0;
 
191
    else
 
192
      x = 2.0 *x/(1.0 - y) - 1.0;
 
193
    y = 2.0*y - 1.0;
 
194
    
 
195
    // Compute number of derivatives
 
196
    unsigned int num_derivatives = 1;
 
197
    
 
198
    for (unsigned int j = 0; j < n; j++)
 
199
      num_derivatives *= 2;
 
200
    
 
201
    
 
202
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
203
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
204
    
 
205
    for (unsigned int j = 0; j < num_derivatives; j++)
 
206
    {
 
207
      combinations[j] = new unsigned int [n];
 
208
      for (unsigned int k = 0; k < n; k++)
 
209
        combinations[j][k] = 0;
 
210
    }
 
211
    
 
212
    // Generate combinations of derivatives
 
213
    for (unsigned int row = 1; row < num_derivatives; row++)
 
214
    {
 
215
      for (unsigned int num = 0; num < row; num++)
 
216
      {
 
217
        for (unsigned int col = n-1; col+1 > 0; col--)
 
218
        {
 
219
          if (combinations[row][col] + 1 > 1)
 
220
            combinations[row][col] = 0;
 
221
          else
 
222
          {
 
223
            combinations[row][col] += 1;
 
224
            break;
 
225
          }
 
226
        }
 
227
      }
 
228
    }
 
229
    
 
230
    // Compute inverse of Jacobian
 
231
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
232
    
 
233
    // Declare transformation matrix
 
234
    // Declare pointer to two dimensional array and initialise
 
235
    double **transform = new double *[num_derivatives];
 
236
    
 
237
    for (unsigned int j = 0; j < num_derivatives; j++)
 
238
    {
 
239
      transform[j] = new double [num_derivatives];
 
240
      for (unsigned int k = 0; k < num_derivatives; k++)
 
241
        transform[j][k] = 1;
 
242
    }
 
243
    
 
244
    // Construct transformation matrix
 
245
    for (unsigned int row = 0; row < num_derivatives; row++)
 
246
    {
 
247
      for (unsigned int col = 0; col < num_derivatives; col++)
 
248
      {
 
249
        for (unsigned int k = 0; k < n; k++)
 
250
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
251
      }
 
252
    }
 
253
    
 
254
    // Reset values
 
255
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
256
      values[j] = 0;
 
257
    
 
258
    // Map degree of freedom to element degree of freedom
 
259
    const unsigned int dof = i;
 
260
    
 
261
    // Generate scalings
 
262
    const double scalings_y_0 = 1;
 
263
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
264
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
 
265
    
 
266
    // Compute psitilde_a
 
267
    const double psitilde_a_0 = 1;
 
268
    const double psitilde_a_1 = x;
 
269
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
 
270
    
 
271
    // Compute psitilde_bs
 
272
    const double psitilde_bs_0_0 = 1;
 
273
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
274
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
 
275
    const double psitilde_bs_1_0 = 1;
 
276
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
 
277
    const double psitilde_bs_2_0 = 1;
 
278
    
 
279
    // Compute basisvalues
 
280
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
281
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
282
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
283
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
 
284
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
 
285
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
 
286
    
 
287
    // Table(s) of coefficients
 
288
    static const double coefficients0[6][6] = \
 
289
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
 
290
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
 
291
    {0, 0, 0.2, 0, 0, 0.163299316185545},
 
292
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
 
293
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
 
294
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
 
295
    
 
296
    // Interesting (new) part
 
297
    // Tables of derivatives of the polynomial base (transpose)
 
298
    static const double dmats0[6][6] = \
 
299
    {{0, 0, 0, 0, 0, 0},
 
300
    {4.89897948556635, 0, 0, 0, 0, 0},
 
301
    {0, 0, 0, 0, 0, 0},
 
302
    {0, 9.48683298050514, 0, 0, 0, 0},
 
303
    {4, 0, 7.07106781186548, 0, 0, 0},
 
304
    {0, 0, 0, 0, 0, 0}};
 
305
    
 
306
    static const double dmats1[6][6] = \
 
307
    {{0, 0, 0, 0, 0, 0},
 
308
    {2.44948974278318, 0, 0, 0, 0, 0},
 
309
    {4.24264068711928, 0, 0, 0, 0, 0},
 
310
    {2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
 
311
    {2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
 
312
    {-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
 
313
    
 
314
    // Compute reference derivatives
 
315
    // Declare pointer to array of derivatives on FIAT element
 
316
    double *derivatives = new double [num_derivatives];
 
317
    
 
318
    // Declare coefficients
 
319
    double coeff0_0 = 0;
 
320
    double coeff0_1 = 0;
 
321
    double coeff0_2 = 0;
 
322
    double coeff0_3 = 0;
 
323
    double coeff0_4 = 0;
 
324
    double coeff0_5 = 0;
 
325
    
 
326
    // Declare new coefficients
 
327
    double new_coeff0_0 = 0;
 
328
    double new_coeff0_1 = 0;
 
329
    double new_coeff0_2 = 0;
 
330
    double new_coeff0_3 = 0;
 
331
    double new_coeff0_4 = 0;
 
332
    double new_coeff0_5 = 0;
 
333
    
 
334
    // Loop possible derivatives
 
335
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
336
    {
 
337
      // Get values from coefficients array
 
338
      new_coeff0_0 = coefficients0[dof][0];
 
339
      new_coeff0_1 = coefficients0[dof][1];
 
340
      new_coeff0_2 = coefficients0[dof][2];
 
341
      new_coeff0_3 = coefficients0[dof][3];
 
342
      new_coeff0_4 = coefficients0[dof][4];
 
343
      new_coeff0_5 = coefficients0[dof][5];
 
344
    
 
345
      // Loop derivative order
 
346
      for (unsigned int j = 0; j < n; j++)
 
347
      {
 
348
        // Update old coefficients
 
349
        coeff0_0 = new_coeff0_0;
 
350
        coeff0_1 = new_coeff0_1;
 
351
        coeff0_2 = new_coeff0_2;
 
352
        coeff0_3 = new_coeff0_3;
 
353
        coeff0_4 = new_coeff0_4;
 
354
        coeff0_5 = new_coeff0_5;
 
355
    
 
356
        if(combinations[deriv_num][j] == 0)
 
357
        {
 
358
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
 
359
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
 
360
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
 
361
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
 
362
          new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
 
363
          new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
 
364
        }
 
365
        if(combinations[deriv_num][j] == 1)
 
366
        {
 
367
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
 
368
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
 
369
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
 
370
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
 
371
          new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
 
372
          new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
 
373
        }
 
374
    
 
375
      }
 
376
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
377
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
 
378
    }
 
379
    
 
380
    // Transform derivatives back to physical element
 
381
    for (unsigned int row = 0; row < num_derivatives; row++)
 
382
    {
 
383
      for (unsigned int col = 0; col < num_derivatives; col++)
 
384
      {
 
385
        values[row] += transform[row][col]*derivatives[col];
 
386
      }
 
387
    }
 
388
    // Delete pointer to array of derivatives on FIAT element
 
389
    delete [] derivatives;
 
390
    
 
391
    // Delete pointer to array of combinations of derivatives and transform
 
392
    for (unsigned int row = 0; row < num_derivatives; row++)
 
393
    {
 
394
      delete [] combinations[row];
 
395
      delete [] transform[row];
 
396
    }
 
397
    
 
398
    delete [] combinations;
 
399
    delete [] transform;
 
400
  }
 
401
 
 
402
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
403
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
404
                                              double* values,
 
405
                                              const double* coordinates,
 
406
                                              const ufc::cell& c) const
 
407
  {
 
408
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
409
  }
 
410
 
 
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
 
415
  {
 
416
    // The reference points, direction and weights:
 
417
    static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
 
418
    static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
 
419
    static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
 
420
    
 
421
    const double * const * x = c.coordinates;
 
422
    double result = 0.0;
 
423
    // Iterate over the points:
 
424
    // Evaluate basis functions for affine mapping
 
425
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
426
    const double w1 = X[i][0][0];
 
427
    const double w2 = X[i][0][1];
 
428
    
 
429
    // Compute affine mapping y = F(X)
 
430
    double y[2];
 
431
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
432
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
433
    
 
434
    // Evaluate function at physical points
 
435
    double values[1];
 
436
    f.evaluate(values, y, c);
 
437
    
 
438
    // Map function values using appropriate mapping
 
439
    // Affine map: Do nothing
 
440
    
 
441
    // Note that we do not map the weights (yet).
 
442
    
 
443
    // Take directional components
 
444
    for(int k = 0; k < 1; k++)
 
445
      result += values[k]*D[i][0][k];
 
446
    // Multiply by weights
 
447
    result *= W[i][0];
 
448
    
 
449
    return result;
 
450
  }
 
451
 
 
452
  /// Evaluate linear functionals for all dofs on the function f
 
453
  virtual void evaluate_dofs(double* values,
 
454
                             const ufc::function& f,
 
455
                             const ufc::cell& c) const
 
456
  {
 
457
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
458
  }
 
459
 
 
460
  /// Interpolate vertex values from dof values
 
461
  virtual void interpolate_vertex_values(double* vertex_values,
 
462
                                         const double* dof_values,
 
463
                                         const ufc::cell& c) const
 
464
  {
 
465
    // Evaluate at vertices and use affine mapping
 
466
    vertex_values[0] = dof_values[0];
 
467
    vertex_values[1] = dof_values[1];
 
468
    vertex_values[2] = dof_values[2];
 
469
  }
 
470
 
 
471
  /// Return the number of sub elements (for a mixed element)
 
472
  virtual unsigned int num_sub_elements() const
 
473
  {
 
474
    return 1;
 
475
  }
 
476
 
 
477
  /// Create a new finite element for sub element i (for a mixed element)
 
478
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
479
  {
 
480
    return new poisson2dp2_0_finite_element_0();
 
481
  }
 
482
 
 
483
};
 
484
 
 
485
/// This class defines the interface for a finite element.
 
486
 
 
487
class poisson2dp2_0_finite_element_1: public ufc::finite_element
 
488
{
 
489
public:
 
490
 
 
491
  /// Constructor
 
492
  poisson2dp2_0_finite_element_1() : ufc::finite_element()
 
493
  {
 
494
    // Do nothing
 
495
  }
 
496
 
 
497
  /// Destructor
 
498
  virtual ~poisson2dp2_0_finite_element_1()
 
499
  {
 
500
    // Do nothing
 
501
  }
 
502
 
 
503
  /// Return a string identifying the finite element
 
504
  virtual const char* signature() const
 
505
  {
 
506
    return "FiniteElement('Lagrange', 'triangle', 2)";
 
507
  }
 
508
 
 
509
  /// Return the cell shape
 
510
  virtual ufc::shape cell_shape() const
 
511
  {
 
512
    return ufc::triangle;
 
513
  }
 
514
 
 
515
  /// Return the dimension of the finite element function space
 
516
  virtual unsigned int space_dimension() const
 
517
  {
 
518
    return 6;
 
519
  }
 
520
 
 
521
  /// Return the rank of the value space
 
522
  virtual unsigned int value_rank() const
 
523
  {
 
524
    return 0;
 
525
  }
 
526
 
 
527
  /// Return the dimension of the value space for axis i
 
528
  virtual unsigned int value_dimension(unsigned int i) const
 
529
  {
 
530
    return 1;
 
531
  }
 
532
 
 
533
  /// Evaluate basis function i at given point in cell
 
534
  virtual void evaluate_basis(unsigned int i,
 
535
                              double* values,
 
536
                              const double* coordinates,
 
537
                              const ufc::cell& c) const
 
538
  {
 
539
    // Extract vertex coordinates
 
540
    const double * const * element_coordinates = c.coordinates;
 
541
    
 
542
    // Compute Jacobian of affine map from reference cell
 
543
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
544
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
545
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
546
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
547
    
 
548
    // Compute determinant of Jacobian
 
549
    const double detJ = J_00*J_11 - J_01*J_10;
 
550
    
 
551
    // Compute inverse of Jacobian
 
552
    
 
553
    // Get coordinates and map to the reference (UFC) element
 
554
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
555
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
556
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
557
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
558
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
559
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
560
    
 
561
    // Map coordinates to the reference square
 
562
    if (std::abs(y - 1.0) < 1e-14)
 
563
      x = -1.0;
 
564
    else
 
565
      x = 2.0 *x/(1.0 - y) - 1.0;
 
566
    y = 2.0*y - 1.0;
 
567
    
 
568
    // Reset values
 
569
    *values = 0;
 
570
    
 
571
    // Map degree of freedom to element degree of freedom
 
572
    const unsigned int dof = i;
 
573
    
 
574
    // Generate scalings
 
575
    const double scalings_y_0 = 1;
 
576
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
577
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
 
578
    
 
579
    // Compute psitilde_a
 
580
    const double psitilde_a_0 = 1;
 
581
    const double psitilde_a_1 = x;
 
582
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
 
583
    
 
584
    // Compute psitilde_bs
 
585
    const double psitilde_bs_0_0 = 1;
 
586
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
587
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
 
588
    const double psitilde_bs_1_0 = 1;
 
589
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
 
590
    const double psitilde_bs_2_0 = 1;
 
591
    
 
592
    // Compute basisvalues
 
593
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
594
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
595
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
596
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
 
597
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
 
598
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
 
599
    
 
600
    // Table(s) of coefficients
 
601
    static const double coefficients0[6][6] = \
 
602
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
 
603
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
 
604
    {0, 0, 0.2, 0, 0, 0.163299316185545},
 
605
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
 
606
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
 
607
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
 
608
    
 
609
    // Extract relevant coefficients
 
610
    const double coeff0_0 = coefficients0[dof][0];
 
611
    const double coeff0_1 = coefficients0[dof][1];
 
612
    const double coeff0_2 = coefficients0[dof][2];
 
613
    const double coeff0_3 = coefficients0[dof][3];
 
614
    const double coeff0_4 = coefficients0[dof][4];
 
615
    const double coeff0_5 = coefficients0[dof][5];
 
616
    
 
617
    // Compute value(s)
 
618
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
 
619
  }
 
620
 
 
621
  /// Evaluate all basis functions at given point in cell
 
622
  virtual void evaluate_basis_all(double* values,
 
623
                                  const double* coordinates,
 
624
                                  const ufc::cell& c) const
 
625
  {
 
626
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
627
  }
 
628
 
 
629
  /// Evaluate order n derivatives of basis function i at given point in cell
 
630
  virtual void evaluate_basis_derivatives(unsigned int i,
 
631
                                          unsigned int n,
 
632
                                          double* values,
 
633
                                          const double* coordinates,
 
634
                                          const ufc::cell& c) const
 
635
  {
 
636
    // Extract vertex coordinates
 
637
    const double * const * element_coordinates = c.coordinates;
 
638
    
 
639
    // Compute Jacobian of affine map from reference cell
 
640
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
641
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
642
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
643
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
644
    
 
645
    // Compute determinant of Jacobian
 
646
    const double detJ = J_00*J_11 - J_01*J_10;
 
647
    
 
648
    // Compute inverse of Jacobian
 
649
    
 
650
    // Get coordinates and map to the reference (UFC) element
 
651
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
652
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
653
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
654
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
655
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
656
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
657
    
 
658
    // Map coordinates to the reference square
 
659
    if (std::abs(y - 1.0) < 1e-14)
 
660
      x = -1.0;
 
661
    else
 
662
      x = 2.0 *x/(1.0 - y) - 1.0;
 
663
    y = 2.0*y - 1.0;
 
664
    
 
665
    // Compute number of derivatives
 
666
    unsigned int num_derivatives = 1;
 
667
    
 
668
    for (unsigned int j = 0; j < n; j++)
 
669
      num_derivatives *= 2;
 
670
    
 
671
    
 
672
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
673
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
674
    
 
675
    for (unsigned int j = 0; j < num_derivatives; j++)
 
676
    {
 
677
      combinations[j] = new unsigned int [n];
 
678
      for (unsigned int k = 0; k < n; k++)
 
679
        combinations[j][k] = 0;
 
680
    }
 
681
    
 
682
    // Generate combinations of derivatives
 
683
    for (unsigned int row = 1; row < num_derivatives; row++)
 
684
    {
 
685
      for (unsigned int num = 0; num < row; num++)
 
686
      {
 
687
        for (unsigned int col = n-1; col+1 > 0; col--)
 
688
        {
 
689
          if (combinations[row][col] + 1 > 1)
 
690
            combinations[row][col] = 0;
 
691
          else
 
692
          {
 
693
            combinations[row][col] += 1;
 
694
            break;
 
695
          }
 
696
        }
 
697
      }
 
698
    }
 
699
    
 
700
    // Compute inverse of Jacobian
 
701
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
702
    
 
703
    // Declare transformation matrix
 
704
    // Declare pointer to two dimensional array and initialise
 
705
    double **transform = new double *[num_derivatives];
 
706
    
 
707
    for (unsigned int j = 0; j < num_derivatives; j++)
 
708
    {
 
709
      transform[j] = new double [num_derivatives];
 
710
      for (unsigned int k = 0; k < num_derivatives; k++)
 
711
        transform[j][k] = 1;
 
712
    }
 
713
    
 
714
    // Construct transformation matrix
 
715
    for (unsigned int row = 0; row < num_derivatives; row++)
 
716
    {
 
717
      for (unsigned int col = 0; col < num_derivatives; col++)
 
718
      {
 
719
        for (unsigned int k = 0; k < n; k++)
 
720
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
721
      }
 
722
    }
 
723
    
 
724
    // Reset values
 
725
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
726
      values[j] = 0;
 
727
    
 
728
    // Map degree of freedom to element degree of freedom
 
729
    const unsigned int dof = i;
 
730
    
 
731
    // Generate scalings
 
732
    const double scalings_y_0 = 1;
 
733
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
734
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
 
735
    
 
736
    // Compute psitilde_a
 
737
    const double psitilde_a_0 = 1;
 
738
    const double psitilde_a_1 = x;
 
739
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
 
740
    
 
741
    // Compute psitilde_bs
 
742
    const double psitilde_bs_0_0 = 1;
 
743
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
744
    const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
 
745
    const double psitilde_bs_1_0 = 1;
 
746
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
 
747
    const double psitilde_bs_2_0 = 1;
 
748
    
 
749
    // Compute basisvalues
 
750
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
751
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
752
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
753
    const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
 
754
    const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
 
755
    const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
 
756
    
 
757
    // Table(s) of coefficients
 
758
    static const double coefficients0[6][6] = \
 
759
    {{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582064, 0.0544331053951817},
 
760
    {0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582063, 0.0544331053951818},
 
761
    {0, 0, 0.2, 0, 0, 0.163299316185545},
 
762
    {0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},
 
763
    {0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},
 
764
    {0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};
 
765
    
 
766
    // Interesting (new) part
 
767
    // Tables of derivatives of the polynomial base (transpose)
 
768
    static const double dmats0[6][6] = \
 
769
    {{0, 0, 0, 0, 0, 0},
 
770
    {4.89897948556635, 0, 0, 0, 0, 0},
 
771
    {0, 0, 0, 0, 0, 0},
 
772
    {0, 9.48683298050514, 0, 0, 0, 0},
 
773
    {4, 0, 7.07106781186548, 0, 0, 0},
 
774
    {0, 0, 0, 0, 0, 0}};
 
775
    
 
776
    static const double dmats1[6][6] = \
 
777
    {{0, 0, 0, 0, 0, 0},
 
778
    {2.44948974278318, 0, 0, 0, 0, 0},
 
779
    {4.24264068711928, 0, 0, 0, 0, 0},
 
780
    {2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},
 
781
    {2, 6.12372435695795, 3.53553390593274, 0, 0, 0},
 
782
    {-2.3094010767585, 0, 8.16496580927726, 0, 0, 0}};
 
783
    
 
784
    // Compute reference derivatives
 
785
    // Declare pointer to array of derivatives on FIAT element
 
786
    double *derivatives = new double [num_derivatives];
 
787
    
 
788
    // Declare coefficients
 
789
    double coeff0_0 = 0;
 
790
    double coeff0_1 = 0;
 
791
    double coeff0_2 = 0;
 
792
    double coeff0_3 = 0;
 
793
    double coeff0_4 = 0;
 
794
    double coeff0_5 = 0;
 
795
    
 
796
    // Declare new coefficients
 
797
    double new_coeff0_0 = 0;
 
798
    double new_coeff0_1 = 0;
 
799
    double new_coeff0_2 = 0;
 
800
    double new_coeff0_3 = 0;
 
801
    double new_coeff0_4 = 0;
 
802
    double new_coeff0_5 = 0;
 
803
    
 
804
    // Loop possible derivatives
 
805
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
806
    {
 
807
      // Get values from coefficients array
 
808
      new_coeff0_0 = coefficients0[dof][0];
 
809
      new_coeff0_1 = coefficients0[dof][1];
 
810
      new_coeff0_2 = coefficients0[dof][2];
 
811
      new_coeff0_3 = coefficients0[dof][3];
 
812
      new_coeff0_4 = coefficients0[dof][4];
 
813
      new_coeff0_5 = coefficients0[dof][5];
 
814
    
 
815
      // Loop derivative order
 
816
      for (unsigned int j = 0; j < n; j++)
 
817
      {
 
818
        // Update old coefficients
 
819
        coeff0_0 = new_coeff0_0;
 
820
        coeff0_1 = new_coeff0_1;
 
821
        coeff0_2 = new_coeff0_2;
 
822
        coeff0_3 = new_coeff0_3;
 
823
        coeff0_4 = new_coeff0_4;
 
824
        coeff0_5 = new_coeff0_5;
 
825
    
 
826
        if(combinations[deriv_num][j] == 0)
 
827
        {
 
828
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
 
829
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
 
830
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
 
831
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
 
832
          new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
 
833
          new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
 
834
        }
 
835
        if(combinations[deriv_num][j] == 1)
 
836
        {
 
837
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
 
838
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
 
839
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
 
840
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
 
841
          new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
 
842
          new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
 
843
        }
 
844
    
 
845
      }
 
846
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
847
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
 
848
    }
 
849
    
 
850
    // Transform derivatives back to physical element
 
851
    for (unsigned int row = 0; row < num_derivatives; row++)
 
852
    {
 
853
      for (unsigned int col = 0; col < num_derivatives; col++)
 
854
      {
 
855
        values[row] += transform[row][col]*derivatives[col];
 
856
      }
 
857
    }
 
858
    // Delete pointer to array of derivatives on FIAT element
 
859
    delete [] derivatives;
 
860
    
 
861
    // Delete pointer to array of combinations of derivatives and transform
 
862
    for (unsigned int row = 0; row < num_derivatives; row++)
 
863
    {
 
864
      delete [] combinations[row];
 
865
      delete [] transform[row];
 
866
    }
 
867
    
 
868
    delete [] combinations;
 
869
    delete [] transform;
 
870
  }
 
871
 
 
872
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
873
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
874
                                              double* values,
 
875
                                              const double* coordinates,
 
876
                                              const ufc::cell& c) const
 
877
  {
 
878
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
879
  }
 
880
 
 
881
  /// Evaluate linear functional for dof i on the function f
 
882
  virtual double evaluate_dof(unsigned int i,
 
883
                              const ufc::function& f,
 
884
                              const ufc::cell& c) const
 
885
  {
 
886
    // The reference points, direction and weights:
 
887
    static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
 
888
    static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
 
889
    static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
 
890
    
 
891
    const double * const * x = c.coordinates;
 
892
    double result = 0.0;
 
893
    // Iterate over the points:
 
894
    // Evaluate basis functions for affine mapping
 
895
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
896
    const double w1 = X[i][0][0];
 
897
    const double w2 = X[i][0][1];
 
898
    
 
899
    // Compute affine mapping y = F(X)
 
900
    double y[2];
 
901
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
902
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
903
    
 
904
    // Evaluate function at physical points
 
905
    double values[1];
 
906
    f.evaluate(values, y, c);
 
907
    
 
908
    // Map function values using appropriate mapping
 
909
    // Affine map: Do nothing
 
910
    
 
911
    // Note that we do not map the weights (yet).
 
912
    
 
913
    // Take directional components
 
914
    for(int k = 0; k < 1; k++)
 
915
      result += values[k]*D[i][0][k];
 
916
    // Multiply by weights
 
917
    result *= W[i][0];
 
918
    
 
919
    return result;
 
920
  }
 
921
 
 
922
  /// Evaluate linear functionals for all dofs on the function f
 
923
  virtual void evaluate_dofs(double* values,
 
924
                             const ufc::function& f,
 
925
                             const ufc::cell& c) const
 
926
  {
 
927
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
928
  }
 
929
 
 
930
  /// Interpolate vertex values from dof values
 
931
  virtual void interpolate_vertex_values(double* vertex_values,
 
932
                                         const double* dof_values,
 
933
                                         const ufc::cell& c) const
 
934
  {
 
935
    // Evaluate at vertices and use affine mapping
 
936
    vertex_values[0] = dof_values[0];
 
937
    vertex_values[1] = dof_values[1];
 
938
    vertex_values[2] = dof_values[2];
 
939
  }
 
940
 
 
941
  /// Return the number of sub elements (for a mixed element)
 
942
  virtual unsigned int num_sub_elements() const
 
943
  {
 
944
    return 1;
 
945
  }
 
946
 
 
947
  /// Create a new finite element for sub element i (for a mixed element)
 
948
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
949
  {
 
950
    return new poisson2dp2_0_finite_element_1();
 
951
  }
 
952
 
 
953
};
 
954
 
 
955
/// This class defines the interface for a local-to-global mapping of
 
956
/// degrees of freedom (dofs).
 
957
 
 
958
class poisson2dp2_0_dof_map_0: public ufc::dof_map
 
959
{
 
960
private:
 
961
 
 
962
  unsigned int __global_dimension;
 
963
 
 
964
public:
 
965
 
 
966
  /// Constructor
 
967
  poisson2dp2_0_dof_map_0() : ufc::dof_map()
 
968
  {
 
969
    __global_dimension = 0;
 
970
  }
 
971
 
 
972
  /// Destructor
 
973
  virtual ~poisson2dp2_0_dof_map_0()
 
974
  {
 
975
    // Do nothing
 
976
  }
 
977
 
 
978
  /// Return a string identifying the dof map
 
979
  virtual const char* signature() const
 
980
  {
 
981
    return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
 
982
  }
 
983
 
 
984
  /// Return true iff mesh entities of topological dimension d are needed
 
985
  virtual bool needs_mesh_entities(unsigned int d) const
 
986
  {
 
987
    switch ( d )
 
988
    {
 
989
    case 0:
 
990
      return true;
 
991
      break;
 
992
    case 1:
 
993
      return true;
 
994
      break;
 
995
    case 2:
 
996
      return false;
 
997
      break;
 
998
    }
 
999
    return false;
 
1000
  }
 
1001
 
 
1002
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1003
  virtual bool init_mesh(const ufc::mesh& m)
 
1004
  {
 
1005
    __global_dimension = m.num_entities[0] + m.num_entities[1];
 
1006
    return false;
 
1007
  }
 
1008
 
 
1009
  /// Initialize dof map for given cell
 
1010
  virtual void init_cell(const ufc::mesh& m,
 
1011
                         const ufc::cell& c)
 
1012
  {
 
1013
    // Do nothing
 
1014
  }
 
1015
 
 
1016
  /// Finish initialization of dof map for cells
 
1017
  virtual void init_cell_finalize()
 
1018
  {
 
1019
    // Do nothing
 
1020
  }
 
1021
 
 
1022
  /// Return the dimension of the global finite element function space
 
1023
  virtual unsigned int global_dimension() const
 
1024
  {
 
1025
    return __global_dimension;
 
1026
  }
 
1027
 
 
1028
  /// Return the dimension of the local finite element function space for a cell
 
1029
  virtual unsigned int local_dimension(const ufc::cell& c) const
 
1030
  {
 
1031
    return 6;
 
1032
  }
 
1033
 
 
1034
  /// Return the maximum dimension of the local finite element function space
 
1035
  virtual unsigned int max_local_dimension() const
 
1036
  {
 
1037
    return 6;
 
1038
  }
 
1039
 
 
1040
  // Return the geometric dimension of the coordinates this dof map provides
 
1041
  virtual unsigned int geometric_dimension() const
 
1042
  {
 
1043
    return 2;
 
1044
  }
 
1045
 
 
1046
  /// Return the number of dofs on each cell facet
 
1047
  virtual unsigned int num_facet_dofs() const
 
1048
  {
 
1049
    return 3;
 
1050
  }
 
1051
 
 
1052
  /// Return the number of dofs associated with each cell entity of dimension d
 
1053
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1054
  {
 
1055
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1056
  }
 
1057
 
 
1058
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1059
  virtual void tabulate_dofs(unsigned int* dofs,
 
1060
                             const ufc::mesh& m,
 
1061
                             const ufc::cell& c) const
 
1062
  {
 
1063
    dofs[0] = c.entity_indices[0][0];
 
1064
    dofs[1] = c.entity_indices[0][1];
 
1065
    dofs[2] = c.entity_indices[0][2];
 
1066
    unsigned int offset = m.num_entities[0];
 
1067
    dofs[3] = offset + c.entity_indices[1][0];
 
1068
    dofs[4] = offset + c.entity_indices[1][1];
 
1069
    dofs[5] = offset + c.entity_indices[1][2];
 
1070
  }
 
1071
 
 
1072
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1073
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1074
                                   unsigned int facet) const
 
1075
  {
 
1076
    switch ( facet )
 
1077
    {
 
1078
    case 0:
 
1079
      dofs[0] = 1;
 
1080
      dofs[1] = 2;
 
1081
      dofs[2] = 3;
 
1082
      break;
 
1083
    case 1:
 
1084
      dofs[0] = 0;
 
1085
      dofs[1] = 2;
 
1086
      dofs[2] = 4;
 
1087
      break;
 
1088
    case 2:
 
1089
      dofs[0] = 0;
 
1090
      dofs[1] = 1;
 
1091
      dofs[2] = 5;
 
1092
      break;
 
1093
    }
 
1094
  }
 
1095
 
 
1096
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1097
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1098
                                    unsigned int d, unsigned int i) const
 
1099
  {
 
1100
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1101
  }
 
1102
 
 
1103
  /// Tabulate the coordinates of all dofs on a cell
 
1104
  virtual void tabulate_coordinates(double** coordinates,
 
1105
                                    const ufc::cell& c) const
 
1106
  {
 
1107
    const double * const * x = c.coordinates;
 
1108
    coordinates[0][0] = x[0][0];
 
1109
    coordinates[0][1] = x[0][1];
 
1110
    coordinates[1][0] = x[1][0];
 
1111
    coordinates[1][1] = x[1][1];
 
1112
    coordinates[2][0] = x[2][0];
 
1113
    coordinates[2][1] = x[2][1];
 
1114
    coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
 
1115
    coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
 
1116
    coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
 
1117
    coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
 
1118
    coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
 
1119
    coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
 
1120
  }
 
1121
 
 
1122
  /// Return the number of sub dof maps (for a mixed element)
 
1123
  virtual unsigned int num_sub_dof_maps() const
 
1124
  {
 
1125
    return 1;
 
1126
  }
 
1127
 
 
1128
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1129
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1130
  {
 
1131
    return new poisson2dp2_0_dof_map_0();
 
1132
  }
 
1133
 
 
1134
};
 
1135
 
 
1136
/// This class defines the interface for a local-to-global mapping of
 
1137
/// degrees of freedom (dofs).
 
1138
 
 
1139
class poisson2dp2_0_dof_map_1: public ufc::dof_map
 
1140
{
 
1141
private:
 
1142
 
 
1143
  unsigned int __global_dimension;
 
1144
 
 
1145
public:
 
1146
 
 
1147
  /// Constructor
 
1148
  poisson2dp2_0_dof_map_1() : ufc::dof_map()
 
1149
  {
 
1150
    __global_dimension = 0;
 
1151
  }
 
1152
 
 
1153
  /// Destructor
 
1154
  virtual ~poisson2dp2_0_dof_map_1()
 
1155
  {
 
1156
    // Do nothing
 
1157
  }
 
1158
 
 
1159
  /// Return a string identifying the dof map
 
1160
  virtual const char* signature() const
 
1161
  {
 
1162
    return "FFC dof map for FiniteElement('Lagrange', 'triangle', 2)";
 
1163
  }
 
1164
 
 
1165
  /// Return true iff mesh entities of topological dimension d are needed
 
1166
  virtual bool needs_mesh_entities(unsigned int d) const
 
1167
  {
 
1168
    switch ( d )
 
1169
    {
 
1170
    case 0:
 
1171
      return true;
 
1172
      break;
 
1173
    case 1:
 
1174
      return true;
 
1175
      break;
 
1176
    case 2:
 
1177
      return false;
 
1178
      break;
 
1179
    }
 
1180
    return false;
 
1181
  }
 
1182
 
 
1183
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1184
  virtual bool init_mesh(const ufc::mesh& m)
 
1185
  {
 
1186
    __global_dimension = m.num_entities[0] + m.num_entities[1];
 
1187
    return false;
 
1188
  }
 
1189
 
 
1190
  /// Initialize dof map for given cell
 
1191
  virtual void init_cell(const ufc::mesh& m,
 
1192
                         const ufc::cell& c)
 
1193
  {
 
1194
    // Do nothing
 
1195
  }
 
1196
 
 
1197
  /// Finish initialization of dof map for cells
 
1198
  virtual void init_cell_finalize()
 
1199
  {
 
1200
    // Do nothing
 
1201
  }
 
1202
 
 
1203
  /// Return the dimension of the global finite element function space
 
1204
  virtual unsigned int global_dimension() const
 
1205
  {
 
1206
    return __global_dimension;
 
1207
  }
 
1208
 
 
1209
  /// Return the dimension of the local finite element function space for a cell
 
1210
  virtual unsigned int local_dimension(const ufc::cell& c) const
 
1211
  {
 
1212
    return 6;
 
1213
  }
 
1214
 
 
1215
  /// Return the maximum dimension of the local finite element function space
 
1216
  virtual unsigned int max_local_dimension() const
 
1217
  {
 
1218
    return 6;
 
1219
  }
 
1220
 
 
1221
  // Return the geometric dimension of the coordinates this dof map provides
 
1222
  virtual unsigned int geometric_dimension() const
 
1223
  {
 
1224
    return 2;
 
1225
  }
 
1226
 
 
1227
  /// Return the number of dofs on each cell facet
 
1228
  virtual unsigned int num_facet_dofs() const
 
1229
  {
 
1230
    return 3;
 
1231
  }
 
1232
 
 
1233
  /// Return the number of dofs associated with each cell entity of dimension d
 
1234
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1235
  {
 
1236
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1237
  }
 
1238
 
 
1239
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1240
  virtual void tabulate_dofs(unsigned int* dofs,
 
1241
                             const ufc::mesh& m,
 
1242
                             const ufc::cell& c) const
 
1243
  {
 
1244
    dofs[0] = c.entity_indices[0][0];
 
1245
    dofs[1] = c.entity_indices[0][1];
 
1246
    dofs[2] = c.entity_indices[0][2];
 
1247
    unsigned int offset = m.num_entities[0];
 
1248
    dofs[3] = offset + c.entity_indices[1][0];
 
1249
    dofs[4] = offset + c.entity_indices[1][1];
 
1250
    dofs[5] = offset + c.entity_indices[1][2];
 
1251
  }
 
1252
 
 
1253
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1254
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1255
                                   unsigned int facet) const
 
1256
  {
 
1257
    switch ( facet )
 
1258
    {
 
1259
    case 0:
 
1260
      dofs[0] = 1;
 
1261
      dofs[1] = 2;
 
1262
      dofs[2] = 3;
 
1263
      break;
 
1264
    case 1:
 
1265
      dofs[0] = 0;
 
1266
      dofs[1] = 2;
 
1267
      dofs[2] = 4;
 
1268
      break;
 
1269
    case 2:
 
1270
      dofs[0] = 0;
 
1271
      dofs[1] = 1;
 
1272
      dofs[2] = 5;
 
1273
      break;
 
1274
    }
 
1275
  }
 
1276
 
 
1277
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1278
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1279
                                    unsigned int d, unsigned int i) const
 
1280
  {
 
1281
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1282
  }
 
1283
 
 
1284
  /// Tabulate the coordinates of all dofs on a cell
 
1285
  virtual void tabulate_coordinates(double** coordinates,
 
1286
                                    const ufc::cell& c) const
 
1287
  {
 
1288
    const double * const * x = c.coordinates;
 
1289
    coordinates[0][0] = x[0][0];
 
1290
    coordinates[0][1] = x[0][1];
 
1291
    coordinates[1][0] = x[1][0];
 
1292
    coordinates[1][1] = x[1][1];
 
1293
    coordinates[2][0] = x[2][0];
 
1294
    coordinates[2][1] = x[2][1];
 
1295
    coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
 
1296
    coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
 
1297
    coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
 
1298
    coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
 
1299
    coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
 
1300
    coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
 
1301
  }
 
1302
 
 
1303
  /// Return the number of sub dof maps (for a mixed element)
 
1304
  virtual unsigned int num_sub_dof_maps() const
 
1305
  {
 
1306
    return 1;
 
1307
  }
 
1308
 
 
1309
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1310
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1311
  {
 
1312
    return new poisson2dp2_0_dof_map_1();
1301
1313
  }
1302
1314
 
1303
1315
};
1306
1318
/// tensor corresponding to the local contribution to a form from
1307
1319
/// the integral over a cell.
1308
1320
 
1309
 
class UFC_Poisson2DP2BilinearForm_cell_integral_0: public ufc::cell_integral
 
1321
class poisson2dp2_0_cell_integral_0_quadrature: public ufc::cell_integral
1310
1322
{
1311
1323
public:
1312
1324
 
1313
1325
  /// Constructor
1314
 
  UFC_Poisson2DP2BilinearForm_cell_integral_0() : ufc::cell_integral()
 
1326
  poisson2dp2_0_cell_integral_0_quadrature() : ufc::cell_integral()
1315
1327
  {
1316
1328
    // Do nothing
1317
1329
  }
1318
1330
 
1319
1331
  /// Destructor
1320
 
  virtual ~UFC_Poisson2DP2BilinearForm_cell_integral_0()
 
1332
  virtual ~poisson2dp2_0_cell_integral_0_quadrature()
1321
1333
  {
1322
1334
    // Do nothing
1323
1335
  }
1335
1347
    const double J_01 = x[2][0] - x[0][0];
1336
1348
    const double J_10 = x[1][1] - x[0][1];
1337
1349
    const double J_11 = x[2][1] - x[0][1];
1338
 
      
 
1350
    
1339
1351
    // Compute determinant of Jacobian
1340
1352
    double detJ = J_00*J_11 - J_01*J_10;
1341
 
      
 
1353
    
1342
1354
    // Compute inverse of Jacobian
1343
1355
    const double Jinv_00 =  J_11 / detJ;
1344
1356
    const double Jinv_01 = -J_01 / detJ;
1348
1360
    // Set scale factor
1349
1361
    const double det = std::abs(detJ);
1350
1362
    
1351
 
    // Number of operations to compute element tensor = 114
1352
 
    // Compute geometry tensors
1353
 
    // Number of operations to compute decalrations = 16
1354
 
    const double G0_0_0 = det*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
1355
 
    const double G0_0_1 = det*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
1356
 
    const double G0_1_0 = det*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
1357
 
    const double G0_1_1 = det*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
1358
 
    
1359
 
    // Compute element tensor
1360
 
    // Number of operations to compute tensor = 98
1361
 
    A[0] = 0.499999999999999*G0_0_0 + 0.499999999999999*G0_0_1 + 0.499999999999999*G0_1_0 + 0.499999999999999*G0_1_1;
1362
 
    A[1] = 0.166666666666666*G0_0_0 + 0.166666666666666*G0_1_0;
1363
 
    A[2] = 0.166666666666666*G0_0_1 + 0.166666666666666*G0_1_1;
1364
 
    A[3] = 0;
1365
 
    A[4] = -0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_1;
1366
 
    A[5] = -0.666666666666666*G0_0_0 - 0.666666666666666*G0_1_0;
1367
 
    A[6] = 0.166666666666666*G0_0_0 + 0.166666666666666*G0_0_1;
1368
 
    A[7] = 0.499999999999999*G0_0_0;
1369
 
    A[8] = -0.166666666666666*G0_0_1;
1370
 
    A[9] = 0.666666666666665*G0_0_1;
1371
 
    A[10] = 0;
1372
 
    A[11] = -0.666666666666666*G0_0_0 - 0.666666666666665*G0_0_1;
1373
 
    A[12] = 0.166666666666666*G0_1_0 + 0.166666666666666*G0_1_1;
1374
 
    A[13] = -0.166666666666666*G0_1_0;
1375
 
    A[14] = 0.499999999999999*G0_1_1;
1376
 
    A[15] = 0.666666666666666*G0_1_0;
1377
 
    A[16] = -0.666666666666665*G0_1_0 - 0.666666666666665*G0_1_1;
1378
 
    A[17] = 0;
1379
 
    A[18] = 0;
1380
 
    A[19] = 0.666666666666665*G0_1_0;
1381
 
    A[20] = 0.666666666666666*G0_0_1;
1382
 
    A[21] = 1.33333333333333*G0_0_0 + 0.666666666666665*G0_0_1 + 0.666666666666665*G0_1_0 + 1.33333333333333*G0_1_1;
1383
 
    A[22] = -1.33333333333333*G0_0_0 - 0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_0;
1384
 
    A[23] = -0.666666666666666*G0_0_1 - 0.666666666666665*G0_1_0 - 1.33333333333333*G0_1_1;
1385
 
    A[24] = -0.666666666666665*G0_1_0 - 0.666666666666665*G0_1_1;
1386
 
    A[25] = 0;
1387
 
    A[26] = -0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_1;
1388
 
    A[27] = -1.33333333333333*G0_0_0 - 0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_0;
1389
 
    A[28] = 1.33333333333333*G0_0_0 + 0.666666666666665*G0_0_1 + 0.666666666666665*G0_1_0 + 1.33333333333333*G0_1_1;
1390
 
    A[29] = 0.666666666666666*G0_0_1 + 0.666666666666665*G0_1_0;
1391
 
    A[30] = -0.666666666666665*G0_0_0 - 0.666666666666666*G0_0_1;
1392
 
    A[31] = -0.666666666666665*G0_0_0 - 0.666666666666665*G0_1_0;
1393
 
    A[32] = 0;
1394
 
    A[33] = -0.666666666666665*G0_0_1 - 0.666666666666666*G0_1_0 - 1.33333333333333*G0_1_1;
1395
 
    A[34] = 0.666666666666665*G0_0_1 + 0.666666666666666*G0_1_0;
1396
 
    A[35] = 1.33333333333333*G0_0_0 + 0.666666666666666*G0_0_1 + 0.666666666666666*G0_1_0 + 1.33333333333333*G0_1_1;
 
1363
    
 
1364
    // Array of quadrature weights
 
1365
    static const double W4[4] = {0.159020690871988, 0.0909793091280113, 0.159020690871988, 0.0909793091280113};
 
1366
    // Quadrature points on the UFC reference element: (0.178558728263616, 0.155051025721682), (0.0750311102226081, 0.644948974278318), (0.666390246014701, 0.155051025721682), (0.280019915499074, 0.644948974278318)
 
1367
    
 
1368
    // Value of basis functions at quadrature points.
 
1369
    static const double FE0_D10[4][5] = \
 
1370
    {{-1.66556098405881, -0.285765086945534, 0.620204102886728, -0.620204102886728, 1.95132607100434},
 
1371
    {-0.120079661996296, -0.699875559109567, 2.57979589711327, -2.57979589711327, 0.819955221105863},
 
1372
    {0.285765086945534, 1.66556098405881, 0.620204102886728, -0.620204102886729, -1.95132607100434},
 
1373
    {0.699875559109568, 0.120079661996296, 2.57979589711327, -2.57979589711327, -0.819955221105864}};
 
1374
    
 
1375
    // Array of non-zero columns
 
1376
    static const unsigned int nzc0[5] = {0, 1, 3, 4, 5};
 
1377
    static const double FE0_D01[4][5] = \
 
1378
    {{-1.66556098405881, -0.379795897113272, 0.714234913054466, 2.04535688117208, -0.714234913054466},
 
1379
    {-0.120079661996296, 1.57979589711327, 0.300124440890432, -1.45971623511698, -0.300124440890433},
 
1380
    {0.285765086945535, -0.379795897113271, 2.66556098405881, 0.0940308101677373, -2.66556098405881},
 
1381
    {0.699875559109568, 1.57979589711327, 1.1200796619963, -2.27967145622284, -1.1200796619963}};
 
1382
    
 
1383
    // Array of non-zero columns
 
1384
    static const unsigned int nzc1[5] = {0, 2, 3, 4, 5};
 
1385
    
 
1386
    // Number of operations to compute geometry constants: 12
 
1387
    const double G0 = det*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1388
    const double G1 = det*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1389
    const double G2 = det*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1390
    
 
1391
    // Compute element tensor using UFL quadrature representation
 
1392
    // Optimisations: ('simplify expressions', True), ('ignore zero tables', True), ('non zero columns', True), ('remove zero terms', True), ('ignore ones', True)
 
1393
    // Total number of operations to compute element tensor: 1224
 
1394
    
 
1395
    // Loop quadrature points for integral
 
1396
    // Number of operations to compute element tensor for following IP loop = 1212
 
1397
    for (unsigned int ip = 0; ip < 4; ip++)
 
1398
    {
 
1399
      
 
1400
      // Number of operations to compute ip constants: 3
 
1401
      // Number of operations: 1
 
1402
      const double Gip0 = G0*W4[ip];
 
1403
      
 
1404
      // Number of operations: 1
 
1405
      const double Gip1 = G1*W4[ip];
 
1406
      
 
1407
      // Number of operations: 1
 
1408
      const double Gip2 = G2*W4[ip];
 
1409
      
 
1410
      
 
1411
      // Number of operations for primary indices = 300
 
1412
      for (unsigned int j = 0; j < 5; j++)
 
1413
      {
 
1414
        for (unsigned int k = 0; k < 5; k++)
 
1415
        {
 
1416
          // Number of operations to compute entry = 3
 
1417
          A[nzc1[j]*6 + nzc1[k]] += FE0_D01[ip][j]*FE0_D01[ip][k]*Gip0;
 
1418
          // Number of operations to compute entry = 3
 
1419
          A[nzc1[j]*6 + nzc0[k]] += FE0_D01[ip][j]*FE0_D10[ip][k]*Gip1;
 
1420
          // Number of operations to compute entry = 3
 
1421
          A[nzc0[j]*6 + nzc1[k]] += FE0_D01[ip][k]*FE0_D10[ip][j]*Gip1;
 
1422
          // Number of operations to compute entry = 3
 
1423
          A[nzc0[j]*6 + nzc0[k]] += FE0_D10[ip][j]*FE0_D10[ip][k]*Gip2;
 
1424
        }// end loop over 'k'
 
1425
      }// end loop over 'j'
 
1426
    }// end loop over 'ip'
 
1427
  }
 
1428
 
 
1429
};
 
1430
 
 
1431
/// This class defines the interface for the tabulation of the cell
 
1432
/// tensor corresponding to the local contribution to a form from
 
1433
/// the integral over a cell.
 
1434
 
 
1435
class poisson2dp2_0_cell_integral_0: public ufc::cell_integral
 
1436
{
 
1437
private:
 
1438
 
 
1439
  poisson2dp2_0_cell_integral_0_quadrature integral_0_quadrature;
 
1440
 
 
1441
public:
 
1442
 
 
1443
  /// Constructor
 
1444
  poisson2dp2_0_cell_integral_0() : ufc::cell_integral()
 
1445
  {
 
1446
    // Do nothing
 
1447
  }
 
1448
 
 
1449
  /// Destructor
 
1450
  virtual ~poisson2dp2_0_cell_integral_0()
 
1451
  {
 
1452
    // Do nothing
 
1453
  }
 
1454
 
 
1455
  /// Tabulate the tensor for the contribution from a local cell
 
1456
  virtual void tabulate_tensor(double* A,
 
1457
                               const double * const * w,
 
1458
                               const ufc::cell& c) const
 
1459
  {
 
1460
    // Reset values of the element tensor block
 
1461
    for (unsigned int j = 0; j < 36; j++)
 
1462
      A[j] = 0;
 
1463
    
 
1464
    // Add all contributions to element tensor
 
1465
    integral_0_quadrature.tabulate_tensor(A, w, c);
1397
1466
  }
1398
1467
 
1399
1468
};
1413
1482
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
1414
1483
/// fixed functions (coefficients).
1415
1484
 
1416
 
class UFC_Poisson2DP2BilinearForm: public ufc::form
 
1485
class poisson2dp2_form_0: public ufc::form
1417
1486
{
1418
1487
public:
1419
1488
 
1420
1489
  /// Constructor
1421
 
  UFC_Poisson2DP2BilinearForm() : ufc::form()
 
1490
  poisson2dp2_form_0() : ufc::form()
1422
1491
  {
1423
1492
    // Do nothing
1424
1493
  }
1425
1494
 
1426
1495
  /// Destructor
1427
 
  virtual ~UFC_Poisson2DP2BilinearForm()
 
1496
  virtual ~poisson2dp2_form_0()
1428
1497
  {
1429
1498
    // Do nothing
1430
1499
  }
1432
1501
  /// Return a string identifying the form
1433
1502
  virtual const char* signature() const
1434
1503
  {
1435
 
    return "(dXa0[0, 1]/dxb0[0, 1])(dXa1[0, 1]/dxb0[0, 1]) | ((d/dXa0[0, 1])vi0[0, 1, 2, 3, 4, 5])*((d/dXa1[0, 1])vi1[0, 1, 2, 3, 4, 5])*dX(0)";
 
1504
    return "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Measure('cell', 0, None))])";
1436
1505
  }
1437
1506
 
1438
1507
  /// Return the rank of the global tensor (r)
1452
1521
  {
1453
1522
    return 1;
1454
1523
  }
1455
 
  
 
1524
 
1456
1525
  /// Return the number of exterior facet integrals
1457
1526
  virtual unsigned int num_exterior_facet_integrals() const
1458
1527
  {
1459
1528
    return 0;
1460
1529
  }
1461
 
  
 
1530
 
1462
1531
  /// Return the number of interior facet integrals
1463
1532
  virtual unsigned int num_interior_facet_integrals() const
1464
1533
  {
1465
1534
    return 0;
1466
1535
  }
1467
 
    
 
1536
 
1468
1537
  /// Create a new finite element for argument function i
1469
1538
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
1470
1539
  {
1471
1540
    switch ( i )
1472
1541
    {
1473
1542
    case 0:
1474
 
      return new UFC_Poisson2DP2BilinearForm_finite_element_0();
 
1543
      return new poisson2dp2_0_finite_element_0();
1475
1544
      break;
1476
1545
    case 1:
1477
 
      return new UFC_Poisson2DP2BilinearForm_finite_element_1();
 
1546
      return new poisson2dp2_0_finite_element_1();
1478
1547
      break;
1479
1548
    }
1480
1549
    return 0;
1481
1550
  }
1482
 
  
 
1551
 
1483
1552
  /// Create a new dof map for argument function i
1484
1553
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
1485
1554
  {
1486
1555
    switch ( i )
1487
1556
    {
1488
1557
    case 0:
1489
 
      return new UFC_Poisson2DP2BilinearForm_dof_map_0();
 
1558
      return new poisson2dp2_0_dof_map_0();
1490
1559
      break;
1491
1560
    case 1:
1492
 
      return new UFC_Poisson2DP2BilinearForm_dof_map_1();
 
1561
      return new poisson2dp2_0_dof_map_1();
1493
1562
      break;
1494
1563
    }
1495
1564
    return 0;
1498
1567
  /// Create a new cell integral on sub domain i
1499
1568
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
1500
1569
  {
1501
 
    return new UFC_Poisson2DP2BilinearForm_cell_integral_0();
 
1570
    return new poisson2dp2_0_cell_integral_0();
1502
1571
  }
1503
1572
 
1504
1573
  /// Create a new exterior facet integral on sub domain i
1517
1586
 
1518
1587
// DOLFIN wrappers
1519
1588
 
1520
 
#include <dolfin/fem/Form.h>
 
1589
// Standard library includes
 
1590
#include <string>
 
1591
 
 
1592
// DOLFIN includes
 
1593
#include <dolfin/common/NoDeleter.h>
1521
1594
#include <dolfin/fem/FiniteElement.h>
1522
1595
#include <dolfin/fem/DofMap.h>
1523
 
#include <dolfin/function/Coefficient.h>
1524
 
#include <dolfin/function/Function.h>
 
1596
#include <dolfin/fem/Form.h>
1525
1597
#include <dolfin/function/FunctionSpace.h>
1526
 
 
1527
 
class Poisson2DP2BilinearFormFunctionSpace0 : public dolfin::FunctionSpace
1528
 
{
1529
 
public:
1530
 
 
1531
 
  Poisson2DP2BilinearFormFunctionSpace0(const dolfin::Mesh& mesh)
1532
 
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
1533
 
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2DP2BilinearForm_finite_element_1()))),
1534
 
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2DP2BilinearForm_dof_map_1()), mesh)))
1535
 
  {
1536
 
    // Do nothing
1537
 
  }
1538
 
 
1539
 
};
1540
 
 
1541
 
class Poisson2DP2BilinearFormFunctionSpace1 : public dolfin::FunctionSpace
1542
 
{
1543
 
public:
1544
 
 
1545
 
  Poisson2DP2BilinearFormFunctionSpace1(const dolfin::Mesh& mesh)
1546
 
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
1547
 
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2DP2BilinearForm_finite_element_1()))),
1548
 
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2DP2BilinearForm_dof_map_1()), mesh)))
1549
 
  {
1550
 
    // Do nothing
1551
 
  }
1552
 
 
1553
 
};
1554
 
 
1555
 
class Poisson2DP2TestSpace : public dolfin::FunctionSpace
1556
 
{
1557
 
public:
1558
 
 
1559
 
  Poisson2DP2TestSpace(const dolfin::Mesh& mesh)
1560
 
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
1561
 
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2DP2BilinearForm_finite_element_1()))),
1562
 
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2DP2BilinearForm_dof_map_1()), mesh)))
1563
 
  {
1564
 
    // Do nothing
1565
 
  }
1566
 
 
1567
 
};
1568
 
 
1569
 
class Poisson2DP2TrialSpace : public dolfin::FunctionSpace
1570
 
{
1571
 
public:
1572
 
 
1573
 
  Poisson2DP2TrialSpace(const dolfin::Mesh& mesh)
1574
 
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
1575
 
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2DP2BilinearForm_finite_element_1()))),
1576
 
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2DP2BilinearForm_dof_map_1()), mesh)))
1577
 
  {
1578
 
    // Do nothing
1579
 
  }
1580
 
 
1581
 
};
1582
 
 
1583
 
class Poisson2DP2FunctionSpace : public dolfin::FunctionSpace
1584
 
{
1585
 
public:
1586
 
 
1587
 
  Poisson2DP2FunctionSpace(const dolfin::Mesh& mesh)
1588
 
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
1589
 
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2DP2BilinearForm_finite_element_1()))),
1590
 
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2DP2BilinearForm_dof_map_1()), mesh)))
1591
 
  {
1592
 
    // Do nothing
1593
 
  }
1594
 
 
1595
 
};
1596
 
 
1597
 
class Poisson2DP2BilinearForm : public dolfin::Form
1598
 
{
1599
 
public:
1600
 
 
1601
 
  // Create form on given function space(s)
1602
 
  Poisson2DP2BilinearForm(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1) : dolfin::Form()
1603
 
  {
1604
 
    boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());
1605
 
    _function_spaces.push_back(_V0);
1606
 
    boost::shared_ptr<const dolfin::FunctionSpace> _V1(&V1, dolfin::NoDeleter<const dolfin::FunctionSpace>());
1607
 
    _function_spaces.push_back(_V1);
1608
 
 
1609
 
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2DP2BilinearForm());
1610
 
  }
1611
 
 
1612
 
  // Create form on given function space(s) (shared data)
1613
 
  Poisson2DP2BilinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::FunctionSpace> V1) : dolfin::Form()
1614
 
  {
1615
 
    _function_spaces.push_back(V0);
1616
 
    _function_spaces.push_back(V1);
1617
 
 
1618
 
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2DP2BilinearForm());
 
1598
#include <dolfin/function/GenericFunction.h>
 
1599
#include <dolfin/function/CoefficientAssigner.h>
 
1600
 
 
1601
namespace Poisson2DP2
 
1602
{
 
1603
 
 
1604
class Form_0_FunctionSpace_0: public dolfin::FunctionSpace
 
1605
{
 
1606
public:
 
1607
 
 
1608
  Form_0_FunctionSpace_0(const dolfin::Mesh& mesh):
 
1609
      dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
 
1610
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_0()))),
 
1611
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_0()), dolfin::reference_to_no_delete_pointer(mesh))))
 
1612
  {
 
1613
    // Do nothing
 
1614
  }
 
1615
 
 
1616
  Form_0_FunctionSpace_0(dolfin::Mesh& mesh):
 
1617
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
 
1618
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_0()))),
 
1619
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_0()), dolfin::reference_to_no_delete_pointer(mesh))))
 
1620
  {
 
1621
    // Do nothing
 
1622
  }
 
1623
 
 
1624
  Form_0_FunctionSpace_0(boost::shared_ptr<dolfin::Mesh> mesh):
 
1625
      dolfin::FunctionSpace(mesh,
 
1626
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_0()))),
 
1627
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_0()), mesh)))
 
1628
  {
 
1629
      // Do nothing
 
1630
  }
 
1631
 
 
1632
  Form_0_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
 
1633
      dolfin::FunctionSpace(mesh,
 
1634
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_0()))),
 
1635
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_0()), mesh)))
 
1636
  {
 
1637
      // Do nothing
 
1638
  }
 
1639
 
 
1640
 
 
1641
  ~Form_0_FunctionSpace_0()
 
1642
  {
 
1643
  }
 
1644
 
 
1645
};
 
1646
 
 
1647
class Form_0_FunctionSpace_1: public dolfin::FunctionSpace
 
1648
{
 
1649
public:
 
1650
 
 
1651
  Form_0_FunctionSpace_1(const dolfin::Mesh& mesh):
 
1652
      dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
 
1653
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_1()))),
 
1654
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_1()), dolfin::reference_to_no_delete_pointer(mesh))))
 
1655
  {
 
1656
    // Do nothing
 
1657
  }
 
1658
 
 
1659
  Form_0_FunctionSpace_1(dolfin::Mesh& mesh):
 
1660
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
 
1661
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_1()))),
 
1662
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_1()), dolfin::reference_to_no_delete_pointer(mesh))))
 
1663
  {
 
1664
    // Do nothing
 
1665
  }
 
1666
 
 
1667
  Form_0_FunctionSpace_1(boost::shared_ptr<dolfin::Mesh> mesh):
 
1668
      dolfin::FunctionSpace(mesh,
 
1669
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_1()))),
 
1670
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_1()), mesh)))
 
1671
  {
 
1672
      // Do nothing
 
1673
  }
 
1674
 
 
1675
  Form_0_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh):
 
1676
      dolfin::FunctionSpace(mesh,
 
1677
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson2dp2_0_finite_element_1()))),
 
1678
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson2dp2_0_dof_map_1()), mesh)))
 
1679
  {
 
1680
      // Do nothing
 
1681
  }
 
1682
 
 
1683
 
 
1684
  ~Form_0_FunctionSpace_1()
 
1685
  {
 
1686
  }
 
1687
 
 
1688
};
 
1689
 
 
1690
class Form_0: public dolfin::Form
 
1691
{
 
1692
public:
 
1693
 
 
1694
  // Constructor
 
1695
  Form_0(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1):
 
1696
    dolfin::Form(2, 0)
 
1697
  {
 
1698
    _function_spaces[0] = reference_to_no_delete_pointer(V0);
 
1699
    _function_spaces[1] = reference_to_no_delete_pointer(V1);
 
1700
 
 
1701
    _ufc_form = boost::shared_ptr<const ufc::form>(new poisson2dp2_form_0());
 
1702
  }
 
1703
 
 
1704
  // Constructor
 
1705
  Form_0(boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::FunctionSpace> V1):
 
1706
    dolfin::Form(2, 0)
 
1707
  {
 
1708
    _function_spaces[0] = V0;
 
1709
    _function_spaces[1] = V1;
 
1710
 
 
1711
    _ufc_form = boost::shared_ptr<const ufc::form>(new poisson2dp2_form_0());
1619
1712
  }
1620
1713
 
1621
1714
  // Destructor
1622
 
  ~Poisson2DP2BilinearForm() {}
1623
 
 
 
1715
  ~Form_0()
 
1716
  {}
 
1717
 
 
1718
  /// Return the number of the coefficient with this name
 
1719
  virtual dolfin::uint coefficient_number(const std::string& name) const
 
1720
  {
 
1721
 
 
1722
    dolfin::error("No coefficients.");
 
1723
    return 0;
 
1724
  }
 
1725
 
 
1726
  /// Return the name of the coefficient with this number
 
1727
  virtual std::string coefficient_name(dolfin::uint i) const
 
1728
  {
 
1729
 
 
1730
    dolfin::error("No coefficients.");
 
1731
    return "unnamed";
 
1732
  }
 
1733
 
 
1734
  // Typedefs
 
1735
  typedef Form_0_FunctionSpace_0 TestSpace;
 
1736
  typedef Form_0_FunctionSpace_1 TrialSpace;
 
1737
 
 
1738
  // Coefficients
1624
1739
};
1625
1740
 
 
1741
// Class typedefs
 
1742
typedef Form_0 BilinearForm;
 
1743
typedef Form_0::TestSpace FunctionSpace;
 
1744
 
 
1745
} // namespace Poisson2DP2
 
1746
 
1626
1747
#endif