← Back to branch summary

~ubuntu-branches/ubuntu/trusty/ffc/trusty

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2010-02-03 20:22:35 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100203202235-fe8d0kajuvgy2sqn
Tags: 0.9.0-1
* New upstream release.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/regression/references/Heat.h
 
1
// This code conforms with the UFC specification version 1.4
 
2
// and was automatically generated by FFC version 0.9.0.
 
3
 
 
4
#ifndef __HEAT_H
 
5
#define __HEAT_H
 
6
 
 
7
#include <cmath>
 
8
#include <stdexcept>
 
9
#include <fstream>
 
10
#include <ufc.h>
 
11
 
 
12
/// This class defines the interface for a finite element.
 
13
 
 
14
class heat_finite_element_0: public ufc::finite_element
 
15
{
 
16
public:
 
17
 
 
18
  /// Constructor
 
19
  heat_finite_element_0() : ufc::finite_element()
 
20
  {
 
21
    // Do nothing
 
22
  }
 
23
 
 
24
  /// Destructor
 
25
  virtual ~heat_finite_element_0()
 
26
  {
 
27
    // Do nothing
 
28
  }
 
29
 
 
30
  /// Return a string identifying the finite element
 
31
  virtual const char* signature() const
 
32
  {
 
33
    return "FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
 
34
  }
 
35
 
 
36
  /// Return the cell shape
 
37
  virtual ufc::shape cell_shape() const
 
38
  {
 
39
    return ufc::triangle;
 
40
  }
 
41
 
 
42
  /// Return the dimension of the finite element function space
 
43
  virtual unsigned int space_dimension() const
 
44
  {
 
45
    return 1;
 
46
  }
 
47
 
 
48
  /// Return the rank of the value space
 
49
  virtual unsigned int value_rank() const
 
50
  {
 
51
    return 0;
 
52
  }
 
53
 
 
54
  /// Return the dimension of the value space for axis i
 
55
  virtual unsigned int value_dimension(unsigned int i) const
 
56
  {
 
57
    return 1;
 
58
  }
 
59
 
 
60
  /// Evaluate basis function i at given point in cell
 
61
  virtual void evaluate_basis(unsigned int i,
 
62
                              double* values,
 
63
                              const double* coordinates,
 
64
                              const ufc::cell& c) const
 
65
  {
 
66
    // Extract vertex coordinates
 
67
    
 
68
    // Compute Jacobian of affine map from reference cell
 
69
    
 
70
    // Compute determinant of Jacobian
 
71
    
 
72
    // Compute inverse of Jacobian
 
73
    
 
74
    // Compute constants
 
75
    
 
76
    // Get coordinates and map to the reference (FIAT) element
 
77
    
 
78
    // Reset values
 
79
    *values = 0.00000000;
 
80
    
 
81
    // Map degree of freedom to element degree of freedom
 
82
    const unsigned int dof = i;
 
83
    
 
84
    // Array of basisvalues
 
85
    double basisvalues[1] = {0.00000000};
 
86
    
 
87
    // Declare helper variables
 
88
    
 
89
    // Compute basisvalues
 
90
    basisvalues[0] = 1.00000000;
 
91
    
 
92
    // Table(s) of coefficients
 
93
    static const double coefficients0[1][1] = \
 
94
    {{1.00000000}};
 
95
    
 
96
    // Compute value(s).
 
97
    for (unsigned int r = 0; r < 1; r++)
 
98
    {
 
99
      *values += coefficients0[dof][r]*basisvalues[r];
 
100
    }// end loop over 'r'
 
101
  }
 
102
 
 
103
  /// Evaluate all basis functions at given point in cell
 
104
  virtual void evaluate_basis_all(double* values,
 
105
                                  const double* coordinates,
 
106
                                  const ufc::cell& c) const
 
107
  {
 
108
    // Element is constant, calling evaluate_basis.
 
109
    evaluate_basis(0, values, coordinates, c);
 
110
  }
 
111
 
 
112
  /// Evaluate order n derivatives of basis function i at given point in cell
 
113
  virtual void evaluate_basis_derivatives(unsigned int i,
 
114
                                          unsigned int n,
 
115
                                          double* values,
 
116
                                          const double* coordinates,
 
117
                                          const ufc::cell& c) const
 
118
  {
 
119
    // Extract vertex coordinates
 
120
    const double * const * x = c.coordinates;
 
121
    
 
122
    // Compute Jacobian of affine map from reference cell
 
123
    const double J_00 = x[1][0] - x[0][0];
 
124
    const double J_01 = x[2][0] - x[0][0];
 
125
    const double J_10 = x[1][1] - x[0][1];
 
126
    const double J_11 = x[2][1] - x[0][1];
 
127
    
 
128
    // Compute determinant of Jacobian
 
129
    double detJ = J_00*J_11 - J_01*J_10;
 
130
    
 
131
    // Compute inverse of Jacobian
 
132
    const double K_00 =  J_11 / detJ;
 
133
    const double K_01 = -J_01 / detJ;
 
134
    const double K_10 = -J_10 / detJ;
 
135
    const double K_11 =  J_00 / detJ;
 
136
    
 
137
    // Compute constants
 
138
    
 
139
    // Get coordinates and map to the reference (FIAT) element
 
140
    
 
141
    // Compute number of derivatives.
 
142
    unsigned int  num_derivatives = 1;
 
143
    
 
144
    for (unsigned int r = 0; r < n; r++)
 
145
    {
 
146
      num_derivatives *= 2;
 
147
    }// end loop over 'r'
 
148
    
 
149
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
150
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
151
    for (unsigned int row = 0; row < num_derivatives; row++)
 
152
    {
 
153
      combinations[row] = new unsigned int [n];
 
154
      for (unsigned int col = 0; col < n; col++)
 
155
        combinations[row][col] = 0;
 
156
    }
 
157
    
 
158
    // Generate combinations of derivatives
 
159
    for (unsigned int row = 1; row < num_derivatives; row++)
 
160
    {
 
161
      for (unsigned int num = 0; num < row; num++)
 
162
      {
 
163
        for (unsigned int col = n-1; col+1 > 0; col--)
 
164
        {
 
165
          if (combinations[row][col] + 1 > 1)
 
166
            combinations[row][col] = 0;
 
167
          else
 
168
          {
 
169
            combinations[row][col] += 1;
 
170
            break;
 
171
          }
 
172
        }
 
173
      }
 
174
    }
 
175
    
 
176
    // Compute inverse of Jacobian
 
177
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
 
178
    
 
179
    // Declare transformation matrix
 
180
    // Declare pointer to two dimensional array and initialise
 
181
    double **transform = new double *[num_derivatives];
 
182
    
 
183
    for (unsigned int j = 0; j < num_derivatives; j++)
 
184
    {
 
185
      transform[j] = new double [num_derivatives];
 
186
      for (unsigned int k = 0; k < num_derivatives; k++)
 
187
        transform[j][k] = 1;
 
188
    }
 
189
    
 
190
    // Construct transformation matrix
 
191
    for (unsigned int row = 0; row < num_derivatives; row++)
 
192
    {
 
193
      for (unsigned int col = 0; col < num_derivatives; col++)
 
194
      {
 
195
        for (unsigned int k = 0; k < n; k++)
 
196
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
197
      }
 
198
    }
 
199
    
 
200
    // Reset values. Assuming that values is always an array.
 
201
    for (unsigned int r = 0; r < num_derivatives; r++)
 
202
    {
 
203
      values[r] = 0.00000000;
 
204
    }// end loop over 'r'
 
205
    
 
206
    // Map degree of freedom to element degree of freedom
 
207
    const unsigned int dof = i;
 
208
    
 
209
    // Array of basisvalues
 
210
    double basisvalues[1] = {0.00000000};
 
211
    
 
212
    // Declare helper variables
 
213
    
 
214
    // Compute basisvalues
 
215
    basisvalues[0] = 1.00000000;
 
216
    
 
217
    // Table(s) of coefficients
 
218
    static const double coefficients0[1][1] = \
 
219
    {{1.00000000}};
 
220
    
 
221
    // Tables of derivatives of the polynomial base (transpose).
 
222
    static const double dmats0[1][1] = \
 
223
    {{0.00000000}};
 
224
    
 
225
    static const double dmats1[1][1] = \
 
226
    {{0.00000000}};
 
227
    
 
228
    // Compute reference derivatives
 
229
    // Declare pointer to array of derivatives on FIAT element
 
230
    double *derivatives = new double [num_derivatives];
 
231
    for (unsigned int r = 0; r < num_derivatives; r++)
 
232
    {
 
233
      derivatives[r] = 0.00000000;
 
234
    }// end loop over 'r'
 
235
    
 
236
    // Declare derivative matrix (of polynomial basis).
 
237
    double dmats[1][1] = \
 
238
    {{1.00000000}};
 
239
    
 
240
    // Declare (auxiliary) derivative matrix (of polynomial basis).
 
241
    double dmats_old[1][1] = \
 
242
    {{1.00000000}};
 
243
    
 
244
    // Loop possible derivatives.
 
245
    for (unsigned int r = 0; r < num_derivatives; r++)
 
246
    {
 
247
      // Resetting dmats values to compute next derivative.
 
248
      for (unsigned int t = 0; t < 1; t++)
 
249
      {
 
250
        for (unsigned int u = 0; u < 1; u++)
 
251
        {
 
252
          dmats[t][u] = 0.00000000;
 
253
          if (t == u)
 
254
          {
 
255
          dmats[t][u] = 1.00000000;
 
256
          }
 
257
          
 
258
        }// end loop over 'u'
 
259
      }// end loop over 't'
 
260
      
 
261
      // Looping derivative order to generate dmats.
 
262
      for (unsigned int s = 0; s < n; s++)
 
263
      {
 
264
        // Updating dmats_old with new values and resetting dmats.
 
265
        for (unsigned int t = 0; t < 1; t++)
 
266
        {
 
267
          for (unsigned int u = 0; u < 1; u++)
 
268
          {
 
269
            dmats_old[t][u] = dmats[t][u];
 
270
            dmats[t][u] = 0.00000000;
 
271
          }// end loop over 'u'
 
272
        }// end loop over 't'
 
273
        
 
274
        // Update dmats using an inner product.
 
275
        if (combinations[r][s] == 0)
 
276
        {
 
277
        for (unsigned int t = 0; t < 1; t++)
 
278
        {
 
279
          for (unsigned int u = 0; u < 1; u++)
 
280
          {
 
281
            for (unsigned int tu = 0; tu < 1; tu++)
 
282
            {
 
283
              dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
284
            }// end loop over 'tu'
 
285
          }// end loop over 'u'
 
286
        }// end loop over 't'
 
287
        }
 
288
        
 
289
        if (combinations[r][s] == 1)
 
290
        {
 
291
        for (unsigned int t = 0; t < 1; t++)
 
292
        {
 
293
          for (unsigned int u = 0; u < 1; u++)
 
294
          {
 
295
            for (unsigned int tu = 0; tu < 1; tu++)
 
296
            {
 
297
              dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
298
            }// end loop over 'tu'
 
299
          }// end loop over 'u'
 
300
        }// end loop over 't'
 
301
        }
 
302
        
 
303
      }// end loop over 's'
 
304
      for (unsigned int s = 0; s < 1; s++)
 
305
      {
 
306
        for (unsigned int t = 0; t < 1; t++)
 
307
        {
 
308
          derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
 
309
        }// end loop over 't'
 
310
      }// end loop over 's'
 
311
    }// end loop over 'r'
 
312
    
 
313
    // Transform derivatives back to physical element
 
314
    for (unsigned int row = 0; row < num_derivatives; row++)
 
315
    {
 
316
      for (unsigned int col = 0; col < num_derivatives; col++)
 
317
      {
 
318
        values[row] += transform[row][col]*derivatives[col];
 
319
      }
 
320
    }
 
321
    
 
322
    // Delete pointer to array of derivatives on FIAT element
 
323
    delete [] derivatives;
 
324
    
 
325
    // Delete pointer to array of combinations of derivatives and transform
 
326
    for (unsigned int r = 0; r < num_derivatives; r++)
 
327
    {
 
328
      delete [] combinations[r];
 
329
      delete [] transform[r];
 
330
    }// end loop over 'r'
 
331
    delete [] combinations;
 
332
    delete [] transform;
 
333
  }
 
334
 
 
335
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
336
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
337
                                              double* values,
 
338
                                              const double* coordinates,
 
339
                                              const ufc::cell& c) const
 
340
  {
 
341
    // Element is constant, calling evaluate_basis_derivatives.
 
342
    evaluate_basis_derivatives(0, n, values, coordinates, c);
 
343
  }
 
344
 
 
345
  /// Evaluate linear functional for dof i on the function f
 
346
  virtual double evaluate_dof(unsigned int i,
 
347
                              const ufc::function& f,
 
348
                              const ufc::cell& c) const
 
349
  {
 
350
    // Declare variables for result of evaluation
 
351
    double vals[1];
 
352
    
 
353
    // Declare variable for physical coordinates
 
354
    double y[2];
 
355
    
 
356
    const double * const * x = c.coordinates;
 
357
    switch (i)
 
358
    {
 
359
    case 0:
 
360
      {
 
361
        y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
 
362
      y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
 
363
      f.evaluate(vals, y, c);
 
364
      return vals[0];
 
365
        break;
 
366
      }
 
367
    }
 
368
    
 
369
    return 0.0;
 
370
  }
 
371
 
 
372
  /// Evaluate linear functionals for all dofs on the function f
 
373
  virtual void evaluate_dofs(double* values,
 
374
                             const ufc::function& f,
 
375
                             const ufc::cell& c) const
 
376
  {
 
377
    // Declare variables for result of evaluation
 
378
    double vals[1];
 
379
    
 
380
    // Declare variable for physical coordinates
 
381
    double y[2];
 
382
    
 
383
    const double * const * x = c.coordinates;
 
384
    y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
 
385
    y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
 
386
    f.evaluate(vals, y, c);
 
387
    values[0] = vals[0];
 
388
  }
 
389
 
 
390
  /// Interpolate vertex values from dof values
 
391
  virtual void interpolate_vertex_values(double* vertex_values,
 
392
                                         const double* dof_values,
 
393
                                         const ufc::cell& c) const
 
394
  {
 
395
    // Evaluate function and change variables
 
396
    vertex_values[0] = dof_values[0];
 
397
    vertex_values[1] = dof_values[0];
 
398
    vertex_values[2] = dof_values[0];
 
399
  }
 
400
 
 
401
  /// Return the number of sub elements (for a mixed element)
 
402
  virtual unsigned int num_sub_elements() const
 
403
  {
 
404
    return 0;
 
405
  }
 
406
 
 
407
  /// Create a new finite element for sub element i (for a mixed element)
 
408
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
409
  {
 
410
    return 0;
 
411
  }
 
412
 
 
413
};
 
414
 
 
415
/// This class defines the interface for a finite element.
 
416
 
 
417
class heat_finite_element_1: public ufc::finite_element
 
418
{
 
419
public:
 
420
 
 
421
  /// Constructor
 
422
  heat_finite_element_1() : ufc::finite_element()
 
423
  {
 
424
    // Do nothing
 
425
  }
 
426
 
 
427
  /// Destructor
 
428
  virtual ~heat_finite_element_1()
 
429
  {
 
430
    // Do nothing
 
431
  }
 
432
 
 
433
  /// Return a string identifying the finite element
 
434
  virtual const char* signature() const
 
435
  {
 
436
    return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
 
437
  }
 
438
 
 
439
  /// Return the cell shape
 
440
  virtual ufc::shape cell_shape() const
 
441
  {
 
442
    return ufc::triangle;
 
443
  }
 
444
 
 
445
  /// Return the dimension of the finite element function space
 
446
  virtual unsigned int space_dimension() const
 
447
  {
 
448
    return 3;
 
449
  }
 
450
 
 
451
  /// Return the rank of the value space
 
452
  virtual unsigned int value_rank() const
 
453
  {
 
454
    return 0;
 
455
  }
 
456
 
 
457
  /// Return the dimension of the value space for axis i
 
458
  virtual unsigned int value_dimension(unsigned int i) const
 
459
  {
 
460
    return 1;
 
461
  }
 
462
 
 
463
  /// Evaluate basis function i at given point in cell
 
464
  virtual void evaluate_basis(unsigned int i,
 
465
                              double* values,
 
466
                              const double* coordinates,
 
467
                              const ufc::cell& c) const
 
468
  {
 
469
    // Extract vertex coordinates
 
470
    const double * const * x = c.coordinates;
 
471
    
 
472
    // Compute Jacobian of affine map from reference cell
 
473
    const double J_00 = x[1][0] - x[0][0];
 
474
    const double J_01 = x[2][0] - x[0][0];
 
475
    const double J_10 = x[1][1] - x[0][1];
 
476
    const double J_11 = x[2][1] - x[0][1];
 
477
    
 
478
    // Compute determinant of Jacobian
 
479
    double detJ = J_00*J_11 - J_01*J_10;
 
480
    
 
481
    // Compute inverse of Jacobian
 
482
    
 
483
    // Compute constants
 
484
    const double C0 = x[1][0] + x[2][0];
 
485
    const double C1 = x[1][1] + x[2][1];
 
486
    
 
487
    // Get coordinates and map to the reference (FIAT) element
 
488
    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
 
489
    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
 
490
    
 
491
    // Reset values
 
492
    *values = 0.00000000;
 
493
    
 
494
    // Map degree of freedom to element degree of freedom
 
495
    const unsigned int dof = i;
 
496
    
 
497
    // Array of basisvalues
 
498
    double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
499
    
 
500
    // Declare helper variables
 
501
    unsigned int rr = 0;
 
502
    unsigned int ss = 0;
 
503
    double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
504
    
 
505
    // Compute basisvalues
 
506
    basisvalues[0] = 1.00000000;
 
507
    basisvalues[1] = tmp0;
 
508
    for (unsigned int r = 0; r < 1; r++)
 
509
    {
 
510
      rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
511
      ss = r*(r + 1)/2;
 
512
      basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
513
    }// end loop over 'r'
 
514
    for (unsigned int r = 0; r < 2; r++)
 
515
    {
 
516
      for (unsigned int s = 0; s < 2 - r; s++)
 
517
      {
 
518
        rr = (r + s)*(r + s + 1)/2 + s;
 
519
        basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
520
      }// end loop over 's'
 
521
    }// end loop over 'r'
 
522
    
 
523
    // Table(s) of coefficients
 
524
    static const double coefficients0[3][3] = \
 
525
    {{0.47140452, -0.28867513, -0.16666667},
 
526
    {0.47140452, 0.28867513, -0.16666667},
 
527
    {0.47140452, 0.00000000, 0.33333333}};
 
528
    
 
529
    // Compute value(s).
 
530
    for (unsigned int r = 0; r < 3; r++)
 
531
    {
 
532
      *values += coefficients0[dof][r]*basisvalues[r];
 
533
    }// end loop over 'r'
 
534
  }
 
535
 
 
536
  /// Evaluate all basis functions at given point in cell
 
537
  virtual void evaluate_basis_all(double* values,
 
538
                                  const double* coordinates,
 
539
                                  const ufc::cell& c) const
 
540
  {
 
541
    // Helper variable to hold values of a single dof.
 
542
    double dof_values = 0.00000000;
 
543
    
 
544
    // Loop dofs and call evaluate_basis.
 
545
    for (unsigned int r = 0; r < 3; r++)
 
546
    {
 
547
      evaluate_basis(r, &dof_values, coordinates, c);
 
548
      values[r] = dof_values;
 
549
    }// end loop over 'r'
 
550
  }
 
551
 
 
552
  /// Evaluate order n derivatives of basis function i at given point in cell
 
553
  virtual void evaluate_basis_derivatives(unsigned int i,
 
554
                                          unsigned int n,
 
555
                                          double* values,
 
556
                                          const double* coordinates,
 
557
                                          const ufc::cell& c) const
 
558
  {
 
559
    // Extract vertex coordinates
 
560
    const double * const * x = c.coordinates;
 
561
    
 
562
    // Compute Jacobian of affine map from reference cell
 
563
    const double J_00 = x[1][0] - x[0][0];
 
564
    const double J_01 = x[2][0] - x[0][0];
 
565
    const double J_10 = x[1][1] - x[0][1];
 
566
    const double J_11 = x[2][1] - x[0][1];
 
567
    
 
568
    // Compute determinant of Jacobian
 
569
    double detJ = J_00*J_11 - J_01*J_10;
 
570
    
 
571
    // Compute inverse of Jacobian
 
572
    const double K_00 =  J_11 / detJ;
 
573
    const double K_01 = -J_01 / detJ;
 
574
    const double K_10 = -J_10 / detJ;
 
575
    const double K_11 =  J_00 / detJ;
 
576
    
 
577
    // Compute constants
 
578
    const double C0 = x[1][0] + x[2][0];
 
579
    const double C1 = x[1][1] + x[2][1];
 
580
    
 
581
    // Get coordinates and map to the reference (FIAT) element
 
582
    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
 
583
    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
 
584
    
 
585
    // Compute number of derivatives.
 
586
    unsigned int  num_derivatives = 1;
 
587
    
 
588
    for (unsigned int r = 0; r < n; r++)
 
589
    {
 
590
      num_derivatives *= 2;
 
591
    }// end loop over 'r'
 
592
    
 
593
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
594
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
595
    for (unsigned int row = 0; row < num_derivatives; row++)
 
596
    {
 
597
      combinations[row] = new unsigned int [n];
 
598
      for (unsigned int col = 0; col < n; col++)
 
599
        combinations[row][col] = 0;
 
600
    }
 
601
    
 
602
    // Generate combinations of derivatives
 
603
    for (unsigned int row = 1; row < num_derivatives; row++)
 
604
    {
 
605
      for (unsigned int num = 0; num < row; num++)
 
606
      {
 
607
        for (unsigned int col = n-1; col+1 > 0; col--)
 
608
        {
 
609
          if (combinations[row][col] + 1 > 1)
 
610
            combinations[row][col] = 0;
 
611
          else
 
612
          {
 
613
            combinations[row][col] += 1;
 
614
            break;
 
615
          }
 
616
        }
 
617
      }
 
618
    }
 
619
    
 
620
    // Compute inverse of Jacobian
 
621
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
 
622
    
 
623
    // Declare transformation matrix
 
624
    // Declare pointer to two dimensional array and initialise
 
625
    double **transform = new double *[num_derivatives];
 
626
    
 
627
    for (unsigned int j = 0; j < num_derivatives; j++)
 
628
    {
 
629
      transform[j] = new double [num_derivatives];
 
630
      for (unsigned int k = 0; k < num_derivatives; k++)
 
631
        transform[j][k] = 1;
 
632
    }
 
633
    
 
634
    // Construct transformation matrix
 
635
    for (unsigned int row = 0; row < num_derivatives; row++)
 
636
    {
 
637
      for (unsigned int col = 0; col < num_derivatives; col++)
 
638
      {
 
639
        for (unsigned int k = 0; k < n; k++)
 
640
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
641
      }
 
642
    }
 
643
    
 
644
    // Reset values. Assuming that values is always an array.
 
645
    for (unsigned int r = 0; r < num_derivatives; r++)
 
646
    {
 
647
      values[r] = 0.00000000;
 
648
    }// end loop over 'r'
 
649
    
 
650
    // Map degree of freedom to element degree of freedom
 
651
    const unsigned int dof = i;
 
652
    
 
653
    // Array of basisvalues
 
654
    double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
655
    
 
656
    // Declare helper variables
 
657
    unsigned int rr = 0;
 
658
    unsigned int ss = 0;
 
659
    double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
660
    
 
661
    // Compute basisvalues
 
662
    basisvalues[0] = 1.00000000;
 
663
    basisvalues[1] = tmp0;
 
664
    for (unsigned int r = 0; r < 1; r++)
 
665
    {
 
666
      rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
667
      ss = r*(r + 1)/2;
 
668
      basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
669
    }// end loop over 'r'
 
670
    for (unsigned int r = 0; r < 2; r++)
 
671
    {
 
672
      for (unsigned int s = 0; s < 2 - r; s++)
 
673
      {
 
674
        rr = (r + s)*(r + s + 1)/2 + s;
 
675
        basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
676
      }// end loop over 's'
 
677
    }// end loop over 'r'
 
678
    
 
679
    // Table(s) of coefficients
 
680
    static const double coefficients0[3][3] = \
 
681
    {{0.47140452, -0.28867513, -0.16666667},
 
682
    {0.47140452, 0.28867513, -0.16666667},
 
683
    {0.47140452, 0.00000000, 0.33333333}};
 
684
    
 
685
    // Tables of derivatives of the polynomial base (transpose).
 
686
    static const double dmats0[3][3] = \
 
687
    {{0.00000000, 0.00000000, 0.00000000},
 
688
    {4.89897949, 0.00000000, 0.00000000},
 
689
    {0.00000000, 0.00000000, 0.00000000}};
 
690
    
 
691
    static const double dmats1[3][3] = \
 
692
    {{0.00000000, 0.00000000, 0.00000000},
 
693
    {2.44948974, 0.00000000, 0.00000000},
 
694
    {4.24264069, 0.00000000, 0.00000000}};
 
695
    
 
696
    // Compute reference derivatives
 
697
    // Declare pointer to array of derivatives on FIAT element
 
698
    double *derivatives = new double [num_derivatives];
 
699
    for (unsigned int r = 0; r < num_derivatives; r++)
 
700
    {
 
701
      derivatives[r] = 0.00000000;
 
702
    }// end loop over 'r'
 
703
    
 
704
    // Declare derivative matrix (of polynomial basis).
 
705
    double dmats[3][3] = \
 
706
    {{1.00000000, 0.00000000, 0.00000000},
 
707
    {0.00000000, 1.00000000, 0.00000000},
 
708
    {0.00000000, 0.00000000, 1.00000000}};
 
709
    
 
710
    // Declare (auxiliary) derivative matrix (of polynomial basis).
 
711
    double dmats_old[3][3] = \
 
712
    {{1.00000000, 0.00000000, 0.00000000},
 
713
    {0.00000000, 1.00000000, 0.00000000},
 
714
    {0.00000000, 0.00000000, 1.00000000}};
 
715
    
 
716
    // Loop possible derivatives.
 
717
    for (unsigned int r = 0; r < num_derivatives; r++)
 
718
    {
 
719
      // Resetting dmats values to compute next derivative.
 
720
      for (unsigned int t = 0; t < 3; t++)
 
721
      {
 
722
        for (unsigned int u = 0; u < 3; u++)
 
723
        {
 
724
          dmats[t][u] = 0.00000000;
 
725
          if (t == u)
 
726
          {
 
727
          dmats[t][u] = 1.00000000;
 
728
          }
 
729
          
 
730
        }// end loop over 'u'
 
731
      }// end loop over 't'
 
732
      
 
733
      // Looping derivative order to generate dmats.
 
734
      for (unsigned int s = 0; s < n; s++)
 
735
      {
 
736
        // Updating dmats_old with new values and resetting dmats.
 
737
        for (unsigned int t = 0; t < 3; t++)
 
738
        {
 
739
          for (unsigned int u = 0; u < 3; u++)
 
740
          {
 
741
            dmats_old[t][u] = dmats[t][u];
 
742
            dmats[t][u] = 0.00000000;
 
743
          }// end loop over 'u'
 
744
        }// end loop over 't'
 
745
        
 
746
        // Update dmats using an inner product.
 
747
        if (combinations[r][s] == 0)
 
748
        {
 
749
        for (unsigned int t = 0; t < 3; t++)
 
750
        {
 
751
          for (unsigned int u = 0; u < 3; u++)
 
752
          {
 
753
            for (unsigned int tu = 0; tu < 3; tu++)
 
754
            {
 
755
              dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
756
            }// end loop over 'tu'
 
757
          }// end loop over 'u'
 
758
        }// end loop over 't'
 
759
        }
 
760
        
 
761
        if (combinations[r][s] == 1)
 
762
        {
 
763
        for (unsigned int t = 0; t < 3; t++)
 
764
        {
 
765
          for (unsigned int u = 0; u < 3; u++)
 
766
          {
 
767
            for (unsigned int tu = 0; tu < 3; tu++)
 
768
            {
 
769
              dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
770
            }// end loop over 'tu'
 
771
          }// end loop over 'u'
 
772
        }// end loop over 't'
 
773
        }
 
774
        
 
775
      }// end loop over 's'
 
776
      for (unsigned int s = 0; s < 3; s++)
 
777
      {
 
778
        for (unsigned int t = 0; t < 3; t++)
 
779
        {
 
780
          derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
 
781
        }// end loop over 't'
 
782
      }// end loop over 's'
 
783
    }// end loop over 'r'
 
784
    
 
785
    // Transform derivatives back to physical element
 
786
    for (unsigned int row = 0; row < num_derivatives; row++)
 
787
    {
 
788
      for (unsigned int col = 0; col < num_derivatives; col++)
 
789
      {
 
790
        values[row] += transform[row][col]*derivatives[col];
 
791
      }
 
792
    }
 
793
    
 
794
    // Delete pointer to array of derivatives on FIAT element
 
795
    delete [] derivatives;
 
796
    
 
797
    // Delete pointer to array of combinations of derivatives and transform
 
798
    for (unsigned int r = 0; r < num_derivatives; r++)
 
799
    {
 
800
      delete [] combinations[r];
 
801
      delete [] transform[r];
 
802
    }// end loop over 'r'
 
803
    delete [] combinations;
 
804
    delete [] transform;
 
805
  }
 
806
 
 
807
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
808
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
809
                                              double* values,
 
810
                                              const double* coordinates,
 
811
                                              const ufc::cell& c) const
 
812
  {
 
813
    // Compute number of derivatives.
 
814
    unsigned int  num_derivatives = 1;
 
815
    
 
816
    for (unsigned int r = 0; r < n; r++)
 
817
    {
 
818
      num_derivatives *= 2;
 
819
    }// end loop over 'r'
 
820
    
 
821
    // Helper variable to hold values of a single dof.
 
822
    double *dof_values = new double [num_derivatives];
 
823
    for (unsigned int r = 0; r < num_derivatives; r++)
 
824
    {
 
825
      dof_values[r] = 0.00000000;
 
826
    }// end loop over 'r'
 
827
    
 
828
    // Loop dofs and call evaluate_basis_derivatives.
 
829
    for (unsigned int r = 0; r < 3; r++)
 
830
    {
 
831
      evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
 
832
      for (unsigned int s = 0; s < num_derivatives; s++)
 
833
      {
 
834
        values[r*num_derivatives + s] = dof_values[s];
 
835
      }// end loop over 's'
 
836
    }// end loop over 'r'
 
837
    
 
838
    // Delete pointer.
 
839
    delete [] dof_values;
 
840
  }
 
841
 
 
842
  /// Evaluate linear functional for dof i on the function f
 
843
  virtual double evaluate_dof(unsigned int i,
 
844
                              const ufc::function& f,
 
845
                              const ufc::cell& c) const
 
846
  {
 
847
    // Declare variables for result of evaluation
 
848
    double vals[1];
 
849
    
 
850
    // Declare variable for physical coordinates
 
851
    double y[2];
 
852
    
 
853
    const double * const * x = c.coordinates;
 
854
    switch (i)
 
855
    {
 
856
    case 0:
 
857
      {
 
858
        y[0] = x[0][0];
 
859
      y[1] = x[0][1];
 
860
      f.evaluate(vals, y, c);
 
861
      return vals[0];
 
862
        break;
 
863
      }
 
864
    case 1:
 
865
      {
 
866
        y[0] = x[1][0];
 
867
      y[1] = x[1][1];
 
868
      f.evaluate(vals, y, c);
 
869
      return vals[0];
 
870
        break;
 
871
      }
 
872
    case 2:
 
873
      {
 
874
        y[0] = x[2][0];
 
875
      y[1] = x[2][1];
 
876
      f.evaluate(vals, y, c);
 
877
      return vals[0];
 
878
        break;
 
879
      }
 
880
    }
 
881
    
 
882
    return 0.0;
 
883
  }
 
884
 
 
885
  /// Evaluate linear functionals for all dofs on the function f
 
886
  virtual void evaluate_dofs(double* values,
 
887
                             const ufc::function& f,
 
888
                             const ufc::cell& c) const
 
889
  {
 
890
    // Declare variables for result of evaluation
 
891
    double vals[1];
 
892
    
 
893
    // Declare variable for physical coordinates
 
894
    double y[2];
 
895
    
 
896
    const double * const * x = c.coordinates;
 
897
    y[0] = x[0][0];
 
898
    y[1] = x[0][1];
 
899
    f.evaluate(vals, y, c);
 
900
    values[0] = vals[0];
 
901
    y[0] = x[1][0];
 
902
    y[1] = x[1][1];
 
903
    f.evaluate(vals, y, c);
 
904
    values[1] = vals[0];
 
905
    y[0] = x[2][0];
 
906
    y[1] = x[2][1];
 
907
    f.evaluate(vals, y, c);
 
908
    values[2] = vals[0];
 
909
  }
 
910
 
 
911
  /// Interpolate vertex values from dof values
 
912
  virtual void interpolate_vertex_values(double* vertex_values,
 
913
                                         const double* dof_values,
 
914
                                         const ufc::cell& c) const
 
915
  {
 
916
    // Evaluate function and change variables
 
917
    vertex_values[0] = dof_values[0];
 
918
    vertex_values[1] = dof_values[1];
 
919
    vertex_values[2] = dof_values[2];
 
920
  }
 
921
 
 
922
  /// Return the number of sub elements (for a mixed element)
 
923
  virtual unsigned int num_sub_elements() const
 
924
  {
 
925
    return 0;
 
926
  }
 
927
 
 
928
  /// Create a new finite element for sub element i (for a mixed element)
 
929
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
930
  {
 
931
    return 0;
 
932
  }
 
933
 
 
934
};
 
935
 
 
936
/// This class defines the interface for a local-to-global mapping of
 
937
/// degrees of freedom (dofs).
 
938
 
 
939
class heat_dof_map_0: public ufc::dof_map
 
940
{
 
941
private:
 
942
 
 
943
  unsigned int _global_dimension;
 
944
 
 
945
public:
 
946
 
 
947
  /// Constructor
 
948
  heat_dof_map_0() : ufc::dof_map()
 
949
  {
 
950
    _global_dimension = 0;
 
951
  }
 
952
 
 
953
  /// Destructor
 
954
  virtual ~heat_dof_map_0()
 
955
  {
 
956
    // Do nothing
 
957
  }
 
958
 
 
959
  /// Return a string identifying the dof map
 
960
  virtual const char* signature() const
 
961
  {
 
962
    return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
 
963
  }
 
964
 
 
965
  /// Return true iff mesh entities of topological dimension d are needed
 
966
  virtual bool needs_mesh_entities(unsigned int d) const
 
967
  {
 
968
    switch (d)
 
969
    {
 
970
    case 0:
 
971
      {
 
972
        return false;
 
973
        break;
 
974
      }
 
975
    case 1:
 
976
      {
 
977
        return false;
 
978
        break;
 
979
      }
 
980
    case 2:
 
981
      {
 
982
        return true;
 
983
        break;
 
984
      }
 
985
    }
 
986
    
 
987
    return false;
 
988
  }
 
989
 
 
990
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
991
  virtual bool init_mesh(const ufc::mesh& m)
 
992
  {
 
993
    _global_dimension = m.num_entities[2];
 
994
    return false;
 
995
  }
 
996
 
 
997
  /// Initialize dof map for given cell
 
998
  virtual void init_cell(const ufc::mesh& m,
 
999
                         const ufc::cell& c)
 
1000
  {
 
1001
    // Do nothing
 
1002
  }
 
1003
 
 
1004
  /// Finish initialization of dof map for cells
 
1005
  virtual void init_cell_finalize()
 
1006
  {
 
1007
    // Do nothing
 
1008
  }
 
1009
 
 
1010
  /// Return the dimension of the global finite element function space
 
1011
  virtual unsigned int global_dimension() const
 
1012
  {
 
1013
    return _global_dimension;
 
1014
  }
 
1015
 
 
1016
  /// Return the dimension of the local finite element function space for a cell
 
1017
  virtual unsigned int local_dimension(const ufc::cell& c) const
 
1018
  {
 
1019
    return 1;
 
1020
  }
 
1021
 
 
1022
  /// Return the maximum dimension of the local finite element function space
 
1023
  virtual unsigned int max_local_dimension() const
 
1024
  {
 
1025
    return 1;
 
1026
  }
 
1027
 
 
1028
  // Return the geometric dimension of the coordinates this dof map provides
 
1029
  virtual unsigned int geometric_dimension() const
 
1030
  {
 
1031
    return 2;
 
1032
  }
 
1033
 
 
1034
  /// Return the number of dofs on each cell facet
 
1035
  virtual unsigned int num_facet_dofs() const
 
1036
  {
 
1037
    return 0;
 
1038
  }
 
1039
 
 
1040
  /// Return the number of dofs associated with each cell entity of dimension d
 
1041
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1042
  {
 
1043
    switch (d)
 
1044
    {
 
1045
    case 0:
 
1046
      {
 
1047
        return 0;
 
1048
        break;
 
1049
      }
 
1050
    case 1:
 
1051
      {
 
1052
        return 0;
 
1053
        break;
 
1054
      }
 
1055
    case 2:
 
1056
      {
 
1057
        return 1;
 
1058
        break;
 
1059
      }
 
1060
    }
 
1061
    
 
1062
    return 0;
 
1063
  }
 
1064
 
 
1065
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1066
  virtual void tabulate_dofs(unsigned int* dofs,
 
1067
                             const ufc::mesh& m,
 
1068
                             const ufc::cell& c) const
 
1069
  {
 
1070
    dofs[0] = c.entity_indices[2][0];
 
1071
  }
 
1072
 
 
1073
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1074
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1075
                                   unsigned int facet) const
 
1076
  {
 
1077
    switch (facet)
 
1078
    {
 
1079
    case 0:
 
1080
      {
 
1081
        
 
1082
        break;
 
1083
      }
 
1084
    case 1:
 
1085
      {
 
1086
        
 
1087
        break;
 
1088
      }
 
1089
    case 2:
 
1090
      {
 
1091
        
 
1092
        break;
 
1093
      }
 
1094
    }
 
1095
    
 
1096
  }
 
1097
 
 
1098
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1099
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1100
                                    unsigned int d, unsigned int i) const
 
1101
  {
 
1102
    if (d > 2)
 
1103
    {
 
1104
    std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
 
1105
    }
 
1106
    
 
1107
    switch (d)
 
1108
    {
 
1109
    case 0:
 
1110
      {
 
1111
        
 
1112
        break;
 
1113
      }
 
1114
    case 1:
 
1115
      {
 
1116
        
 
1117
        break;
 
1118
      }
 
1119
    case 2:
 
1120
      {
 
1121
        if (i > 0)
 
1122
      {
 
1123
      std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
 
1124
      }
 
1125
      
 
1126
      dofs[0] = 0;
 
1127
        break;
 
1128
      }
 
1129
    }
 
1130
    
 
1131
  }
 
1132
 
 
1133
  /// Tabulate the coordinates of all dofs on a cell
 
1134
  virtual void tabulate_coordinates(double** coordinates,
 
1135
                                    const ufc::cell& c) const
 
1136
  {
 
1137
    const double * const * x = c.coordinates;
 
1138
    
 
1139
    coordinates[0][0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
 
1140
    coordinates[0][1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
 
1141
  }
 
1142
 
 
1143
  /// Return the number of sub dof maps (for a mixed element)
 
1144
  virtual unsigned int num_sub_dof_maps() const
 
1145
  {
 
1146
    return 0;
 
1147
  }
 
1148
 
 
1149
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1150
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1151
  {
 
1152
    return 0;
 
1153
  }
 
1154
 
 
1155
};
 
1156
 
 
1157
/// This class defines the interface for a local-to-global mapping of
 
1158
/// degrees of freedom (dofs).
 
1159
 
 
1160
class heat_dof_map_1: public ufc::dof_map
 
1161
{
 
1162
private:
 
1163
 
 
1164
  unsigned int _global_dimension;
 
1165
 
 
1166
public:
 
1167
 
 
1168
  /// Constructor
 
1169
  heat_dof_map_1() : ufc::dof_map()
 
1170
  {
 
1171
    _global_dimension = 0;
 
1172
  }
 
1173
 
 
1174
  /// Destructor
 
1175
  virtual ~heat_dof_map_1()
 
1176
  {
 
1177
    // Do nothing
 
1178
  }
 
1179
 
 
1180
  /// Return a string identifying the dof map
 
1181
  virtual const char* signature() const
 
1182
  {
 
1183
    return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
 
1184
  }
 
1185
 
 
1186
  /// Return true iff mesh entities of topological dimension d are needed
 
1187
  virtual bool needs_mesh_entities(unsigned int d) const
 
1188
  {
 
1189
    switch (d)
 
1190
    {
 
1191
    case 0:
 
1192
      {
 
1193
        return true;
 
1194
        break;
 
1195
      }
 
1196
    case 1:
 
1197
      {
 
1198
        return false;
 
1199
        break;
 
1200
      }
 
1201
    case 2:
 
1202
      {
 
1203
        return false;
 
1204
        break;
 
1205
      }
 
1206
    }
 
1207
    
 
1208
    return false;
 
1209
  }
 
1210
 
 
1211
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1212
  virtual bool init_mesh(const ufc::mesh& m)
 
1213
  {
 
1214
    _global_dimension = m.num_entities[0];
 
1215
    return false;
 
1216
  }
 
1217
 
 
1218
  /// Initialize dof map for given cell
 
1219
  virtual void init_cell(const ufc::mesh& m,
 
1220
                         const ufc::cell& c)
 
1221
  {
 
1222
    // Do nothing
 
1223
  }
 
1224
 
 
1225
  /// Finish initialization of dof map for cells
 
1226
  virtual void init_cell_finalize()
 
1227
  {
 
1228
    // Do nothing
 
1229
  }
 
1230
 
 
1231
  /// Return the dimension of the global finite element function space
 
1232
  virtual unsigned int global_dimension() const
 
1233
  {
 
1234
    return _global_dimension;
 
1235
  }
 
1236
 
 
1237
  /// Return the dimension of the local finite element function space for a cell
 
1238
  virtual unsigned int local_dimension(const ufc::cell& c) const
 
1239
  {
 
1240
    return 3;
 
1241
  }
 
1242
 
 
1243
  /// Return the maximum dimension of the local finite element function space
 
1244
  virtual unsigned int max_local_dimension() const
 
1245
  {
 
1246
    return 3;
 
1247
  }
 
1248
 
 
1249
  // Return the geometric dimension of the coordinates this dof map provides
 
1250
  virtual unsigned int geometric_dimension() const
 
1251
  {
 
1252
    return 2;
 
1253
  }
 
1254
 
 
1255
  /// Return the number of dofs on each cell facet
 
1256
  virtual unsigned int num_facet_dofs() const
 
1257
  {
 
1258
    return 2;
 
1259
  }
 
1260
 
 
1261
  /// Return the number of dofs associated with each cell entity of dimension d
 
1262
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1263
  {
 
1264
    switch (d)
 
1265
    {
 
1266
    case 0:
 
1267
      {
 
1268
        return 1;
 
1269
        break;
 
1270
      }
 
1271
    case 1:
 
1272
      {
 
1273
        return 0;
 
1274
        break;
 
1275
      }
 
1276
    case 2:
 
1277
      {
 
1278
        return 0;
 
1279
        break;
 
1280
      }
 
1281
    }
 
1282
    
 
1283
    return 0;
 
1284
  }
 
1285
 
 
1286
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1287
  virtual void tabulate_dofs(unsigned int* dofs,
 
1288
                             const ufc::mesh& m,
 
1289
                             const ufc::cell& c) const
 
1290
  {
 
1291
    dofs[0] = c.entity_indices[0][0];
 
1292
    dofs[1] = c.entity_indices[0][1];
 
1293
    dofs[2] = c.entity_indices[0][2];
 
1294
  }
 
1295
 
 
1296
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1297
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1298
                                   unsigned int facet) const
 
1299
  {
 
1300
    switch (facet)
 
1301
    {
 
1302
    case 0:
 
1303
      {
 
1304
        dofs[0] = 1;
 
1305
      dofs[1] = 2;
 
1306
        break;
 
1307
      }
 
1308
    case 1:
 
1309
      {
 
1310
        dofs[0] = 0;
 
1311
      dofs[1] = 2;
 
1312
        break;
 
1313
      }
 
1314
    case 2:
 
1315
      {
 
1316
        dofs[0] = 0;
 
1317
      dofs[1] = 1;
 
1318
        break;
 
1319
      }
 
1320
    }
 
1321
    
 
1322
  }
 
1323
 
 
1324
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1325
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1326
                                    unsigned int d, unsigned int i) const
 
1327
  {
 
1328
    if (d > 2)
 
1329
    {
 
1330
    std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
 
1331
    }
 
1332
    
 
1333
    switch (d)
 
1334
    {
 
1335
    case 0:
 
1336
      {
 
1337
        if (i > 2)
 
1338
      {
 
1339
      std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
 
1340
      }
 
1341
      
 
1342
      switch (i)
 
1343
      {
 
1344
      case 0:
 
1345
        {
 
1346
          dofs[0] = 0;
 
1347
          break;
 
1348
        }
 
1349
      case 1:
 
1350
        {
 
1351
          dofs[0] = 1;
 
1352
          break;
 
1353
        }
 
1354
      case 2:
 
1355
        {
 
1356
          dofs[0] = 2;
 
1357
          break;
 
1358
        }
 
1359
      }
 
1360
      
 
1361
        break;
 
1362
      }
 
1363
    case 1:
 
1364
      {
 
1365
        
 
1366
        break;
 
1367
      }
 
1368
    case 2:
 
1369
      {
 
1370
        
 
1371
        break;
 
1372
      }
 
1373
    }
 
1374
    
 
1375
  }
 
1376
 
 
1377
  /// Tabulate the coordinates of all dofs on a cell
 
1378
  virtual void tabulate_coordinates(double** coordinates,
 
1379
                                    const ufc::cell& c) const
 
1380
  {
 
1381
    const double * const * x = c.coordinates;
 
1382
    
 
1383
    coordinates[0][0] = x[0][0];
 
1384
    coordinates[0][1] = x[0][1];
 
1385
    coordinates[1][0] = x[1][0];
 
1386
    coordinates[1][1] = x[1][1];
 
1387
    coordinates[2][0] = x[2][0];
 
1388
    coordinates[2][1] = x[2][1];
 
1389
  }
 
1390
 
 
1391
  /// Return the number of sub dof maps (for a mixed element)
 
1392
  virtual unsigned int num_sub_dof_maps() const
 
1393
  {
 
1394
    return 0;
 
1395
  }
 
1396
 
 
1397
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1398
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1399
  {
 
1400
    return 0;
 
1401
  }
 
1402
 
 
1403
};
 
1404
 
 
1405
/// This class defines the interface for the tabulation of the cell
 
1406
/// tensor corresponding to the local contribution to a form from
 
1407
/// the integral over a cell.
 
1408
 
 
1409
class heat_cell_integral_0_0: public ufc::cell_integral
 
1410
{
 
1411
public:
 
1412
 
 
1413
  /// Constructor
 
1414
  heat_cell_integral_0_0() : ufc::cell_integral()
 
1415
  {
 
1416
    // Do nothing
 
1417
  }
 
1418
 
 
1419
  /// Destructor
 
1420
  virtual ~heat_cell_integral_0_0()
 
1421
  {
 
1422
    // Do nothing
 
1423
  }
 
1424
 
 
1425
  /// Tabulate the tensor for the contribution from a local cell
 
1426
  virtual void tabulate_tensor(double* A,
 
1427
                               const double * const * w,
 
1428
                               const ufc::cell& c) const
 
1429
  {
 
1430
    // Extract vertex coordinates
 
1431
    const double * const * x = c.coordinates;
 
1432
    
 
1433
    // Compute Jacobian of affine map from reference cell
 
1434
    const double J_00 = x[1][0] - x[0][0];
 
1435
    const double J_01 = x[2][0] - x[0][0];
 
1436
    const double J_10 = x[1][1] - x[0][1];
 
1437
    const double J_11 = x[2][1] - x[0][1];
 
1438
    
 
1439
    // Compute determinant of Jacobian
 
1440
    double detJ = J_00*J_11 - J_01*J_10;
 
1441
    
 
1442
    // Compute inverse of Jacobian
 
1443
    const double K_00 =  J_11 / detJ;
 
1444
    const double K_01 = -J_01 / detJ;
 
1445
    const double K_10 = -J_10 / detJ;
 
1446
    const double K_11 =  J_00 / detJ;
 
1447
    
 
1448
    // Set scale factor
 
1449
    const double det = std::abs(detJ);
 
1450
    
 
1451
    // Array of quadrature weights
 
1452
    static const double W4[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
 
1453
    // Quadrature points on the UFC reference element: (0.17855873, 0.15505103), (0.07503111, 0.64494897), (0.66639025, 0.15505103), (0.28001992, 0.64494897)
 
1454
    
 
1455
    // Value of basis functions at quadrature points.
 
1456
    static const double FE0[4][3] = \
 
1457
    {{0.66639025, 0.17855873, 0.15505103},
 
1458
    {0.28001992, 0.07503111, 0.64494897},
 
1459
    {0.17855873, 0.66639025, 0.15505103},
 
1460
    {0.07503111, 0.28001992, 0.64494897}};
 
1461
    
 
1462
    static const double FE0_D01[4][3] = \
 
1463
    {{-1.00000000, 0.00000000, 1.00000000},
 
1464
    {-1.00000000, 0.00000000, 1.00000000},
 
1465
    {-1.00000000, 0.00000000, 1.00000000},
 
1466
    {-1.00000000, 0.00000000, 1.00000000}};
 
1467
    
 
1468
    static const double FE0_D10[4][3] = \
 
1469
    {{-1.00000000, 1.00000000, 0.00000000},
 
1470
    {-1.00000000, 1.00000000, 0.00000000},
 
1471
    {-1.00000000, 1.00000000, 0.00000000},
 
1472
    {-1.00000000, 1.00000000, 0.00000000}};
 
1473
    
 
1474
    for (unsigned int r = 0; r < 9; r++)
 
1475
    {
 
1476
      A[r] = 0.00000000;
 
1477
    }// end loop over 'r'
 
1478
    
 
1479
    // Compute element tensor using UFL quadrature representation
 
1480
    // Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
 
1481
    
 
1482
    // Loop quadrature points for integral
 
1483
    // Number of operations to compute element tensor for following IP loop = 816
 
1484
    for (unsigned int ip = 0; ip < 4; ip++)
 
1485
    {
 
1486
      
 
1487
      // Coefficient declarations
 
1488
      double F0 = 0.00000000;
 
1489
      
 
1490
      // Total number of operations to compute function values = 6
 
1491
      for (unsigned int r = 0; r < 3; r++)
 
1492
      {
 
1493
        F0 += FE0[ip][r]*w[0][r];
 
1494
      }// end loop over 'r'
 
1495
      
 
1496
      // Number of operations for primary indices: 198
 
1497
      for (unsigned int j = 0; j < 3; j++)
 
1498
      {
 
1499
        for (unsigned int k = 0; k < 3; k++)
 
1500
        {
 
1501
          // Number of operations to compute entry: 22
 
1502
          A[j*3 + k] += ((FE0[ip][j]*FE0[ip][k] + ((((K_00*FE0_D10[ip][j] + K_10*FE0_D01[ip][j]))*((K_00*FE0_D10[ip][k] + K_10*FE0_D01[ip][k])) + ((K_01*FE0_D10[ip][j] + K_11*FE0_D01[ip][j]))*((K_01*FE0_D10[ip][k] + K_11*FE0_D01[ip][k]))))*F0*w[1][0]))*W4[ip]*det;
 
1503
        }// end loop over 'k'
 
1504
      }// end loop over 'j'
 
1505
    }// end loop over 'ip'
 
1506
  }
 
1507
 
 
1508
};
 
1509
 
 
1510
/// This class defines the interface for the tabulation of the cell
 
1511
/// tensor corresponding to the local contribution to a form from
 
1512
/// the integral over a cell.
 
1513
 
 
1514
class heat_cell_integral_1_0: public ufc::cell_integral
 
1515
{
 
1516
public:
 
1517
 
 
1518
  /// Constructor
 
1519
  heat_cell_integral_1_0() : ufc::cell_integral()
 
1520
  {
 
1521
    // Do nothing
 
1522
  }
 
1523
 
 
1524
  /// Destructor
 
1525
  virtual ~heat_cell_integral_1_0()
 
1526
  {
 
1527
    // Do nothing
 
1528
  }
 
1529
 
 
1530
  /// Tabulate the tensor for the contribution from a local cell
 
1531
  virtual void tabulate_tensor(double* A,
 
1532
                               const double * const * w,
 
1533
                               const ufc::cell& c) const
 
1534
  {
 
1535
    // Number of operations (multiply-add pairs) for Jacobian data:      12
 
1536
    // Number of operations (multiply-add pairs) for geometry tensor:    7
 
1537
    // Number of operations (multiply-add pairs) for tensor contraction: 16
 
1538
    // Total number of operations (multiply-add pairs):                  35
 
1539
    
 
1540
    // Extract vertex coordinates
 
1541
    const double * const * x = c.coordinates;
 
1542
    
 
1543
    // Compute Jacobian of affine map from reference cell
 
1544
    const double J_00 = x[1][0] - x[0][0];
 
1545
    const double J_01 = x[2][0] - x[0][0];
 
1546
    const double J_10 = x[1][1] - x[0][1];
 
1547
    const double J_11 = x[2][1] - x[0][1];
 
1548
    
 
1549
    // Compute determinant of Jacobian
 
1550
    double detJ = J_00*J_11 - J_01*J_10;
 
1551
    
 
1552
    // Compute inverse of Jacobian
 
1553
    
 
1554
    // Set scale factor
 
1555
    const double det = std::abs(detJ);
 
1556
    
 
1557
    // Compute geometry tensor
 
1558
    const double G0_0 = det*w[0][0]*(1.0);
 
1559
    const double G0_1 = det*w[0][1]*(1.0);
 
1560
    const double G0_2 = det*w[0][2]*(1.0);
 
1561
    const double G1_0_0 = det*w[1][0]*w[2][0]*(1.0);
 
1562
    const double G1_1_0 = det*w[1][1]*w[2][0]*(1.0);
 
1563
    const double G1_2_0 = det*w[1][2]*w[2][0]*(1.0);
 
1564
    
 
1565
    // Compute element tensor
 
1566
    A[0] = 0.08333333*G0_0 + 0.04166667*G0_1 + 0.04166667*G0_2 + 0.08333333*G1_0_0 + 0.04166667*G1_1_0 + 0.04166667*G1_2_0;
 
1567
    A[1] = 0.04166667*G0_0 + 0.08333333*G0_1 + 0.04166667*G0_2 + 0.04166667*G1_0_0 + 0.08333333*G1_1_0 + 0.04166667*G1_2_0;
 
1568
    A[2] = 0.04166667*G0_0 + 0.04166667*G0_1 + 0.08333333*G0_2 + 0.04166667*G1_0_0 + 0.04166667*G1_1_0 + 0.08333333*G1_2_0;
 
1569
  }
 
1570
 
 
1571
};
 
1572
 
 
1573
/// This class defines the interface for the assembly of the global
 
1574
/// tensor corresponding to a form with r + n arguments, that is, a
 
1575
/// mapping
 
1576
///
 
1577
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
1578
///
 
1579
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
1580
/// global tensor A is defined by
 
1581
///
 
1582
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
1583
///
 
1584
/// where each argument Vj represents the application to the
 
1585
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
1586
/// fixed functions (coefficients).
 
1587
 
 
1588
class heat_form_0: public ufc::form
 
1589
{
 
1590
public:
 
1591
 
 
1592
  /// Constructor
 
1593
  heat_form_0() : ufc::form()
 
1594
  {
 
1595
    // Do nothing
 
1596
  }
 
1597
 
 
1598
  /// Destructor
 
1599
  virtual ~heat_form_0()
 
1600
  {
 
1601
    // Do nothing
 
1602
  }
 
1603
 
 
1604
  /// Return a string identifying the form
 
1605
  virtual const char* signature() const
 
1606
  {
 
1607
    return "Form([Integral(Sum(Product(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1)), Product(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 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', 1, Space(2)), 1), 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', 1, Space(2)), 1), 0), Constant(Cell('triangle', 1, Space(2)), 1)))), Measure('cell', 0, None))])";
 
1608
  }
 
1609
 
 
1610
  /// Return the rank of the global tensor (r)
 
1611
  virtual unsigned int rank() const
 
1612
  {
 
1613
    return 2;
 
1614
  }
 
1615
 
 
1616
  /// Return the number of coefficients (n)
 
1617
  virtual unsigned int num_coefficients() const
 
1618
  {
 
1619
    return 2;
 
1620
  }
 
1621
 
 
1622
  /// Return the number of cell integrals
 
1623
  virtual unsigned int num_cell_integrals() const
 
1624
  {
 
1625
    return 1;
 
1626
  }
 
1627
 
 
1628
  /// Return the number of exterior facet integrals
 
1629
  virtual unsigned int num_exterior_facet_integrals() const
 
1630
  {
 
1631
    return 0;
 
1632
  }
 
1633
 
 
1634
  /// Return the number of interior facet integrals
 
1635
  virtual unsigned int num_interior_facet_integrals() const
 
1636
  {
 
1637
    return 0;
 
1638
  }
 
1639
 
 
1640
  /// Create a new finite element for argument function i
 
1641
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
1642
  {
 
1643
    switch (i)
 
1644
    {
 
1645
    case 0:
 
1646
      {
 
1647
        return new heat_finite_element_1();
 
1648
        break;
 
1649
      }
 
1650
    case 1:
 
1651
      {
 
1652
        return new heat_finite_element_1();
 
1653
        break;
 
1654
      }
 
1655
    case 2:
 
1656
      {
 
1657
        return new heat_finite_element_1();
 
1658
        break;
 
1659
      }
 
1660
    case 3:
 
1661
      {
 
1662
        return new heat_finite_element_0();
 
1663
        break;
 
1664
      }
 
1665
    }
 
1666
    
 
1667
    return 0;
 
1668
  }
 
1669
 
 
1670
  /// Create a new dof map for argument function i
 
1671
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
1672
  {
 
1673
    switch (i)
 
1674
    {
 
1675
    case 0:
 
1676
      {
 
1677
        return new heat_dof_map_1();
 
1678
        break;
 
1679
      }
 
1680
    case 1:
 
1681
      {
 
1682
        return new heat_dof_map_1();
 
1683
        break;
 
1684
      }
 
1685
    case 2:
 
1686
      {
 
1687
        return new heat_dof_map_1();
 
1688
        break;
 
1689
      }
 
1690
    case 3:
 
1691
      {
 
1692
        return new heat_dof_map_0();
 
1693
        break;
 
1694
      }
 
1695
    }
 
1696
    
 
1697
    return 0;
 
1698
  }
 
1699
 
 
1700
  /// Create a new cell integral on sub domain i
 
1701
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
1702
  {
 
1703
    switch (i)
 
1704
    {
 
1705
    case 0:
 
1706
      {
 
1707
        return new heat_cell_integral_0_0();
 
1708
        break;
 
1709
      }
 
1710
    }
 
1711
    
 
1712
    return 0;
 
1713
  }
 
1714
 
 
1715
  /// Create a new exterior facet integral on sub domain i
 
1716
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
1717
  {
 
1718
    return 0;
 
1719
  }
 
1720
 
 
1721
  /// Create a new interior facet integral on sub domain i
 
1722
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
1723
  {
 
1724
    return 0;
 
1725
  }
 
1726
 
 
1727
};
 
1728
 
 
1729
/// This class defines the interface for the assembly of the global
 
1730
/// tensor corresponding to a form with r + n arguments, that is, a
 
1731
/// mapping
 
1732
///
 
1733
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
1734
///
 
1735
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
1736
/// global tensor A is defined by
 
1737
///
 
1738
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
1739
///
 
1740
/// where each argument Vj represents the application to the
 
1741
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
1742
/// fixed functions (coefficients).
 
1743
 
 
1744
class heat_form_1: public ufc::form
 
1745
{
 
1746
public:
 
1747
 
 
1748
  /// Constructor
 
1749
  heat_form_1() : ufc::form()
 
1750
  {
 
1751
    // Do nothing
 
1752
  }
 
1753
 
 
1754
  /// Destructor
 
1755
  virtual ~heat_form_1()
 
1756
  {
 
1757
    // Do nothing
 
1758
  }
 
1759
 
 
1760
  /// Return a string identifying the form
 
1761
  virtual const char* signature() const
 
1762
  {
 
1763
    return "Form([Integral(Sum(Product(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0)), Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1), Product(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Constant(Cell('triangle', 1, Space(2)), 2)))), Measure('cell', 0, None))])";
 
1764
  }
 
1765
 
 
1766
  /// Return the rank of the global tensor (r)
 
1767
  virtual unsigned int rank() const
 
1768
  {
 
1769
    return 1;
 
1770
  }
 
1771
 
 
1772
  /// Return the number of coefficients (n)
 
1773
  virtual unsigned int num_coefficients() const
 
1774
  {
 
1775
    return 3;
 
1776
  }
 
1777
 
 
1778
  /// Return the number of cell integrals
 
1779
  virtual unsigned int num_cell_integrals() const
 
1780
  {
 
1781
    return 1;
 
1782
  }
 
1783
 
 
1784
  /// Return the number of exterior facet integrals
 
1785
  virtual unsigned int num_exterior_facet_integrals() const
 
1786
  {
 
1787
    return 0;
 
1788
  }
 
1789
 
 
1790
  /// Return the number of interior facet integrals
 
1791
  virtual unsigned int num_interior_facet_integrals() const
 
1792
  {
 
1793
    return 0;
 
1794
  }
 
1795
 
 
1796
  /// Create a new finite element for argument function i
 
1797
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
1798
  {
 
1799
    switch (i)
 
1800
    {
 
1801
    case 0:
 
1802
      {
 
1803
        return new heat_finite_element_1();
 
1804
        break;
 
1805
      }
 
1806
    case 1:
 
1807
      {
 
1808
        return new heat_finite_element_1();
 
1809
        break;
 
1810
      }
 
1811
    case 2:
 
1812
      {
 
1813
        return new heat_finite_element_1();
 
1814
        break;
 
1815
      }
 
1816
    case 3:
 
1817
      {
 
1818
        return new heat_finite_element_0();
 
1819
        break;
 
1820
      }
 
1821
    }
 
1822
    
 
1823
    return 0;
 
1824
  }
 
1825
 
 
1826
  /// Create a new dof map for argument function i
 
1827
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
1828
  {
 
1829
    switch (i)
 
1830
    {
 
1831
    case 0:
 
1832
      {
 
1833
        return new heat_dof_map_1();
 
1834
        break;
 
1835
      }
 
1836
    case 1:
 
1837
      {
 
1838
        return new heat_dof_map_1();
 
1839
        break;
 
1840
      }
 
1841
    case 2:
 
1842
      {
 
1843
        return new heat_dof_map_1();
 
1844
        break;
 
1845
      }
 
1846
    case 3:
 
1847
      {
 
1848
        return new heat_dof_map_0();
 
1849
        break;
 
1850
      }
 
1851
    }
 
1852
    
 
1853
    return 0;
 
1854
  }
 
1855
 
 
1856
  /// Create a new cell integral on sub domain i
 
1857
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
1858
  {
 
1859
    switch (i)
 
1860
    {
 
1861
    case 0:
 
1862
      {
 
1863
        return new heat_cell_integral_1_0();
 
1864
        break;
 
1865
      }
 
1866
    }
 
1867
    
 
1868
    return 0;
 
1869
  }
 
1870
 
 
1871
  /// Create a new exterior facet integral on sub domain i
 
1872
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
1873
  {
 
1874
    return 0;
 
1875
  }
 
1876
 
 
1877
  /// Create a new interior facet integral on sub domain i
 
1878
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
1879
  {
 
1880
    return 0;
 
1881
  }
 
1882
 
 
1883
};
 
1884
 
 
1885
#endif