← Back to branch summary

~ubuntu-branches/ubuntu/vivid/ffc/vivid-proposed

« back to all changes in this revision

Viewing changes to test/regression/references/r_quadrature_O/Heat.h

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2011-11-29 11:38:38 UTC
  • mfrom: (1.1.11)
  • Revision ID: package-import@ubuntu.com-20111129113838-vx815cxjdf50zuq3
Tags: 1.0-rc1-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/regression/references/r_quadrature_O/Heat.h
 
1
// This code conforms with the UFC specification version 2.0.3
 
2
// and was automatically generated by FFC version 1.0-beta2+.
 
3
// 
 
4
// This code was generated with the following parameters:
 
5
// 
 
6
//   cache_dir:                      ''
 
7
//   convert_exceptions_to_warnings: True
 
8
//   cpp_optimize:                   False
 
9
//   cpp_optimize_flags:             '-O2'
 
10
//   epsilon:                        1e-14
 
11
//   error_control:                  False
 
12
//   form_postfix:                   True
 
13
//   format:                         'ufc'
 
14
//   log_level:                      20
 
15
//   log_prefix:                     ''
 
16
//   optimize:                       True
 
17
//   output_dir:                     '.'
 
18
//   precision:                      '8'
 
19
//   quadrature_degree:              'auto'
 
20
//   quadrature_rule:                'auto'
 
21
//   representation:                 'quadrature'
 
22
//   split:                          False
 
23
//   swig_binary:                    'swig'
 
24
//   swig_path:                      ''
 
25
 
 
26
#ifndef __HEAT_H
 
27
#define __HEAT_H
 
28
 
 
29
#include <cmath>
 
30
#include <stdexcept>
 
31
#include <fstream>
 
32
#include <ufc.h>
 
33
 
 
34
/// This class defines the interface for a finite element.
 
35
 
 
36
class heat_finite_element_0: public ufc::finite_element
 
37
{
 
38
public:
 
39
 
 
40
  /// Constructor
 
41
  heat_finite_element_0() : ufc::finite_element()
 
42
  {
 
43
    // Do nothing
 
44
  }
 
45
 
 
46
  /// Destructor
 
47
  virtual ~heat_finite_element_0()
 
48
  {
 
49
    // Do nothing
 
50
  }
 
51
 
 
52
  /// Return a string identifying the finite element
 
53
  virtual const char* signature() const
 
54
  {
 
55
    return "FiniteElement('Real', Cell('triangle', Space(2)), 0, None)";
 
56
  }
 
57
 
 
58
  /// Return the cell shape
 
59
  virtual ufc::shape cell_shape() const
 
60
  {
 
61
    return ufc::triangle;
 
62
  }
 
63
 
 
64
  /// Return the topological dimension of the cell shape
 
65
  virtual unsigned int topological_dimension() const
 
66
  {
 
67
    return 2;
 
68
  }
 
69
 
 
70
  /// Return the geometric dimension of the cell shape
 
71
  virtual unsigned int geometric_dimension() const
 
72
  {
 
73
    return 2;
 
74
  }
 
75
 
 
76
  /// Return the dimension of the finite element function space
 
77
  virtual unsigned int space_dimension() const
 
78
  {
 
79
    return 1;
 
80
  }
 
81
 
 
82
  /// Return the rank of the value space
 
83
  virtual unsigned int value_rank() const
 
84
  {
 
85
    return 0;
 
86
  }
 
87
 
 
88
  /// Return the dimension of the value space for axis i
 
89
  virtual unsigned int value_dimension(unsigned int i) const
 
90
  {
 
91
    return 1;
 
92
  }
 
93
 
 
94
  /// Evaluate basis function i at given point in cell
 
95
  virtual void evaluate_basis(unsigned int i,
 
96
                              double* values,
 
97
                              const double* coordinates,
 
98
                              const ufc::cell& c) const
 
99
  {
 
100
    // Extract vertex coordinates
 
101
    
 
102
    // Compute Jacobian of affine map from reference cell
 
103
    
 
104
    // Compute determinant of Jacobian
 
105
    
 
106
    // Compute inverse of Jacobian
 
107
    
 
108
    // Compute constants
 
109
    
 
110
    // Get coordinates and map to the reference (FIAT) element
 
111
    
 
112
    // Reset values.
 
113
    *values = 0.0;
 
114
    
 
115
    // Array of basisvalues.
 
116
    double basisvalues[1] = {0.0};
 
117
    
 
118
    // Declare helper variables.
 
119
    
 
120
    // Compute basisvalues.
 
121
    basisvalues[0] = 1.0;
 
122
    
 
123
    // Table(s) of coefficients.
 
124
    static const double coefficients0[1] = \
 
125
    {1.0};
 
126
    
 
127
    // Compute value(s).
 
128
    for (unsigned int r = 0; r < 1; r++)
 
129
    {
 
130
      *values += coefficients0[r]*basisvalues[r];
 
131
    }// end loop over 'r'
 
132
  }
 
133
 
 
134
  /// Evaluate all basis functions at given point in cell
 
135
  virtual void evaluate_basis_all(double* values,
 
136
                                  const double* coordinates,
 
137
                                  const ufc::cell& c) const
 
138
  {
 
139
    // Element is constant, calling evaluate_basis.
 
140
    evaluate_basis(0, values, coordinates, c);
 
141
  }
 
142
 
 
143
  /// Evaluate order n derivatives of basis function i at given point in cell
 
144
  virtual void evaluate_basis_derivatives(unsigned int i,
 
145
                                          unsigned int n,
 
146
                                          double* values,
 
147
                                          const double* coordinates,
 
148
                                          const ufc::cell& c) const
 
149
  {
 
150
    // Extract vertex coordinates
 
151
    const double * const * x = c.coordinates;
 
152
    
 
153
    // Compute Jacobian of affine map from reference cell
 
154
    const double J_00 = x[1][0] - x[0][0];
 
155
    const double J_01 = x[2][0] - x[0][0];
 
156
    const double J_10 = x[1][1] - x[0][1];
 
157
    const double J_11 = x[2][1] - x[0][1];
 
158
    
 
159
    // Compute determinant of Jacobian
 
160
    double detJ = J_00*J_11 - J_01*J_10;
 
161
    
 
162
    // Compute inverse of Jacobian
 
163
    const double K_00 =  J_11 / detJ;
 
164
    const double K_01 = -J_01 / detJ;
 
165
    const double K_10 = -J_10 / detJ;
 
166
    const double K_11 =  J_00 / detJ;
 
167
    
 
168
    // Compute constants
 
169
    
 
170
    // Get coordinates and map to the reference (FIAT) element
 
171
    
 
172
    // Compute number of derivatives.
 
173
    unsigned int num_derivatives = 1;
 
174
    for (unsigned int r = 0; r < n; r++)
 
175
    {
 
176
      num_derivatives *= 2;
 
177
    }// end loop over 'r'
 
178
    
 
179
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
180
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
181
    for (unsigned int row = 0; row < num_derivatives; row++)
 
182
    {
 
183
      combinations[row] = new unsigned int [n];
 
184
      for (unsigned int col = 0; col < n; col++)
 
185
        combinations[row][col] = 0;
 
186
    }
 
187
    
 
188
    // Generate combinations of derivatives
 
189
    for (unsigned int row = 1; row < num_derivatives; row++)
 
190
    {
 
191
      for (unsigned int num = 0; num < row; num++)
 
192
      {
 
193
        for (unsigned int col = n-1; col+1 > 0; col--)
 
194
        {
 
195
          if (combinations[row][col] + 1 > 1)
 
196
            combinations[row][col] = 0;
 
197
          else
 
198
          {
 
199
            combinations[row][col] += 1;
 
200
            break;
 
201
          }
 
202
        }
 
203
      }
 
204
    }
 
205
    
 
206
    // Compute inverse of Jacobian
 
207
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
 
208
    
 
209
    // Declare transformation matrix
 
210
    // Declare pointer to two dimensional array and initialise
 
211
    double **transform = new double *[num_derivatives];
 
212
    
 
213
    for (unsigned int j = 0; j < num_derivatives; j++)
 
214
    {
 
215
      transform[j] = new double [num_derivatives];
 
216
      for (unsigned int k = 0; k < num_derivatives; k++)
 
217
        transform[j][k] = 1;
 
218
    }
 
219
    
 
220
    // Construct transformation matrix
 
221
    for (unsigned int row = 0; row < num_derivatives; row++)
 
222
    {
 
223
      for (unsigned int col = 0; col < num_derivatives; col++)
 
224
      {
 
225
        for (unsigned int k = 0; k < n; k++)
 
226
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
227
      }
 
228
    }
 
229
    
 
230
    // Reset values. Assuming that values is always an array.
 
231
    for (unsigned int r = 0; r < num_derivatives; r++)
 
232
    {
 
233
      values[r] = 0.0;
 
234
    }// end loop over 'r'
 
235
    
 
236
    
 
237
    // Array of basisvalues.
 
238
    double basisvalues[1] = {0.0};
 
239
    
 
240
    // Declare helper variables.
 
241
    
 
242
    // Compute basisvalues.
 
243
    basisvalues[0] = 1.0;
 
244
    
 
245
    // Table(s) of coefficients.
 
246
    static const double coefficients0[1] = \
 
247
    {1.0};
 
248
    
 
249
    // Tables of derivatives of the polynomial base (transpose).
 
250
    static const double dmats0[1][1] = \
 
251
    {{0.0}};
 
252
    
 
253
    static const double dmats1[1][1] = \
 
254
    {{0.0}};
 
255
    
 
256
    // Compute reference derivatives.
 
257
    // Declare pointer to array of derivatives on FIAT element.
 
258
    double *derivatives = new double[num_derivatives];
 
259
    for (unsigned int r = 0; r < num_derivatives; r++)
 
260
    {
 
261
      derivatives[r] = 0.0;
 
262
    }// end loop over 'r'
 
263
    
 
264
    // Declare derivative matrix (of polynomial basis).
 
265
    double dmats[1][1] = \
 
266
    {{1.0}};
 
267
    
 
268
    // Declare (auxiliary) derivative matrix (of polynomial basis).
 
269
    double dmats_old[1][1] = \
 
270
    {{1.0}};
 
271
    
 
272
    // Loop possible derivatives.
 
273
    for (unsigned int r = 0; r < num_derivatives; r++)
 
274
    {
 
275
      // Resetting dmats values to compute next derivative.
 
276
      for (unsigned int t = 0; t < 1; t++)
 
277
      {
 
278
        for (unsigned int u = 0; u < 1; u++)
 
279
        {
 
280
          dmats[t][u] = 0.0;
 
281
          if (t == u)
 
282
          {
 
283
          dmats[t][u] = 1.0;
 
284
          }
 
285
          
 
286
        }// end loop over 'u'
 
287
      }// end loop over 't'
 
288
      
 
289
      // Looping derivative order to generate dmats.
 
290
      for (unsigned int s = 0; s < n; s++)
 
291
      {
 
292
        // Updating dmats_old with new values and resetting dmats.
 
293
        for (unsigned int t = 0; t < 1; t++)
 
294
        {
 
295
          for (unsigned int u = 0; u < 1; u++)
 
296
          {
 
297
            dmats_old[t][u] = dmats[t][u];
 
298
            dmats[t][u] = 0.0;
 
299
          }// end loop over 'u'
 
300
        }// end loop over 't'
 
301
        
 
302
        // Update dmats using an inner product.
 
303
        if (combinations[r][s] == 0)
 
304
        {
 
305
        for (unsigned int t = 0; t < 1; t++)
 
306
        {
 
307
          for (unsigned int u = 0; u < 1; u++)
 
308
          {
 
309
            for (unsigned int tu = 0; tu < 1; tu++)
 
310
            {
 
311
              dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
312
            }// end loop over 'tu'
 
313
          }// end loop over 'u'
 
314
        }// end loop over 't'
 
315
        }
 
316
        
 
317
        if (combinations[r][s] == 1)
 
318
        {
 
319
        for (unsigned int t = 0; t < 1; t++)
 
320
        {
 
321
          for (unsigned int u = 0; u < 1; u++)
 
322
          {
 
323
            for (unsigned int tu = 0; tu < 1; tu++)
 
324
            {
 
325
              dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
326
            }// end loop over 'tu'
 
327
          }// end loop over 'u'
 
328
        }// end loop over 't'
 
329
        }
 
330
        
 
331
      }// end loop over 's'
 
332
      for (unsigned int s = 0; s < 1; s++)
 
333
      {
 
334
        for (unsigned int t = 0; t < 1; t++)
 
335
        {
 
336
          derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
337
        }// end loop over 't'
 
338
      }// end loop over 's'
 
339
    }// end loop over 'r'
 
340
    
 
341
    // Transform derivatives back to physical element
 
342
    for (unsigned int r = 0; r < num_derivatives; r++)
 
343
    {
 
344
      for (unsigned int s = 0; s < num_derivatives; s++)
 
345
      {
 
346
        values[r] += transform[r][s]*derivatives[s];
 
347
      }// end loop over 's'
 
348
    }// end loop over 'r'
 
349
    
 
350
    // Delete pointer to array of derivatives on FIAT element
 
351
    delete [] derivatives;
 
352
    
 
353
    // Delete pointer to array of combinations of derivatives and transform
 
354
    for (unsigned int r = 0; r < num_derivatives; r++)
 
355
    {
 
356
      delete [] combinations[r];
 
357
    }// end loop over 'r'
 
358
    delete [] combinations;
 
359
    for (unsigned int r = 0; r < num_derivatives; r++)
 
360
    {
 
361
      delete [] transform[r];
 
362
    }// end loop over 'r'
 
363
    delete [] transform;
 
364
  }
 
365
 
 
366
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
367
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
368
                                              double* values,
 
369
                                              const double* coordinates,
 
370
                                              const ufc::cell& c) const
 
371
  {
 
372
    // Element is constant, calling evaluate_basis_derivatives.
 
373
    evaluate_basis_derivatives(0, n, values, coordinates, c);
 
374
  }
 
375
 
 
376
  /// Evaluate linear functional for dof i on the function f
 
377
  virtual double evaluate_dof(unsigned int i,
 
378
                              const ufc::function& f,
 
379
                              const ufc::cell& c) const
 
380
  {
 
381
    // Declare variables for result of evaluation.
 
382
    double vals[1];
 
383
    
 
384
    // Declare variable for physical coordinates.
 
385
    double y[2];
 
386
    const double * const * x = c.coordinates;
 
387
    switch (i)
 
388
    {
 
389
    case 0:
 
390
      {
 
391
        y[0] = 0.33333333*x[0][0] + 0.33333333*x[1][0] + 0.33333333*x[2][0];
 
392
      y[1] = 0.33333333*x[0][1] + 0.33333333*x[1][1] + 0.33333333*x[2][1];
 
393
      f.evaluate(vals, y, c);
 
394
      return vals[0];
 
395
        break;
 
396
      }
 
397
    }
 
398
    
 
399
    return 0.0;
 
400
  }
 
401
 
 
402
  /// Evaluate linear functionals for all dofs on the function f
 
403
  virtual void evaluate_dofs(double* values,
 
404
                             const ufc::function& f,
 
405
                             const ufc::cell& c) const
 
406
  {
 
407
    // Declare variables for result of evaluation.
 
408
    double vals[1];
 
409
    
 
410
    // Declare variable for physical coordinates.
 
411
    double y[2];
 
412
    const double * const * x = c.coordinates;
 
413
    y[0] = 0.33333333*x[0][0] + 0.33333333*x[1][0] + 0.33333333*x[2][0];
 
414
    y[1] = 0.33333333*x[0][1] + 0.33333333*x[1][1] + 0.33333333*x[2][1];
 
415
    f.evaluate(vals, y, c);
 
416
    values[0] = vals[0];
 
417
  }
 
418
 
 
419
  /// Interpolate vertex values from dof values
 
420
  virtual void interpolate_vertex_values(double* vertex_values,
 
421
                                         const double* dof_values,
 
422
                                         const ufc::cell& c) const
 
423
  {
 
424
    // Evaluate function and change variables
 
425
    vertex_values[0] = dof_values[0];
 
426
    vertex_values[1] = dof_values[0];
 
427
    vertex_values[2] = dof_values[0];
 
428
  }
 
429
 
 
430
  /// Map coordinate xhat from reference cell to coordinate x in cell
 
431
  virtual void map_from_reference_cell(double* x,
 
432
                                       const double* xhat,
 
433
                                       const ufc::cell& c) const
 
434
  {
 
435
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
436
  }
 
437
 
 
438
  /// Map from coordinate x in cell to coordinate xhat in reference cell
 
439
  virtual void map_to_reference_cell(double* xhat,
 
440
                                     const double* x,
 
441
                                     const ufc::cell& c) const
 
442
  {
 
443
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
444
  }
 
445
 
 
446
  /// Return the number of sub elements (for a mixed element)
 
447
  virtual unsigned int num_sub_elements() const
 
448
  {
 
449
    return 0;
 
450
  }
 
451
 
 
452
  /// Create a new finite element for sub element i (for a mixed element)
 
453
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
454
  {
 
455
    return 0;
 
456
  }
 
457
 
 
458
  /// Create a new class instance
 
459
  virtual ufc::finite_element* create() const
 
460
  {
 
461
    return new heat_finite_element_0();
 
462
  }
 
463
 
 
464
};
 
465
 
 
466
/// This class defines the interface for a finite element.
 
467
 
 
468
class heat_finite_element_1: public ufc::finite_element
 
469
{
 
470
public:
 
471
 
 
472
  /// Constructor
 
473
  heat_finite_element_1() : ufc::finite_element()
 
474
  {
 
475
    // Do nothing
 
476
  }
 
477
 
 
478
  /// Destructor
 
479
  virtual ~heat_finite_element_1()
 
480
  {
 
481
    // Do nothing
 
482
  }
 
483
 
 
484
  /// Return a string identifying the finite element
 
485
  virtual const char* signature() const
 
486
  {
 
487
    return "FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)";
 
488
  }
 
489
 
 
490
  /// Return the cell shape
 
491
  virtual ufc::shape cell_shape() const
 
492
  {
 
493
    return ufc::triangle;
 
494
  }
 
495
 
 
496
  /// Return the topological dimension of the cell shape
 
497
  virtual unsigned int topological_dimension() const
 
498
  {
 
499
    return 2;
 
500
  }
 
501
 
 
502
  /// Return the geometric dimension of the cell shape
 
503
  virtual unsigned int geometric_dimension() const
 
504
  {
 
505
    return 2;
 
506
  }
 
507
 
 
508
  /// Return the dimension of the finite element function space
 
509
  virtual unsigned int space_dimension() const
 
510
  {
 
511
    return 3;
 
512
  }
 
513
 
 
514
  /// Return the rank of the value space
 
515
  virtual unsigned int value_rank() const
 
516
  {
 
517
    return 0;
 
518
  }
 
519
 
 
520
  /// Return the dimension of the value space for axis i
 
521
  virtual unsigned int value_dimension(unsigned int i) const
 
522
  {
 
523
    return 1;
 
524
  }
 
525
 
 
526
  /// Evaluate basis function i at given point in cell
 
527
  virtual void evaluate_basis(unsigned int i,
 
528
                              double* values,
 
529
                              const double* coordinates,
 
530
                              const ufc::cell& c) const
 
531
  {
 
532
    // Extract vertex coordinates
 
533
    const double * const * x = c.coordinates;
 
534
    
 
535
    // Compute Jacobian of affine map from reference cell
 
536
    const double J_00 = x[1][0] - x[0][0];
 
537
    const double J_01 = x[2][0] - x[0][0];
 
538
    const double J_10 = x[1][1] - x[0][1];
 
539
    const double J_11 = x[2][1] - x[0][1];
 
540
    
 
541
    // Compute determinant of Jacobian
 
542
    double detJ = J_00*J_11 - J_01*J_10;
 
543
    
 
544
    // Compute inverse of Jacobian
 
545
    
 
546
    // Compute constants
 
547
    const double C0 = x[1][0] + x[2][0];
 
548
    const double C1 = x[1][1] + x[2][1];
 
549
    
 
550
    // Get coordinates and map to the reference (FIAT) element
 
551
    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
 
552
    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
 
553
    
 
554
    // Reset values.
 
555
    *values = 0.0;
 
556
    switch (i)
 
557
    {
 
558
    case 0:
 
559
      {
 
560
        
 
561
      // Array of basisvalues.
 
562
      double basisvalues[3] = {0.0, 0.0, 0.0};
 
563
      
 
564
      // Declare helper variables.
 
565
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
566
      
 
567
      // Compute basisvalues.
 
568
      basisvalues[0] = 1.0;
 
569
      basisvalues[1] = tmp0;
 
570
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
571
      basisvalues[0] *= std::sqrt(0.5);
 
572
      basisvalues[2] *= std::sqrt(1.0);
 
573
      basisvalues[1] *= std::sqrt(3.0);
 
574
      
 
575
      // Table(s) of coefficients.
 
576
      static const double coefficients0[3] = \
 
577
      {0.47140452, -0.28867513, -0.16666667};
 
578
      
 
579
      // Compute value(s).
 
580
      for (unsigned int r = 0; r < 3; r++)
 
581
      {
 
582
        *values += coefficients0[r]*basisvalues[r];
 
583
      }// end loop over 'r'
 
584
        break;
 
585
      }
 
586
    case 1:
 
587
      {
 
588
        
 
589
      // Array of basisvalues.
 
590
      double basisvalues[3] = {0.0, 0.0, 0.0};
 
591
      
 
592
      // Declare helper variables.
 
593
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
594
      
 
595
      // Compute basisvalues.
 
596
      basisvalues[0] = 1.0;
 
597
      basisvalues[1] = tmp0;
 
598
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
599
      basisvalues[0] *= std::sqrt(0.5);
 
600
      basisvalues[2] *= std::sqrt(1.0);
 
601
      basisvalues[1] *= std::sqrt(3.0);
 
602
      
 
603
      // Table(s) of coefficients.
 
604
      static const double coefficients0[3] = \
 
605
      {0.47140452, 0.28867513, -0.16666667};
 
606
      
 
607
      // Compute value(s).
 
608
      for (unsigned int r = 0; r < 3; r++)
 
609
      {
 
610
        *values += coefficients0[r]*basisvalues[r];
 
611
      }// end loop over 'r'
 
612
        break;
 
613
      }
 
614
    case 2:
 
615
      {
 
616
        
 
617
      // Array of basisvalues.
 
618
      double basisvalues[3] = {0.0, 0.0, 0.0};
 
619
      
 
620
      // Declare helper variables.
 
621
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
622
      
 
623
      // Compute basisvalues.
 
624
      basisvalues[0] = 1.0;
 
625
      basisvalues[1] = tmp0;
 
626
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
627
      basisvalues[0] *= std::sqrt(0.5);
 
628
      basisvalues[2] *= std::sqrt(1.0);
 
629
      basisvalues[1] *= std::sqrt(3.0);
 
630
      
 
631
      // Table(s) of coefficients.
 
632
      static const double coefficients0[3] = \
 
633
      {0.47140452, 0.0, 0.33333333};
 
634
      
 
635
      // Compute value(s).
 
636
      for (unsigned int r = 0; r < 3; r++)
 
637
      {
 
638
        *values += coefficients0[r]*basisvalues[r];
 
639
      }// end loop over 'r'
 
640
        break;
 
641
      }
 
642
    }
 
643
    
 
644
  }
 
645
 
 
646
  /// Evaluate all basis functions at given point in cell
 
647
  virtual void evaluate_basis_all(double* values,
 
648
                                  const double* coordinates,
 
649
                                  const ufc::cell& c) const
 
650
  {
 
651
    // Helper variable to hold values of a single dof.
 
652
    double dof_values = 0.0;
 
653
    
 
654
    // Loop dofs and call evaluate_basis.
 
655
    for (unsigned int r = 0; r < 3; r++)
 
656
    {
 
657
      evaluate_basis(r, &dof_values, coordinates, c);
 
658
      values[r] = dof_values;
 
659
    }// end loop over 'r'
 
660
  }
 
661
 
 
662
  /// Evaluate order n derivatives of basis function i at given point in cell
 
663
  virtual void evaluate_basis_derivatives(unsigned int i,
 
664
                                          unsigned int n,
 
665
                                          double* values,
 
666
                                          const double* coordinates,
 
667
                                          const ufc::cell& c) const
 
668
  {
 
669
    // Extract vertex coordinates
 
670
    const double * const * x = c.coordinates;
 
671
    
 
672
    // Compute Jacobian of affine map from reference cell
 
673
    const double J_00 = x[1][0] - x[0][0];
 
674
    const double J_01 = x[2][0] - x[0][0];
 
675
    const double J_10 = x[1][1] - x[0][1];
 
676
    const double J_11 = x[2][1] - x[0][1];
 
677
    
 
678
    // Compute determinant of Jacobian
 
679
    double detJ = J_00*J_11 - J_01*J_10;
 
680
    
 
681
    // Compute inverse of Jacobian
 
682
    const double K_00 =  J_11 / detJ;
 
683
    const double K_01 = -J_01 / detJ;
 
684
    const double K_10 = -J_10 / detJ;
 
685
    const double K_11 =  J_00 / detJ;
 
686
    
 
687
    // Compute constants
 
688
    const double C0 = x[1][0] + x[2][0];
 
689
    const double C1 = x[1][1] + x[2][1];
 
690
    
 
691
    // Get coordinates and map to the reference (FIAT) element
 
692
    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
 
693
    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
 
694
    
 
695
    // Compute number of derivatives.
 
696
    unsigned int num_derivatives = 1;
 
697
    for (unsigned int r = 0; r < n; r++)
 
698
    {
 
699
      num_derivatives *= 2;
 
700
    }// end loop over 'r'
 
701
    
 
702
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
703
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
704
    for (unsigned int row = 0; row < num_derivatives; row++)
 
705
    {
 
706
      combinations[row] = new unsigned int [n];
 
707
      for (unsigned int col = 0; col < n; col++)
 
708
        combinations[row][col] = 0;
 
709
    }
 
710
    
 
711
    // Generate combinations of derivatives
 
712
    for (unsigned int row = 1; row < num_derivatives; row++)
 
713
    {
 
714
      for (unsigned int num = 0; num < row; num++)
 
715
      {
 
716
        for (unsigned int col = n-1; col+1 > 0; col--)
 
717
        {
 
718
          if (combinations[row][col] + 1 > 1)
 
719
            combinations[row][col] = 0;
 
720
          else
 
721
          {
 
722
            combinations[row][col] += 1;
 
723
            break;
 
724
          }
 
725
        }
 
726
      }
 
727
    }
 
728
    
 
729
    // Compute inverse of Jacobian
 
730
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
 
731
    
 
732
    // Declare transformation matrix
 
733
    // Declare pointer to two dimensional array and initialise
 
734
    double **transform = new double *[num_derivatives];
 
735
    
 
736
    for (unsigned int j = 0; j < num_derivatives; j++)
 
737
    {
 
738
      transform[j] = new double [num_derivatives];
 
739
      for (unsigned int k = 0; k < num_derivatives; k++)
 
740
        transform[j][k] = 1;
 
741
    }
 
742
    
 
743
    // Construct transformation matrix
 
744
    for (unsigned int row = 0; row < num_derivatives; row++)
 
745
    {
 
746
      for (unsigned int col = 0; col < num_derivatives; col++)
 
747
      {
 
748
        for (unsigned int k = 0; k < n; k++)
 
749
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
750
      }
 
751
    }
 
752
    
 
753
    // Reset values. Assuming that values is always an array.
 
754
    for (unsigned int r = 0; r < num_derivatives; r++)
 
755
    {
 
756
      values[r] = 0.0;
 
757
    }// end loop over 'r'
 
758
    
 
759
    switch (i)
 
760
    {
 
761
    case 0:
 
762
      {
 
763
        
 
764
      // Array of basisvalues.
 
765
      double basisvalues[3] = {0.0, 0.0, 0.0};
 
766
      
 
767
      // Declare helper variables.
 
768
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
769
      
 
770
      // Compute basisvalues.
 
771
      basisvalues[0] = 1.0;
 
772
      basisvalues[1] = tmp0;
 
773
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
774
      basisvalues[0] *= std::sqrt(0.5);
 
775
      basisvalues[2] *= std::sqrt(1.0);
 
776
      basisvalues[1] *= std::sqrt(3.0);
 
777
      
 
778
      // Table(s) of coefficients.
 
779
      static const double coefficients0[3] = \
 
780
      {0.47140452, -0.28867513, -0.16666667};
 
781
      
 
782
      // Tables of derivatives of the polynomial base (transpose).
 
783
      static const double dmats0[3][3] = \
 
784
      {{0.0, 0.0, 0.0},
 
785
      {4.8989795, 0.0, 0.0},
 
786
      {0.0, 0.0, 0.0}};
 
787
      
 
788
      static const double dmats1[3][3] = \
 
789
      {{0.0, 0.0, 0.0},
 
790
      {2.4494897, 0.0, 0.0},
 
791
      {4.2426407, 0.0, 0.0}};
 
792
      
 
793
      // Compute reference derivatives.
 
794
      // Declare pointer to array of derivatives on FIAT element.
 
795
      double *derivatives = new double[num_derivatives];
 
796
      for (unsigned int r = 0; r < num_derivatives; r++)
 
797
      {
 
798
        derivatives[r] = 0.0;
 
799
      }// end loop over 'r'
 
800
      
 
801
      // Declare derivative matrix (of polynomial basis).
 
802
      double dmats[3][3] = \
 
803
      {{1.0, 0.0, 0.0},
 
804
      {0.0, 1.0, 0.0},
 
805
      {0.0, 0.0, 1.0}};
 
806
      
 
807
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
808
      double dmats_old[3][3] = \
 
809
      {{1.0, 0.0, 0.0},
 
810
      {0.0, 1.0, 0.0},
 
811
      {0.0, 0.0, 1.0}};
 
812
      
 
813
      // Loop possible derivatives.
 
814
      for (unsigned int r = 0; r < num_derivatives; r++)
 
815
      {
 
816
        // Resetting dmats values to compute next derivative.
 
817
        for (unsigned int t = 0; t < 3; t++)
 
818
        {
 
819
          for (unsigned int u = 0; u < 3; u++)
 
820
          {
 
821
            dmats[t][u] = 0.0;
 
822
            if (t == u)
 
823
            {
 
824
            dmats[t][u] = 1.0;
 
825
            }
 
826
            
 
827
          }// end loop over 'u'
 
828
        }// end loop over 't'
 
829
        
 
830
        // Looping derivative order to generate dmats.
 
831
        for (unsigned int s = 0; s < n; s++)
 
832
        {
 
833
          // Updating dmats_old with new values and resetting dmats.
 
834
          for (unsigned int t = 0; t < 3; t++)
 
835
          {
 
836
            for (unsigned int u = 0; u < 3; u++)
 
837
            {
 
838
              dmats_old[t][u] = dmats[t][u];
 
839
              dmats[t][u] = 0.0;
 
840
            }// end loop over 'u'
 
841
          }// end loop over 't'
 
842
          
 
843
          // Update dmats using an inner product.
 
844
          if (combinations[r][s] == 0)
 
845
          {
 
846
          for (unsigned int t = 0; t < 3; t++)
 
847
          {
 
848
            for (unsigned int u = 0; u < 3; u++)
 
849
            {
 
850
              for (unsigned int tu = 0; tu < 3; tu++)
 
851
              {
 
852
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
853
              }// end loop over 'tu'
 
854
            }// end loop over 'u'
 
855
          }// end loop over 't'
 
856
          }
 
857
          
 
858
          if (combinations[r][s] == 1)
 
859
          {
 
860
          for (unsigned int t = 0; t < 3; t++)
 
861
          {
 
862
            for (unsigned int u = 0; u < 3; u++)
 
863
            {
 
864
              for (unsigned int tu = 0; tu < 3; tu++)
 
865
              {
 
866
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
867
              }// end loop over 'tu'
 
868
            }// end loop over 'u'
 
869
          }// end loop over 't'
 
870
          }
 
871
          
 
872
        }// end loop over 's'
 
873
        for (unsigned int s = 0; s < 3; s++)
 
874
        {
 
875
          for (unsigned int t = 0; t < 3; t++)
 
876
          {
 
877
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
878
          }// end loop over 't'
 
879
        }// end loop over 's'
 
880
      }// end loop over 'r'
 
881
      
 
882
      // Transform derivatives back to physical element
 
883
      for (unsigned int r = 0; r < num_derivatives; r++)
 
884
      {
 
885
        for (unsigned int s = 0; s < num_derivatives; s++)
 
886
        {
 
887
          values[r] += transform[r][s]*derivatives[s];
 
888
        }// end loop over 's'
 
889
      }// end loop over 'r'
 
890
      
 
891
      // Delete pointer to array of derivatives on FIAT element
 
892
      delete [] derivatives;
 
893
      
 
894
      // Delete pointer to array of combinations of derivatives and transform
 
895
      for (unsigned int r = 0; r < num_derivatives; r++)
 
896
      {
 
897
        delete [] combinations[r];
 
898
      }// end loop over 'r'
 
899
      delete [] combinations;
 
900
      for (unsigned int r = 0; r < num_derivatives; r++)
 
901
      {
 
902
        delete [] transform[r];
 
903
      }// end loop over 'r'
 
904
      delete [] transform;
 
905
        break;
 
906
      }
 
907
    case 1:
 
908
      {
 
909
        
 
910
      // Array of basisvalues.
 
911
      double basisvalues[3] = {0.0, 0.0, 0.0};
 
912
      
 
913
      // Declare helper variables.
 
914
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
915
      
 
916
      // Compute basisvalues.
 
917
      basisvalues[0] = 1.0;
 
918
      basisvalues[1] = tmp0;
 
919
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
920
      basisvalues[0] *= std::sqrt(0.5);
 
921
      basisvalues[2] *= std::sqrt(1.0);
 
922
      basisvalues[1] *= std::sqrt(3.0);
 
923
      
 
924
      // Table(s) of coefficients.
 
925
      static const double coefficients0[3] = \
 
926
      {0.47140452, 0.28867513, -0.16666667};
 
927
      
 
928
      // Tables of derivatives of the polynomial base (transpose).
 
929
      static const double dmats0[3][3] = \
 
930
      {{0.0, 0.0, 0.0},
 
931
      {4.8989795, 0.0, 0.0},
 
932
      {0.0, 0.0, 0.0}};
 
933
      
 
934
      static const double dmats1[3][3] = \
 
935
      {{0.0, 0.0, 0.0},
 
936
      {2.4494897, 0.0, 0.0},
 
937
      {4.2426407, 0.0, 0.0}};
 
938
      
 
939
      // Compute reference derivatives.
 
940
      // Declare pointer to array of derivatives on FIAT element.
 
941
      double *derivatives = new double[num_derivatives];
 
942
      for (unsigned int r = 0; r < num_derivatives; r++)
 
943
      {
 
944
        derivatives[r] = 0.0;
 
945
      }// end loop over 'r'
 
946
      
 
947
      // Declare derivative matrix (of polynomial basis).
 
948
      double dmats[3][3] = \
 
949
      {{1.0, 0.0, 0.0},
 
950
      {0.0, 1.0, 0.0},
 
951
      {0.0, 0.0, 1.0}};
 
952
      
 
953
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
954
      double dmats_old[3][3] = \
 
955
      {{1.0, 0.0, 0.0},
 
956
      {0.0, 1.0, 0.0},
 
957
      {0.0, 0.0, 1.0}};
 
958
      
 
959
      // Loop possible derivatives.
 
960
      for (unsigned int r = 0; r < num_derivatives; r++)
 
961
      {
 
962
        // Resetting dmats values to compute next derivative.
 
963
        for (unsigned int t = 0; t < 3; t++)
 
964
        {
 
965
          for (unsigned int u = 0; u < 3; u++)
 
966
          {
 
967
            dmats[t][u] = 0.0;
 
968
            if (t == u)
 
969
            {
 
970
            dmats[t][u] = 1.0;
 
971
            }
 
972
            
 
973
          }// end loop over 'u'
 
974
        }// end loop over 't'
 
975
        
 
976
        // Looping derivative order to generate dmats.
 
977
        for (unsigned int s = 0; s < n; s++)
 
978
        {
 
979
          // Updating dmats_old with new values and resetting dmats.
 
980
          for (unsigned int t = 0; t < 3; t++)
 
981
          {
 
982
            for (unsigned int u = 0; u < 3; u++)
 
983
            {
 
984
              dmats_old[t][u] = dmats[t][u];
 
985
              dmats[t][u] = 0.0;
 
986
            }// end loop over 'u'
 
987
          }// end loop over 't'
 
988
          
 
989
          // Update dmats using an inner product.
 
990
          if (combinations[r][s] == 0)
 
991
          {
 
992
          for (unsigned int t = 0; t < 3; t++)
 
993
          {
 
994
            for (unsigned int u = 0; u < 3; u++)
 
995
            {
 
996
              for (unsigned int tu = 0; tu < 3; tu++)
 
997
              {
 
998
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
999
              }// end loop over 'tu'
 
1000
            }// end loop over 'u'
 
1001
          }// end loop over 't'
 
1002
          }
 
1003
          
 
1004
          if (combinations[r][s] == 1)
 
1005
          {
 
1006
          for (unsigned int t = 0; t < 3; t++)
 
1007
          {
 
1008
            for (unsigned int u = 0; u < 3; u++)
 
1009
            {
 
1010
              for (unsigned int tu = 0; tu < 3; tu++)
 
1011
              {
 
1012
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
1013
              }// end loop over 'tu'
 
1014
            }// end loop over 'u'
 
1015
          }// end loop over 't'
 
1016
          }
 
1017
          
 
1018
        }// end loop over 's'
 
1019
        for (unsigned int s = 0; s < 3; s++)
 
1020
        {
 
1021
          for (unsigned int t = 0; t < 3; t++)
 
1022
          {
 
1023
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
1024
          }// end loop over 't'
 
1025
        }// end loop over 's'
 
1026
      }// end loop over 'r'
 
1027
      
 
1028
      // Transform derivatives back to physical element
 
1029
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1030
      {
 
1031
        for (unsigned int s = 0; s < num_derivatives; s++)
 
1032
        {
 
1033
          values[r] += transform[r][s]*derivatives[s];
 
1034
        }// end loop over 's'
 
1035
      }// end loop over 'r'
 
1036
      
 
1037
      // Delete pointer to array of derivatives on FIAT element
 
1038
      delete [] derivatives;
 
1039
      
 
1040
      // Delete pointer to array of combinations of derivatives and transform
 
1041
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1042
      {
 
1043
        delete [] combinations[r];
 
1044
      }// end loop over 'r'
 
1045
      delete [] combinations;
 
1046
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1047
      {
 
1048
        delete [] transform[r];
 
1049
      }// end loop over 'r'
 
1050
      delete [] transform;
 
1051
        break;
 
1052
      }
 
1053
    case 2:
 
1054
      {
 
1055
        
 
1056
      // Array of basisvalues.
 
1057
      double basisvalues[3] = {0.0, 0.0, 0.0};
 
1058
      
 
1059
      // Declare helper variables.
 
1060
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
1061
      
 
1062
      // Compute basisvalues.
 
1063
      basisvalues[0] = 1.0;
 
1064
      basisvalues[1] = tmp0;
 
1065
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
1066
      basisvalues[0] *= std::sqrt(0.5);
 
1067
      basisvalues[2] *= std::sqrt(1.0);
 
1068
      basisvalues[1] *= std::sqrt(3.0);
 
1069
      
 
1070
      // Table(s) of coefficients.
 
1071
      static const double coefficients0[3] = \
 
1072
      {0.47140452, 0.0, 0.33333333};
 
1073
      
 
1074
      // Tables of derivatives of the polynomial base (transpose).
 
1075
      static const double dmats0[3][3] = \
 
1076
      {{0.0, 0.0, 0.0},
 
1077
      {4.8989795, 0.0, 0.0},
 
1078
      {0.0, 0.0, 0.0}};
 
1079
      
 
1080
      static const double dmats1[3][3] = \
 
1081
      {{0.0, 0.0, 0.0},
 
1082
      {2.4494897, 0.0, 0.0},
 
1083
      {4.2426407, 0.0, 0.0}};
 
1084
      
 
1085
      // Compute reference derivatives.
 
1086
      // Declare pointer to array of derivatives on FIAT element.
 
1087
      double *derivatives = new double[num_derivatives];
 
1088
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1089
      {
 
1090
        derivatives[r] = 0.0;
 
1091
      }// end loop over 'r'
 
1092
      
 
1093
      // Declare derivative matrix (of polynomial basis).
 
1094
      double dmats[3][3] = \
 
1095
      {{1.0, 0.0, 0.0},
 
1096
      {0.0, 1.0, 0.0},
 
1097
      {0.0, 0.0, 1.0}};
 
1098
      
 
1099
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
1100
      double dmats_old[3][3] = \
 
1101
      {{1.0, 0.0, 0.0},
 
1102
      {0.0, 1.0, 0.0},
 
1103
      {0.0, 0.0, 1.0}};
 
1104
      
 
1105
      // Loop possible derivatives.
 
1106
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1107
      {
 
1108
        // Resetting dmats values to compute next derivative.
 
1109
        for (unsigned int t = 0; t < 3; t++)
 
1110
        {
 
1111
          for (unsigned int u = 0; u < 3; u++)
 
1112
          {
 
1113
            dmats[t][u] = 0.0;
 
1114
            if (t == u)
 
1115
            {
 
1116
            dmats[t][u] = 1.0;
 
1117
            }
 
1118
            
 
1119
          }// end loop over 'u'
 
1120
        }// end loop over 't'
 
1121
        
 
1122
        // Looping derivative order to generate dmats.
 
1123
        for (unsigned int s = 0; s < n; s++)
 
1124
        {
 
1125
          // Updating dmats_old with new values and resetting dmats.
 
1126
          for (unsigned int t = 0; t < 3; t++)
 
1127
          {
 
1128
            for (unsigned int u = 0; u < 3; u++)
 
1129
            {
 
1130
              dmats_old[t][u] = dmats[t][u];
 
1131
              dmats[t][u] = 0.0;
 
1132
            }// end loop over 'u'
 
1133
          }// end loop over 't'
 
1134
          
 
1135
          // Update dmats using an inner product.
 
1136
          if (combinations[r][s] == 0)
 
1137
          {
 
1138
          for (unsigned int t = 0; t < 3; t++)
 
1139
          {
 
1140
            for (unsigned int u = 0; u < 3; u++)
 
1141
            {
 
1142
              for (unsigned int tu = 0; tu < 3; tu++)
 
1143
              {
 
1144
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
1145
              }// end loop over 'tu'
 
1146
            }// end loop over 'u'
 
1147
          }// end loop over 't'
 
1148
          }
 
1149
          
 
1150
          if (combinations[r][s] == 1)
 
1151
          {
 
1152
          for (unsigned int t = 0; t < 3; t++)
 
1153
          {
 
1154
            for (unsigned int u = 0; u < 3; u++)
 
1155
            {
 
1156
              for (unsigned int tu = 0; tu < 3; tu++)
 
1157
              {
 
1158
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
1159
              }// end loop over 'tu'
 
1160
            }// end loop over 'u'
 
1161
          }// end loop over 't'
 
1162
          }
 
1163
          
 
1164
        }// end loop over 's'
 
1165
        for (unsigned int s = 0; s < 3; s++)
 
1166
        {
 
1167
          for (unsigned int t = 0; t < 3; t++)
 
1168
          {
 
1169
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
1170
          }// end loop over 't'
 
1171
        }// end loop over 's'
 
1172
      }// end loop over 'r'
 
1173
      
 
1174
      // Transform derivatives back to physical element
 
1175
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1176
      {
 
1177
        for (unsigned int s = 0; s < num_derivatives; s++)
 
1178
        {
 
1179
          values[r] += transform[r][s]*derivatives[s];
 
1180
        }// end loop over 's'
 
1181
      }// end loop over 'r'
 
1182
      
 
1183
      // Delete pointer to array of derivatives on FIAT element
 
1184
      delete [] derivatives;
 
1185
      
 
1186
      // Delete pointer to array of combinations of derivatives and transform
 
1187
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1188
      {
 
1189
        delete [] combinations[r];
 
1190
      }// end loop over 'r'
 
1191
      delete [] combinations;
 
1192
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1193
      {
 
1194
        delete [] transform[r];
 
1195
      }// end loop over 'r'
 
1196
      delete [] transform;
 
1197
        break;
 
1198
      }
 
1199
    }
 
1200
    
 
1201
  }
 
1202
 
 
1203
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
1204
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
1205
                                              double* values,
 
1206
                                              const double* coordinates,
 
1207
                                              const ufc::cell& c) const
 
1208
  {
 
1209
    // Compute number of derivatives.
 
1210
    unsigned int num_derivatives = 1;
 
1211
    for (unsigned int r = 0; r < n; r++)
 
1212
    {
 
1213
      num_derivatives *= 2;
 
1214
    }// end loop over 'r'
 
1215
    
 
1216
    // Helper variable to hold values of a single dof.
 
1217
    double *dof_values = new double[num_derivatives];
 
1218
    for (unsigned int r = 0; r < num_derivatives; r++)
 
1219
    {
 
1220
      dof_values[r] = 0.0;
 
1221
    }// end loop over 'r'
 
1222
    
 
1223
    // Loop dofs and call evaluate_basis_derivatives.
 
1224
    for (unsigned int r = 0; r < 3; r++)
 
1225
    {
 
1226
      evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
 
1227
      for (unsigned int s = 0; s < num_derivatives; s++)
 
1228
      {
 
1229
        values[r*num_derivatives + s] = dof_values[s];
 
1230
      }// end loop over 's'
 
1231
    }// end loop over 'r'
 
1232
    
 
1233
    // Delete pointer.
 
1234
    delete [] dof_values;
 
1235
  }
 
1236
 
 
1237
  /// Evaluate linear functional for dof i on the function f
 
1238
  virtual double evaluate_dof(unsigned int i,
 
1239
                              const ufc::function& f,
 
1240
                              const ufc::cell& c) const
 
1241
  {
 
1242
    // Declare variables for result of evaluation.
 
1243
    double vals[1];
 
1244
    
 
1245
    // Declare variable for physical coordinates.
 
1246
    double y[2];
 
1247
    const double * const * x = c.coordinates;
 
1248
    switch (i)
 
1249
    {
 
1250
    case 0:
 
1251
      {
 
1252
        y[0] = x[0][0];
 
1253
      y[1] = x[0][1];
 
1254
      f.evaluate(vals, y, c);
 
1255
      return vals[0];
 
1256
        break;
 
1257
      }
 
1258
    case 1:
 
1259
      {
 
1260
        y[0] = x[1][0];
 
1261
      y[1] = x[1][1];
 
1262
      f.evaluate(vals, y, c);
 
1263
      return vals[0];
 
1264
        break;
 
1265
      }
 
1266
    case 2:
 
1267
      {
 
1268
        y[0] = x[2][0];
 
1269
      y[1] = x[2][1];
 
1270
      f.evaluate(vals, y, c);
 
1271
      return vals[0];
 
1272
        break;
 
1273
      }
 
1274
    }
 
1275
    
 
1276
    return 0.0;
 
1277
  }
 
1278
 
 
1279
  /// Evaluate linear functionals for all dofs on the function f
 
1280
  virtual void evaluate_dofs(double* values,
 
1281
                             const ufc::function& f,
 
1282
                             const ufc::cell& c) const
 
1283
  {
 
1284
    // Declare variables for result of evaluation.
 
1285
    double vals[1];
 
1286
    
 
1287
    // Declare variable for physical coordinates.
 
1288
    double y[2];
 
1289
    const double * const * x = c.coordinates;
 
1290
    y[0] = x[0][0];
 
1291
    y[1] = x[0][1];
 
1292
    f.evaluate(vals, y, c);
 
1293
    values[0] = vals[0];
 
1294
    y[0] = x[1][0];
 
1295
    y[1] = x[1][1];
 
1296
    f.evaluate(vals, y, c);
 
1297
    values[1] = vals[0];
 
1298
    y[0] = x[2][0];
 
1299
    y[1] = x[2][1];
 
1300
    f.evaluate(vals, y, c);
 
1301
    values[2] = vals[0];
 
1302
  }
 
1303
 
 
1304
  /// Interpolate vertex values from dof values
 
1305
  virtual void interpolate_vertex_values(double* vertex_values,
 
1306
                                         const double* dof_values,
 
1307
                                         const ufc::cell& c) const
 
1308
  {
 
1309
    // Evaluate function and change variables
 
1310
    vertex_values[0] = dof_values[0];
 
1311
    vertex_values[1] = dof_values[1];
 
1312
    vertex_values[2] = dof_values[2];
 
1313
  }
 
1314
 
 
1315
  /// Map coordinate xhat from reference cell to coordinate x in cell
 
1316
  virtual void map_from_reference_cell(double* x,
 
1317
                                       const double* xhat,
 
1318
                                       const ufc::cell& c) const
 
1319
  {
 
1320
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
1321
  }
 
1322
 
 
1323
  /// Map from coordinate x in cell to coordinate xhat in reference cell
 
1324
  virtual void map_to_reference_cell(double* xhat,
 
1325
                                     const double* x,
 
1326
                                     const ufc::cell& c) const
 
1327
  {
 
1328
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
1329
  }
 
1330
 
 
1331
  /// Return the number of sub elements (for a mixed element)
 
1332
  virtual unsigned int num_sub_elements() const
 
1333
  {
 
1334
    return 0;
 
1335
  }
 
1336
 
 
1337
  /// Create a new finite element for sub element i (for a mixed element)
 
1338
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
1339
  {
 
1340
    return 0;
 
1341
  }
 
1342
 
 
1343
  /// Create a new class instance
 
1344
  virtual ufc::finite_element* create() const
 
1345
  {
 
1346
    return new heat_finite_element_1();
 
1347
  }
 
1348
 
 
1349
};
 
1350
 
 
1351
/// This class defines the interface for a local-to-global mapping of
 
1352
/// degrees of freedom (dofs).
 
1353
 
 
1354
class heat_dofmap_0: public ufc::dofmap
 
1355
{
 
1356
private:
 
1357
 
 
1358
  unsigned int _global_dimension;
 
1359
public:
 
1360
 
 
1361
  /// Constructor
 
1362
  heat_dofmap_0() : ufc::dofmap()
 
1363
  {
 
1364
    _global_dimension = 0;
 
1365
  }
 
1366
 
 
1367
  /// Destructor
 
1368
  virtual ~heat_dofmap_0()
 
1369
  {
 
1370
    // Do nothing
 
1371
  }
 
1372
 
 
1373
  /// Return a string identifying the dofmap
 
1374
  virtual const char* signature() const
 
1375
  {
 
1376
    return "FFC dofmap for FiniteElement('Real', Cell('triangle', Space(2)), 0, None)";
 
1377
  }
 
1378
 
 
1379
  /// Return true iff mesh entities of topological dimension d are needed
 
1380
  virtual bool needs_mesh_entities(unsigned int d) const
 
1381
  {
 
1382
    switch (d)
 
1383
    {
 
1384
    case 0:
 
1385
      {
 
1386
        return false;
 
1387
        break;
 
1388
      }
 
1389
    case 1:
 
1390
      {
 
1391
        return false;
 
1392
        break;
 
1393
      }
 
1394
    case 2:
 
1395
      {
 
1396
        return false;
 
1397
        break;
 
1398
      }
 
1399
    }
 
1400
    
 
1401
    return false;
 
1402
  }
 
1403
 
 
1404
  /// Initialize dofmap for mesh (return true iff init_cell() is needed)
 
1405
  virtual bool init_mesh(const ufc::mesh& m)
 
1406
  {
 
1407
    _global_dimension = 1;
 
1408
    return false;
 
1409
  }
 
1410
 
 
1411
  /// Initialize dofmap for given cell
 
1412
  virtual void init_cell(const ufc::mesh& m,
 
1413
                         const ufc::cell& c)
 
1414
  {
 
1415
    // Do nothing
 
1416
  }
 
1417
 
 
1418
  /// Finish initialization of dofmap for cells
 
1419
  virtual void init_cell_finalize()
 
1420
  {
 
1421
    // Do nothing
 
1422
  }
 
1423
 
 
1424
  /// Return the topological dimension of the associated cell shape
 
1425
  virtual unsigned int topological_dimension() const
 
1426
  {
 
1427
    return 2;
 
1428
  }
 
1429
 
 
1430
  /// Return the geometric dimension of the associated cell shape
 
1431
  virtual unsigned int geometric_dimension() const
 
1432
  {
 
1433
    return 2;
 
1434
  }
 
1435
 
 
1436
  /// Return the dimension of the global finite element function space
 
1437
  virtual unsigned int global_dimension() const
 
1438
  {
 
1439
    return _global_dimension;
 
1440
  }
 
1441
 
 
1442
  /// Return the dimension of the local finite element function space for a cell
 
1443
  virtual unsigned int local_dimension(const ufc::cell& c) const
 
1444
  {
 
1445
    return 1;
 
1446
  }
 
1447
 
 
1448
  /// Return the maximum dimension of the local finite element function space
 
1449
  virtual unsigned int max_local_dimension() const
 
1450
  {
 
1451
    return 1;
 
1452
  }
 
1453
 
 
1454
  /// Return the number of dofs on each cell facet
 
1455
  virtual unsigned int num_facet_dofs() const
 
1456
  {
 
1457
    return 0;
 
1458
  }
 
1459
 
 
1460
  /// Return the number of dofs associated with each cell entity of dimension d
 
1461
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1462
  {
 
1463
    switch (d)
 
1464
    {
 
1465
    case 0:
 
1466
      {
 
1467
        return 0;
 
1468
        break;
 
1469
      }
 
1470
    case 1:
 
1471
      {
 
1472
        return 0;
 
1473
        break;
 
1474
      }
 
1475
    case 2:
 
1476
      {
 
1477
        return 1;
 
1478
        break;
 
1479
      }
 
1480
    }
 
1481
    
 
1482
    return 0;
 
1483
  }
 
1484
 
 
1485
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1486
  virtual void tabulate_dofs(unsigned int* dofs,
 
1487
                             const ufc::mesh& m,
 
1488
                             const ufc::cell& c) const
 
1489
  {
 
1490
    dofs[0] = 0;
 
1491
  }
 
1492
 
 
1493
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1494
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1495
                                   unsigned int facet) const
 
1496
  {
 
1497
    switch (facet)
 
1498
    {
 
1499
    case 0:
 
1500
      {
 
1501
        
 
1502
        break;
 
1503
      }
 
1504
    case 1:
 
1505
      {
 
1506
        
 
1507
        break;
 
1508
      }
 
1509
    case 2:
 
1510
      {
 
1511
        
 
1512
        break;
 
1513
      }
 
1514
    }
 
1515
    
 
1516
  }
 
1517
 
 
1518
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1519
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1520
                                    unsigned int d, unsigned int i) const
 
1521
  {
 
1522
    if (d > 2)
 
1523
    {
 
1524
    std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
 
1525
    }
 
1526
    
 
1527
    switch (d)
 
1528
    {
 
1529
    case 0:
 
1530
      {
 
1531
        
 
1532
        break;
 
1533
      }
 
1534
    case 1:
 
1535
      {
 
1536
        
 
1537
        break;
 
1538
      }
 
1539
    case 2:
 
1540
      {
 
1541
        if (i > 0)
 
1542
      {
 
1543
      std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
 
1544
      }
 
1545
      
 
1546
      dofs[0] = 0;
 
1547
        break;
 
1548
      }
 
1549
    }
 
1550
    
 
1551
  }
 
1552
 
 
1553
  /// Tabulate the coordinates of all dofs on a cell
 
1554
  virtual void tabulate_coordinates(double** coordinates,
 
1555
                                    const ufc::cell& c) const
 
1556
  {
 
1557
    const double * const * x = c.coordinates;
 
1558
    
 
1559
    coordinates[0][0] = 0.33333333*x[0][0] + 0.33333333*x[1][0] + 0.33333333*x[2][0];
 
1560
    coordinates[0][1] = 0.33333333*x[0][1] + 0.33333333*x[1][1] + 0.33333333*x[2][1];
 
1561
  }
 
1562
 
 
1563
  /// Return the number of sub dofmaps (for a mixed element)
 
1564
  virtual unsigned int num_sub_dofmaps() const
 
1565
  {
 
1566
    return 0;
 
1567
  }
 
1568
 
 
1569
  /// Create a new dofmap for sub dofmap i (for a mixed element)
 
1570
  virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
 
1571
  {
 
1572
    return 0;
 
1573
  }
 
1574
 
 
1575
  /// Create a new class instance
 
1576
  virtual ufc::dofmap* create() const
 
1577
  {
 
1578
    return new heat_dofmap_0();
 
1579
  }
 
1580
 
 
1581
};
 
1582
 
 
1583
/// This class defines the interface for a local-to-global mapping of
 
1584
/// degrees of freedom (dofs).
 
1585
 
 
1586
class heat_dofmap_1: public ufc::dofmap
 
1587
{
 
1588
private:
 
1589
 
 
1590
  unsigned int _global_dimension;
 
1591
public:
 
1592
 
 
1593
  /// Constructor
 
1594
  heat_dofmap_1() : ufc::dofmap()
 
1595
  {
 
1596
    _global_dimension = 0;
 
1597
  }
 
1598
 
 
1599
  /// Destructor
 
1600
  virtual ~heat_dofmap_1()
 
1601
  {
 
1602
    // Do nothing
 
1603
  }
 
1604
 
 
1605
  /// Return a string identifying the dofmap
 
1606
  virtual const char* signature() const
 
1607
  {
 
1608
    return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None)";
 
1609
  }
 
1610
 
 
1611
  /// Return true iff mesh entities of topological dimension d are needed
 
1612
  virtual bool needs_mesh_entities(unsigned int d) const
 
1613
  {
 
1614
    switch (d)
 
1615
    {
 
1616
    case 0:
 
1617
      {
 
1618
        return true;
 
1619
        break;
 
1620
      }
 
1621
    case 1:
 
1622
      {
 
1623
        return false;
 
1624
        break;
 
1625
      }
 
1626
    case 2:
 
1627
      {
 
1628
        return false;
 
1629
        break;
 
1630
      }
 
1631
    }
 
1632
    
 
1633
    return false;
 
1634
  }
 
1635
 
 
1636
  /// Initialize dofmap for mesh (return true iff init_cell() is needed)
 
1637
  virtual bool init_mesh(const ufc::mesh& m)
 
1638
  {
 
1639
    _global_dimension = m.num_entities[0];
 
1640
    return false;
 
1641
  }
 
1642
 
 
1643
  /// Initialize dofmap for given cell
 
1644
  virtual void init_cell(const ufc::mesh& m,
 
1645
                         const ufc::cell& c)
 
1646
  {
 
1647
    // Do nothing
 
1648
  }
 
1649
 
 
1650
  /// Finish initialization of dofmap for cells
 
1651
  virtual void init_cell_finalize()
 
1652
  {
 
1653
    // Do nothing
 
1654
  }
 
1655
 
 
1656
  /// Return the topological dimension of the associated cell shape
 
1657
  virtual unsigned int topological_dimension() const
 
1658
  {
 
1659
    return 2;
 
1660
  }
 
1661
 
 
1662
  /// Return the geometric dimension of the associated cell shape
 
1663
  virtual unsigned int geometric_dimension() const
 
1664
  {
 
1665
    return 2;
 
1666
  }
 
1667
 
 
1668
  /// Return the dimension of the global finite element function space
 
1669
  virtual unsigned int global_dimension() const
 
1670
  {
 
1671
    return _global_dimension;
 
1672
  }
 
1673
 
 
1674
  /// Return the dimension of the local finite element function space for a cell
 
1675
  virtual unsigned int local_dimension(const ufc::cell& c) const
 
1676
  {
 
1677
    return 3;
 
1678
  }
 
1679
 
 
1680
  /// Return the maximum dimension of the local finite element function space
 
1681
  virtual unsigned int max_local_dimension() const
 
1682
  {
 
1683
    return 3;
 
1684
  }
 
1685
 
 
1686
  /// Return the number of dofs on each cell facet
 
1687
  virtual unsigned int num_facet_dofs() const
 
1688
  {
 
1689
    return 2;
 
1690
  }
 
1691
 
 
1692
  /// Return the number of dofs associated with each cell entity of dimension d
 
1693
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1694
  {
 
1695
    switch (d)
 
1696
    {
 
1697
    case 0:
 
1698
      {
 
1699
        return 1;
 
1700
        break;
 
1701
      }
 
1702
    case 1:
 
1703
      {
 
1704
        return 0;
 
1705
        break;
 
1706
      }
 
1707
    case 2:
 
1708
      {
 
1709
        return 0;
 
1710
        break;
 
1711
      }
 
1712
    }
 
1713
    
 
1714
    return 0;
 
1715
  }
 
1716
 
 
1717
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1718
  virtual void tabulate_dofs(unsigned int* dofs,
 
1719
                             const ufc::mesh& m,
 
1720
                             const ufc::cell& c) const
 
1721
  {
 
1722
    dofs[0] = c.entity_indices[0][0];
 
1723
    dofs[1] = c.entity_indices[0][1];
 
1724
    dofs[2] = c.entity_indices[0][2];
 
1725
  }
 
1726
 
 
1727
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1728
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1729
                                   unsigned int facet) const
 
1730
  {
 
1731
    switch (facet)
 
1732
    {
 
1733
    case 0:
 
1734
      {
 
1735
        dofs[0] = 1;
 
1736
      dofs[1] = 2;
 
1737
        break;
 
1738
      }
 
1739
    case 1:
 
1740
      {
 
1741
        dofs[0] = 0;
 
1742
      dofs[1] = 2;
 
1743
        break;
 
1744
      }
 
1745
    case 2:
 
1746
      {
 
1747
        dofs[0] = 0;
 
1748
      dofs[1] = 1;
 
1749
        break;
 
1750
      }
 
1751
    }
 
1752
    
 
1753
  }
 
1754
 
 
1755
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1756
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1757
                                    unsigned int d, unsigned int i) const
 
1758
  {
 
1759
    if (d > 2)
 
1760
    {
 
1761
    std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
 
1762
    }
 
1763
    
 
1764
    switch (d)
 
1765
    {
 
1766
    case 0:
 
1767
      {
 
1768
        if (i > 2)
 
1769
      {
 
1770
      std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
 
1771
      }
 
1772
      
 
1773
      switch (i)
 
1774
      {
 
1775
      case 0:
 
1776
        {
 
1777
          dofs[0] = 0;
 
1778
          break;
 
1779
        }
 
1780
      case 1:
 
1781
        {
 
1782
          dofs[0] = 1;
 
1783
          break;
 
1784
        }
 
1785
      case 2:
 
1786
        {
 
1787
          dofs[0] = 2;
 
1788
          break;
 
1789
        }
 
1790
      }
 
1791
      
 
1792
        break;
 
1793
      }
 
1794
    case 1:
 
1795
      {
 
1796
        
 
1797
        break;
 
1798
      }
 
1799
    case 2:
 
1800
      {
 
1801
        
 
1802
        break;
 
1803
      }
 
1804
    }
 
1805
    
 
1806
  }
 
1807
 
 
1808
  /// Tabulate the coordinates of all dofs on a cell
 
1809
  virtual void tabulate_coordinates(double** coordinates,
 
1810
                                    const ufc::cell& c) const
 
1811
  {
 
1812
    const double * const * x = c.coordinates;
 
1813
    
 
1814
    coordinates[0][0] = x[0][0];
 
1815
    coordinates[0][1] = x[0][1];
 
1816
    coordinates[1][0] = x[1][0];
 
1817
    coordinates[1][1] = x[1][1];
 
1818
    coordinates[2][0] = x[2][0];
 
1819
    coordinates[2][1] = x[2][1];
 
1820
  }
 
1821
 
 
1822
  /// Return the number of sub dofmaps (for a mixed element)
 
1823
  virtual unsigned int num_sub_dofmaps() const
 
1824
  {
 
1825
    return 0;
 
1826
  }
 
1827
 
 
1828
  /// Create a new dofmap for sub dofmap i (for a mixed element)
 
1829
  virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
 
1830
  {
 
1831
    return 0;
 
1832
  }
 
1833
 
 
1834
  /// Create a new class instance
 
1835
  virtual ufc::dofmap* create() const
 
1836
  {
 
1837
    return new heat_dofmap_1();
 
1838
  }
 
1839
 
 
1840
};
 
1841
 
 
1842
/// This class defines the interface for the tabulation of the cell
 
1843
/// tensor corresponding to the local contribution to a form from
 
1844
/// the integral over a cell.
 
1845
 
 
1846
class heat_cell_integral_0_0: public ufc::cell_integral
 
1847
{
 
1848
public:
 
1849
 
 
1850
  /// Constructor
 
1851
  heat_cell_integral_0_0() : ufc::cell_integral()
 
1852
  {
 
1853
    // Do nothing
 
1854
  }
 
1855
 
 
1856
  /// Destructor
 
1857
  virtual ~heat_cell_integral_0_0()
 
1858
  {
 
1859
    // Do nothing
 
1860
  }
 
1861
 
 
1862
  /// Tabulate the tensor for the contribution from a local cell
 
1863
  virtual void tabulate_tensor(double* A,
 
1864
                               const double * const * w,
 
1865
                               const ufc::cell& c) const
 
1866
  {
 
1867
    // Extract vertex coordinates
 
1868
    const double * const * x = c.coordinates;
 
1869
    
 
1870
    // Compute Jacobian of affine map from reference cell
 
1871
    const double J_00 = x[1][0] - x[0][0];
 
1872
    const double J_01 = x[2][0] - x[0][0];
 
1873
    const double J_10 = x[1][1] - x[0][1];
 
1874
    const double J_11 = x[2][1] - x[0][1];
 
1875
    
 
1876
    // Compute determinant of Jacobian
 
1877
    double detJ = J_00*J_11 - J_01*J_10;
 
1878
    
 
1879
    // Compute inverse of Jacobian
 
1880
    const double K_00 =  J_11 / detJ;
 
1881
    const double K_01 = -J_01 / detJ;
 
1882
    const double K_10 = -J_10 / detJ;
 
1883
    const double K_11 =  J_00 / detJ;
 
1884
    
 
1885
    // Set scale factor
 
1886
    const double det = std::abs(detJ);
 
1887
    
 
1888
    // Cell Volume.
 
1889
    
 
1890
    // Compute circumradius, assuming triangle is embedded in 2D.
 
1891
    
 
1892
    
 
1893
    // Facet Area.
 
1894
    
 
1895
    // Array of quadrature weights.
 
1896
    static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
 
1897
    // Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
 
1898
    
 
1899
    // Value of basis functions at quadrature points.
 
1900
    static const double FE0[3][3] = \
 
1901
    {{0.66666667, 0.16666667, 0.16666667},
 
1902
    {0.16666667, 0.16666667, 0.66666667},
 
1903
    {0.16666667, 0.66666667, 0.16666667}};
 
1904
    
 
1905
    static const double FE0_D01[3][2] = \
 
1906
    {{-1.0, 1.0},
 
1907
    {-1.0, 1.0},
 
1908
    {-1.0, 1.0}};
 
1909
    
 
1910
    // Array of non-zero columns
 
1911
    static const unsigned int nzc1[2] = {0, 1};
 
1912
    
 
1913
    // Array of non-zero columns
 
1914
    static const unsigned int nzc0[2] = {0, 2};
 
1915
    
 
1916
    // Reset values in the element tensor.
 
1917
    for (unsigned int r = 0; r < 9; r++)
 
1918
    {
 
1919
      A[r] = 0.0;
 
1920
    }// end loop over 'r'
 
1921
    // Number of operations to compute geometry constants: 15.
 
1922
    double G[3];
 
1923
    G[0] = det*w[1][0]*(K_10*K_10 + K_11*K_11);
 
1924
    G[1] = det*w[1][0]*(K_00*K_10 + K_01*K_11);
 
1925
    G[2] = det*w[1][0]*(K_00*K_00 + K_01*K_01);
 
1926
    
 
1927
    // Compute element tensor using UFL quadrature representation
 
1928
    // Optimisations: ('eliminate zeros', True), ('ignore ones', True), ('ignore zero tables', True), ('optimisation', 'simplify_expressions'), ('remove zero terms', True)
 
1929
    
 
1930
    // Loop quadrature points for integral.
 
1931
    // Number of operations to compute element tensor for following IP loop = 264
 
1932
    for (unsigned int ip = 0; ip < 3; ip++)
 
1933
    {
 
1934
      
 
1935
      // Coefficient declarations.
 
1936
      double F0 = 0.0;
 
1937
      
 
1938
      // Total number of operations to compute function values = 6
 
1939
      for (unsigned int r = 0; r < 3; r++)
 
1940
      {
 
1941
        F0 += FE0[ip][r]*w[0][r];
 
1942
      }// end loop over 'r'
 
1943
      
 
1944
      // Number of operations to compute ip constants: 7
 
1945
      double I[4];
 
1946
      // Number of operations: 2
 
1947
      I[0] = F0*G[0]*W3[ip];
 
1948
      
 
1949
      // Number of operations: 2
 
1950
      I[1] = F0*G[1]*W3[ip];
 
1951
      
 
1952
      // Number of operations: 2
 
1953
      I[2] = F0*G[2]*W3[ip];
 
1954
      
 
1955
      // Number of operations: 1
 
1956
      I[3] = W3[ip]*det;
 
1957
      
 
1958
      
 
1959
      // Number of operations for primary indices: 27
 
1960
      for (unsigned int j = 0; j < 3; j++)
 
1961
      {
 
1962
        for (unsigned int k = 0; k < 3; k++)
 
1963
        {
 
1964
          // Number of operations to compute entry: 3
 
1965
          A[j*3 + k] += FE0[ip][j]*FE0[ip][k]*I[3];
 
1966
        }// end loop over 'k'
 
1967
      }// end loop over 'j'
 
1968
      
 
1969
      // Number of operations for primary indices: 48
 
1970
      for (unsigned int j = 0; j < 2; j++)
 
1971
      {
 
1972
        for (unsigned int k = 0; k < 2; k++)
 
1973
        {
 
1974
          // Number of operations to compute entry: 3
 
1975
          A[nzc0[j]*3 + nzc0[k]] += FE0_D01[ip][j]*FE0_D01[ip][k]*I[0];
 
1976
          // Number of operations to compute entry: 3
 
1977
          A[nzc0[j]*3 + nzc1[k]] += FE0_D01[ip][j]*FE0_D01[ip][k]*I[1];
 
1978
          // Number of operations to compute entry: 3
 
1979
          A[nzc1[j]*3 + nzc0[k]] += FE0_D01[ip][j]*FE0_D01[ip][k]*I[1];
 
1980
          // Number of operations to compute entry: 3
 
1981
          A[nzc1[j]*3 + nzc1[k]] += FE0_D01[ip][j]*FE0_D01[ip][k]*I[2];
 
1982
        }// end loop over 'k'
 
1983
      }// end loop over 'j'
 
1984
    }// end loop over 'ip'
 
1985
  }
 
1986
 
 
1987
  /// Tabulate the tensor for the contribution from a local cell
 
1988
  /// using the specified reference cell quadrature points/weights
 
1989
  virtual void tabulate_tensor(double* A,
 
1990
                               const double * const * w,
 
1991
                               const ufc::cell& c,
 
1992
                               unsigned int num_quadrature_points,
 
1993
                               const double * const * quadrature_points,
 
1994
                               const double* quadrature_weights) const
 
1995
  {
 
1996
    std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
 
1997
  }
 
1998
 
 
1999
};
 
2000
 
 
2001
/// This class defines the interface for the tabulation of the cell
 
2002
/// tensor corresponding to the local contribution to a form from
 
2003
/// the integral over a cell.
 
2004
 
 
2005
class heat_cell_integral_1_0: public ufc::cell_integral
 
2006
{
 
2007
public:
 
2008
 
 
2009
  /// Constructor
 
2010
  heat_cell_integral_1_0() : ufc::cell_integral()
 
2011
  {
 
2012
    // Do nothing
 
2013
  }
 
2014
 
 
2015
  /// Destructor
 
2016
  virtual ~heat_cell_integral_1_0()
 
2017
  {
 
2018
    // Do nothing
 
2019
  }
 
2020
 
 
2021
  /// Tabulate the tensor for the contribution from a local cell
 
2022
  virtual void tabulate_tensor(double* A,
 
2023
                               const double * const * w,
 
2024
                               const ufc::cell& c) const
 
2025
  {
 
2026
    // Extract vertex coordinates
 
2027
    const double * const * x = c.coordinates;
 
2028
    
 
2029
    // Compute Jacobian of affine map from reference cell
 
2030
    const double J_00 = x[1][0] - x[0][0];
 
2031
    const double J_01 = x[2][0] - x[0][0];
 
2032
    const double J_10 = x[1][1] - x[0][1];
 
2033
    const double J_11 = x[2][1] - x[0][1];
 
2034
    
 
2035
    // Compute determinant of Jacobian
 
2036
    double detJ = J_00*J_11 - J_01*J_10;
 
2037
    
 
2038
    // Compute inverse of Jacobian
 
2039
    
 
2040
    // Set scale factor
 
2041
    const double det = std::abs(detJ);
 
2042
    
 
2043
    // Cell Volume.
 
2044
    
 
2045
    // Compute circumradius, assuming triangle is embedded in 2D.
 
2046
    
 
2047
    
 
2048
    // Facet Area.
 
2049
    
 
2050
    // Array of quadrature weights.
 
2051
    static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
 
2052
    // Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
 
2053
    
 
2054
    // Value of basis functions at quadrature points.
 
2055
    static const double FE0[3][3] = \
 
2056
    {{0.66666667, 0.16666667, 0.16666667},
 
2057
    {0.16666667, 0.16666667, 0.66666667},
 
2058
    {0.16666667, 0.66666667, 0.16666667}};
 
2059
    
 
2060
    // Reset values in the element tensor.
 
2061
    for (unsigned int r = 0; r < 3; r++)
 
2062
    {
 
2063
      A[r] = 0.0;
 
2064
    }// end loop over 'r'
 
2065
    // Number of operations to compute geometry constants: 1.
 
2066
    double G[1];
 
2067
    G[0] = det*w[2][0];
 
2068
    
 
2069
    // Compute element tensor using UFL quadrature representation
 
2070
    // Optimisations: ('eliminate zeros', True), ('ignore ones', True), ('ignore zero tables', True), ('optimisation', 'simplify_expressions'), ('remove zero terms', True)
 
2071
    
 
2072
    // Loop quadrature points for integral.
 
2073
    // Number of operations to compute element tensor for following IP loop = 66
 
2074
    for (unsigned int ip = 0; ip < 3; ip++)
 
2075
    {
 
2076
      
 
2077
      // Coefficient declarations.
 
2078
      double F0 = 0.0;
 
2079
      double F1 = 0.0;
 
2080
      
 
2081
      // Total number of operations to compute function values = 12
 
2082
      for (unsigned int r = 0; r < 3; r++)
 
2083
      {
 
2084
        F0 += FE0[ip][r]*w[0][r];
 
2085
        F1 += FE0[ip][r]*w[1][r];
 
2086
      }// end loop over 'r'
 
2087
      
 
2088
      // Number of operations to compute ip constants: 4
 
2089
      double I[1];
 
2090
      // Number of operations: 4
 
2091
      I[0] = W3[ip]*(F0*det + F1*G[0]);
 
2092
      
 
2093
      
 
2094
      // Number of operations for primary indices: 6
 
2095
      for (unsigned int j = 0; j < 3; j++)
 
2096
      {
 
2097
        // Number of operations to compute entry: 2
 
2098
        A[j] += FE0[ip][j]*I[0];
 
2099
      }// end loop over 'j'
 
2100
    }// end loop over 'ip'
 
2101
  }
 
2102
 
 
2103
  /// Tabulate the tensor for the contribution from a local cell
 
2104
  /// using the specified reference cell quadrature points/weights
 
2105
  virtual void tabulate_tensor(double* A,
 
2106
                               const double * const * w,
 
2107
                               const ufc::cell& c,
 
2108
                               unsigned int num_quadrature_points,
 
2109
                               const double * const * quadrature_points,
 
2110
                               const double* quadrature_weights) const
 
2111
  {
 
2112
    std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
 
2113
  }
 
2114
 
 
2115
};
 
2116
 
 
2117
/// This class defines the interface for the assembly of the global
 
2118
/// tensor corresponding to a form with r + n arguments, that is, a
 
2119
/// mapping
 
2120
///
 
2121
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
2122
///
 
2123
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
2124
/// global tensor A is defined by
 
2125
///
 
2126
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
2127
///
 
2128
/// where each argument Vj represents the application to the
 
2129
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
2130
/// fixed functions (coefficients).
 
2131
 
 
2132
class heat_form_0: public ufc::form
 
2133
{
 
2134
public:
 
2135
 
 
2136
  /// Constructor
 
2137
  heat_form_0() : ufc::form()
 
2138
  {
 
2139
    // Do nothing
 
2140
  }
 
2141
 
 
2142
  /// Destructor
 
2143
  virtual ~heat_form_0()
 
2144
  {
 
2145
    // Do nothing
 
2146
  }
 
2147
 
 
2148
  /// Return a string identifying the form
 
2149
  virtual const char* signature() const
 
2150
  {
 
2151
    return "Form([Integral(Sum(Product(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0), Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 1)), Product(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0), Constant(Cell('triangle', Space(2)), 1)))), Measure('cell', 0, None))])";
 
2152
  }
 
2153
 
 
2154
  /// Return the rank of the global tensor (r)
 
2155
  virtual unsigned int rank() const
 
2156
  {
 
2157
    return 2;
 
2158
  }
 
2159
 
 
2160
  /// Return the number of coefficients (n)
 
2161
  virtual unsigned int num_coefficients() const
 
2162
  {
 
2163
    return 2;
 
2164
  }
 
2165
 
 
2166
  /// Return the number of cell domains
 
2167
  virtual unsigned int num_cell_domains() const
 
2168
  {
 
2169
    return 1;
 
2170
  }
 
2171
 
 
2172
  /// Return the number of exterior facet domains
 
2173
  virtual unsigned int num_exterior_facet_domains() const
 
2174
  {
 
2175
    return 0;
 
2176
  }
 
2177
 
 
2178
  /// Return the number of interior facet domains
 
2179
  virtual unsigned int num_interior_facet_domains() const
 
2180
  {
 
2181
    return 0;
 
2182
  }
 
2183
 
 
2184
  /// Create a new finite element for argument function i
 
2185
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
2186
  {
 
2187
    switch (i)
 
2188
    {
 
2189
    case 0:
 
2190
      {
 
2191
        return new heat_finite_element_1();
 
2192
        break;
 
2193
      }
 
2194
    case 1:
 
2195
      {
 
2196
        return new heat_finite_element_1();
 
2197
        break;
 
2198
      }
 
2199
    case 2:
 
2200
      {
 
2201
        return new heat_finite_element_1();
 
2202
        break;
 
2203
      }
 
2204
    case 3:
 
2205
      {
 
2206
        return new heat_finite_element_0();
 
2207
        break;
 
2208
      }
 
2209
    }
 
2210
    
 
2211
    return 0;
 
2212
  }
 
2213
 
 
2214
  /// Create a new dofmap for argument function i
 
2215
  virtual ufc::dofmap* create_dofmap(unsigned int i) const
 
2216
  {
 
2217
    switch (i)
 
2218
    {
 
2219
    case 0:
 
2220
      {
 
2221
        return new heat_dofmap_1();
 
2222
        break;
 
2223
      }
 
2224
    case 1:
 
2225
      {
 
2226
        return new heat_dofmap_1();
 
2227
        break;
 
2228
      }
 
2229
    case 2:
 
2230
      {
 
2231
        return new heat_dofmap_1();
 
2232
        break;
 
2233
      }
 
2234
    case 3:
 
2235
      {
 
2236
        return new heat_dofmap_0();
 
2237
        break;
 
2238
      }
 
2239
    }
 
2240
    
 
2241
    return 0;
 
2242
  }
 
2243
 
 
2244
  /// Create a new cell integral on sub domain i
 
2245
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
2246
  {
 
2247
    switch (i)
 
2248
    {
 
2249
    case 0:
 
2250
      {
 
2251
        return new heat_cell_integral_0_0();
 
2252
        break;
 
2253
      }
 
2254
    }
 
2255
    
 
2256
    return 0;
 
2257
  }
 
2258
 
 
2259
  /// Create a new exterior facet integral on sub domain i
 
2260
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
2261
  {
 
2262
    return 0;
 
2263
  }
 
2264
 
 
2265
  /// Create a new interior facet integral on sub domain i
 
2266
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
2267
  {
 
2268
    return 0;
 
2269
  }
 
2270
 
 
2271
};
 
2272
 
 
2273
/// This class defines the interface for the assembly of the global
 
2274
/// tensor corresponding to a form with r + n arguments, that is, a
 
2275
/// mapping
 
2276
///
 
2277
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
2278
///
 
2279
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
2280
/// global tensor A is defined by
 
2281
///
 
2282
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
2283
///
 
2284
/// where each argument Vj represents the application to the
 
2285
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
2286
/// fixed functions (coefficients).
 
2287
 
 
2288
class heat_form_1: public ufc::form
 
2289
{
 
2290
public:
 
2291
 
 
2292
  /// Constructor
 
2293
  heat_form_1() : ufc::form()
 
2294
  {
 
2295
    // Do nothing
 
2296
  }
 
2297
 
 
2298
  /// Destructor
 
2299
  virtual ~heat_form_1()
 
2300
  {
 
2301
    // Do nothing
 
2302
  }
 
2303
 
 
2304
  /// Return a string identifying the form
 
2305
  virtual const char* signature() const
 
2306
  {
 
2307
    return "Form([Integral(Sum(Product(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0), Coefficient(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0)), Product(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 0), Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', Space(2)), 1, None), 1), Constant(Cell('triangle', Space(2)), 2)))), Measure('cell', 0, None))])";
 
2308
  }
 
2309
 
 
2310
  /// Return the rank of the global tensor (r)
 
2311
  virtual unsigned int rank() const
 
2312
  {
 
2313
    return 1;
 
2314
  }
 
2315
 
 
2316
  /// Return the number of coefficients (n)
 
2317
  virtual unsigned int num_coefficients() const
 
2318
  {
 
2319
    return 3;
 
2320
  }
 
2321
 
 
2322
  /// Return the number of cell domains
 
2323
  virtual unsigned int num_cell_domains() const
 
2324
  {
 
2325
    return 1;
 
2326
  }
 
2327
 
 
2328
  /// Return the number of exterior facet domains
 
2329
  virtual unsigned int num_exterior_facet_domains() const
 
2330
  {
 
2331
    return 0;
 
2332
  }
 
2333
 
 
2334
  /// Return the number of interior facet domains
 
2335
  virtual unsigned int num_interior_facet_domains() const
 
2336
  {
 
2337
    return 0;
 
2338
  }
 
2339
 
 
2340
  /// Create a new finite element for argument function i
 
2341
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
2342
  {
 
2343
    switch (i)
 
2344
    {
 
2345
    case 0:
 
2346
      {
 
2347
        return new heat_finite_element_1();
 
2348
        break;
 
2349
      }
 
2350
    case 1:
 
2351
      {
 
2352
        return new heat_finite_element_1();
 
2353
        break;
 
2354
      }
 
2355
    case 2:
 
2356
      {
 
2357
        return new heat_finite_element_1();
 
2358
        break;
 
2359
      }
 
2360
    case 3:
 
2361
      {
 
2362
        return new heat_finite_element_0();
 
2363
        break;
 
2364
      }
 
2365
    }
 
2366
    
 
2367
    return 0;
 
2368
  }
 
2369
 
 
2370
  /// Create a new dofmap for argument function i
 
2371
  virtual ufc::dofmap* create_dofmap(unsigned int i) const
 
2372
  {
 
2373
    switch (i)
 
2374
    {
 
2375
    case 0:
 
2376
      {
 
2377
        return new heat_dofmap_1();
 
2378
        break;
 
2379
      }
 
2380
    case 1:
 
2381
      {
 
2382
        return new heat_dofmap_1();
 
2383
        break;
 
2384
      }
 
2385
    case 2:
 
2386
      {
 
2387
        return new heat_dofmap_1();
 
2388
        break;
 
2389
      }
 
2390
    case 3:
 
2391
      {
 
2392
        return new heat_dofmap_0();
 
2393
        break;
 
2394
      }
 
2395
    }
 
2396
    
 
2397
    return 0;
 
2398
  }
 
2399
 
 
2400
  /// Create a new cell integral on sub domain i
 
2401
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
2402
  {
 
2403
    switch (i)
 
2404
    {
 
2405
    case 0:
 
2406
      {
 
2407
        return new heat_cell_integral_1_0();
 
2408
        break;
 
2409
      }
 
2410
    }
 
2411
    
 
2412
    return 0;
 
2413
  }
 
2414
 
 
2415
  /// Create a new exterior facet integral on sub domain i
 
2416
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
2417
  {
 
2418
    return 0;
 
2419
  }
 
2420
 
 
2421
  /// Create a new interior facet integral on sub domain i
 
2422
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
2423
  {
 
2424
    return 0;
 
2425
  }
 
2426
 
 
2427
};
 
2428
 
 
2429
#endif