← Back to branch summary

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

« back to all changes in this revision

Viewing changes to demo/pde/dg/poisson/cpp/P1Projection.h

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
Import upstream version 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/pde/dg/poisson/cpp/P1Projection.h
 
1
// This code conforms with the UFC specification version 1.0
 
2
// and was automatically generated by FFC version 0.6.0.
 
3
//
 
4
// Warning: This code was generated with the option '-l dolfin'
 
5
// and contains DOLFIN-specific wrappers that depend on DOLFIN.
 
6
 
 
7
#ifndef __P1PROJECTION_H
 
8
#define __P1PROJECTION_H
 
9
 
 
10
#include <cmath>
 
11
#include <stdexcept>
 
12
#include <fstream>
 
13
#include <ufc.h>
 
14
 
 
15
/// This class defines the interface for a finite element.
 
16
 
 
17
class UFC_P1ProjectionBilinearForm_finite_element_0: public ufc::finite_element
 
18
{
 
19
public:
 
20
 
 
21
  /// Constructor
 
22
  UFC_P1ProjectionBilinearForm_finite_element_0() : ufc::finite_element()
 
23
  {
 
24
    // Do nothing
 
25
  }
 
26
 
 
27
  /// Destructor
 
28
  virtual ~UFC_P1ProjectionBilinearForm_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', 1)";
 
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 3;
 
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
    
 
108
    // Compute psitilde_a
 
109
    const double psitilde_a_0 = 1;
 
110
    const double psitilde_a_1 = x;
 
111
    
 
112
    // Compute psitilde_bs
 
113
    const double psitilde_bs_0_0 = 1;
 
114
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
115
    const double psitilde_bs_1_0 = 1;
 
116
    
 
117
    // Compute basisvalues
 
118
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
119
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
120
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
121
    
 
122
    // Table(s) of coefficients
 
123
    const static double coefficients0[3][3] = \
 
124
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
125
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
126
    {0.471404520791032, 0, 0.333333333333333}};
 
127
    
 
128
    // Extract relevant coefficients
 
129
    const double coeff0_0 = coefficients0[dof][0];
 
130
    const double coeff0_1 = coefficients0[dof][1];
 
131
    const double coeff0_2 = coefficients0[dof][2];
 
132
    
 
133
    // Compute value(s)
 
134
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
135
  }
 
136
 
 
137
  /// Evaluate all basis functions at given point in cell
 
138
  virtual void evaluate_basis_all(double* values,
 
139
                                  const double* coordinates,
 
140
                                  const ufc::cell& c) const
 
141
  {
 
142
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
143
  }
 
144
 
 
145
  /// Evaluate order n derivatives of basis function i at given point in cell
 
146
  virtual void evaluate_basis_derivatives(unsigned int i,
 
147
                                          unsigned int n,
 
148
                                          double* values,
 
149
                                          const double* coordinates,
 
150
                                          const ufc::cell& c) const
 
151
  {
 
152
    // Extract vertex coordinates
 
153
    const double * const * element_coordinates = c.coordinates;
 
154
    
 
155
    // Compute Jacobian of affine map from reference cell
 
156
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
157
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
158
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
159
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
160
      
 
161
    // Compute determinant of Jacobian
 
162
    const double detJ = J_00*J_11 - J_01*J_10;
 
163
    
 
164
    // Compute inverse of Jacobian
 
165
    
 
166
    // Get coordinates and map to the reference (UFC) element
 
167
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
168
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
169
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
170
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
171
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
172
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
173
    
 
174
    // Map coordinates to the reference square
 
175
    if (std::abs(y - 1.0) < 1e-14)
 
176
      x = -1.0;
 
177
    else
 
178
      x = 2.0 *x/(1.0 - y) - 1.0;
 
179
    y = 2.0*y - 1.0;
 
180
    
 
181
    // Compute number of derivatives
 
182
    unsigned int num_derivatives = 1;
 
183
    
 
184
    for (unsigned int j = 0; j < n; j++)
 
185
      num_derivatives *= 2;
 
186
    
 
187
    
 
188
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
189
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
190
        
 
191
    for (unsigned int j = 0; j < num_derivatives; j++)
 
192
    {
 
193
      combinations[j] = new unsigned int [n];
 
194
      for (unsigned int k = 0; k < n; k++)
 
195
        combinations[j][k] = 0;
 
196
    }
 
197
        
 
198
    // Generate combinations of derivatives
 
199
    for (unsigned int row = 1; row < num_derivatives; row++)
 
200
    {
 
201
      for (unsigned int num = 0; num < row; num++)
 
202
      {
 
203
        for (unsigned int col = n-1; col+1 > 0; col--)
 
204
        {
 
205
          if (combinations[row][col] + 1 > 1)
 
206
            combinations[row][col] = 0;
 
207
          else
 
208
          {
 
209
            combinations[row][col] += 1;
 
210
            break;
 
211
          }
 
212
        }
 
213
      }
 
214
    }
 
215
    
 
216
    // Compute inverse of Jacobian
 
217
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
218
    
 
219
    // Declare transformation matrix
 
220
    // Declare pointer to two dimensional array and initialise
 
221
    double **transform = new double *[num_derivatives];
 
222
        
 
223
    for (unsigned int j = 0; j < num_derivatives; j++)
 
224
    {
 
225
      transform[j] = new double [num_derivatives];
 
226
      for (unsigned int k = 0; k < num_derivatives; k++)
 
227
        transform[j][k] = 1;
 
228
    }
 
229
    
 
230
    // Construct transformation matrix
 
231
    for (unsigned int row = 0; row < num_derivatives; row++)
 
232
    {
 
233
      for (unsigned int col = 0; col < num_derivatives; col++)
 
234
      {
 
235
        for (unsigned int k = 0; k < n; k++)
 
236
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
237
      }
 
238
    }
 
239
    
 
240
    // Reset values
 
241
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
242
      values[j] = 0;
 
243
    
 
244
    // Map degree of freedom to element degree of freedom
 
245
    const unsigned int dof = i;
 
246
    
 
247
    // Generate scalings
 
248
    const double scalings_y_0 = 1;
 
249
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
250
    
 
251
    // Compute psitilde_a
 
252
    const double psitilde_a_0 = 1;
 
253
    const double psitilde_a_1 = x;
 
254
    
 
255
    // Compute psitilde_bs
 
256
    const double psitilde_bs_0_0 = 1;
 
257
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
258
    const double psitilde_bs_1_0 = 1;
 
259
    
 
260
    // Compute basisvalues
 
261
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
262
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
263
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
264
    
 
265
    // Table(s) of coefficients
 
266
    const static double coefficients0[3][3] = \
 
267
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
268
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
269
    {0.471404520791032, 0, 0.333333333333333}};
 
270
    
 
271
    // Interesting (new) part
 
272
    // Tables of derivatives of the polynomial base (transpose)
 
273
    const static double dmats0[3][3] = \
 
274
    {{0, 0, 0},
 
275
    {4.89897948556636, 0, 0},
 
276
    {0, 0, 0}};
 
277
    
 
278
    const static double dmats1[3][3] = \
 
279
    {{0, 0, 0},
 
280
    {2.44948974278318, 0, 0},
 
281
    {4.24264068711928, 0, 0}};
 
282
    
 
283
    // Compute reference derivatives
 
284
    // Declare pointer to array of derivatives on FIAT element
 
285
    double *derivatives = new double [num_derivatives];
 
286
    
 
287
    // Declare coefficients
 
288
    double coeff0_0 = 0;
 
289
    double coeff0_1 = 0;
 
290
    double coeff0_2 = 0;
 
291
    
 
292
    // Declare new coefficients
 
293
    double new_coeff0_0 = 0;
 
294
    double new_coeff0_1 = 0;
 
295
    double new_coeff0_2 = 0;
 
296
    
 
297
    // Loop possible derivatives
 
298
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
299
    {
 
300
      // Get values from coefficients array
 
301
      new_coeff0_0 = coefficients0[dof][0];
 
302
      new_coeff0_1 = coefficients0[dof][1];
 
303
      new_coeff0_2 = coefficients0[dof][2];
 
304
    
 
305
      // Loop derivative order
 
306
      for (unsigned int j = 0; j < n; j++)
 
307
      {
 
308
        // Update old coefficients
 
309
        coeff0_0 = new_coeff0_0;
 
310
        coeff0_1 = new_coeff0_1;
 
311
        coeff0_2 = new_coeff0_2;
 
312
    
 
313
        if(combinations[deriv_num][j] == 0)
 
314
        {
 
315
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
316
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
317
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
318
        }
 
319
        if(combinations[deriv_num][j] == 1)
 
320
        {
 
321
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
322
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
323
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
324
        }
 
325
    
 
326
      }
 
327
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
328
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
329
    }
 
330
    
 
331
    // Transform derivatives back to physical element
 
332
    for (unsigned int row = 0; row < num_derivatives; row++)
 
333
    {
 
334
      for (unsigned int col = 0; col < num_derivatives; col++)
 
335
      {
 
336
        values[row] += transform[row][col]*derivatives[col];
 
337
      }
 
338
    }
 
339
    // Delete pointer to array of derivatives on FIAT element
 
340
    delete [] derivatives;
 
341
    
 
342
    // Delete pointer to array of combinations of derivatives and transform
 
343
    for (unsigned int row = 0; row < num_derivatives; row++)
 
344
    {
 
345
      delete [] combinations[row];
 
346
      delete [] transform[row];
 
347
    }
 
348
    
 
349
    delete [] combinations;
 
350
    delete [] transform;
 
351
  }
 
352
 
 
353
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
354
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
355
                                              double* values,
 
356
                                              const double* coordinates,
 
357
                                              const ufc::cell& c) const
 
358
  {
 
359
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
360
  }
 
361
 
 
362
  /// Evaluate linear functional for dof i on the function f
 
363
  virtual double evaluate_dof(unsigned int i,
 
364
                              const ufc::function& f,
 
365
                              const ufc::cell& c) const
 
366
  {
 
367
    // The reference points, direction and weights:
 
368
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
369
    const static double W[3][1] = {{1}, {1}, {1}};
 
370
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
371
    
 
372
    const double * const * x = c.coordinates;
 
373
    double result = 0.0;
 
374
    // Iterate over the points:
 
375
    // Evaluate basis functions for affine mapping
 
376
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
377
    const double w1 = X[i][0][0];
 
378
    const double w2 = X[i][0][1];
 
379
    
 
380
    // Compute affine mapping y = F(X)
 
381
    double y[2];
 
382
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
383
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
384
    
 
385
    // Evaluate function at physical points
 
386
    double values[1];
 
387
    f.evaluate(values, y, c);
 
388
    
 
389
    // Map function values using appropriate mapping
 
390
    // Affine map: Do nothing
 
391
    
 
392
    // Note that we do not map the weights (yet).
 
393
    
 
394
    // Take directional components
 
395
    for(int k = 0; k < 1; k++)
 
396
      result += values[k]*D[i][0][k];
 
397
    // Multiply by weights 
 
398
    result *= W[i][0];
 
399
    
 
400
    return result;
 
401
  }
 
402
 
 
403
  /// Evaluate linear functionals for all dofs on the function f
 
404
  virtual void evaluate_dofs(double* values,
 
405
                             const ufc::function& f,
 
406
                             const ufc::cell& c) const
 
407
  {
 
408
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
409
  }
 
410
 
 
411
  /// Interpolate vertex values from dof values
 
412
  virtual void interpolate_vertex_values(double* vertex_values,
 
413
                                         const double* dof_values,
 
414
                                         const ufc::cell& c) const
 
415
  {
 
416
    // Evaluate at vertices and use affine mapping
 
417
    vertex_values[0] = dof_values[0];
 
418
    vertex_values[1] = dof_values[1];
 
419
    vertex_values[2] = dof_values[2];
 
420
  }
 
421
 
 
422
  /// Return the number of sub elements (for a mixed element)
 
423
  virtual unsigned int num_sub_elements() const
 
424
  {
 
425
    return 1;
 
426
  }
 
427
 
 
428
  /// Create a new finite element for sub element i (for a mixed element)
 
429
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
430
  {
 
431
    return new UFC_P1ProjectionBilinearForm_finite_element_0();
 
432
  }
 
433
 
 
434
};
 
435
 
 
436
/// This class defines the interface for a finite element.
 
437
 
 
438
class UFC_P1ProjectionBilinearForm_finite_element_1: public ufc::finite_element
 
439
{
 
440
public:
 
441
 
 
442
  /// Constructor
 
443
  UFC_P1ProjectionBilinearForm_finite_element_1() : ufc::finite_element()
 
444
  {
 
445
    // Do nothing
 
446
  }
 
447
 
 
448
  /// Destructor
 
449
  virtual ~UFC_P1ProjectionBilinearForm_finite_element_1()
 
450
  {
 
451
    // Do nothing
 
452
  }
 
453
 
 
454
  /// Return a string identifying the finite element
 
455
  virtual const char* signature() const
 
456
  {
 
457
    return "FiniteElement('Lagrange', 'triangle', 1)";
 
458
  }
 
459
 
 
460
  /// Return the cell shape
 
461
  virtual ufc::shape cell_shape() const
 
462
  {
 
463
    return ufc::triangle;
 
464
  }
 
465
 
 
466
  /// Return the dimension of the finite element function space
 
467
  virtual unsigned int space_dimension() const
 
468
  {
 
469
    return 3;
 
470
  }
 
471
 
 
472
  /// Return the rank of the value space
 
473
  virtual unsigned int value_rank() const
 
474
  {
 
475
    return 0;
 
476
  }
 
477
 
 
478
  /// Return the dimension of the value space for axis i
 
479
  virtual unsigned int value_dimension(unsigned int i) const
 
480
  {
 
481
    return 1;
 
482
  }
 
483
 
 
484
  /// Evaluate basis function i at given point in cell
 
485
  virtual void evaluate_basis(unsigned int i,
 
486
                              double* values,
 
487
                              const double* coordinates,
 
488
                              const ufc::cell& c) const
 
489
  {
 
490
    // Extract vertex coordinates
 
491
    const double * const * element_coordinates = c.coordinates;
 
492
    
 
493
    // Compute Jacobian of affine map from reference cell
 
494
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
495
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
496
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
497
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
498
      
 
499
    // Compute determinant of Jacobian
 
500
    const double detJ = J_00*J_11 - J_01*J_10;
 
501
    
 
502
    // Compute inverse of Jacobian
 
503
    
 
504
    // Get coordinates and map to the reference (UFC) element
 
505
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
506
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
507
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
508
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
509
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
510
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
511
    
 
512
    // Map coordinates to the reference square
 
513
    if (std::abs(y - 1.0) < 1e-14)
 
514
      x = -1.0;
 
515
    else
 
516
      x = 2.0 *x/(1.0 - y) - 1.0;
 
517
    y = 2.0*y - 1.0;
 
518
    
 
519
    // Reset values
 
520
    *values = 0;
 
521
    
 
522
    // Map degree of freedom to element degree of freedom
 
523
    const unsigned int dof = i;
 
524
    
 
525
    // Generate scalings
 
526
    const double scalings_y_0 = 1;
 
527
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
528
    
 
529
    // Compute psitilde_a
 
530
    const double psitilde_a_0 = 1;
 
531
    const double psitilde_a_1 = x;
 
532
    
 
533
    // Compute psitilde_bs
 
534
    const double psitilde_bs_0_0 = 1;
 
535
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
536
    const double psitilde_bs_1_0 = 1;
 
537
    
 
538
    // Compute basisvalues
 
539
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
540
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
541
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
542
    
 
543
    // Table(s) of coefficients
 
544
    const static double coefficients0[3][3] = \
 
545
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
546
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
547
    {0.471404520791032, 0, 0.333333333333333}};
 
548
    
 
549
    // Extract relevant coefficients
 
550
    const double coeff0_0 = coefficients0[dof][0];
 
551
    const double coeff0_1 = coefficients0[dof][1];
 
552
    const double coeff0_2 = coefficients0[dof][2];
 
553
    
 
554
    // Compute value(s)
 
555
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
556
  }
 
557
 
 
558
  /// Evaluate all basis functions at given point in cell
 
559
  virtual void evaluate_basis_all(double* values,
 
560
                                  const double* coordinates,
 
561
                                  const ufc::cell& c) const
 
562
  {
 
563
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
564
  }
 
565
 
 
566
  /// Evaluate order n derivatives of basis function i at given point in cell
 
567
  virtual void evaluate_basis_derivatives(unsigned int i,
 
568
                                          unsigned int n,
 
569
                                          double* values,
 
570
                                          const double* coordinates,
 
571
                                          const ufc::cell& c) const
 
572
  {
 
573
    // Extract vertex coordinates
 
574
    const double * const * element_coordinates = c.coordinates;
 
575
    
 
576
    // Compute Jacobian of affine map from reference cell
 
577
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
578
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
579
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
580
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
581
      
 
582
    // Compute determinant of Jacobian
 
583
    const double detJ = J_00*J_11 - J_01*J_10;
 
584
    
 
585
    // Compute inverse of Jacobian
 
586
    
 
587
    // Get coordinates and map to the reference (UFC) element
 
588
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
589
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
590
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
591
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
592
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
593
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
594
    
 
595
    // Map coordinates to the reference square
 
596
    if (std::abs(y - 1.0) < 1e-14)
 
597
      x = -1.0;
 
598
    else
 
599
      x = 2.0 *x/(1.0 - y) - 1.0;
 
600
    y = 2.0*y - 1.0;
 
601
    
 
602
    // Compute number of derivatives
 
603
    unsigned int num_derivatives = 1;
 
604
    
 
605
    for (unsigned int j = 0; j < n; j++)
 
606
      num_derivatives *= 2;
 
607
    
 
608
    
 
609
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
610
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
611
        
 
612
    for (unsigned int j = 0; j < num_derivatives; j++)
 
613
    {
 
614
      combinations[j] = new unsigned int [n];
 
615
      for (unsigned int k = 0; k < n; k++)
 
616
        combinations[j][k] = 0;
 
617
    }
 
618
        
 
619
    // Generate combinations of derivatives
 
620
    for (unsigned int row = 1; row < num_derivatives; row++)
 
621
    {
 
622
      for (unsigned int num = 0; num < row; num++)
 
623
      {
 
624
        for (unsigned int col = n-1; col+1 > 0; col--)
 
625
        {
 
626
          if (combinations[row][col] + 1 > 1)
 
627
            combinations[row][col] = 0;
 
628
          else
 
629
          {
 
630
            combinations[row][col] += 1;
 
631
            break;
 
632
          }
 
633
        }
 
634
      }
 
635
    }
 
636
    
 
637
    // Compute inverse of Jacobian
 
638
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
639
    
 
640
    // Declare transformation matrix
 
641
    // Declare pointer to two dimensional array and initialise
 
642
    double **transform = new double *[num_derivatives];
 
643
        
 
644
    for (unsigned int j = 0; j < num_derivatives; j++)
 
645
    {
 
646
      transform[j] = new double [num_derivatives];
 
647
      for (unsigned int k = 0; k < num_derivatives; k++)
 
648
        transform[j][k] = 1;
 
649
    }
 
650
    
 
651
    // Construct transformation matrix
 
652
    for (unsigned int row = 0; row < num_derivatives; row++)
 
653
    {
 
654
      for (unsigned int col = 0; col < num_derivatives; col++)
 
655
      {
 
656
        for (unsigned int k = 0; k < n; k++)
 
657
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
658
      }
 
659
    }
 
660
    
 
661
    // Reset values
 
662
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
663
      values[j] = 0;
 
664
    
 
665
    // Map degree of freedom to element degree of freedom
 
666
    const unsigned int dof = i;
 
667
    
 
668
    // Generate scalings
 
669
    const double scalings_y_0 = 1;
 
670
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
671
    
 
672
    // Compute psitilde_a
 
673
    const double psitilde_a_0 = 1;
 
674
    const double psitilde_a_1 = x;
 
675
    
 
676
    // Compute psitilde_bs
 
677
    const double psitilde_bs_0_0 = 1;
 
678
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
679
    const double psitilde_bs_1_0 = 1;
 
680
    
 
681
    // Compute basisvalues
 
682
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
683
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
684
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
685
    
 
686
    // Table(s) of coefficients
 
687
    const static double coefficients0[3][3] = \
 
688
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
689
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
690
    {0.471404520791032, 0, 0.333333333333333}};
 
691
    
 
692
    // Interesting (new) part
 
693
    // Tables of derivatives of the polynomial base (transpose)
 
694
    const static double dmats0[3][3] = \
 
695
    {{0, 0, 0},
 
696
    {4.89897948556636, 0, 0},
 
697
    {0, 0, 0}};
 
698
    
 
699
    const static double dmats1[3][3] = \
 
700
    {{0, 0, 0},
 
701
    {2.44948974278318, 0, 0},
 
702
    {4.24264068711928, 0, 0}};
 
703
    
 
704
    // Compute reference derivatives
 
705
    // Declare pointer to array of derivatives on FIAT element
 
706
    double *derivatives = new double [num_derivatives];
 
707
    
 
708
    // Declare coefficients
 
709
    double coeff0_0 = 0;
 
710
    double coeff0_1 = 0;
 
711
    double coeff0_2 = 0;
 
712
    
 
713
    // Declare new coefficients
 
714
    double new_coeff0_0 = 0;
 
715
    double new_coeff0_1 = 0;
 
716
    double new_coeff0_2 = 0;
 
717
    
 
718
    // Loop possible derivatives
 
719
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
720
    {
 
721
      // Get values from coefficients array
 
722
      new_coeff0_0 = coefficients0[dof][0];
 
723
      new_coeff0_1 = coefficients0[dof][1];
 
724
      new_coeff0_2 = coefficients0[dof][2];
 
725
    
 
726
      // Loop derivative order
 
727
      for (unsigned int j = 0; j < n; j++)
 
728
      {
 
729
        // Update old coefficients
 
730
        coeff0_0 = new_coeff0_0;
 
731
        coeff0_1 = new_coeff0_1;
 
732
        coeff0_2 = new_coeff0_2;
 
733
    
 
734
        if(combinations[deriv_num][j] == 0)
 
735
        {
 
736
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
737
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
738
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
739
        }
 
740
        if(combinations[deriv_num][j] == 1)
 
741
        {
 
742
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
743
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
744
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
745
        }
 
746
    
 
747
      }
 
748
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
749
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
750
    }
 
751
    
 
752
    // Transform derivatives back to physical element
 
753
    for (unsigned int row = 0; row < num_derivatives; row++)
 
754
    {
 
755
      for (unsigned int col = 0; col < num_derivatives; col++)
 
756
      {
 
757
        values[row] += transform[row][col]*derivatives[col];
 
758
      }
 
759
    }
 
760
    // Delete pointer to array of derivatives on FIAT element
 
761
    delete [] derivatives;
 
762
    
 
763
    // Delete pointer to array of combinations of derivatives and transform
 
764
    for (unsigned int row = 0; row < num_derivatives; row++)
 
765
    {
 
766
      delete [] combinations[row];
 
767
      delete [] transform[row];
 
768
    }
 
769
    
 
770
    delete [] combinations;
 
771
    delete [] transform;
 
772
  }
 
773
 
 
774
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
775
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
776
                                              double* values,
 
777
                                              const double* coordinates,
 
778
                                              const ufc::cell& c) const
 
779
  {
 
780
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
781
  }
 
782
 
 
783
  /// Evaluate linear functional for dof i on the function f
 
784
  virtual double evaluate_dof(unsigned int i,
 
785
                              const ufc::function& f,
 
786
                              const ufc::cell& c) const
 
787
  {
 
788
    // The reference points, direction and weights:
 
789
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
790
    const static double W[3][1] = {{1}, {1}, {1}};
 
791
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
792
    
 
793
    const double * const * x = c.coordinates;
 
794
    double result = 0.0;
 
795
    // Iterate over the points:
 
796
    // Evaluate basis functions for affine mapping
 
797
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
798
    const double w1 = X[i][0][0];
 
799
    const double w2 = X[i][0][1];
 
800
    
 
801
    // Compute affine mapping y = F(X)
 
802
    double y[2];
 
803
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
804
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
805
    
 
806
    // Evaluate function at physical points
 
807
    double values[1];
 
808
    f.evaluate(values, y, c);
 
809
    
 
810
    // Map function values using appropriate mapping
 
811
    // Affine map: Do nothing
 
812
    
 
813
    // Note that we do not map the weights (yet).
 
814
    
 
815
    // Take directional components
 
816
    for(int k = 0; k < 1; k++)
 
817
      result += values[k]*D[i][0][k];
 
818
    // Multiply by weights 
 
819
    result *= W[i][0];
 
820
    
 
821
    return result;
 
822
  }
 
823
 
 
824
  /// Evaluate linear functionals for all dofs on the function f
 
825
  virtual void evaluate_dofs(double* values,
 
826
                             const ufc::function& f,
 
827
                             const ufc::cell& c) const
 
828
  {
 
829
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
830
  }
 
831
 
 
832
  /// Interpolate vertex values from dof values
 
833
  virtual void interpolate_vertex_values(double* vertex_values,
 
834
                                         const double* dof_values,
 
835
                                         const ufc::cell& c) const
 
836
  {
 
837
    // Evaluate at vertices and use affine mapping
 
838
    vertex_values[0] = dof_values[0];
 
839
    vertex_values[1] = dof_values[1];
 
840
    vertex_values[2] = dof_values[2];
 
841
  }
 
842
 
 
843
  /// Return the number of sub elements (for a mixed element)
 
844
  virtual unsigned int num_sub_elements() const
 
845
  {
 
846
    return 1;
 
847
  }
 
848
 
 
849
  /// Create a new finite element for sub element i (for a mixed element)
 
850
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
851
  {
 
852
    return new UFC_P1ProjectionBilinearForm_finite_element_1();
 
853
  }
 
854
 
 
855
};
 
856
 
 
857
/// This class defines the interface for a local-to-global mapping of
 
858
/// degrees of freedom (dofs).
 
859
 
 
860
class UFC_P1ProjectionBilinearForm_dof_map_0: public ufc::dof_map
 
861
{
 
862
private:
 
863
 
 
864
  unsigned int __global_dimension;
 
865
 
 
866
public:
 
867
 
 
868
  /// Constructor
 
869
  UFC_P1ProjectionBilinearForm_dof_map_0() : ufc::dof_map()
 
870
  {
 
871
    __global_dimension = 0;
 
872
  }
 
873
 
 
874
  /// Destructor
 
875
  virtual ~UFC_P1ProjectionBilinearForm_dof_map_0()
 
876
  {
 
877
    // Do nothing
 
878
  }
 
879
 
 
880
  /// Return a string identifying the dof map
 
881
  virtual const char* signature() const
 
882
  {
 
883
    return "FFC dof map for FiniteElement('Lagrange', 'triangle', 1)";
 
884
  }
 
885
 
 
886
  /// Return true iff mesh entities of topological dimension d are needed
 
887
  virtual bool needs_mesh_entities(unsigned int d) const
 
888
  {
 
889
    switch ( d )
 
890
    {
 
891
    case 0:
 
892
      return true;
 
893
      break;
 
894
    case 1:
 
895
      return false;
 
896
      break;
 
897
    case 2:
 
898
      return false;
 
899
      break;
 
900
    }
 
901
    return false;
 
902
  }
 
903
 
 
904
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
905
  virtual bool init_mesh(const ufc::mesh& m)
 
906
  {
 
907
    __global_dimension = m.num_entities[0];
 
908
    return false;
 
909
  }
 
910
 
 
911
  /// Initialize dof map for given cell
 
912
  virtual void init_cell(const ufc::mesh& m,
 
913
                         const ufc::cell& c)
 
914
  {
 
915
    // Do nothing
 
916
  }
 
917
 
 
918
  /// Finish initialization of dof map for cells
 
919
  virtual void init_cell_finalize()
 
920
  {
 
921
    // Do nothing
 
922
  }
 
923
 
 
924
  /// Return the dimension of the global finite element function space
 
925
  virtual unsigned int global_dimension() const
 
926
  {
 
927
    return __global_dimension;
 
928
  }
 
929
 
 
930
  /// Return the dimension of the local finite element function space
 
931
  virtual unsigned int local_dimension() const
 
932
  {
 
933
    return 3;
 
934
  }
 
935
 
 
936
  // Return the geometric dimension of the coordinates this dof map provides
 
937
  virtual unsigned int geometric_dimension() const
 
938
  {
 
939
    return 2;
 
940
  }
 
941
 
 
942
  /// Return the number of dofs on each cell facet
 
943
  virtual unsigned int num_facet_dofs() const
 
944
  {
 
945
    return 2;
 
946
  }
 
947
 
 
948
  /// Return the number of dofs associated with each cell entity of dimension d
 
949
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
950
  {
 
951
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
952
  }
 
953
 
 
954
  /// Tabulate the local-to-global mapping of dofs on a cell
 
955
  virtual void tabulate_dofs(unsigned int* dofs,
 
956
                             const ufc::mesh& m,
 
957
                             const ufc::cell& c) const
 
958
  {
 
959
    dofs[0] = c.entity_indices[0][0];
 
960
    dofs[1] = c.entity_indices[0][1];
 
961
    dofs[2] = c.entity_indices[0][2];
 
962
  }
 
963
 
 
964
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
965
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
966
                                   unsigned int facet) const
 
967
  {
 
968
    switch ( facet )
 
969
    {
 
970
    case 0:
 
971
      dofs[0] = 1;
 
972
      dofs[1] = 2;
 
973
      break;
 
974
    case 1:
 
975
      dofs[0] = 0;
 
976
      dofs[1] = 2;
 
977
      break;
 
978
    case 2:
 
979
      dofs[0] = 0;
 
980
      dofs[1] = 1;
 
981
      break;
 
982
    }
 
983
  }
 
984
 
 
985
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
986
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
987
                                    unsigned int d, unsigned int i) const
 
988
  {
 
989
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
990
  }
 
991
 
 
992
  /// Tabulate the coordinates of all dofs on a cell
 
993
  virtual void tabulate_coordinates(double** coordinates,
 
994
                                    const ufc::cell& c) const
 
995
  {
 
996
    const double * const * x = c.coordinates;
 
997
    coordinates[0][0] = x[0][0];
 
998
    coordinates[0][1] = x[0][1];
 
999
    coordinates[1][0] = x[1][0];
 
1000
    coordinates[1][1] = x[1][1];
 
1001
    coordinates[2][0] = x[2][0];
 
1002
    coordinates[2][1] = x[2][1];
 
1003
  }
 
1004
 
 
1005
  /// Return the number of sub dof maps (for a mixed element)
 
1006
  virtual unsigned int num_sub_dof_maps() const
 
1007
  {
 
1008
    return 1;
 
1009
  }
 
1010
 
 
1011
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1012
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1013
  {
 
1014
    return new UFC_P1ProjectionBilinearForm_dof_map_0();
 
1015
  }
 
1016
 
 
1017
};
 
1018
 
 
1019
/// This class defines the interface for a local-to-global mapping of
 
1020
/// degrees of freedom (dofs).
 
1021
 
 
1022
class UFC_P1ProjectionBilinearForm_dof_map_1: public ufc::dof_map
 
1023
{
 
1024
private:
 
1025
 
 
1026
  unsigned int __global_dimension;
 
1027
 
 
1028
public:
 
1029
 
 
1030
  /// Constructor
 
1031
  UFC_P1ProjectionBilinearForm_dof_map_1() : ufc::dof_map()
 
1032
  {
 
1033
    __global_dimension = 0;
 
1034
  }
 
1035
 
 
1036
  /// Destructor
 
1037
  virtual ~UFC_P1ProjectionBilinearForm_dof_map_1()
 
1038
  {
 
1039
    // Do nothing
 
1040
  }
 
1041
 
 
1042
  /// Return a string identifying the dof map
 
1043
  virtual const char* signature() const
 
1044
  {
 
1045
    return "FFC dof map for FiniteElement('Lagrange', 'triangle', 1)";
 
1046
  }
 
1047
 
 
1048
  /// Return true iff mesh entities of topological dimension d are needed
 
1049
  virtual bool needs_mesh_entities(unsigned int d) const
 
1050
  {
 
1051
    switch ( d )
 
1052
    {
 
1053
    case 0:
 
1054
      return true;
 
1055
      break;
 
1056
    case 1:
 
1057
      return false;
 
1058
      break;
 
1059
    case 2:
 
1060
      return false;
 
1061
      break;
 
1062
    }
 
1063
    return false;
 
1064
  }
 
1065
 
 
1066
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1067
  virtual bool init_mesh(const ufc::mesh& m)
 
1068
  {
 
1069
    __global_dimension = m.num_entities[0];
 
1070
    return false;
 
1071
  }
 
1072
 
 
1073
  /// Initialize dof map for given cell
 
1074
  virtual void init_cell(const ufc::mesh& m,
 
1075
                         const ufc::cell& c)
 
1076
  {
 
1077
    // Do nothing
 
1078
  }
 
1079
 
 
1080
  /// Finish initialization of dof map for cells
 
1081
  virtual void init_cell_finalize()
 
1082
  {
 
1083
    // Do nothing
 
1084
  }
 
1085
 
 
1086
  /// Return the dimension of the global finite element function space
 
1087
  virtual unsigned int global_dimension() const
 
1088
  {
 
1089
    return __global_dimension;
 
1090
  }
 
1091
 
 
1092
  /// Return the dimension of the local finite element function space
 
1093
  virtual unsigned int local_dimension() const
 
1094
  {
 
1095
    return 3;
 
1096
  }
 
1097
 
 
1098
  // Return the geometric dimension of the coordinates this dof map provides
 
1099
  virtual unsigned int geometric_dimension() const
 
1100
  {
 
1101
    return 2;
 
1102
  }
 
1103
 
 
1104
  /// Return the number of dofs on each cell facet
 
1105
  virtual unsigned int num_facet_dofs() const
 
1106
  {
 
1107
    return 2;
 
1108
  }
 
1109
 
 
1110
  /// Return the number of dofs associated with each cell entity of dimension d
 
1111
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1112
  {
 
1113
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1114
  }
 
1115
 
 
1116
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1117
  virtual void tabulate_dofs(unsigned int* dofs,
 
1118
                             const ufc::mesh& m,
 
1119
                             const ufc::cell& c) const
 
1120
  {
 
1121
    dofs[0] = c.entity_indices[0][0];
 
1122
    dofs[1] = c.entity_indices[0][1];
 
1123
    dofs[2] = c.entity_indices[0][2];
 
1124
  }
 
1125
 
 
1126
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1127
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1128
                                   unsigned int facet) const
 
1129
  {
 
1130
    switch ( facet )
 
1131
    {
 
1132
    case 0:
 
1133
      dofs[0] = 1;
 
1134
      dofs[1] = 2;
 
1135
      break;
 
1136
    case 1:
 
1137
      dofs[0] = 0;
 
1138
      dofs[1] = 2;
 
1139
      break;
 
1140
    case 2:
 
1141
      dofs[0] = 0;
 
1142
      dofs[1] = 1;
 
1143
      break;
 
1144
    }
 
1145
  }
 
1146
 
 
1147
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1148
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1149
                                    unsigned int d, unsigned int i) const
 
1150
  {
 
1151
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1152
  }
 
1153
 
 
1154
  /// Tabulate the coordinates of all dofs on a cell
 
1155
  virtual void tabulate_coordinates(double** coordinates,
 
1156
                                    const ufc::cell& c) const
 
1157
  {
 
1158
    const double * const * x = c.coordinates;
 
1159
    coordinates[0][0] = x[0][0];
 
1160
    coordinates[0][1] = x[0][1];
 
1161
    coordinates[1][0] = x[1][0];
 
1162
    coordinates[1][1] = x[1][1];
 
1163
    coordinates[2][0] = x[2][0];
 
1164
    coordinates[2][1] = x[2][1];
 
1165
  }
 
1166
 
 
1167
  /// Return the number of sub dof maps (for a mixed element)
 
1168
  virtual unsigned int num_sub_dof_maps() const
 
1169
  {
 
1170
    return 1;
 
1171
  }
 
1172
 
 
1173
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1174
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1175
  {
 
1176
    return new UFC_P1ProjectionBilinearForm_dof_map_1();
 
1177
  }
 
1178
 
 
1179
};
 
1180
 
 
1181
/// This class defines the interface for the tabulation of the cell
 
1182
/// tensor corresponding to the local contribution to a form from
 
1183
/// the integral over a cell.
 
1184
 
 
1185
class UFC_P1ProjectionBilinearForm_cell_integral_0: public ufc::cell_integral
 
1186
{
 
1187
public:
 
1188
 
 
1189
  /// Constructor
 
1190
  UFC_P1ProjectionBilinearForm_cell_integral_0() : ufc::cell_integral()
 
1191
  {
 
1192
    // Do nothing
 
1193
  }
 
1194
 
 
1195
  /// Destructor
 
1196
  virtual ~UFC_P1ProjectionBilinearForm_cell_integral_0()
 
1197
  {
 
1198
    // Do nothing
 
1199
  }
 
1200
 
 
1201
  /// Tabulate the tensor for the contribution from a local cell
 
1202
  virtual void tabulate_tensor(double* A,
 
1203
                               const double * const * w,
 
1204
                               const ufc::cell& c) const
 
1205
  {
 
1206
    // Extract vertex coordinates
 
1207
    const double * const * x = c.coordinates;
 
1208
    
 
1209
    // Compute Jacobian of affine map from reference cell
 
1210
    const double J_00 = x[1][0] - x[0][0];
 
1211
    const double J_01 = x[2][0] - x[0][0];
 
1212
    const double J_10 = x[1][1] - x[0][1];
 
1213
    const double J_11 = x[2][1] - x[0][1];
 
1214
      
 
1215
    // Compute determinant of Jacobian
 
1216
    double detJ = J_00*J_11 - J_01*J_10;
 
1217
      
 
1218
    // Compute inverse of Jacobian
 
1219
    
 
1220
    // Set scale factor
 
1221
    const double det = std::abs(detJ);
 
1222
    
 
1223
    // Number of operations to compute element tensor = 10
 
1224
    // Compute geometry tensors
 
1225
    // Number of operations to compute decalrations = 1
 
1226
    const double G0_ = det;
 
1227
    
 
1228
    // Compute element tensor
 
1229
    // Number of operations to compute tensor = 9
 
1230
    A[0] = 0.0833333333333332*G0_;
 
1231
    A[1] = 0.0416666666666666*G0_;
 
1232
    A[2] = 0.0416666666666666*G0_;
 
1233
    A[3] = 0.0416666666666666*G0_;
 
1234
    A[4] = 0.0833333333333332*G0_;
 
1235
    A[5] = 0.0416666666666666*G0_;
 
1236
    A[6] = 0.0416666666666666*G0_;
 
1237
    A[7] = 0.0416666666666666*G0_;
 
1238
    A[8] = 0.0833333333333332*G0_;
 
1239
  }
 
1240
 
 
1241
};
 
1242
 
 
1243
/// This class defines the interface for the assembly of the global
 
1244
/// tensor corresponding to a form with r + n arguments, that is, a
 
1245
/// mapping
 
1246
///
 
1247
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
1248
///
 
1249
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
1250
/// global tensor A is defined by
 
1251
///
 
1252
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
1253
///
 
1254
/// where each argument Vj represents the application to the
 
1255
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
1256
/// fixed functions (coefficients).
 
1257
 
 
1258
class UFC_P1ProjectionBilinearForm: public ufc::form
 
1259
{
 
1260
public:
 
1261
 
 
1262
  /// Constructor
 
1263
  UFC_P1ProjectionBilinearForm() : ufc::form()
 
1264
  {
 
1265
    // Do nothing
 
1266
  }
 
1267
 
 
1268
  /// Destructor
 
1269
  virtual ~UFC_P1ProjectionBilinearForm()
 
1270
  {
 
1271
    // Do nothing
 
1272
  }
 
1273
 
 
1274
  /// Return a string identifying the form
 
1275
  virtual const char* signature() const
 
1276
  {
 
1277
    return " | vi0[0, 1, 2]*vi1[0, 1, 2]*dX(0)";
 
1278
  }
 
1279
 
 
1280
  /// Return the rank of the global tensor (r)
 
1281
  virtual unsigned int rank() const
 
1282
  {
 
1283
    return 2;
 
1284
  }
 
1285
 
 
1286
  /// Return the number of coefficients (n)
 
1287
  virtual unsigned int num_coefficients() const
 
1288
  {
 
1289
    return 0;
 
1290
  }
 
1291
 
 
1292
  /// Return the number of cell integrals
 
1293
  virtual unsigned int num_cell_integrals() const
 
1294
  {
 
1295
    return 1;
 
1296
  }
 
1297
  
 
1298
  /// Return the number of exterior facet integrals
 
1299
  virtual unsigned int num_exterior_facet_integrals() const
 
1300
  {
 
1301
    return 0;
 
1302
  }
 
1303
  
 
1304
  /// Return the number of interior facet integrals
 
1305
  virtual unsigned int num_interior_facet_integrals() const
 
1306
  {
 
1307
    return 0;
 
1308
  }
 
1309
    
 
1310
  /// Create a new finite element for argument function i
 
1311
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
1312
  {
 
1313
    switch ( i )
 
1314
    {
 
1315
    case 0:
 
1316
      return new UFC_P1ProjectionBilinearForm_finite_element_0();
 
1317
      break;
 
1318
    case 1:
 
1319
      return new UFC_P1ProjectionBilinearForm_finite_element_1();
 
1320
      break;
 
1321
    }
 
1322
    return 0;
 
1323
  }
 
1324
  
 
1325
  /// Create a new dof map for argument function i
 
1326
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
1327
  {
 
1328
    switch ( i )
 
1329
    {
 
1330
    case 0:
 
1331
      return new UFC_P1ProjectionBilinearForm_dof_map_0();
 
1332
      break;
 
1333
    case 1:
 
1334
      return new UFC_P1ProjectionBilinearForm_dof_map_1();
 
1335
      break;
 
1336
    }
 
1337
    return 0;
 
1338
  }
 
1339
 
 
1340
  /// Create a new cell integral on sub domain i
 
1341
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
1342
  {
 
1343
    return new UFC_P1ProjectionBilinearForm_cell_integral_0();
 
1344
  }
 
1345
 
 
1346
  /// Create a new exterior facet integral on sub domain i
 
1347
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
1348
  {
 
1349
    return 0;
 
1350
  }
 
1351
 
 
1352
  /// Create a new interior facet integral on sub domain i
 
1353
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
1354
  {
 
1355
    return 0;
 
1356
  }
 
1357
 
 
1358
};
 
1359
 
 
1360
/// This class defines the interface for a finite element.
 
1361
 
 
1362
class UFC_P1ProjectionLinearForm_finite_element_0: public ufc::finite_element
 
1363
{
 
1364
public:
 
1365
 
 
1366
  /// Constructor
 
1367
  UFC_P1ProjectionLinearForm_finite_element_0() : ufc::finite_element()
 
1368
  {
 
1369
    // Do nothing
 
1370
  }
 
1371
 
 
1372
  /// Destructor
 
1373
  virtual ~UFC_P1ProjectionLinearForm_finite_element_0()
 
1374
  {
 
1375
    // Do nothing
 
1376
  }
 
1377
 
 
1378
  /// Return a string identifying the finite element
 
1379
  virtual const char* signature() const
 
1380
  {
 
1381
    return "FiniteElement('Lagrange', 'triangle', 1)";
 
1382
  }
 
1383
 
 
1384
  /// Return the cell shape
 
1385
  virtual ufc::shape cell_shape() const
 
1386
  {
 
1387
    return ufc::triangle;
 
1388
  }
 
1389
 
 
1390
  /// Return the dimension of the finite element function space
 
1391
  virtual unsigned int space_dimension() const
 
1392
  {
 
1393
    return 3;
 
1394
  }
 
1395
 
 
1396
  /// Return the rank of the value space
 
1397
  virtual unsigned int value_rank() const
 
1398
  {
 
1399
    return 0;
 
1400
  }
 
1401
 
 
1402
  /// Return the dimension of the value space for axis i
 
1403
  virtual unsigned int value_dimension(unsigned int i) const
 
1404
  {
 
1405
    return 1;
 
1406
  }
 
1407
 
 
1408
  /// Evaluate basis function i at given point in cell
 
1409
  virtual void evaluate_basis(unsigned int i,
 
1410
                              double* values,
 
1411
                              const double* coordinates,
 
1412
                              const ufc::cell& c) const
 
1413
  {
 
1414
    // Extract vertex coordinates
 
1415
    const double * const * element_coordinates = c.coordinates;
 
1416
    
 
1417
    // Compute Jacobian of affine map from reference cell
 
1418
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
1419
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
1420
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
1421
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
1422
      
 
1423
    // Compute determinant of Jacobian
 
1424
    const double detJ = J_00*J_11 - J_01*J_10;
 
1425
    
 
1426
    // Compute inverse of Jacobian
 
1427
    
 
1428
    // Get coordinates and map to the reference (UFC) element
 
1429
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
1430
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
1431
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
1432
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
1433
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
1434
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
1435
    
 
1436
    // Map coordinates to the reference square
 
1437
    if (std::abs(y - 1.0) < 1e-14)
 
1438
      x = -1.0;
 
1439
    else
 
1440
      x = 2.0 *x/(1.0 - y) - 1.0;
 
1441
    y = 2.0*y - 1.0;
 
1442
    
 
1443
    // Reset values
 
1444
    *values = 0;
 
1445
    
 
1446
    // Map degree of freedom to element degree of freedom
 
1447
    const unsigned int dof = i;
 
1448
    
 
1449
    // Generate scalings
 
1450
    const double scalings_y_0 = 1;
 
1451
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
1452
    
 
1453
    // Compute psitilde_a
 
1454
    const double psitilde_a_0 = 1;
 
1455
    const double psitilde_a_1 = x;
 
1456
    
 
1457
    // Compute psitilde_bs
 
1458
    const double psitilde_bs_0_0 = 1;
 
1459
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
1460
    const double psitilde_bs_1_0 = 1;
 
1461
    
 
1462
    // Compute basisvalues
 
1463
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
1464
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
1465
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
1466
    
 
1467
    // Table(s) of coefficients
 
1468
    const static double coefficients0[3][3] = \
 
1469
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
1470
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
1471
    {0.471404520791032, 0, 0.333333333333333}};
 
1472
    
 
1473
    // Extract relevant coefficients
 
1474
    const double coeff0_0 = coefficients0[dof][0];
 
1475
    const double coeff0_1 = coefficients0[dof][1];
 
1476
    const double coeff0_2 = coefficients0[dof][2];
 
1477
    
 
1478
    // Compute value(s)
 
1479
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
1480
  }
 
1481
 
 
1482
  /// Evaluate all basis functions at given point in cell
 
1483
  virtual void evaluate_basis_all(double* values,
 
1484
                                  const double* coordinates,
 
1485
                                  const ufc::cell& c) const
 
1486
  {
 
1487
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
1488
  }
 
1489
 
 
1490
  /// Evaluate order n derivatives of basis function i at given point in cell
 
1491
  virtual void evaluate_basis_derivatives(unsigned int i,
 
1492
                                          unsigned int n,
 
1493
                                          double* values,
 
1494
                                          const double* coordinates,
 
1495
                                          const ufc::cell& c) const
 
1496
  {
 
1497
    // Extract vertex coordinates
 
1498
    const double * const * element_coordinates = c.coordinates;
 
1499
    
 
1500
    // Compute Jacobian of affine map from reference cell
 
1501
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
1502
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
1503
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
1504
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
1505
      
 
1506
    // Compute determinant of Jacobian
 
1507
    const double detJ = J_00*J_11 - J_01*J_10;
 
1508
    
 
1509
    // Compute inverse of Jacobian
 
1510
    
 
1511
    // Get coordinates and map to the reference (UFC) element
 
1512
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
1513
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
1514
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
1515
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
1516
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
1517
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
1518
    
 
1519
    // Map coordinates to the reference square
 
1520
    if (std::abs(y - 1.0) < 1e-14)
 
1521
      x = -1.0;
 
1522
    else
 
1523
      x = 2.0 *x/(1.0 - y) - 1.0;
 
1524
    y = 2.0*y - 1.0;
 
1525
    
 
1526
    // Compute number of derivatives
 
1527
    unsigned int num_derivatives = 1;
 
1528
    
 
1529
    for (unsigned int j = 0; j < n; j++)
 
1530
      num_derivatives *= 2;
 
1531
    
 
1532
    
 
1533
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
1534
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
1535
        
 
1536
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1537
    {
 
1538
      combinations[j] = new unsigned int [n];
 
1539
      for (unsigned int k = 0; k < n; k++)
 
1540
        combinations[j][k] = 0;
 
1541
    }
 
1542
        
 
1543
    // Generate combinations of derivatives
 
1544
    for (unsigned int row = 1; row < num_derivatives; row++)
 
1545
    {
 
1546
      for (unsigned int num = 0; num < row; num++)
 
1547
      {
 
1548
        for (unsigned int col = n-1; col+1 > 0; col--)
 
1549
        {
 
1550
          if (combinations[row][col] + 1 > 1)
 
1551
            combinations[row][col] = 0;
 
1552
          else
 
1553
          {
 
1554
            combinations[row][col] += 1;
 
1555
            break;
 
1556
          }
 
1557
        }
 
1558
      }
 
1559
    }
 
1560
    
 
1561
    // Compute inverse of Jacobian
 
1562
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
1563
    
 
1564
    // Declare transformation matrix
 
1565
    // Declare pointer to two dimensional array and initialise
 
1566
    double **transform = new double *[num_derivatives];
 
1567
        
 
1568
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1569
    {
 
1570
      transform[j] = new double [num_derivatives];
 
1571
      for (unsigned int k = 0; k < num_derivatives; k++)
 
1572
        transform[j][k] = 1;
 
1573
    }
 
1574
    
 
1575
    // Construct transformation matrix
 
1576
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1577
    {
 
1578
      for (unsigned int col = 0; col < num_derivatives; col++)
 
1579
      {
 
1580
        for (unsigned int k = 0; k < n; k++)
 
1581
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
1582
      }
 
1583
    }
 
1584
    
 
1585
    // Reset values
 
1586
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
1587
      values[j] = 0;
 
1588
    
 
1589
    // Map degree of freedom to element degree of freedom
 
1590
    const unsigned int dof = i;
 
1591
    
 
1592
    // Generate scalings
 
1593
    const double scalings_y_0 = 1;
 
1594
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
1595
    
 
1596
    // Compute psitilde_a
 
1597
    const double psitilde_a_0 = 1;
 
1598
    const double psitilde_a_1 = x;
 
1599
    
 
1600
    // Compute psitilde_bs
 
1601
    const double psitilde_bs_0_0 = 1;
 
1602
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
1603
    const double psitilde_bs_1_0 = 1;
 
1604
    
 
1605
    // Compute basisvalues
 
1606
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
1607
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
1608
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
1609
    
 
1610
    // Table(s) of coefficients
 
1611
    const static double coefficients0[3][3] = \
 
1612
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
1613
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
1614
    {0.471404520791032, 0, 0.333333333333333}};
 
1615
    
 
1616
    // Interesting (new) part
 
1617
    // Tables of derivatives of the polynomial base (transpose)
 
1618
    const static double dmats0[3][3] = \
 
1619
    {{0, 0, 0},
 
1620
    {4.89897948556636, 0, 0},
 
1621
    {0, 0, 0}};
 
1622
    
 
1623
    const static double dmats1[3][3] = \
 
1624
    {{0, 0, 0},
 
1625
    {2.44948974278318, 0, 0},
 
1626
    {4.24264068711928, 0, 0}};
 
1627
    
 
1628
    // Compute reference derivatives
 
1629
    // Declare pointer to array of derivatives on FIAT element
 
1630
    double *derivatives = new double [num_derivatives];
 
1631
    
 
1632
    // Declare coefficients
 
1633
    double coeff0_0 = 0;
 
1634
    double coeff0_1 = 0;
 
1635
    double coeff0_2 = 0;
 
1636
    
 
1637
    // Declare new coefficients
 
1638
    double new_coeff0_0 = 0;
 
1639
    double new_coeff0_1 = 0;
 
1640
    double new_coeff0_2 = 0;
 
1641
    
 
1642
    // Loop possible derivatives
 
1643
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
1644
    {
 
1645
      // Get values from coefficients array
 
1646
      new_coeff0_0 = coefficients0[dof][0];
 
1647
      new_coeff0_1 = coefficients0[dof][1];
 
1648
      new_coeff0_2 = coefficients0[dof][2];
 
1649
    
 
1650
      // Loop derivative order
 
1651
      for (unsigned int j = 0; j < n; j++)
 
1652
      {
 
1653
        // Update old coefficients
 
1654
        coeff0_0 = new_coeff0_0;
 
1655
        coeff0_1 = new_coeff0_1;
 
1656
        coeff0_2 = new_coeff0_2;
 
1657
    
 
1658
        if(combinations[deriv_num][j] == 0)
 
1659
        {
 
1660
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
1661
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
1662
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
1663
        }
 
1664
        if(combinations[deriv_num][j] == 1)
 
1665
        {
 
1666
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
1667
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
1668
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
1669
        }
 
1670
    
 
1671
      }
 
1672
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
1673
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
1674
    }
 
1675
    
 
1676
    // Transform derivatives back to physical element
 
1677
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1678
    {
 
1679
      for (unsigned int col = 0; col < num_derivatives; col++)
 
1680
      {
 
1681
        values[row] += transform[row][col]*derivatives[col];
 
1682
      }
 
1683
    }
 
1684
    // Delete pointer to array of derivatives on FIAT element
 
1685
    delete [] derivatives;
 
1686
    
 
1687
    // Delete pointer to array of combinations of derivatives and transform
 
1688
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1689
    {
 
1690
      delete [] combinations[row];
 
1691
      delete [] transform[row];
 
1692
    }
 
1693
    
 
1694
    delete [] combinations;
 
1695
    delete [] transform;
 
1696
  }
 
1697
 
 
1698
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
1699
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
1700
                                              double* values,
 
1701
                                              const double* coordinates,
 
1702
                                              const ufc::cell& c) const
 
1703
  {
 
1704
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
1705
  }
 
1706
 
 
1707
  /// Evaluate linear functional for dof i on the function f
 
1708
  virtual double evaluate_dof(unsigned int i,
 
1709
                              const ufc::function& f,
 
1710
                              const ufc::cell& c) const
 
1711
  {
 
1712
    // The reference points, direction and weights:
 
1713
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
1714
    const static double W[3][1] = {{1}, {1}, {1}};
 
1715
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
1716
    
 
1717
    const double * const * x = c.coordinates;
 
1718
    double result = 0.0;
 
1719
    // Iterate over the points:
 
1720
    // Evaluate basis functions for affine mapping
 
1721
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
1722
    const double w1 = X[i][0][0];
 
1723
    const double w2 = X[i][0][1];
 
1724
    
 
1725
    // Compute affine mapping y = F(X)
 
1726
    double y[2];
 
1727
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
1728
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
1729
    
 
1730
    // Evaluate function at physical points
 
1731
    double values[1];
 
1732
    f.evaluate(values, y, c);
 
1733
    
 
1734
    // Map function values using appropriate mapping
 
1735
    // Affine map: Do nothing
 
1736
    
 
1737
    // Note that we do not map the weights (yet).
 
1738
    
 
1739
    // Take directional components
 
1740
    for(int k = 0; k < 1; k++)
 
1741
      result += values[k]*D[i][0][k];
 
1742
    // Multiply by weights 
 
1743
    result *= W[i][0];
 
1744
    
 
1745
    return result;
 
1746
  }
 
1747
 
 
1748
  /// Evaluate linear functionals for all dofs on the function f
 
1749
  virtual void evaluate_dofs(double* values,
 
1750
                             const ufc::function& f,
 
1751
                             const ufc::cell& c) const
 
1752
  {
 
1753
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1754
  }
 
1755
 
 
1756
  /// Interpolate vertex values from dof values
 
1757
  virtual void interpolate_vertex_values(double* vertex_values,
 
1758
                                         const double* dof_values,
 
1759
                                         const ufc::cell& c) const
 
1760
  {
 
1761
    // Evaluate at vertices and use affine mapping
 
1762
    vertex_values[0] = dof_values[0];
 
1763
    vertex_values[1] = dof_values[1];
 
1764
    vertex_values[2] = dof_values[2];
 
1765
  }
 
1766
 
 
1767
  /// Return the number of sub elements (for a mixed element)
 
1768
  virtual unsigned int num_sub_elements() const
 
1769
  {
 
1770
    return 1;
 
1771
  }
 
1772
 
 
1773
  /// Create a new finite element for sub element i (for a mixed element)
 
1774
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
1775
  {
 
1776
    return new UFC_P1ProjectionLinearForm_finite_element_0();
 
1777
  }
 
1778
 
 
1779
};
 
1780
 
 
1781
/// This class defines the interface for a finite element.
 
1782
 
 
1783
class UFC_P1ProjectionLinearForm_finite_element_1: public ufc::finite_element
 
1784
{
 
1785
public:
 
1786
 
 
1787
  /// Constructor
 
1788
  UFC_P1ProjectionLinearForm_finite_element_1() : ufc::finite_element()
 
1789
  {
 
1790
    // Do nothing
 
1791
  }
 
1792
 
 
1793
  /// Destructor
 
1794
  virtual ~UFC_P1ProjectionLinearForm_finite_element_1()
 
1795
  {
 
1796
    // Do nothing
 
1797
  }
 
1798
 
 
1799
  /// Return a string identifying the finite element
 
1800
  virtual const char* signature() const
 
1801
  {
 
1802
    return "FiniteElement('Discontinuous Lagrange', 'triangle', 1)";
 
1803
  }
 
1804
 
 
1805
  /// Return the cell shape
 
1806
  virtual ufc::shape cell_shape() const
 
1807
  {
 
1808
    return ufc::triangle;
 
1809
  }
 
1810
 
 
1811
  /// Return the dimension of the finite element function space
 
1812
  virtual unsigned int space_dimension() const
 
1813
  {
 
1814
    return 3;
 
1815
  }
 
1816
 
 
1817
  /// Return the rank of the value space
 
1818
  virtual unsigned int value_rank() const
 
1819
  {
 
1820
    return 0;
 
1821
  }
 
1822
 
 
1823
  /// Return the dimension of the value space for axis i
 
1824
  virtual unsigned int value_dimension(unsigned int i) const
 
1825
  {
 
1826
    return 1;
 
1827
  }
 
1828
 
 
1829
  /// Evaluate basis function i at given point in cell
 
1830
  virtual void evaluate_basis(unsigned int i,
 
1831
                              double* values,
 
1832
                              const double* coordinates,
 
1833
                              const ufc::cell& c) const
 
1834
  {
 
1835
    // Extract vertex coordinates
 
1836
    const double * const * element_coordinates = c.coordinates;
 
1837
    
 
1838
    // Compute Jacobian of affine map from reference cell
 
1839
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
1840
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
1841
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
1842
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
1843
      
 
1844
    // Compute determinant of Jacobian
 
1845
    const double detJ = J_00*J_11 - J_01*J_10;
 
1846
    
 
1847
    // Compute inverse of Jacobian
 
1848
    
 
1849
    // Get coordinates and map to the reference (UFC) element
 
1850
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
1851
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
1852
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
1853
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
1854
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
1855
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
1856
    
 
1857
    // Map coordinates to the reference square
 
1858
    if (std::abs(y - 1.0) < 1e-14)
 
1859
      x = -1.0;
 
1860
    else
 
1861
      x = 2.0 *x/(1.0 - y) - 1.0;
 
1862
    y = 2.0*y - 1.0;
 
1863
    
 
1864
    // Reset values
 
1865
    *values = 0;
 
1866
    
 
1867
    // Map degree of freedom to element degree of freedom
 
1868
    const unsigned int dof = i;
 
1869
    
 
1870
    // Generate scalings
 
1871
    const double scalings_y_0 = 1;
 
1872
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
1873
    
 
1874
    // Compute psitilde_a
 
1875
    const double psitilde_a_0 = 1;
 
1876
    const double psitilde_a_1 = x;
 
1877
    
 
1878
    // Compute psitilde_bs
 
1879
    const double psitilde_bs_0_0 = 1;
 
1880
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
1881
    const double psitilde_bs_1_0 = 1;
 
1882
    
 
1883
    // Compute basisvalues
 
1884
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
1885
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
1886
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
1887
    
 
1888
    // Table(s) of coefficients
 
1889
    const static double coefficients0[3][3] = \
 
1890
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
1891
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
1892
    {0.471404520791032, 0, 0.333333333333333}};
 
1893
    
 
1894
    // Extract relevant coefficients
 
1895
    const double coeff0_0 = coefficients0[dof][0];
 
1896
    const double coeff0_1 = coefficients0[dof][1];
 
1897
    const double coeff0_2 = coefficients0[dof][2];
 
1898
    
 
1899
    // Compute value(s)
 
1900
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
1901
  }
 
1902
 
 
1903
  /// Evaluate all basis functions at given point in cell
 
1904
  virtual void evaluate_basis_all(double* values,
 
1905
                                  const double* coordinates,
 
1906
                                  const ufc::cell& c) const
 
1907
  {
 
1908
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
1909
  }
 
1910
 
 
1911
  /// Evaluate order n derivatives of basis function i at given point in cell
 
1912
  virtual void evaluate_basis_derivatives(unsigned int i,
 
1913
                                          unsigned int n,
 
1914
                                          double* values,
 
1915
                                          const double* coordinates,
 
1916
                                          const ufc::cell& c) const
 
1917
  {
 
1918
    // Extract vertex coordinates
 
1919
    const double * const * element_coordinates = c.coordinates;
 
1920
    
 
1921
    // Compute Jacobian of affine map from reference cell
 
1922
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
1923
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
1924
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
1925
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
1926
      
 
1927
    // Compute determinant of Jacobian
 
1928
    const double detJ = J_00*J_11 - J_01*J_10;
 
1929
    
 
1930
    // Compute inverse of Jacobian
 
1931
    
 
1932
    // Get coordinates and map to the reference (UFC) element
 
1933
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
1934
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
1935
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
1936
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
1937
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
1938
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
1939
    
 
1940
    // Map coordinates to the reference square
 
1941
    if (std::abs(y - 1.0) < 1e-14)
 
1942
      x = -1.0;
 
1943
    else
 
1944
      x = 2.0 *x/(1.0 - y) - 1.0;
 
1945
    y = 2.0*y - 1.0;
 
1946
    
 
1947
    // Compute number of derivatives
 
1948
    unsigned int num_derivatives = 1;
 
1949
    
 
1950
    for (unsigned int j = 0; j < n; j++)
 
1951
      num_derivatives *= 2;
 
1952
    
 
1953
    
 
1954
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
1955
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
1956
        
 
1957
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1958
    {
 
1959
      combinations[j] = new unsigned int [n];
 
1960
      for (unsigned int k = 0; k < n; k++)
 
1961
        combinations[j][k] = 0;
 
1962
    }
 
1963
        
 
1964
    // Generate combinations of derivatives
 
1965
    for (unsigned int row = 1; row < num_derivatives; row++)
 
1966
    {
 
1967
      for (unsigned int num = 0; num < row; num++)
 
1968
      {
 
1969
        for (unsigned int col = n-1; col+1 > 0; col--)
 
1970
        {
 
1971
          if (combinations[row][col] + 1 > 1)
 
1972
            combinations[row][col] = 0;
 
1973
          else
 
1974
          {
 
1975
            combinations[row][col] += 1;
 
1976
            break;
 
1977
          }
 
1978
        }
 
1979
      }
 
1980
    }
 
1981
    
 
1982
    // Compute inverse of Jacobian
 
1983
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
1984
    
 
1985
    // Declare transformation matrix
 
1986
    // Declare pointer to two dimensional array and initialise
 
1987
    double **transform = new double *[num_derivatives];
 
1988
        
 
1989
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1990
    {
 
1991
      transform[j] = new double [num_derivatives];
 
1992
      for (unsigned int k = 0; k < num_derivatives; k++)
 
1993
        transform[j][k] = 1;
 
1994
    }
 
1995
    
 
1996
    // Construct transformation matrix
 
1997
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1998
    {
 
1999
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2000
      {
 
2001
        for (unsigned int k = 0; k < n; k++)
 
2002
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
2003
      }
 
2004
    }
 
2005
    
 
2006
    // Reset values
 
2007
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
2008
      values[j] = 0;
 
2009
    
 
2010
    // Map degree of freedom to element degree of freedom
 
2011
    const unsigned int dof = i;
 
2012
    
 
2013
    // Generate scalings
 
2014
    const double scalings_y_0 = 1;
 
2015
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2016
    
 
2017
    // Compute psitilde_a
 
2018
    const double psitilde_a_0 = 1;
 
2019
    const double psitilde_a_1 = x;
 
2020
    
 
2021
    // Compute psitilde_bs
 
2022
    const double psitilde_bs_0_0 = 1;
 
2023
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2024
    const double psitilde_bs_1_0 = 1;
 
2025
    
 
2026
    // Compute basisvalues
 
2027
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
2028
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
2029
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
2030
    
 
2031
    // Table(s) of coefficients
 
2032
    const static double coefficients0[3][3] = \
 
2033
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
2034
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
2035
    {0.471404520791032, 0, 0.333333333333333}};
 
2036
    
 
2037
    // Interesting (new) part
 
2038
    // Tables of derivatives of the polynomial base (transpose)
 
2039
    const static double dmats0[3][3] = \
 
2040
    {{0, 0, 0},
 
2041
    {4.89897948556636, 0, 0},
 
2042
    {0, 0, 0}};
 
2043
    
 
2044
    const static double dmats1[3][3] = \
 
2045
    {{0, 0, 0},
 
2046
    {2.44948974278318, 0, 0},
 
2047
    {4.24264068711928, 0, 0}};
 
2048
    
 
2049
    // Compute reference derivatives
 
2050
    // Declare pointer to array of derivatives on FIAT element
 
2051
    double *derivatives = new double [num_derivatives];
 
2052
    
 
2053
    // Declare coefficients
 
2054
    double coeff0_0 = 0;
 
2055
    double coeff0_1 = 0;
 
2056
    double coeff0_2 = 0;
 
2057
    
 
2058
    // Declare new coefficients
 
2059
    double new_coeff0_0 = 0;
 
2060
    double new_coeff0_1 = 0;
 
2061
    double new_coeff0_2 = 0;
 
2062
    
 
2063
    // Loop possible derivatives
 
2064
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
2065
    {
 
2066
      // Get values from coefficients array
 
2067
      new_coeff0_0 = coefficients0[dof][0];
 
2068
      new_coeff0_1 = coefficients0[dof][1];
 
2069
      new_coeff0_2 = coefficients0[dof][2];
 
2070
    
 
2071
      // Loop derivative order
 
2072
      for (unsigned int j = 0; j < n; j++)
 
2073
      {
 
2074
        // Update old coefficients
 
2075
        coeff0_0 = new_coeff0_0;
 
2076
        coeff0_1 = new_coeff0_1;
 
2077
        coeff0_2 = new_coeff0_2;
 
2078
    
 
2079
        if(combinations[deriv_num][j] == 0)
 
2080
        {
 
2081
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
2082
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
2083
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
2084
        }
 
2085
        if(combinations[deriv_num][j] == 1)
 
2086
        {
 
2087
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
2088
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
2089
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
2090
        }
 
2091
    
 
2092
      }
 
2093
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
2094
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
2095
    }
 
2096
    
 
2097
    // Transform derivatives back to physical element
 
2098
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2099
    {
 
2100
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2101
      {
 
2102
        values[row] += transform[row][col]*derivatives[col];
 
2103
      }
 
2104
    }
 
2105
    // Delete pointer to array of derivatives on FIAT element
 
2106
    delete [] derivatives;
 
2107
    
 
2108
    // Delete pointer to array of combinations of derivatives and transform
 
2109
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2110
    {
 
2111
      delete [] combinations[row];
 
2112
      delete [] transform[row];
 
2113
    }
 
2114
    
 
2115
    delete [] combinations;
 
2116
    delete [] transform;
 
2117
  }
 
2118
 
 
2119
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
2120
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
2121
                                              double* values,
 
2122
                                              const double* coordinates,
 
2123
                                              const ufc::cell& c) const
 
2124
  {
 
2125
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
2126
  }
 
2127
 
 
2128
  /// Evaluate linear functional for dof i on the function f
 
2129
  virtual double evaluate_dof(unsigned int i,
 
2130
                              const ufc::function& f,
 
2131
                              const ufc::cell& c) const
 
2132
  {
 
2133
    // The reference points, direction and weights:
 
2134
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
2135
    const static double W[3][1] = {{1}, {1}, {1}};
 
2136
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
2137
    
 
2138
    const double * const * x = c.coordinates;
 
2139
    double result = 0.0;
 
2140
    // Iterate over the points:
 
2141
    // Evaluate basis functions for affine mapping
 
2142
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
2143
    const double w1 = X[i][0][0];
 
2144
    const double w2 = X[i][0][1];
 
2145
    
 
2146
    // Compute affine mapping y = F(X)
 
2147
    double y[2];
 
2148
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
2149
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
2150
    
 
2151
    // Evaluate function at physical points
 
2152
    double values[1];
 
2153
    f.evaluate(values, y, c);
 
2154
    
 
2155
    // Map function values using appropriate mapping
 
2156
    // Affine map: Do nothing
 
2157
    
 
2158
    // Note that we do not map the weights (yet).
 
2159
    
 
2160
    // Take directional components
 
2161
    for(int k = 0; k < 1; k++)
 
2162
      result += values[k]*D[i][0][k];
 
2163
    // Multiply by weights 
 
2164
    result *= W[i][0];
 
2165
    
 
2166
    return result;
 
2167
  }
 
2168
 
 
2169
  /// Evaluate linear functionals for all dofs on the function f
 
2170
  virtual void evaluate_dofs(double* values,
 
2171
                             const ufc::function& f,
 
2172
                             const ufc::cell& c) const
 
2173
  {
 
2174
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2175
  }
 
2176
 
 
2177
  /// Interpolate vertex values from dof values
 
2178
  virtual void interpolate_vertex_values(double* vertex_values,
 
2179
                                         const double* dof_values,
 
2180
                                         const ufc::cell& c) const
 
2181
  {
 
2182
    // Evaluate at vertices and use affine mapping
 
2183
    vertex_values[0] = dof_values[0];
 
2184
    vertex_values[1] = dof_values[1];
 
2185
    vertex_values[2] = dof_values[2];
 
2186
  }
 
2187
 
 
2188
  /// Return the number of sub elements (for a mixed element)
 
2189
  virtual unsigned int num_sub_elements() const
 
2190
  {
 
2191
    return 1;
 
2192
  }
 
2193
 
 
2194
  /// Create a new finite element for sub element i (for a mixed element)
 
2195
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
2196
  {
 
2197
    return new UFC_P1ProjectionLinearForm_finite_element_1();
 
2198
  }
 
2199
 
 
2200
};
 
2201
 
 
2202
/// This class defines the interface for a local-to-global mapping of
 
2203
/// degrees of freedom (dofs).
 
2204
 
 
2205
class UFC_P1ProjectionLinearForm_dof_map_0: public ufc::dof_map
 
2206
{
 
2207
private:
 
2208
 
 
2209
  unsigned int __global_dimension;
 
2210
 
 
2211
public:
 
2212
 
 
2213
  /// Constructor
 
2214
  UFC_P1ProjectionLinearForm_dof_map_0() : ufc::dof_map()
 
2215
  {
 
2216
    __global_dimension = 0;
 
2217
  }
 
2218
 
 
2219
  /// Destructor
 
2220
  virtual ~UFC_P1ProjectionLinearForm_dof_map_0()
 
2221
  {
 
2222
    // Do nothing
 
2223
  }
 
2224
 
 
2225
  /// Return a string identifying the dof map
 
2226
  virtual const char* signature() const
 
2227
  {
 
2228
    return "FFC dof map for FiniteElement('Lagrange', 'triangle', 1)";
 
2229
  }
 
2230
 
 
2231
  /// Return true iff mesh entities of topological dimension d are needed
 
2232
  virtual bool needs_mesh_entities(unsigned int d) const
 
2233
  {
 
2234
    switch ( d )
 
2235
    {
 
2236
    case 0:
 
2237
      return true;
 
2238
      break;
 
2239
    case 1:
 
2240
      return false;
 
2241
      break;
 
2242
    case 2:
 
2243
      return false;
 
2244
      break;
 
2245
    }
 
2246
    return false;
 
2247
  }
 
2248
 
 
2249
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
2250
  virtual bool init_mesh(const ufc::mesh& m)
 
2251
  {
 
2252
    __global_dimension = m.num_entities[0];
 
2253
    return false;
 
2254
  }
 
2255
 
 
2256
  /// Initialize dof map for given cell
 
2257
  virtual void init_cell(const ufc::mesh& m,
 
2258
                         const ufc::cell& c)
 
2259
  {
 
2260
    // Do nothing
 
2261
  }
 
2262
 
 
2263
  /// Finish initialization of dof map for cells
 
2264
  virtual void init_cell_finalize()
 
2265
  {
 
2266
    // Do nothing
 
2267
  }
 
2268
 
 
2269
  /// Return the dimension of the global finite element function space
 
2270
  virtual unsigned int global_dimension() const
 
2271
  {
 
2272
    return __global_dimension;
 
2273
  }
 
2274
 
 
2275
  /// Return the dimension of the local finite element function space
 
2276
  virtual unsigned int local_dimension() const
 
2277
  {
 
2278
    return 3;
 
2279
  }
 
2280
 
 
2281
  // Return the geometric dimension of the coordinates this dof map provides
 
2282
  virtual unsigned int geometric_dimension() const
 
2283
  {
 
2284
    return 2;
 
2285
  }
 
2286
 
 
2287
  /// Return the number of dofs on each cell facet
 
2288
  virtual unsigned int num_facet_dofs() const
 
2289
  {
 
2290
    return 2;
 
2291
  }
 
2292
 
 
2293
  /// Return the number of dofs associated with each cell entity of dimension d
 
2294
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
2295
  {
 
2296
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2297
  }
 
2298
 
 
2299
  /// Tabulate the local-to-global mapping of dofs on a cell
 
2300
  virtual void tabulate_dofs(unsigned int* dofs,
 
2301
                             const ufc::mesh& m,
 
2302
                             const ufc::cell& c) const
 
2303
  {
 
2304
    dofs[0] = c.entity_indices[0][0];
 
2305
    dofs[1] = c.entity_indices[0][1];
 
2306
    dofs[2] = c.entity_indices[0][2];
 
2307
  }
 
2308
 
 
2309
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
2310
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
2311
                                   unsigned int facet) const
 
2312
  {
 
2313
    switch ( facet )
 
2314
    {
 
2315
    case 0:
 
2316
      dofs[0] = 1;
 
2317
      dofs[1] = 2;
 
2318
      break;
 
2319
    case 1:
 
2320
      dofs[0] = 0;
 
2321
      dofs[1] = 2;
 
2322
      break;
 
2323
    case 2:
 
2324
      dofs[0] = 0;
 
2325
      dofs[1] = 1;
 
2326
      break;
 
2327
    }
 
2328
  }
 
2329
 
 
2330
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
2331
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
2332
                                    unsigned int d, unsigned int i) const
 
2333
  {
 
2334
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2335
  }
 
2336
 
 
2337
  /// Tabulate the coordinates of all dofs on a cell
 
2338
  virtual void tabulate_coordinates(double** coordinates,
 
2339
                                    const ufc::cell& c) const
 
2340
  {
 
2341
    const double * const * x = c.coordinates;
 
2342
    coordinates[0][0] = x[0][0];
 
2343
    coordinates[0][1] = x[0][1];
 
2344
    coordinates[1][0] = x[1][0];
 
2345
    coordinates[1][1] = x[1][1];
 
2346
    coordinates[2][0] = x[2][0];
 
2347
    coordinates[2][1] = x[2][1];
 
2348
  }
 
2349
 
 
2350
  /// Return the number of sub dof maps (for a mixed element)
 
2351
  virtual unsigned int num_sub_dof_maps() const
 
2352
  {
 
2353
    return 1;
 
2354
  }
 
2355
 
 
2356
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
2357
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
2358
  {
 
2359
    return new UFC_P1ProjectionLinearForm_dof_map_0();
 
2360
  }
 
2361
 
 
2362
};
 
2363
 
 
2364
/// This class defines the interface for a local-to-global mapping of
 
2365
/// degrees of freedom (dofs).
 
2366
 
 
2367
class UFC_P1ProjectionLinearForm_dof_map_1: public ufc::dof_map
 
2368
{
 
2369
private:
 
2370
 
 
2371
  unsigned int __global_dimension;
 
2372
 
 
2373
public:
 
2374
 
 
2375
  /// Constructor
 
2376
  UFC_P1ProjectionLinearForm_dof_map_1() : ufc::dof_map()
 
2377
  {
 
2378
    __global_dimension = 0;
 
2379
  }
 
2380
 
 
2381
  /// Destructor
 
2382
  virtual ~UFC_P1ProjectionLinearForm_dof_map_1()
 
2383
  {
 
2384
    // Do nothing
 
2385
  }
 
2386
 
 
2387
  /// Return a string identifying the dof map
 
2388
  virtual const char* signature() const
 
2389
  {
 
2390
    return "FFC dof map for FiniteElement('Discontinuous Lagrange', 'triangle', 1)";
 
2391
  }
 
2392
 
 
2393
  /// Return true iff mesh entities of topological dimension d are needed
 
2394
  virtual bool needs_mesh_entities(unsigned int d) const
 
2395
  {
 
2396
    switch ( d )
 
2397
    {
 
2398
    case 0:
 
2399
      return false;
 
2400
      break;
 
2401
    case 1:
 
2402
      return false;
 
2403
      break;
 
2404
    case 2:
 
2405
      return true;
 
2406
      break;
 
2407
    }
 
2408
    return false;
 
2409
  }
 
2410
 
 
2411
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
2412
  virtual bool init_mesh(const ufc::mesh& m)
 
2413
  {
 
2414
    __global_dimension = 3*m.num_entities[2];
 
2415
    return false;
 
2416
  }
 
2417
 
 
2418
  /// Initialize dof map for given cell
 
2419
  virtual void init_cell(const ufc::mesh& m,
 
2420
                         const ufc::cell& c)
 
2421
  {
 
2422
    // Do nothing
 
2423
  }
 
2424
 
 
2425
  /// Finish initialization of dof map for cells
 
2426
  virtual void init_cell_finalize()
 
2427
  {
 
2428
    // Do nothing
 
2429
  }
 
2430
 
 
2431
  /// Return the dimension of the global finite element function space
 
2432
  virtual unsigned int global_dimension() const
 
2433
  {
 
2434
    return __global_dimension;
 
2435
  }
 
2436
 
 
2437
  /// Return the dimension of the local finite element function space
 
2438
  virtual unsigned int local_dimension() const
 
2439
  {
 
2440
    return 3;
 
2441
  }
 
2442
 
 
2443
  // Return the geometric dimension of the coordinates this dof map provides
 
2444
  virtual unsigned int geometric_dimension() const
 
2445
  {
 
2446
    return 2;
 
2447
  }
 
2448
 
 
2449
  /// Return the number of dofs on each cell facet
 
2450
  virtual unsigned int num_facet_dofs() const
 
2451
  {
 
2452
    return 0;
 
2453
  }
 
2454
 
 
2455
  /// Return the number of dofs associated with each cell entity of dimension d
 
2456
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
2457
  {
 
2458
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2459
  }
 
2460
 
 
2461
  /// Tabulate the local-to-global mapping of dofs on a cell
 
2462
  virtual void tabulate_dofs(unsigned int* dofs,
 
2463
                             const ufc::mesh& m,
 
2464
                             const ufc::cell& c) const
 
2465
  {
 
2466
    dofs[0] = 3*c.entity_indices[2][0];
 
2467
    dofs[1] = 3*c.entity_indices[2][0] + 1;
 
2468
    dofs[2] = 3*c.entity_indices[2][0] + 2;
 
2469
  }
 
2470
 
 
2471
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
2472
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
2473
                                   unsigned int facet) const
 
2474
  {
 
2475
    switch ( facet )
 
2476
    {
 
2477
    case 0:
 
2478
      
 
2479
      break;
 
2480
    case 1:
 
2481
      
 
2482
      break;
 
2483
    case 2:
 
2484
      
 
2485
      break;
 
2486
    }
 
2487
  }
 
2488
 
 
2489
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
2490
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
2491
                                    unsigned int d, unsigned int i) const
 
2492
  {
 
2493
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2494
  }
 
2495
 
 
2496
  /// Tabulate the coordinates of all dofs on a cell
 
2497
  virtual void tabulate_coordinates(double** coordinates,
 
2498
                                    const ufc::cell& c) const
 
2499
  {
 
2500
    const double * const * x = c.coordinates;
 
2501
    coordinates[0][0] = x[0][0];
 
2502
    coordinates[0][1] = x[0][1];
 
2503
    coordinates[1][0] = x[1][0];
 
2504
    coordinates[1][1] = x[1][1];
 
2505
    coordinates[2][0] = x[2][0];
 
2506
    coordinates[2][1] = x[2][1];
 
2507
  }
 
2508
 
 
2509
  /// Return the number of sub dof maps (for a mixed element)
 
2510
  virtual unsigned int num_sub_dof_maps() const
 
2511
  {
 
2512
    return 1;
 
2513
  }
 
2514
 
 
2515
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
2516
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
2517
  {
 
2518
    return new UFC_P1ProjectionLinearForm_dof_map_1();
 
2519
  }
 
2520
 
 
2521
};
 
2522
 
 
2523
/// This class defines the interface for the tabulation of the cell
 
2524
/// tensor corresponding to the local contribution to a form from
 
2525
/// the integral over a cell.
 
2526
 
 
2527
class UFC_P1ProjectionLinearForm_cell_integral_0: public ufc::cell_integral
 
2528
{
 
2529
public:
 
2530
 
 
2531
  /// Constructor
 
2532
  UFC_P1ProjectionLinearForm_cell_integral_0() : ufc::cell_integral()
 
2533
  {
 
2534
    // Do nothing
 
2535
  }
 
2536
 
 
2537
  /// Destructor
 
2538
  virtual ~UFC_P1ProjectionLinearForm_cell_integral_0()
 
2539
  {
 
2540
    // Do nothing
 
2541
  }
 
2542
 
 
2543
  /// Tabulate the tensor for the contribution from a local cell
 
2544
  virtual void tabulate_tensor(double* A,
 
2545
                               const double * const * w,
 
2546
                               const ufc::cell& c) const
 
2547
  {
 
2548
    // Extract vertex coordinates
 
2549
    const double * const * x = c.coordinates;
 
2550
    
 
2551
    // Compute Jacobian of affine map from reference cell
 
2552
    const double J_00 = x[1][0] - x[0][0];
 
2553
    const double J_01 = x[2][0] - x[0][0];
 
2554
    const double J_10 = x[1][1] - x[0][1];
 
2555
    const double J_11 = x[2][1] - x[0][1];
 
2556
      
 
2557
    // Compute determinant of Jacobian
 
2558
    double detJ = J_00*J_11 - J_01*J_10;
 
2559
      
 
2560
    // Compute inverse of Jacobian
 
2561
    
 
2562
    // Set scale factor
 
2563
    const double det = std::abs(detJ);
 
2564
    
 
2565
    // Number of operations to compute element tensor = 18
 
2566
    // Compute coefficients
 
2567
    const double c0_0_0_0 = w[0][0];
 
2568
    const double c0_0_0_1 = w[0][1];
 
2569
    const double c0_0_0_2 = w[0][2];
 
2570
    
 
2571
    // Compute geometry tensors
 
2572
    // Number of operations to compute decalrations = 3
 
2573
    const double G0_0 = det*c0_0_0_0;
 
2574
    const double G0_1 = det*c0_0_0_1;
 
2575
    const double G0_2 = det*c0_0_0_2;
 
2576
    
 
2577
    // Compute element tensor
 
2578
    // Number of operations to compute tensor = 15
 
2579
    A[0] = 0.0833333333333332*G0_0 + 0.0416666666666666*G0_1 + 0.0416666666666666*G0_2;
 
2580
    A[1] = 0.0416666666666666*G0_0 + 0.0833333333333332*G0_1 + 0.0416666666666666*G0_2;
 
2581
    A[2] = 0.0416666666666666*G0_0 + 0.0416666666666666*G0_1 + 0.0833333333333332*G0_2;
 
2582
  }
 
2583
 
 
2584
};
 
2585
 
 
2586
/// This class defines the interface for the assembly of the global
 
2587
/// tensor corresponding to a form with r + n arguments, that is, a
 
2588
/// mapping
 
2589
///
 
2590
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
2591
///
 
2592
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
2593
/// global tensor A is defined by
 
2594
///
 
2595
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
2596
///
 
2597
/// where each argument Vj represents the application to the
 
2598
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
2599
/// fixed functions (coefficients).
 
2600
 
 
2601
class UFC_P1ProjectionLinearForm: public ufc::form
 
2602
{
 
2603
public:
 
2604
 
 
2605
  /// Constructor
 
2606
  UFC_P1ProjectionLinearForm() : ufc::form()
 
2607
  {
 
2608
    // Do nothing
 
2609
  }
 
2610
 
 
2611
  /// Destructor
 
2612
  virtual ~UFC_P1ProjectionLinearForm()
 
2613
  {
 
2614
    // Do nothing
 
2615
  }
 
2616
 
 
2617
  /// Return a string identifying the form
 
2618
  virtual const char* signature() const
 
2619
  {
 
2620
    return "w0_a0[0, 1, 2] | vi0[0, 1, 2]*va0[0, 1, 2]*dX(0)";
 
2621
  }
 
2622
 
 
2623
  /// Return the rank of the global tensor (r)
 
2624
  virtual unsigned int rank() const
 
2625
  {
 
2626
    return 1;
 
2627
  }
 
2628
 
 
2629
  /// Return the number of coefficients (n)
 
2630
  virtual unsigned int num_coefficients() const
 
2631
  {
 
2632
    return 1;
 
2633
  }
 
2634
 
 
2635
  /// Return the number of cell integrals
 
2636
  virtual unsigned int num_cell_integrals() const
 
2637
  {
 
2638
    return 1;
 
2639
  }
 
2640
  
 
2641
  /// Return the number of exterior facet integrals
 
2642
  virtual unsigned int num_exterior_facet_integrals() const
 
2643
  {
 
2644
    return 0;
 
2645
  }
 
2646
  
 
2647
  /// Return the number of interior facet integrals
 
2648
  virtual unsigned int num_interior_facet_integrals() const
 
2649
  {
 
2650
    return 0;
 
2651
  }
 
2652
    
 
2653
  /// Create a new finite element for argument function i
 
2654
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
2655
  {
 
2656
    switch ( i )
 
2657
    {
 
2658
    case 0:
 
2659
      return new UFC_P1ProjectionLinearForm_finite_element_0();
 
2660
      break;
 
2661
    case 1:
 
2662
      return new UFC_P1ProjectionLinearForm_finite_element_1();
 
2663
      break;
 
2664
    }
 
2665
    return 0;
 
2666
  }
 
2667
  
 
2668
  /// Create a new dof map for argument function i
 
2669
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
2670
  {
 
2671
    switch ( i )
 
2672
    {
 
2673
    case 0:
 
2674
      return new UFC_P1ProjectionLinearForm_dof_map_0();
 
2675
      break;
 
2676
    case 1:
 
2677
      return new UFC_P1ProjectionLinearForm_dof_map_1();
 
2678
      break;
 
2679
    }
 
2680
    return 0;
 
2681
  }
 
2682
 
 
2683
  /// Create a new cell integral on sub domain i
 
2684
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
2685
  {
 
2686
    return new UFC_P1ProjectionLinearForm_cell_integral_0();
 
2687
  }
 
2688
 
 
2689
  /// Create a new exterior facet integral on sub domain i
 
2690
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
2691
  {
 
2692
    return 0;
 
2693
  }
 
2694
 
 
2695
  /// Create a new interior facet integral on sub domain i
 
2696
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
2697
  {
 
2698
    return 0;
 
2699
  }
 
2700
 
 
2701
};
 
2702
 
 
2703
// DOLFIN wrappers
 
2704
 
 
2705
#include <dolfin/fem/Form.h>
 
2706
#include <dolfin/fem/FiniteElement.h>
 
2707
#include <dolfin/fem/DofMap.h>
 
2708
#include <dolfin/function/Coefficient.h>
 
2709
#include <dolfin/function/Function.h>
 
2710
#include <dolfin/function/FunctionSpace.h>
 
2711
 
 
2712
class P1ProjectionBilinearFormFunctionSpace0 : public dolfin::FunctionSpace
 
2713
{
 
2714
public:
 
2715
 
 
2716
  P1ProjectionBilinearFormFunctionSpace0(const dolfin::Mesh& mesh)
 
2717
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2718
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_0()))),
 
2719
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_0()), mesh)))
 
2720
  {
 
2721
    // Do nothing
 
2722
  }
 
2723
 
 
2724
};
 
2725
 
 
2726
class P1ProjectionBilinearFormFunctionSpace1 : public dolfin::FunctionSpace
 
2727
{
 
2728
public:
 
2729
 
 
2730
  P1ProjectionBilinearFormFunctionSpace1(const dolfin::Mesh& mesh)
 
2731
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2732
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_0()))),
 
2733
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_0()), mesh)))
 
2734
  {
 
2735
    // Do nothing
 
2736
  }
 
2737
 
 
2738
};
 
2739
 
 
2740
class P1ProjectionLinearFormFunctionSpace0 : public dolfin::FunctionSpace
 
2741
{
 
2742
public:
 
2743
 
 
2744
  P1ProjectionLinearFormFunctionSpace0(const dolfin::Mesh& mesh)
 
2745
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2746
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_0()))),
 
2747
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_0()), mesh)))
 
2748
  {
 
2749
    // Do nothing
 
2750
  }
 
2751
 
 
2752
};
 
2753
 
 
2754
class P1ProjectionLinearFormCoefficientSpace0 : public dolfin::FunctionSpace
 
2755
{
 
2756
public:
 
2757
 
 
2758
  P1ProjectionLinearFormCoefficientSpace0(const dolfin::Mesh& mesh)
 
2759
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2760
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_1()))),
 
2761
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_1()), mesh)))
 
2762
  {
 
2763
    // Do nothing
 
2764
  }
 
2765
 
 
2766
};
 
2767
 
 
2768
class P1ProjectionTestSpace : public dolfin::FunctionSpace
 
2769
{
 
2770
public:
 
2771
 
 
2772
  P1ProjectionTestSpace(const dolfin::Mesh& mesh)
 
2773
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2774
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_0()))),
 
2775
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_0()), mesh)))
 
2776
  {
 
2777
    // Do nothing
 
2778
  }
 
2779
 
 
2780
};
 
2781
 
 
2782
class P1ProjectionTrialSpace : public dolfin::FunctionSpace
 
2783
{
 
2784
public:
 
2785
 
 
2786
  P1ProjectionTrialSpace(const dolfin::Mesh& mesh)
 
2787
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2788
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_0()))),
 
2789
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_0()), mesh)))
 
2790
  {
 
2791
    // Do nothing
 
2792
  }
 
2793
 
 
2794
};
 
2795
 
 
2796
class P1ProjectionCoefficientSpace : public dolfin::FunctionSpace
 
2797
{
 
2798
public:
 
2799
 
 
2800
  P1ProjectionCoefficientSpace(const dolfin::Mesh& mesh)
 
2801
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2802
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_1()))),
 
2803
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_1()), mesh)))
 
2804
  {
 
2805
    // Do nothing
 
2806
  }
 
2807
 
 
2808
};
 
2809
 
 
2810
class P1ProjectionFunctionSpace : public dolfin::FunctionSpace
 
2811
{
 
2812
public:
 
2813
 
 
2814
  P1ProjectionFunctionSpace(const dolfin::Mesh& mesh)
 
2815
    : dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),
 
2816
                            boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_P1ProjectionLinearForm_finite_element_0()))),
 
2817
                            boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_P1ProjectionLinearForm_dof_map_0()), mesh)))
 
2818
  {
 
2819
    // Do nothing
 
2820
  }
 
2821
 
 
2822
};
 
2823
 
 
2824
class P1ProjectionBilinearForm : public dolfin::Form
 
2825
{
 
2826
public:
 
2827
 
 
2828
  // Create form on given function space(s)
 
2829
  P1ProjectionBilinearForm(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1) : dolfin::Form()
 
2830
  {
 
2831
    boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());
 
2832
    _function_spaces.push_back(_V0);
 
2833
    boost::shared_ptr<const dolfin::FunctionSpace> _V1(&V1, dolfin::NoDeleter<const dolfin::FunctionSpace>());
 
2834
    _function_spaces.push_back(_V1);
 
2835
 
 
2836
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_P1ProjectionBilinearForm());
 
2837
  }
 
2838
 
 
2839
  // Create form on given function space(s) (shared data)
 
2840
  P1ProjectionBilinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::FunctionSpace> V1) : dolfin::Form()
 
2841
  {
 
2842
    _function_spaces.push_back(V0);
 
2843
    _function_spaces.push_back(V1);
 
2844
 
 
2845
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_P1ProjectionBilinearForm());
 
2846
  }
 
2847
 
 
2848
  // Destructor
 
2849
  ~P1ProjectionBilinearForm() {}
 
2850
 
 
2851
};
 
2852
 
 
2853
class P1ProjectionLinearFormCoefficient0 : public dolfin::Coefficient
 
2854
{
 
2855
public:
 
2856
 
 
2857
  // Constructor
 
2858
  P1ProjectionLinearFormCoefficient0(dolfin::Form& form) : dolfin::Coefficient(form) {}
 
2859
 
 
2860
  // Destructor  
 
2861
  ~P1ProjectionLinearFormCoefficient0() {}
 
2862
 
 
2863
  // Attach function to coefficient
 
2864
  const P1ProjectionLinearFormCoefficient0& operator= (dolfin::Function& v)
 
2865
  {
 
2866
    attach(v);
 
2867
    return *this;
 
2868
  }
 
2869
 
 
2870
  /// Create function space for coefficient
 
2871
  const dolfin::FunctionSpace* create_function_space() const
 
2872
  {
 
2873
    return new P1ProjectionLinearFormCoefficientSpace0(form.mesh());
 
2874
  }
 
2875
  
 
2876
  /// Return coefficient number
 
2877
  dolfin::uint number() const
 
2878
  {
 
2879
    return 0;
 
2880
  }
 
2881
  
 
2882
  /// Return coefficient name
 
2883
  virtual std::string name() const
 
2884
  {
 
2885
    return "u";
 
2886
  }
 
2887
  
 
2888
};
 
2889
class P1ProjectionLinearForm : public dolfin::Form
 
2890
{
 
2891
public:
 
2892
 
 
2893
  // Create form on given function space(s)
 
2894
  P1ProjectionLinearForm(const dolfin::FunctionSpace& V0) : dolfin::Form(), u(*this)
 
2895
  {
 
2896
    boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());
 
2897
    _function_spaces.push_back(_V0);
 
2898
 
 
2899
    _coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
 
2900
 
 
2901
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_P1ProjectionLinearForm());
 
2902
  }
 
2903
 
 
2904
  // Create form on given function space(s) (shared data)
 
2905
  P1ProjectionLinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0) : dolfin::Form(), u(*this)
 
2906
  {
 
2907
    _function_spaces.push_back(V0);
 
2908
 
 
2909
    _coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
 
2910
 
 
2911
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_P1ProjectionLinearForm());
 
2912
  }
 
2913
 
 
2914
  // Create form on given function space(s) with given coefficient(s)
 
2915
  P1ProjectionLinearForm(const dolfin::FunctionSpace& V0, dolfin::Function& w0) : dolfin::Form(), u(*this)
 
2916
  {
 
2917
    boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());
 
2918
    _function_spaces.push_back(_V0);
 
2919
 
 
2920
    _coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
 
2921
 
 
2922
    this->u = w0;
 
2923
 
 
2924
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_P1ProjectionLinearForm());
 
2925
  }
 
2926
 
 
2927
  // Create form on given function space(s) with given coefficient(s) (shared data)
 
2928
  P1ProjectionLinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, dolfin::Function& w0) : dolfin::Form(), u(*this)
 
2929
  {
 
2930
    _function_spaces.push_back(V0);
 
2931
 
 
2932
    _coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));
 
2933
 
 
2934
    this->u = w0;
 
2935
 
 
2936
    _ufc_form = boost::shared_ptr<const ufc::form>(new UFC_P1ProjectionLinearForm());
 
2937
  }
 
2938
 
 
2939
  // Destructor
 
2940
  ~P1ProjectionLinearForm() {}
 
2941
 
 
2942
  // Coefficients
 
2943
  P1ProjectionLinearFormCoefficient0 u;
 
2944
 
 
2945
};
 
2946
 
 
2947
#endif