← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/mf/ffc-forms/ConvectionMatrix2D.h

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/mf/ffc-forms/ConvectionMatrix2D.h
1
 
// This code conforms with the UFC specification version 1.0
2
 
// and was automatically generated by FFC version 0.4.3.
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 __CONVECTIONMATRIX2D_H
8
 
#define __CONVECTIONMATRIX2D_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_ConvectionMatrix2DBilinearForm_finite_element_0: public ufc::finite_element
18
 
{
19
 
public:
20
 
 
21
 
  /// Constructor
22
 
  UFC_ConvectionMatrix2DBilinearForm_finite_element_0() : ufc::finite_element()
23
 
  {
24
 
    // Do nothing
25
 
  }
26
 
 
27
 
  /// Destructor
28
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_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 "Lagrange finite element of degree 1 on a triangle";
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_ConvectionMatrix2DBilinearForm_finite_element_0();
432
 
  }
433
 
 
434
 
};
435
 
 
436
 
/// This class defines the interface for a finite element.
437
 
 
438
 
class UFC_ConvectionMatrix2DBilinearForm_finite_element_1: public ufc::finite_element
439
 
{
440
 
public:
441
 
 
442
 
  /// Constructor
443
 
  UFC_ConvectionMatrix2DBilinearForm_finite_element_1() : ufc::finite_element()
444
 
  {
445
 
    // Do nothing
446
 
  }
447
 
 
448
 
  /// Destructor
449
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_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 "Lagrange finite element of degree 1 on a triangle";
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_ConvectionMatrix2DBilinearForm_finite_element_1();
853
 
  }
854
 
 
855
 
};
856
 
 
857
 
/// This class defines the interface for a finite element.
858
 
 
859
 
class UFC_ConvectionMatrix2DBilinearForm_finite_element_2: public ufc::finite_element
860
 
{
861
 
public:
862
 
 
863
 
  /// Constructor
864
 
  UFC_ConvectionMatrix2DBilinearForm_finite_element_2() : ufc::finite_element()
865
 
  {
866
 
    // Do nothing
867
 
  }
868
 
 
869
 
  /// Destructor
870
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_finite_element_2()
871
 
  {
872
 
    // Do nothing
873
 
  }
874
 
 
875
 
  /// Return a string identifying the finite element
876
 
  virtual const char* signature() const
877
 
  {
878
 
    return "Discontinuous Lagrange finite element of degree 0 on a triangle";
879
 
  }
880
 
 
881
 
  /// Return the cell shape
882
 
  virtual ufc::shape cell_shape() const
883
 
  {
884
 
    return ufc::triangle;
885
 
  }
886
 
 
887
 
  /// Return the dimension of the finite element function space
888
 
  virtual unsigned int space_dimension() const
889
 
  {
890
 
    return 1;
891
 
  }
892
 
 
893
 
  /// Return the rank of the value space
894
 
  virtual unsigned int value_rank() const
895
 
  {
896
 
    return 0;
897
 
  }
898
 
 
899
 
  /// Return the dimension of the value space for axis i
900
 
  virtual unsigned int value_dimension(unsigned int i) const
901
 
  {
902
 
    return 1;
903
 
  }
904
 
 
905
 
  /// Evaluate basis function i at given point in cell
906
 
  virtual void evaluate_basis(unsigned int i,
907
 
                              double* values,
908
 
                              const double* coordinates,
909
 
                              const ufc::cell& c) const
910
 
  {
911
 
    // Extract vertex coordinates
912
 
    const double * const * element_coordinates = c.coordinates;
913
 
    
914
 
    // Compute Jacobian of affine map from reference cell
915
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
916
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
917
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
918
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
919
 
      
920
 
    // Compute determinant of Jacobian
921
 
    const double detJ = J_00*J_11 - J_01*J_10;
922
 
    
923
 
    // Compute inverse of Jacobian
924
 
    
925
 
    // Get coordinates and map to the reference (UFC) element
926
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
927
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
928
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
929
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
930
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
931
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
932
 
    
933
 
    // Map coordinates to the reference square
934
 
    if (std::abs(y - 1.0) < 1e-14)
935
 
      x = -1.0;
936
 
    else
937
 
      x = 2.0 *x/(1.0 - y) - 1.0;
938
 
    y = 2.0*y - 1.0;
939
 
    
940
 
    // Reset values
941
 
    *values = 0;
942
 
    
943
 
    // Map degree of freedom to element degree of freedom
944
 
    const unsigned int dof = i;
945
 
    
946
 
    // Generate scalings
947
 
    const double scalings_y_0 = 1;
948
 
    
949
 
    // Compute psitilde_a
950
 
    const double psitilde_a_0 = 1;
951
 
    
952
 
    // Compute psitilde_bs
953
 
    const double psitilde_bs_0_0 = 1;
954
 
    
955
 
    // Compute basisvalues
956
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
957
 
    
958
 
    // Table(s) of coefficients
959
 
    const static double coefficients0[1][1] = \
960
 
    {{1.41421356237309}};
961
 
    
962
 
    // Extract relevant coefficients
963
 
    const double coeff0_0 = coefficients0[dof][0];
964
 
    
965
 
    // Compute value(s)
966
 
    *values = coeff0_0*basisvalue0;
967
 
  }
968
 
 
969
 
  /// Evaluate all basis functions at given point in cell
970
 
  virtual void evaluate_basis_all(double* values,
971
 
                                  const double* coordinates,
972
 
                                  const ufc::cell& c) const
973
 
  {
974
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
975
 
  }
976
 
 
977
 
  /// Evaluate order n derivatives of basis function i at given point in cell
978
 
  virtual void evaluate_basis_derivatives(unsigned int i,
979
 
                                          unsigned int n,
980
 
                                          double* values,
981
 
                                          const double* coordinates,
982
 
                                          const ufc::cell& c) const
983
 
  {
984
 
    // Extract vertex coordinates
985
 
    const double * const * element_coordinates = c.coordinates;
986
 
    
987
 
    // Compute Jacobian of affine map from reference cell
988
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
989
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
990
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
991
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
992
 
      
993
 
    // Compute determinant of Jacobian
994
 
    const double detJ = J_00*J_11 - J_01*J_10;
995
 
    
996
 
    // Compute inverse of Jacobian
997
 
    
998
 
    // Get coordinates and map to the reference (UFC) element
999
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1000
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
1001
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1002
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1003
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
1004
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1005
 
    
1006
 
    // Map coordinates to the reference square
1007
 
    if (std::abs(y - 1.0) < 1e-14)
1008
 
      x = -1.0;
1009
 
    else
1010
 
      x = 2.0 *x/(1.0 - y) - 1.0;
1011
 
    y = 2.0*y - 1.0;
1012
 
    
1013
 
    // Compute number of derivatives
1014
 
    unsigned int num_derivatives = 1;
1015
 
    
1016
 
    for (unsigned int j = 0; j < n; j++)
1017
 
      num_derivatives *= 2;
1018
 
    
1019
 
    
1020
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1021
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
1022
 
        
1023
 
    for (unsigned int j = 0; j < num_derivatives; j++)
1024
 
    {
1025
 
      combinations[j] = new unsigned int [n];
1026
 
      for (unsigned int k = 0; k < n; k++)
1027
 
        combinations[j][k] = 0;
1028
 
    }
1029
 
        
1030
 
    // Generate combinations of derivatives
1031
 
    for (unsigned int row = 1; row < num_derivatives; row++)
1032
 
    {
1033
 
      for (unsigned int num = 0; num < row; num++)
1034
 
      {
1035
 
        for (unsigned int col = n-1; col+1 > 0; col--)
1036
 
        {
1037
 
          if (combinations[row][col] + 1 > 1)
1038
 
            combinations[row][col] = 0;
1039
 
          else
1040
 
          {
1041
 
            combinations[row][col] += 1;
1042
 
            break;
1043
 
          }
1044
 
        }
1045
 
      }
1046
 
    }
1047
 
    
1048
 
    // Compute inverse of Jacobian
1049
 
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1050
 
    
1051
 
    // Declare transformation matrix
1052
 
    // Declare pointer to two dimensional array and initialise
1053
 
    double **transform = new double *[num_derivatives];
1054
 
        
1055
 
    for (unsigned int j = 0; j < num_derivatives; j++)
1056
 
    {
1057
 
      transform[j] = new double [num_derivatives];
1058
 
      for (unsigned int k = 0; k < num_derivatives; k++)
1059
 
        transform[j][k] = 1;
1060
 
    }
1061
 
    
1062
 
    // Construct transformation matrix
1063
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1064
 
    {
1065
 
      for (unsigned int col = 0; col < num_derivatives; col++)
1066
 
      {
1067
 
        for (unsigned int k = 0; k < n; k++)
1068
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1069
 
      }
1070
 
    }
1071
 
    
1072
 
    // Reset values
1073
 
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
1074
 
      values[j] = 0;
1075
 
    
1076
 
    // Map degree of freedom to element degree of freedom
1077
 
    const unsigned int dof = i;
1078
 
    
1079
 
    // Generate scalings
1080
 
    const double scalings_y_0 = 1;
1081
 
    
1082
 
    // Compute psitilde_a
1083
 
    const double psitilde_a_0 = 1;
1084
 
    
1085
 
    // Compute psitilde_bs
1086
 
    const double psitilde_bs_0_0 = 1;
1087
 
    
1088
 
    // Compute basisvalues
1089
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1090
 
    
1091
 
    // Table(s) of coefficients
1092
 
    const static double coefficients0[1][1] = \
1093
 
    {{1.41421356237309}};
1094
 
    
1095
 
    // Interesting (new) part
1096
 
    // Tables of derivatives of the polynomial base (transpose)
1097
 
    const static double dmats0[1][1] = \
1098
 
    {{0}};
1099
 
    
1100
 
    const static double dmats1[1][1] = \
1101
 
    {{0}};
1102
 
    
1103
 
    // Compute reference derivatives
1104
 
    // Declare pointer to array of derivatives on FIAT element
1105
 
    double *derivatives = new double [num_derivatives];
1106
 
    
1107
 
    // Declare coefficients
1108
 
    double coeff0_0 = 0;
1109
 
    
1110
 
    // Declare new coefficients
1111
 
    double new_coeff0_0 = 0;
1112
 
    
1113
 
    // Loop possible derivatives
1114
 
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1115
 
    {
1116
 
      // Get values from coefficients array
1117
 
      new_coeff0_0 = coefficients0[dof][0];
1118
 
    
1119
 
      // Loop derivative order
1120
 
      for (unsigned int j = 0; j < n; j++)
1121
 
      {
1122
 
        // Update old coefficients
1123
 
        coeff0_0 = new_coeff0_0;
1124
 
    
1125
 
        if(combinations[deriv_num][j] == 0)
1126
 
        {
1127
 
          new_coeff0_0 = coeff0_0*dmats0[0][0];
1128
 
        }
1129
 
        if(combinations[deriv_num][j] == 1)
1130
 
        {
1131
 
          new_coeff0_0 = coeff0_0*dmats1[0][0];
1132
 
        }
1133
 
    
1134
 
      }
1135
 
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
1136
 
      derivatives[deriv_num] = new_coeff0_0*basisvalue0;
1137
 
    }
1138
 
    
1139
 
    // Transform derivatives back to physical element
1140
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1141
 
    {
1142
 
      for (unsigned int col = 0; col < num_derivatives; col++)
1143
 
      {
1144
 
        values[row] += transform[row][col]*derivatives[col];
1145
 
      }
1146
 
    }
1147
 
    // Delete pointer to array of derivatives on FIAT element
1148
 
    delete [] derivatives;
1149
 
    
1150
 
    // Delete pointer to array of combinations of derivatives and transform
1151
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1152
 
    {
1153
 
      delete [] combinations[row];
1154
 
      delete [] transform[row];
1155
 
    }
1156
 
    
1157
 
    delete [] combinations;
1158
 
    delete [] transform;
1159
 
  }
1160
 
 
1161
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
1162
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
1163
 
                                              double* values,
1164
 
                                              const double* coordinates,
1165
 
                                              const ufc::cell& c) const
1166
 
  {
1167
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1168
 
  }
1169
 
 
1170
 
  /// Evaluate linear functional for dof i on the function f
1171
 
  virtual double evaluate_dof(unsigned int i,
1172
 
                              const ufc::function& f,
1173
 
                              const ufc::cell& c) const
1174
 
  {
1175
 
    // The reference points, direction and weights:
1176
 
    const static double X[1][1][2] = {{{0.333333333333333, 0.333333333333333}}};
1177
 
    const static double W[1][1] = {{1}};
1178
 
    const static double D[1][1][1] = {{{1}}};
1179
 
    
1180
 
    const double * const * x = c.coordinates;
1181
 
    double result = 0.0;
1182
 
    // Iterate over the points:
1183
 
    // Evaluate basis functions for affine mapping
1184
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1185
 
    const double w1 = X[i][0][0];
1186
 
    const double w2 = X[i][0][1];
1187
 
    
1188
 
    // Compute affine mapping y = F(X)
1189
 
    double y[2];
1190
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1191
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1192
 
    
1193
 
    // Evaluate function at physical points
1194
 
    double values[1];
1195
 
    f.evaluate(values, y, c);
1196
 
    
1197
 
    // Map function values using appropriate mapping
1198
 
    // Affine map: Do nothing
1199
 
    
1200
 
    // Note that we do not map the weights (yet).
1201
 
    
1202
 
    // Take directional components
1203
 
    for(int k = 0; k < 1; k++)
1204
 
      result += values[k]*D[i][0][k];
1205
 
    // Multiply by weights 
1206
 
    result *= W[i][0];
1207
 
    
1208
 
    return result;
1209
 
  }
1210
 
 
1211
 
  /// Evaluate linear functionals for all dofs on the function f
1212
 
  virtual void evaluate_dofs(double* values,
1213
 
                             const ufc::function& f,
1214
 
                             const ufc::cell& c) const
1215
 
  {
1216
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1217
 
  }
1218
 
 
1219
 
  /// Interpolate vertex values from dof values
1220
 
  virtual void interpolate_vertex_values(double* vertex_values,
1221
 
                                         const double* dof_values,
1222
 
                                         const ufc::cell& c) const
1223
 
  {
1224
 
    // Evaluate at vertices and use affine mapping
1225
 
    vertex_values[0] = dof_values[0];
1226
 
    vertex_values[1] = dof_values[0];
1227
 
    vertex_values[2] = dof_values[0];
1228
 
  }
1229
 
 
1230
 
  /// Return the number of sub elements (for a mixed element)
1231
 
  virtual unsigned int num_sub_elements() const
1232
 
  {
1233
 
    return 1;
1234
 
  }
1235
 
 
1236
 
  /// Create a new finite element for sub element i (for a mixed element)
1237
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
1238
 
  {
1239
 
    return new UFC_ConvectionMatrix2DBilinearForm_finite_element_2();
1240
 
  }
1241
 
 
1242
 
};
1243
 
 
1244
 
/// This class defines the interface for a finite element.
1245
 
 
1246
 
class UFC_ConvectionMatrix2DBilinearForm_finite_element_3: public ufc::finite_element
1247
 
{
1248
 
public:
1249
 
 
1250
 
  /// Constructor
1251
 
  UFC_ConvectionMatrix2DBilinearForm_finite_element_3() : ufc::finite_element()
1252
 
  {
1253
 
    // Do nothing
1254
 
  }
1255
 
 
1256
 
  /// Destructor
1257
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_finite_element_3()
1258
 
  {
1259
 
    // Do nothing
1260
 
  }
1261
 
 
1262
 
  /// Return a string identifying the finite element
1263
 
  virtual const char* signature() const
1264
 
  {
1265
 
    return "Discontinuous Lagrange finite element of degree 0 on a triangle";
1266
 
  }
1267
 
 
1268
 
  /// Return the cell shape
1269
 
  virtual ufc::shape cell_shape() const
1270
 
  {
1271
 
    return ufc::triangle;
1272
 
  }
1273
 
 
1274
 
  /// Return the dimension of the finite element function space
1275
 
  virtual unsigned int space_dimension() const
1276
 
  {
1277
 
    return 1;
1278
 
  }
1279
 
 
1280
 
  /// Return the rank of the value space
1281
 
  virtual unsigned int value_rank() const
1282
 
  {
1283
 
    return 0;
1284
 
  }
1285
 
 
1286
 
  /// Return the dimension of the value space for axis i
1287
 
  virtual unsigned int value_dimension(unsigned int i) const
1288
 
  {
1289
 
    return 1;
1290
 
  }
1291
 
 
1292
 
  /// Evaluate basis function i at given point in cell
1293
 
  virtual void evaluate_basis(unsigned int i,
1294
 
                              double* values,
1295
 
                              const double* coordinates,
1296
 
                              const ufc::cell& c) const
1297
 
  {
1298
 
    // Extract vertex coordinates
1299
 
    const double * const * element_coordinates = c.coordinates;
1300
 
    
1301
 
    // Compute Jacobian of affine map from reference cell
1302
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1303
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1304
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1305
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1306
 
      
1307
 
    // Compute determinant of Jacobian
1308
 
    const double detJ = J_00*J_11 - J_01*J_10;
1309
 
    
1310
 
    // Compute inverse of Jacobian
1311
 
    
1312
 
    // Get coordinates and map to the reference (UFC) element
1313
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1314
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
1315
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1316
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1317
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
1318
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1319
 
    
1320
 
    // Map coordinates to the reference square
1321
 
    if (std::abs(y - 1.0) < 1e-14)
1322
 
      x = -1.0;
1323
 
    else
1324
 
      x = 2.0 *x/(1.0 - y) - 1.0;
1325
 
    y = 2.0*y - 1.0;
1326
 
    
1327
 
    // Reset values
1328
 
    *values = 0;
1329
 
    
1330
 
    // Map degree of freedom to element degree of freedom
1331
 
    const unsigned int dof = i;
1332
 
    
1333
 
    // Generate scalings
1334
 
    const double scalings_y_0 = 1;
1335
 
    
1336
 
    // Compute psitilde_a
1337
 
    const double psitilde_a_0 = 1;
1338
 
    
1339
 
    // Compute psitilde_bs
1340
 
    const double psitilde_bs_0_0 = 1;
1341
 
    
1342
 
    // Compute basisvalues
1343
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1344
 
    
1345
 
    // Table(s) of coefficients
1346
 
    const static double coefficients0[1][1] = \
1347
 
    {{1.41421356237309}};
1348
 
    
1349
 
    // Extract relevant coefficients
1350
 
    const double coeff0_0 = coefficients0[dof][0];
1351
 
    
1352
 
    // Compute value(s)
1353
 
    *values = coeff0_0*basisvalue0;
1354
 
  }
1355
 
 
1356
 
  /// Evaluate all basis functions at given point in cell
1357
 
  virtual void evaluate_basis_all(double* values,
1358
 
                                  const double* coordinates,
1359
 
                                  const ufc::cell& c) const
1360
 
  {
1361
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1362
 
  }
1363
 
 
1364
 
  /// Evaluate order n derivatives of basis function i at given point in cell
1365
 
  virtual void evaluate_basis_derivatives(unsigned int i,
1366
 
                                          unsigned int n,
1367
 
                                          double* values,
1368
 
                                          const double* coordinates,
1369
 
                                          const ufc::cell& c) const
1370
 
  {
1371
 
    // Extract vertex coordinates
1372
 
    const double * const * element_coordinates = c.coordinates;
1373
 
    
1374
 
    // Compute Jacobian of affine map from reference cell
1375
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1376
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1377
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1378
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1379
 
      
1380
 
    // Compute determinant of Jacobian
1381
 
    const double detJ = J_00*J_11 - J_01*J_10;
1382
 
    
1383
 
    // Compute inverse of Jacobian
1384
 
    
1385
 
    // Get coordinates and map to the reference (UFC) element
1386
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1387
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
1388
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1389
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1390
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
1391
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1392
 
    
1393
 
    // Map coordinates to the reference square
1394
 
    if (std::abs(y - 1.0) < 1e-14)
1395
 
      x = -1.0;
1396
 
    else
1397
 
      x = 2.0 *x/(1.0 - y) - 1.0;
1398
 
    y = 2.0*y - 1.0;
1399
 
    
1400
 
    // Compute number of derivatives
1401
 
    unsigned int num_derivatives = 1;
1402
 
    
1403
 
    for (unsigned int j = 0; j < n; j++)
1404
 
      num_derivatives *= 2;
1405
 
    
1406
 
    
1407
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1408
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
1409
 
        
1410
 
    for (unsigned int j = 0; j < num_derivatives; j++)
1411
 
    {
1412
 
      combinations[j] = new unsigned int [n];
1413
 
      for (unsigned int k = 0; k < n; k++)
1414
 
        combinations[j][k] = 0;
1415
 
    }
1416
 
        
1417
 
    // Generate combinations of derivatives
1418
 
    for (unsigned int row = 1; row < num_derivatives; row++)
1419
 
    {
1420
 
      for (unsigned int num = 0; num < row; num++)
1421
 
      {
1422
 
        for (unsigned int col = n-1; col+1 > 0; col--)
1423
 
        {
1424
 
          if (combinations[row][col] + 1 > 1)
1425
 
            combinations[row][col] = 0;
1426
 
          else
1427
 
          {
1428
 
            combinations[row][col] += 1;
1429
 
            break;
1430
 
          }
1431
 
        }
1432
 
      }
1433
 
    }
1434
 
    
1435
 
    // Compute inverse of Jacobian
1436
 
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1437
 
    
1438
 
    // Declare transformation matrix
1439
 
    // Declare pointer to two dimensional array and initialise
1440
 
    double **transform = new double *[num_derivatives];
1441
 
        
1442
 
    for (unsigned int j = 0; j < num_derivatives; j++)
1443
 
    {
1444
 
      transform[j] = new double [num_derivatives];
1445
 
      for (unsigned int k = 0; k < num_derivatives; k++)
1446
 
        transform[j][k] = 1;
1447
 
    }
1448
 
    
1449
 
    // Construct transformation matrix
1450
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1451
 
    {
1452
 
      for (unsigned int col = 0; col < num_derivatives; col++)
1453
 
      {
1454
 
        for (unsigned int k = 0; k < n; k++)
1455
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1456
 
      }
1457
 
    }
1458
 
    
1459
 
    // Reset values
1460
 
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
1461
 
      values[j] = 0;
1462
 
    
1463
 
    // Map degree of freedom to element degree of freedom
1464
 
    const unsigned int dof = i;
1465
 
    
1466
 
    // Generate scalings
1467
 
    const double scalings_y_0 = 1;
1468
 
    
1469
 
    // Compute psitilde_a
1470
 
    const double psitilde_a_0 = 1;
1471
 
    
1472
 
    // Compute psitilde_bs
1473
 
    const double psitilde_bs_0_0 = 1;
1474
 
    
1475
 
    // Compute basisvalues
1476
 
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1477
 
    
1478
 
    // Table(s) of coefficients
1479
 
    const static double coefficients0[1][1] = \
1480
 
    {{1.41421356237309}};
1481
 
    
1482
 
    // Interesting (new) part
1483
 
    // Tables of derivatives of the polynomial base (transpose)
1484
 
    const static double dmats0[1][1] = \
1485
 
    {{0}};
1486
 
    
1487
 
    const static double dmats1[1][1] = \
1488
 
    {{0}};
1489
 
    
1490
 
    // Compute reference derivatives
1491
 
    // Declare pointer to array of derivatives on FIAT element
1492
 
    double *derivatives = new double [num_derivatives];
1493
 
    
1494
 
    // Declare coefficients
1495
 
    double coeff0_0 = 0;
1496
 
    
1497
 
    // Declare new coefficients
1498
 
    double new_coeff0_0 = 0;
1499
 
    
1500
 
    // Loop possible derivatives
1501
 
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1502
 
    {
1503
 
      // Get values from coefficients array
1504
 
      new_coeff0_0 = coefficients0[dof][0];
1505
 
    
1506
 
      // Loop derivative order
1507
 
      for (unsigned int j = 0; j < n; j++)
1508
 
      {
1509
 
        // Update old coefficients
1510
 
        coeff0_0 = new_coeff0_0;
1511
 
    
1512
 
        if(combinations[deriv_num][j] == 0)
1513
 
        {
1514
 
          new_coeff0_0 = coeff0_0*dmats0[0][0];
1515
 
        }
1516
 
        if(combinations[deriv_num][j] == 1)
1517
 
        {
1518
 
          new_coeff0_0 = coeff0_0*dmats1[0][0];
1519
 
        }
1520
 
    
1521
 
      }
1522
 
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
1523
 
      derivatives[deriv_num] = new_coeff0_0*basisvalue0;
1524
 
    }
1525
 
    
1526
 
    // Transform derivatives back to physical element
1527
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1528
 
    {
1529
 
      for (unsigned int col = 0; col < num_derivatives; col++)
1530
 
      {
1531
 
        values[row] += transform[row][col]*derivatives[col];
1532
 
      }
1533
 
    }
1534
 
    // Delete pointer to array of derivatives on FIAT element
1535
 
    delete [] derivatives;
1536
 
    
1537
 
    // Delete pointer to array of combinations of derivatives and transform
1538
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1539
 
    {
1540
 
      delete [] combinations[row];
1541
 
      delete [] transform[row];
1542
 
    }
1543
 
    
1544
 
    delete [] combinations;
1545
 
    delete [] transform;
1546
 
  }
1547
 
 
1548
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
1549
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
1550
 
                                              double* values,
1551
 
                                              const double* coordinates,
1552
 
                                              const ufc::cell& c) const
1553
 
  {
1554
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1555
 
  }
1556
 
 
1557
 
  /// Evaluate linear functional for dof i on the function f
1558
 
  virtual double evaluate_dof(unsigned int i,
1559
 
                              const ufc::function& f,
1560
 
                              const ufc::cell& c) const
1561
 
  {
1562
 
    // The reference points, direction and weights:
1563
 
    const static double X[1][1][2] = {{{0.333333333333333, 0.333333333333333}}};
1564
 
    const static double W[1][1] = {{1}};
1565
 
    const static double D[1][1][1] = {{{1}}};
1566
 
    
1567
 
    const double * const * x = c.coordinates;
1568
 
    double result = 0.0;
1569
 
    // Iterate over the points:
1570
 
    // Evaluate basis functions for affine mapping
1571
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1572
 
    const double w1 = X[i][0][0];
1573
 
    const double w2 = X[i][0][1];
1574
 
    
1575
 
    // Compute affine mapping y = F(X)
1576
 
    double y[2];
1577
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1578
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1579
 
    
1580
 
    // Evaluate function at physical points
1581
 
    double values[1];
1582
 
    f.evaluate(values, y, c);
1583
 
    
1584
 
    // Map function values using appropriate mapping
1585
 
    // Affine map: Do nothing
1586
 
    
1587
 
    // Note that we do not map the weights (yet).
1588
 
    
1589
 
    // Take directional components
1590
 
    for(int k = 0; k < 1; k++)
1591
 
      result += values[k]*D[i][0][k];
1592
 
    // Multiply by weights 
1593
 
    result *= W[i][0];
1594
 
    
1595
 
    return result;
1596
 
  }
1597
 
 
1598
 
  /// Evaluate linear functionals for all dofs on the function f
1599
 
  virtual void evaluate_dofs(double* values,
1600
 
                             const ufc::function& f,
1601
 
                             const ufc::cell& c) const
1602
 
  {
1603
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1604
 
  }
1605
 
 
1606
 
  /// Interpolate vertex values from dof values
1607
 
  virtual void interpolate_vertex_values(double* vertex_values,
1608
 
                                         const double* dof_values,
1609
 
                                         const ufc::cell& c) const
1610
 
  {
1611
 
    // Evaluate at vertices and use affine mapping
1612
 
    vertex_values[0] = dof_values[0];
1613
 
    vertex_values[1] = dof_values[0];
1614
 
    vertex_values[2] = dof_values[0];
1615
 
  }
1616
 
 
1617
 
  /// Return the number of sub elements (for a mixed element)
1618
 
  virtual unsigned int num_sub_elements() const
1619
 
  {
1620
 
    return 1;
1621
 
  }
1622
 
 
1623
 
  /// Create a new finite element for sub element i (for a mixed element)
1624
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
1625
 
  {
1626
 
    return new UFC_ConvectionMatrix2DBilinearForm_finite_element_3();
1627
 
  }
1628
 
 
1629
 
};
1630
 
 
1631
 
/// This class defines the interface for a local-to-global mapping of
1632
 
/// degrees of freedom (dofs).
1633
 
 
1634
 
class UFC_ConvectionMatrix2DBilinearForm_dof_map_0: public ufc::dof_map
1635
 
{
1636
 
private:
1637
 
 
1638
 
  unsigned int __global_dimension;
1639
 
 
1640
 
public:
1641
 
 
1642
 
  /// Constructor
1643
 
  UFC_ConvectionMatrix2DBilinearForm_dof_map_0() : ufc::dof_map()
1644
 
  {
1645
 
    __global_dimension = 0;
1646
 
  }
1647
 
 
1648
 
  /// Destructor
1649
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_dof_map_0()
1650
 
  {
1651
 
    // Do nothing
1652
 
  }
1653
 
 
1654
 
  /// Return a string identifying the dof map
1655
 
  virtual const char* signature() const
1656
 
  {
1657
 
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
1658
 
  }
1659
 
 
1660
 
  /// Return true iff mesh entities of topological dimension d are needed
1661
 
  virtual bool needs_mesh_entities(unsigned int d) const
1662
 
  {
1663
 
    switch ( d )
1664
 
    {
1665
 
    case 0:
1666
 
      return true;
1667
 
      break;
1668
 
    case 1:
1669
 
      return false;
1670
 
      break;
1671
 
    case 2:
1672
 
      return false;
1673
 
      break;
1674
 
    }
1675
 
    return false;
1676
 
  }
1677
 
 
1678
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
1679
 
  virtual bool init_mesh(const ufc::mesh& m)
1680
 
  {
1681
 
    __global_dimension = m.num_entities[0];
1682
 
    return false;
1683
 
  }
1684
 
 
1685
 
  /// Initialize dof map for given cell
1686
 
  virtual void init_cell(const ufc::mesh& m,
1687
 
                         const ufc::cell& c)
1688
 
  {
1689
 
    // Do nothing
1690
 
  }
1691
 
 
1692
 
  /// Finish initialization of dof map for cells
1693
 
  virtual void init_cell_finalize()
1694
 
  {
1695
 
    // Do nothing
1696
 
  }
1697
 
 
1698
 
  /// Return the dimension of the global finite element function space
1699
 
  virtual unsigned int global_dimension() const
1700
 
  {
1701
 
    return __global_dimension;
1702
 
  }
1703
 
 
1704
 
  /// Return the dimension of the local finite element function space
1705
 
  virtual unsigned int local_dimension() const
1706
 
  {
1707
 
    return 3;
1708
 
  }
1709
 
 
1710
 
  // Return the geometric dimension of the coordinates this dof map provides
1711
 
  virtual unsigned int geometric_dimension() const
1712
 
  {
1713
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1714
 
  }
1715
 
 
1716
 
  /// Return the number of dofs on each cell facet
1717
 
  virtual unsigned int num_facet_dofs() const
1718
 
  {
1719
 
    return 2;
1720
 
  }
1721
 
 
1722
 
  /// Return the number of dofs associated with each cell entity of dimension d
1723
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
1724
 
  {
1725
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1726
 
  }
1727
 
 
1728
 
  /// Tabulate the local-to-global mapping of dofs on a cell
1729
 
  virtual void tabulate_dofs(unsigned int* dofs,
1730
 
                             const ufc::mesh& m,
1731
 
                             const ufc::cell& c) const
1732
 
  {
1733
 
    dofs[0] = c.entity_indices[0][0];
1734
 
    dofs[1] = c.entity_indices[0][1];
1735
 
    dofs[2] = c.entity_indices[0][2];
1736
 
  }
1737
 
 
1738
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
1739
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
1740
 
                                   unsigned int facet) const
1741
 
  {
1742
 
    switch ( facet )
1743
 
    {
1744
 
    case 0:
1745
 
      dofs[0] = 1;
1746
 
      dofs[1] = 2;
1747
 
      break;
1748
 
    case 1:
1749
 
      dofs[0] = 0;
1750
 
      dofs[1] = 2;
1751
 
      break;
1752
 
    case 2:
1753
 
      dofs[0] = 0;
1754
 
      dofs[1] = 1;
1755
 
      break;
1756
 
    }
1757
 
  }
1758
 
 
1759
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
1760
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
1761
 
                                    unsigned int d, unsigned int i) const
1762
 
  {
1763
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1764
 
  }
1765
 
 
1766
 
  /// Tabulate the coordinates of all dofs on a cell
1767
 
  virtual void tabulate_coordinates(double** coordinates,
1768
 
                                    const ufc::cell& c) const
1769
 
  {
1770
 
    const double * const * x = c.coordinates;
1771
 
    coordinates[0][0] = x[0][0];
1772
 
    coordinates[0][1] = x[0][1];
1773
 
    coordinates[1][0] = x[1][0];
1774
 
    coordinates[1][1] = x[1][1];
1775
 
    coordinates[2][0] = x[2][0];
1776
 
    coordinates[2][1] = x[2][1];
1777
 
  }
1778
 
 
1779
 
  /// Return the number of sub dof maps (for a mixed element)
1780
 
  virtual unsigned int num_sub_dof_maps() const
1781
 
  {
1782
 
    return 1;
1783
 
  }
1784
 
 
1785
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
1786
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1787
 
  {
1788
 
    return new UFC_ConvectionMatrix2DBilinearForm_dof_map_0();
1789
 
  }
1790
 
 
1791
 
};
1792
 
 
1793
 
/// This class defines the interface for a local-to-global mapping of
1794
 
/// degrees of freedom (dofs).
1795
 
 
1796
 
class UFC_ConvectionMatrix2DBilinearForm_dof_map_1: public ufc::dof_map
1797
 
{
1798
 
private:
1799
 
 
1800
 
  unsigned int __global_dimension;
1801
 
 
1802
 
public:
1803
 
 
1804
 
  /// Constructor
1805
 
  UFC_ConvectionMatrix2DBilinearForm_dof_map_1() : ufc::dof_map()
1806
 
  {
1807
 
    __global_dimension = 0;
1808
 
  }
1809
 
 
1810
 
  /// Destructor
1811
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_dof_map_1()
1812
 
  {
1813
 
    // Do nothing
1814
 
  }
1815
 
 
1816
 
  /// Return a string identifying the dof map
1817
 
  virtual const char* signature() const
1818
 
  {
1819
 
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
1820
 
  }
1821
 
 
1822
 
  /// Return true iff mesh entities of topological dimension d are needed
1823
 
  virtual bool needs_mesh_entities(unsigned int d) const
1824
 
  {
1825
 
    switch ( d )
1826
 
    {
1827
 
    case 0:
1828
 
      return true;
1829
 
      break;
1830
 
    case 1:
1831
 
      return false;
1832
 
      break;
1833
 
    case 2:
1834
 
      return false;
1835
 
      break;
1836
 
    }
1837
 
    return false;
1838
 
  }
1839
 
 
1840
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
1841
 
  virtual bool init_mesh(const ufc::mesh& m)
1842
 
  {
1843
 
    __global_dimension = m.num_entities[0];
1844
 
    return false;
1845
 
  }
1846
 
 
1847
 
  /// Initialize dof map for given cell
1848
 
  virtual void init_cell(const ufc::mesh& m,
1849
 
                         const ufc::cell& c)
1850
 
  {
1851
 
    // Do nothing
1852
 
  }
1853
 
 
1854
 
  /// Finish initialization of dof map for cells
1855
 
  virtual void init_cell_finalize()
1856
 
  {
1857
 
    // Do nothing
1858
 
  }
1859
 
 
1860
 
  /// Return the dimension of the global finite element function space
1861
 
  virtual unsigned int global_dimension() const
1862
 
  {
1863
 
    return __global_dimension;
1864
 
  }
1865
 
 
1866
 
  /// Return the dimension of the local finite element function space
1867
 
  virtual unsigned int local_dimension() const
1868
 
  {
1869
 
    return 3;
1870
 
  }
1871
 
 
1872
 
  // Return the geometric dimension of the coordinates this dof map provides
1873
 
  virtual unsigned int geometric_dimension() const
1874
 
  {
1875
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1876
 
  }
1877
 
 
1878
 
  /// Return the number of dofs on each cell facet
1879
 
  virtual unsigned int num_facet_dofs() const
1880
 
  {
1881
 
    return 2;
1882
 
  }
1883
 
 
1884
 
  /// Return the number of dofs associated with each cell entity of dimension d
1885
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
1886
 
  {
1887
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1888
 
  }
1889
 
 
1890
 
  /// Tabulate the local-to-global mapping of dofs on a cell
1891
 
  virtual void tabulate_dofs(unsigned int* dofs,
1892
 
                             const ufc::mesh& m,
1893
 
                             const ufc::cell& c) const
1894
 
  {
1895
 
    dofs[0] = c.entity_indices[0][0];
1896
 
    dofs[1] = c.entity_indices[0][1];
1897
 
    dofs[2] = c.entity_indices[0][2];
1898
 
  }
1899
 
 
1900
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
1901
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
1902
 
                                   unsigned int facet) const
1903
 
  {
1904
 
    switch ( facet )
1905
 
    {
1906
 
    case 0:
1907
 
      dofs[0] = 1;
1908
 
      dofs[1] = 2;
1909
 
      break;
1910
 
    case 1:
1911
 
      dofs[0] = 0;
1912
 
      dofs[1] = 2;
1913
 
      break;
1914
 
    case 2:
1915
 
      dofs[0] = 0;
1916
 
      dofs[1] = 1;
1917
 
      break;
1918
 
    }
1919
 
  }
1920
 
 
1921
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
1922
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
1923
 
                                    unsigned int d, unsigned int i) const
1924
 
  {
1925
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1926
 
  }
1927
 
 
1928
 
  /// Tabulate the coordinates of all dofs on a cell
1929
 
  virtual void tabulate_coordinates(double** coordinates,
1930
 
                                    const ufc::cell& c) const
1931
 
  {
1932
 
    const double * const * x = c.coordinates;
1933
 
    coordinates[0][0] = x[0][0];
1934
 
    coordinates[0][1] = x[0][1];
1935
 
    coordinates[1][0] = x[1][0];
1936
 
    coordinates[1][1] = x[1][1];
1937
 
    coordinates[2][0] = x[2][0];
1938
 
    coordinates[2][1] = x[2][1];
1939
 
  }
1940
 
 
1941
 
  /// Return the number of sub dof maps (for a mixed element)
1942
 
  virtual unsigned int num_sub_dof_maps() const
1943
 
  {
1944
 
    return 1;
1945
 
  }
1946
 
 
1947
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
1948
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1949
 
  {
1950
 
    return new UFC_ConvectionMatrix2DBilinearForm_dof_map_1();
1951
 
  }
1952
 
 
1953
 
};
1954
 
 
1955
 
/// This class defines the interface for a local-to-global mapping of
1956
 
/// degrees of freedom (dofs).
1957
 
 
1958
 
class UFC_ConvectionMatrix2DBilinearForm_dof_map_2: public ufc::dof_map
1959
 
{
1960
 
private:
1961
 
 
1962
 
  unsigned int __global_dimension;
1963
 
 
1964
 
public:
1965
 
 
1966
 
  /// Constructor
1967
 
  UFC_ConvectionMatrix2DBilinearForm_dof_map_2() : ufc::dof_map()
1968
 
  {
1969
 
    __global_dimension = 0;
1970
 
  }
1971
 
 
1972
 
  /// Destructor
1973
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_dof_map_2()
1974
 
  {
1975
 
    // Do nothing
1976
 
  }
1977
 
 
1978
 
  /// Return a string identifying the dof map
1979
 
  virtual const char* signature() const
1980
 
  {
1981
 
    return "FFC dof map for Discontinuous Lagrange finite element of degree 0 on a triangle";
1982
 
  }
1983
 
 
1984
 
  /// Return true iff mesh entities of topological dimension d are needed
1985
 
  virtual bool needs_mesh_entities(unsigned int d) const
1986
 
  {
1987
 
    switch ( d )
1988
 
    {
1989
 
    case 0:
1990
 
      return false;
1991
 
      break;
1992
 
    case 1:
1993
 
      return false;
1994
 
      break;
1995
 
    case 2:
1996
 
      return true;
1997
 
      break;
1998
 
    }
1999
 
    return false;
2000
 
  }
2001
 
 
2002
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
2003
 
  virtual bool init_mesh(const ufc::mesh& m)
2004
 
  {
2005
 
    __global_dimension = m.num_entities[2];
2006
 
    return false;
2007
 
  }
2008
 
 
2009
 
  /// Initialize dof map for given cell
2010
 
  virtual void init_cell(const ufc::mesh& m,
2011
 
                         const ufc::cell& c)
2012
 
  {
2013
 
    // Do nothing
2014
 
  }
2015
 
 
2016
 
  /// Finish initialization of dof map for cells
2017
 
  virtual void init_cell_finalize()
2018
 
  {
2019
 
    // Do nothing
2020
 
  }
2021
 
 
2022
 
  /// Return the dimension of the global finite element function space
2023
 
  virtual unsigned int global_dimension() const
2024
 
  {
2025
 
    return __global_dimension;
2026
 
  }
2027
 
 
2028
 
  /// Return the dimension of the local finite element function space
2029
 
  virtual unsigned int local_dimension() const
2030
 
  {
2031
 
    return 1;
2032
 
  }
2033
 
 
2034
 
  // Return the geometric dimension of the coordinates this dof map provides
2035
 
  virtual unsigned int geometric_dimension() const
2036
 
  {
2037
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2038
 
  }
2039
 
 
2040
 
  /// Return the number of dofs on each cell facet
2041
 
  virtual unsigned int num_facet_dofs() const
2042
 
  {
2043
 
    return 0;
2044
 
  }
2045
 
 
2046
 
  /// Return the number of dofs associated with each cell entity of dimension d
2047
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
2048
 
  {
2049
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2050
 
  }
2051
 
 
2052
 
  /// Tabulate the local-to-global mapping of dofs on a cell
2053
 
  virtual void tabulate_dofs(unsigned int* dofs,
2054
 
                             const ufc::mesh& m,
2055
 
                             const ufc::cell& c) const
2056
 
  {
2057
 
    dofs[0] = c.entity_indices[2][0];
2058
 
  }
2059
 
 
2060
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
2061
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
2062
 
                                   unsigned int facet) const
2063
 
  {
2064
 
    switch ( facet )
2065
 
    {
2066
 
    case 0:
2067
 
      
2068
 
      break;
2069
 
    case 1:
2070
 
      
2071
 
      break;
2072
 
    case 2:
2073
 
      
2074
 
      break;
2075
 
    }
2076
 
  }
2077
 
 
2078
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
2079
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
2080
 
                                    unsigned int d, unsigned int i) const
2081
 
  {
2082
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2083
 
  }
2084
 
 
2085
 
  /// Tabulate the coordinates of all dofs on a cell
2086
 
  virtual void tabulate_coordinates(double** coordinates,
2087
 
                                    const ufc::cell& c) const
2088
 
  {
2089
 
    const double * const * x = c.coordinates;
2090
 
    coordinates[0][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
2091
 
    coordinates[0][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
2092
 
  }
2093
 
 
2094
 
  /// Return the number of sub dof maps (for a mixed element)
2095
 
  virtual unsigned int num_sub_dof_maps() const
2096
 
  {
2097
 
    return 1;
2098
 
  }
2099
 
 
2100
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
2101
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2102
 
  {
2103
 
    return new UFC_ConvectionMatrix2DBilinearForm_dof_map_2();
2104
 
  }
2105
 
 
2106
 
};
2107
 
 
2108
 
/// This class defines the interface for a local-to-global mapping of
2109
 
/// degrees of freedom (dofs).
2110
 
 
2111
 
class UFC_ConvectionMatrix2DBilinearForm_dof_map_3: public ufc::dof_map
2112
 
{
2113
 
private:
2114
 
 
2115
 
  unsigned int __global_dimension;
2116
 
 
2117
 
public:
2118
 
 
2119
 
  /// Constructor
2120
 
  UFC_ConvectionMatrix2DBilinearForm_dof_map_3() : ufc::dof_map()
2121
 
  {
2122
 
    __global_dimension = 0;
2123
 
  }
2124
 
 
2125
 
  /// Destructor
2126
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_dof_map_3()
2127
 
  {
2128
 
    // Do nothing
2129
 
  }
2130
 
 
2131
 
  /// Return a string identifying the dof map
2132
 
  virtual const char* signature() const
2133
 
  {
2134
 
    return "FFC dof map for Discontinuous Lagrange finite element of degree 0 on a triangle";
2135
 
  }
2136
 
 
2137
 
  /// Return true iff mesh entities of topological dimension d are needed
2138
 
  virtual bool needs_mesh_entities(unsigned int d) const
2139
 
  {
2140
 
    switch ( d )
2141
 
    {
2142
 
    case 0:
2143
 
      return false;
2144
 
      break;
2145
 
    case 1:
2146
 
      return false;
2147
 
      break;
2148
 
    case 2:
2149
 
      return true;
2150
 
      break;
2151
 
    }
2152
 
    return false;
2153
 
  }
2154
 
 
2155
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
2156
 
  virtual bool init_mesh(const ufc::mesh& m)
2157
 
  {
2158
 
    __global_dimension = m.num_entities[2];
2159
 
    return false;
2160
 
  }
2161
 
 
2162
 
  /// Initialize dof map for given cell
2163
 
  virtual void init_cell(const ufc::mesh& m,
2164
 
                         const ufc::cell& c)
2165
 
  {
2166
 
    // Do nothing
2167
 
  }
2168
 
 
2169
 
  /// Finish initialization of dof map for cells
2170
 
  virtual void init_cell_finalize()
2171
 
  {
2172
 
    // Do nothing
2173
 
  }
2174
 
 
2175
 
  /// Return the dimension of the global finite element function space
2176
 
  virtual unsigned int global_dimension() const
2177
 
  {
2178
 
    return __global_dimension;
2179
 
  }
2180
 
 
2181
 
  /// Return the dimension of the local finite element function space
2182
 
  virtual unsigned int local_dimension() const
2183
 
  {
2184
 
    return 1;
2185
 
  }
2186
 
 
2187
 
  // Return the geometric dimension of the coordinates this dof map provides
2188
 
  virtual unsigned int geometric_dimension() const
2189
 
  {
2190
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2191
 
  }
2192
 
 
2193
 
  /// Return the number of dofs on each cell facet
2194
 
  virtual unsigned int num_facet_dofs() const
2195
 
  {
2196
 
    return 0;
2197
 
  }
2198
 
 
2199
 
  /// Return the number of dofs associated with each cell entity of dimension d
2200
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
2201
 
  {
2202
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2203
 
  }
2204
 
 
2205
 
  /// Tabulate the local-to-global mapping of dofs on a cell
2206
 
  virtual void tabulate_dofs(unsigned int* dofs,
2207
 
                             const ufc::mesh& m,
2208
 
                             const ufc::cell& c) const
2209
 
  {
2210
 
    dofs[0] = c.entity_indices[2][0];
2211
 
  }
2212
 
 
2213
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
2214
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
2215
 
                                   unsigned int facet) const
2216
 
  {
2217
 
    switch ( facet )
2218
 
    {
2219
 
    case 0:
2220
 
      
2221
 
      break;
2222
 
    case 1:
2223
 
      
2224
 
      break;
2225
 
    case 2:
2226
 
      
2227
 
      break;
2228
 
    }
2229
 
  }
2230
 
 
2231
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
2232
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
2233
 
                                    unsigned int d, unsigned int i) const
2234
 
  {
2235
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2236
 
  }
2237
 
 
2238
 
  /// Tabulate the coordinates of all dofs on a cell
2239
 
  virtual void tabulate_coordinates(double** coordinates,
2240
 
                                    const ufc::cell& c) const
2241
 
  {
2242
 
    const double * const * x = c.coordinates;
2243
 
    coordinates[0][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
2244
 
    coordinates[0][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
2245
 
  }
2246
 
 
2247
 
  /// Return the number of sub dof maps (for a mixed element)
2248
 
  virtual unsigned int num_sub_dof_maps() const
2249
 
  {
2250
 
    return 1;
2251
 
  }
2252
 
 
2253
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
2254
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2255
 
  {
2256
 
    return new UFC_ConvectionMatrix2DBilinearForm_dof_map_3();
2257
 
  }
2258
 
 
2259
 
};
2260
 
 
2261
 
/// This class defines the interface for the tabulation of the cell
2262
 
/// tensor corresponding to the local contribution to a form from
2263
 
/// the integral over a cell.
2264
 
 
2265
 
class UFC_ConvectionMatrix2DBilinearForm_cell_integral_0: public ufc::cell_integral
2266
 
{
2267
 
public:
2268
 
 
2269
 
  /// Constructor
2270
 
  UFC_ConvectionMatrix2DBilinearForm_cell_integral_0() : ufc::cell_integral()
2271
 
  {
2272
 
    // Do nothing
2273
 
  }
2274
 
 
2275
 
  /// Destructor
2276
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm_cell_integral_0()
2277
 
  {
2278
 
    // Do nothing
2279
 
  }
2280
 
 
2281
 
  /// Tabulate the tensor for the contribution from a local cell
2282
 
  virtual void tabulate_tensor(double* A,
2283
 
                               const double * const * w,
2284
 
                               const ufc::cell& c) const
2285
 
  {
2286
 
    // Extract vertex coordinates
2287
 
    const double * const * x = c.coordinates;
2288
 
    
2289
 
    // Compute Jacobian of affine map from reference cell
2290
 
    const double J_00 = x[1][0] - x[0][0];
2291
 
    const double J_01 = x[2][0] - x[0][0];
2292
 
    const double J_10 = x[1][1] - x[0][1];
2293
 
    const double J_11 = x[2][1] - x[0][1];
2294
 
      
2295
 
    // Compute determinant of Jacobian
2296
 
    double detJ = J_00*J_11 - J_01*J_10;
2297
 
      
2298
 
    // Compute inverse of Jacobian
2299
 
    const double Jinv_00 =  J_11 / detJ;
2300
 
    const double Jinv_01 = -J_01 / detJ;
2301
 
    const double Jinv_10 = -J_10 / detJ;
2302
 
    const double Jinv_11 =  J_00 / detJ;
2303
 
    
2304
 
    // Set scale factor
2305
 
    const double det = std::abs(detJ);
2306
 
    
2307
 
    // Compute coefficients
2308
 
    const double c0_0_0_0 = w[0][0];
2309
 
    const double c1_1_0_0 = w[1][0];
2310
 
    
2311
 
    // Compute geometry tensors
2312
 
    const double G0_0_0 = det*(c0_0_0_0*Jinv_00 + c1_1_0_0*Jinv_01);
2313
 
    const double G0_0_1 = det*(c0_0_0_0*Jinv_10 + c1_1_0_0*Jinv_11);
2314
 
    
2315
 
    // Compute element tensor
2316
 
    A[0] = -0.166666666666667*G0_0_0 - 0.166666666666667*G0_0_1;
2317
 
    A[1] = 0.166666666666667*G0_0_0;
2318
 
    A[2] = 0.166666666666667*G0_0_1;
2319
 
    A[3] = -0.166666666666667*G0_0_0 - 0.166666666666667*G0_0_1;
2320
 
    A[4] = 0.166666666666667*G0_0_0;
2321
 
    A[5] = 0.166666666666667*G0_0_1;
2322
 
    A[6] = -0.166666666666667*G0_0_0 - 0.166666666666667*G0_0_1;
2323
 
    A[7] = 0.166666666666667*G0_0_0;
2324
 
    A[8] = 0.166666666666667*G0_0_1;
2325
 
  }
2326
 
 
2327
 
};
2328
 
 
2329
 
/// This class defines the interface for the assembly of the global
2330
 
/// tensor corresponding to a form with r + n arguments, that is, a
2331
 
/// mapping
2332
 
///
2333
 
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
2334
 
///
2335
 
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
2336
 
/// global tensor A is defined by
2337
 
///
2338
 
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
2339
 
///
2340
 
/// where each argument Vj represents the application to the
2341
 
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
2342
 
/// fixed functions (coefficients).
2343
 
 
2344
 
class UFC_ConvectionMatrix2DBilinearForm: public ufc::form
2345
 
{
2346
 
public:
2347
 
 
2348
 
  /// Constructor
2349
 
  UFC_ConvectionMatrix2DBilinearForm() : ufc::form()
2350
 
  {
2351
 
    // Do nothing
2352
 
  }
2353
 
 
2354
 
  /// Destructor
2355
 
  virtual ~UFC_ConvectionMatrix2DBilinearForm()
2356
 
  {
2357
 
    // Do nothing
2358
 
  }
2359
 
 
2360
 
  /// Return a string identifying the form
2361
 
  virtual const char* signature() const
2362
 
  {
2363
 
    return "w0_a0(dXa1/dx0) | vi0*va0*((d/dXa1)vi1)*dX(0) + w1_a0(dXa1/dx1) | vi0*va0*((d/dXa1)vi1)*dX(0)";
2364
 
  }
2365
 
 
2366
 
  /// Return the rank of the global tensor (r)
2367
 
  virtual unsigned int rank() const
2368
 
  {
2369
 
    return 2;
2370
 
  }
2371
 
 
2372
 
  /// Return the number of coefficients (n)
2373
 
  virtual unsigned int num_coefficients() const
2374
 
  {
2375
 
    return 2;
2376
 
  }
2377
 
 
2378
 
  /// Return the number of cell integrals
2379
 
  virtual unsigned int num_cell_integrals() const
2380
 
  {
2381
 
    return 1;
2382
 
  }
2383
 
  
2384
 
  /// Return the number of exterior facet integrals
2385
 
  virtual unsigned int num_exterior_facet_integrals() const
2386
 
  {
2387
 
    return 0;
2388
 
  }
2389
 
  
2390
 
  /// Return the number of interior facet integrals
2391
 
  virtual unsigned int num_interior_facet_integrals() const
2392
 
  {
2393
 
    return 0;
2394
 
  }
2395
 
    
2396
 
  /// Create a new finite element for argument function i
2397
 
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
2398
 
  {
2399
 
    switch ( i )
2400
 
    {
2401
 
    case 0:
2402
 
      return new UFC_ConvectionMatrix2DBilinearForm_finite_element_0();
2403
 
      break;
2404
 
    case 1:
2405
 
      return new UFC_ConvectionMatrix2DBilinearForm_finite_element_1();
2406
 
      break;
2407
 
    case 2:
2408
 
      return new UFC_ConvectionMatrix2DBilinearForm_finite_element_2();
2409
 
      break;
2410
 
    case 3:
2411
 
      return new UFC_ConvectionMatrix2DBilinearForm_finite_element_3();
2412
 
      break;
2413
 
    }
2414
 
    return 0;
2415
 
  }
2416
 
  
2417
 
  /// Create a new dof map for argument function i
2418
 
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
2419
 
  {
2420
 
    switch ( i )
2421
 
    {
2422
 
    case 0:
2423
 
      return new UFC_ConvectionMatrix2DBilinearForm_dof_map_0();
2424
 
      break;
2425
 
    case 1:
2426
 
      return new UFC_ConvectionMatrix2DBilinearForm_dof_map_1();
2427
 
      break;
2428
 
    case 2:
2429
 
      return new UFC_ConvectionMatrix2DBilinearForm_dof_map_2();
2430
 
      break;
2431
 
    case 3:
2432
 
      return new UFC_ConvectionMatrix2DBilinearForm_dof_map_3();
2433
 
      break;
2434
 
    }
2435
 
    return 0;
2436
 
  }
2437
 
 
2438
 
  /// Create a new cell integral on sub domain i
2439
 
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
2440
 
  {
2441
 
    return new UFC_ConvectionMatrix2DBilinearForm_cell_integral_0();
2442
 
  }
2443
 
 
2444
 
  /// Create a new exterior facet integral on sub domain i
2445
 
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
2446
 
  {
2447
 
    return 0;
2448
 
  }
2449
 
 
2450
 
  /// Create a new interior facet integral on sub domain i
2451
 
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
2452
 
  {
2453
 
    return 0;
2454
 
  }
2455
 
 
2456
 
};
2457
 
 
2458
 
// DOLFIN wrappers
2459
 
 
2460
 
#include <dolfin/Form.h>
2461
 
 
2462
 
class ConvectionMatrix2DBilinearForm : public dolfin::Form
2463
 
{
2464
 
public:
2465
 
 
2466
 
  ConvectionMatrix2DBilinearForm(dolfin::Function& w0, dolfin::Function& w1) : dolfin::Form()
2467
 
  {
2468
 
    __coefficients.push_back(&w0);
2469
 
    __coefficients.push_back(&w1);
2470
 
  }
2471
 
 
2472
 
  /// Return UFC form
2473
 
  virtual const ufc::form& form() const
2474
 
  {
2475
 
    return __form;
2476
 
  }
2477
 
  
2478
 
  /// Return array of coefficients
2479
 
  virtual const dolfin::Array<dolfin::Function*>& coefficients() const
2480
 
  {
2481
 
    return __coefficients;
2482
 
  }
2483
 
 
2484
 
private:
2485
 
 
2486
 
  // UFC form
2487
 
  UFC_ConvectionMatrix2DBilinearForm __form;
2488
 
 
2489
 
  /// Array of coefficients
2490
 
  dolfin::Array<dolfin::Function*> __coefficients;
2491
 
 
2492
 
};
2493
 
 
2494
 
#endif