← Back to branch summary

~ubuntu-branches/ubuntu/utopic/ffc/utopic

« back to all changes in this revision

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

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2013-06-26 14:48:32 UTC
  • mfrom: (1.1.13)
  • Revision ID: package-import@ubuntu.com-20130626144832-1xd8htax4s3utybz
Tags: 1.2.0-1
* New upstream release.
* debian/control:
  - Bump required version for python-ufc, python-fiat, python-instant
    and python-ufl in Depends field.
  - Bump Standards-Version to 3.9.4.
  - Remove DM-Upload-Allowed field.
  - Bump required debhelper version in Build-Depends.
  - Remove cdbs from Build-Depends.
  - Use canonical URIs for Vcs-* fields.
* debian/compat: Bump to compatibility level 9.
* debian/rules: Rewrite for debhelper (drop cdbs).

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/regression/references/r_quadrature_O/P5tri.h
1
 
// This code conforms with the UFC specification version 2.0.5
2
 
// and was automatically generated by FFC version 1.0.0.
 
1
// This code conforms with the UFC specification version 2.2.0
 
2
// and was automatically generated by FFC version 1.2.0.
3
3
// 
4
4
// This code was generated with the following parameters:
5
5
// 
52
52
  /// Return a string identifying the finite element
53
53
  virtual const char* signature() const
54
54
  {
55
 
    return "FiniteElement('Lagrange', Cell('triangle', Space(2)), 5, None)";
 
55
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 5, None)";
56
56
  }
57
57
 
58
58
  /// Return the cell shape
62
62
  }
63
63
 
64
64
  /// Return the topological dimension of the cell shape
65
 
  virtual unsigned int topological_dimension() const
 
65
  virtual std::size_t topological_dimension() const
66
66
  {
67
67
    return 2;
68
68
  }
69
69
 
70
70
  /// Return the geometric dimension of the cell shape
71
 
  virtual unsigned int geometric_dimension() const
 
71
  virtual std::size_t geometric_dimension() const
72
72
  {
73
73
    return 2;
74
74
  }
75
75
 
76
76
  /// Return the dimension of the finite element function space
77
 
  virtual unsigned int space_dimension() const
 
77
  virtual std::size_t space_dimension() const
78
78
  {
79
79
    return 21;
80
80
  }
81
81
 
82
82
  /// Return the rank of the value space
83
 
  virtual unsigned int value_rank() const
 
83
  virtual std::size_t value_rank() const
84
84
  {
85
85
    return 0;
86
86
  }
87
87
 
88
88
  /// Return the dimension of the value space for axis i
89
 
  virtual unsigned int value_dimension(unsigned int i) const
 
89
  virtual std::size_t value_dimension(std::size_t i) const
90
90
  {
91
91
    return 1;
92
92
  }
93
93
 
94
 
  /// Evaluate basis function i at given point in cell
95
 
  virtual void evaluate_basis(unsigned int i,
 
94
  /// Evaluate basis function i at given point x in cell
 
95
  virtual void evaluate_basis(std::size_t i,
96
96
                              double* values,
97
 
                              const double* coordinates,
98
 
                              const ufc::cell& c) const
 
97
                              const double* x,
 
98
                              const double* vertex_coordinates,
 
99
                              int cell_orientation) const
99
100
  {
100
 
    // Extract vertex coordinates
101
 
    const double * const * x = c.coordinates;
102
 
    
103
 
    // Compute Jacobian of affine map from reference cell
104
 
    const double J_00 = x[1][0] - x[0][0];
105
 
    const double J_01 = x[2][0] - x[0][0];
106
 
    const double J_10 = x[1][1] - x[0][1];
107
 
    const double J_11 = x[2][1] - x[0][1];
108
 
    
109
 
    // Compute determinant of Jacobian
110
 
    double detJ = J_00*J_11 - J_01*J_10;
111
 
    
112
 
    // Compute inverse of Jacobian
 
101
    // Compute Jacobian
 
102
    double J[4];
 
103
    compute_jacobian_triangle_2d(J, vertex_coordinates);
 
104
    
 
105
    // Compute Jacobian inverse and determinant
 
106
    double K[4];
 
107
    double detJ;
 
108
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
 
109
    
113
110
    
114
111
    // Compute constants
115
 
    const double C0 = x[1][0] + x[2][0];
116
 
    const double C1 = x[1][1] + x[2][1];
 
112
    const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
 
113
    const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
117
114
    
118
115
    // Get coordinates and map to the reference (FIAT) element
119
 
    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
120
 
    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
 
116
    double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
 
117
    double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
121
118
    
122
 
    // Reset values.
 
119
    // Reset values
123
120
    *values = 0.0;
124
121
    switch (i)
125
122
    {
126
123
    case 0:
127
124
      {
128
125
        
129
 
      // Array of basisvalues.
 
126
      // Array of basisvalues
130
127
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
131
128
      
132
 
      // Declare helper variables.
 
129
      // Declare helper variables
133
130
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
134
131
      double tmp1 = (1.0 - Y)/2.0;
135
132
      double tmp2 = tmp1*tmp1;
136
133
      
137
 
      // Compute basisvalues.
 
134
      // Compute basisvalues
138
135
      basisvalues[0] = 1.0;
139
136
      basisvalues[1] = tmp0;
140
137
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
178
175
      basisvalues[16] *= std::sqrt(27.0);
179
176
      basisvalues[15] *= std::sqrt(33.0);
180
177
      
181
 
      // Table(s) of coefficients.
 
178
      // Table(s) of coefficients
182
179
      static const double coefficients0[21] = \
183
180
      {0.015432886, -0.009450674, -0.0054563492, 0.030187531, 0.023383161, 0.013500274, -0.023199761, -0.019607376, -0.015187808, -0.0087686853, 0.021389865, 0.018864088, 0.015943064, 0.012349444, 0.007129955, -0.017989176, -0.016271822, -0.014350398, -0.0121283, -0.0093945406, -0.0054239406};
184
181
      
185
 
      // Compute value(s).
 
182
      // Compute value(s)
186
183
      for (unsigned int r = 0; r < 21; r++)
187
184
      {
188
185
        *values += coefficients0[r]*basisvalues[r];
192
189
    case 1:
193
190
      {
194
191
        
195
 
      // Array of basisvalues.
 
192
      // Array of basisvalues
196
193
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
197
194
      
198
 
      // Declare helper variables.
 
195
      // Declare helper variables
199
196
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
200
197
      double tmp1 = (1.0 - Y)/2.0;
201
198
      double tmp2 = tmp1*tmp1;
202
199
      
203
 
      // Compute basisvalues.
 
200
      // Compute basisvalues
204
201
      basisvalues[0] = 1.0;
205
202
      basisvalues[1] = tmp0;
206
203
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
244
241
      basisvalues[16] *= std::sqrt(27.0);
245
242
      basisvalues[15] *= std::sqrt(33.0);
246
243
      
247
 
      // Table(s) of coefficients.
 
244
      // Table(s) of coefficients
248
245
      static const double coefficients0[21] = \
249
246
      {0.015432886, 0.009450674, -0.0054563492, 0.030187531, -0.023383161, 0.013500274, 0.023199761, -0.019607376, 0.015187808, -0.0087686853, 0.021389865, -0.018864088, 0.015943064, -0.012349444, 0.007129955, 0.017989176, -0.016271822, 0.014350398, -0.0121283, 0.0093945406, -0.0054239406};
250
247
      
251
 
      // Compute value(s).
 
248
      // Compute value(s)
252
249
      for (unsigned int r = 0; r < 21; r++)
253
250
      {
254
251
        *values += coefficients0[r]*basisvalues[r];
258
255
    case 2:
259
256
      {
260
257
        
261
 
      // Array of basisvalues.
 
258
      // Array of basisvalues
262
259
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
263
260
      
264
 
      // Declare helper variables.
 
261
      // Declare helper variables
265
262
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
266
263
      double tmp1 = (1.0 - Y)/2.0;
267
264
      double tmp2 = tmp1*tmp1;
268
265
      
269
 
      // Compute basisvalues.
 
266
      // Compute basisvalues
270
267
      basisvalues[0] = 1.0;
271
268
      basisvalues[1] = tmp0;
272
269
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
310
307
      basisvalues[16] *= std::sqrt(27.0);
311
308
      basisvalues[15] *= std::sqrt(33.0);
312
309
      
313
 
      // Table(s) of coefficients.
 
310
      // Table(s) of coefficients
314
311
      static const double coefficients0[21] = \
315
312
      {0.015432886, 0.0, 0.010912698, 0.0, 0.0, 0.040500822, 0.0, 0.0, 0.0, 0.035074741, 0.0, 0.0, 0.0, 0.0, 0.035649775, 0.0, 0.0, 0.0, 0.0, 0.0, 0.032543643};
316
313
      
317
 
      // Compute value(s).
 
314
      // Compute value(s)
318
315
      for (unsigned int r = 0; r < 21; r++)
319
316
      {
320
317
        *values += coefficients0[r]*basisvalues[r];
324
321
    case 3:
325
322
      {
326
323
        
327
 
      // Array of basisvalues.
 
324
      // Array of basisvalues
328
325
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
329
326
      
330
 
      // Declare helper variables.
 
327
      // Declare helper variables
331
328
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
332
329
      double tmp1 = (1.0 - Y)/2.0;
333
330
      double tmp2 = tmp1*tmp1;
334
331
      
335
 
      // Compute basisvalues.
 
332
      // Compute basisvalues
336
333
      basisvalues[0] = 1.0;
337
334
      basisvalues[1] = tmp0;
338
335
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
376
373
      basisvalues[16] *= std::sqrt(27.0);
377
374
      basisvalues[15] *= std::sqrt(33.0);
378
375
      
379
 
      // Table(s) of coefficients.
 
376
      // Table(s) of coefficients
380
377
      static const double coefficients0[21] = \
381
378
      {0.035074741, 0.12600899, -0.084325397, 0.10188292, -0.028644372, 0.0020250411, 0.077332535, -0.026143168, 0.0, 0.0058457902, 0.035649775, 0.012576059, -0.038529072, 0.044252176, -0.029708146, 0.0, 0.032543643, -0.051661432, 0.058215839, -0.052609428, 0.032543643};
382
379
      
383
 
      // Compute value(s).
 
380
      // Compute value(s)
384
381
      for (unsigned int r = 0; r < 21; r++)
385
382
      {
386
383
        *values += coefficients0[r]*basisvalues[r];
390
387
    case 4:
391
388
      {
392
389
        
393
 
      // Array of basisvalues.
 
390
      // Array of basisvalues
394
391
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
395
392
      
396
 
      // Declare helper variables.
 
393
      // Declare helper variables
397
394
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
398
395
      double tmp1 = (1.0 - Y)/2.0;
399
396
      double tmp2 = tmp1*tmp1;
400
397
      
401
 
      // Compute basisvalues.
 
398
      // Compute basisvalues
402
399
      basisvalues[0] = 1.0;
403
400
      basisvalues[1] = tmp0;
404
401
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
442
439
      basisvalues[16] *= std::sqrt(27.0);
443
440
      basisvalues[15] *= std::sqrt(33.0);
444
441
      
445
 
      // Table(s) of coefficients.
 
442
      // Table(s) of coefficients
446
443
      static const double coefficients0[21] = \
447
444
      {0.035074741, -0.032934167, 0.12152778, -0.011320324, 0.075410693, -0.022275452, 0.0, 0.10457267, -0.10125206, 0.058457902, 0.0, 0.056592264, -0.0093001208, -0.040135694, 0.041591404, 0.0, 0.0, 0.051661432, -0.10187772, 0.11837121, -0.081359109};
448
445
      
449
 
      // Compute value(s).
 
446
      // Compute value(s)
450
447
      for (unsigned int r = 0; r < 21; r++)
451
448
      {
452
449
        *values += coefficients0[r]*basisvalues[r];
456
453
    case 5:
457
454
      {
458
455
        
459
 
      // Array of basisvalues.
 
456
      // Array of basisvalues
460
457
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
461
458
      
462
 
      // Declare helper variables.
 
459
      // Declare helper variables
463
460
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
464
461
      double tmp1 = (1.0 - Y)/2.0;
465
462
      double tmp2 = tmp1*tmp1;
466
463
      
467
 
      // Compute basisvalues.
 
464
      // Compute basisvalues
468
465
      basisvalues[0] = 1.0;
469
466
      basisvalues[1] = tmp0;
470
467
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
508
505
      basisvalues[16] *= std::sqrt(27.0);
509
506
      basisvalues[15] *= std::sqrt(33.0);
510
507
      
511
 
      // Table(s) of coefficients.
 
508
      // Table(s) of coefficients
512
509
      static const double coefficients0[21] = \
513
510
      {0.035074741, 0.088779059, -0.089285714, 0.030187531, 0.043258847, -0.059401206, 0.0, 0.026143168, 0.10125206, -0.1169158, 0.0, 0.0, 0.074400966, -0.030873611, -0.011883258, 0.0, 0.0, 0.0, 0.067918479, -0.13152357, 0.10847881};
514
511
      
515
 
      // Compute value(s).
 
512
      // Compute value(s)
516
513
      for (unsigned int r = 0; r < 21; r++)
517
514
      {
518
515
        *values += coefficients0[r]*basisvalues[r];
522
519
    case 6:
523
520
      {
524
521
        
525
 
      // Array of basisvalues.
 
522
      // Array of basisvalues
526
523
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
527
524
      
528
 
      // Declare helper variables.
 
525
      // Declare helper variables
529
526
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
530
527
      double tmp1 = (1.0 - Y)/2.0;
531
528
      double tmp2 = tmp1*tmp1;
532
529
      
533
 
      // Compute basisvalues.
 
530
      // Compute basisvalues
534
531
      basisvalues[0] = 1.0;
535
532
      basisvalues[1] = tmp0;
536
533
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
574
571
      basisvalues[16] *= std::sqrt(27.0);
575
572
      basisvalues[15] *= std::sqrt(33.0);
576
573
      
577
 
      // Table(s) of coefficients.
 
574
      // Table(s) of coefficients
578
575
      static const double coefficients0[21] = \
579
576
      {0.035074741, -0.010023442, 0.15128968, 0.0, 0.050273796, 0.093151892, 0.0, 0.0, 0.050626028, 0.064303692, 0.0, 0.0, 0.0, 0.072038426, -0.023766517, 0.0, 0.0, 0.0, 0.0, 0.065761785, -0.081359109};
580
577
      
581
 
      // Compute value(s).
 
578
      // Compute value(s)
582
579
      for (unsigned int r = 0; r < 21; r++)
583
580
      {
584
581
        *values += coefficients0[r]*basisvalues[r];
588
585
    case 7:
589
586
      {
590
587
        
591
 
      // Array of basisvalues.
 
588
      // Array of basisvalues
592
589
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
593
590
      
594
 
      // Declare helper variables.
 
591
      // Declare helper variables
595
592
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
596
593
      double tmp1 = (1.0 - Y)/2.0;
597
594
      double tmp2 = tmp1*tmp1;
598
595
      
599
 
      // Compute basisvalues.
 
596
      // Compute basisvalues
600
597
      basisvalues[0] = 1.0;
601
598
      basisvalues[1] = tmp0;
602
599
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
640
637
      basisvalues[16] *= std::sqrt(27.0);
641
638
      basisvalues[15] *= std::sqrt(33.0);
642
639
      
643
 
      // Table(s) of coefficients.
 
640
      // Table(s) of coefficients
644
641
      static const double coefficients0[21] = \
645
642
      {0.035074741, -0.12600899, -0.084325397, 0.10188292, 0.028644372, 0.0020250411, -0.077332535, -0.026143168, 0.0, 0.0058457902, 0.035649775, -0.012576059, -0.038529072, -0.044252176, -0.029708146, 0.0, 0.032543643, 0.051661432, 0.058215839, 0.052609428, 0.032543643};
646
643
      
647
 
      // Compute value(s).
 
644
      // Compute value(s)
648
645
      for (unsigned int r = 0; r < 21; r++)
649
646
      {
650
647
        *values += coefficients0[r]*basisvalues[r];
654
651
    case 8:
655
652
      {
656
653
        
657
 
      // Array of basisvalues.
 
654
      // Array of basisvalues
658
655
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
659
656
      
660
 
      // Declare helper variables.
 
657
      // Declare helper variables
661
658
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
662
659
      double tmp1 = (1.0 - Y)/2.0;
663
660
      double tmp2 = tmp1*tmp1;
664
661
      
665
 
      // Compute basisvalues.
 
662
      // Compute basisvalues
666
663
      basisvalues[0] = 1.0;
667
664
      basisvalues[1] = tmp0;
668
665
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
706
703
      basisvalues[16] *= std::sqrt(27.0);
707
704
      basisvalues[15] *= std::sqrt(33.0);
708
705
      
709
 
      // Table(s) of coefficients.
 
706
      // Table(s) of coefficients
710
707
      static const double coefficients0[21] = \
711
708
      {0.035074741, 0.032934167, 0.12152778, -0.011320324, -0.075410693, -0.022275452, 0.0, 0.10457267, 0.10125206, 0.058457902, 0.0, -0.056592264, -0.0093001208, 0.040135694, 0.041591404, 0.0, 0.0, -0.051661432, -0.10187772, -0.11837121, -0.081359109};
712
709
      
713
 
      // Compute value(s).
 
710
      // Compute value(s)
714
711
      for (unsigned int r = 0; r < 21; r++)
715
712
      {
716
713
        *values += coefficients0[r]*basisvalues[r];
720
717
    case 9:
721
718
      {
722
719
        
723
 
      // Array of basisvalues.
 
720
      // Array of basisvalues
724
721
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
725
722
      
726
 
      // Declare helper variables.
 
723
      // Declare helper variables
727
724
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
728
725
      double tmp1 = (1.0 - Y)/2.0;
729
726
      double tmp2 = tmp1*tmp1;
730
727
      
731
 
      // Compute basisvalues.
 
728
      // Compute basisvalues
732
729
      basisvalues[0] = 1.0;
733
730
      basisvalues[1] = tmp0;
734
731
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
772
769
      basisvalues[16] *= std::sqrt(27.0);
773
770
      basisvalues[15] *= std::sqrt(33.0);
774
771
      
775
 
      // Table(s) of coefficients.
 
772
      // Table(s) of coefficients
776
773
      static const double coefficients0[21] = \
777
774
      {0.035074741, -0.088779059, -0.089285714, 0.030187531, -0.043258847, -0.059401206, 0.0, 0.026143168, -0.10125206, -0.1169158, 0.0, 0.0, 0.074400966, 0.030873611, -0.011883258, 0.0, 0.0, 0.0, 0.067918479, 0.13152357, 0.10847881};
778
775
      
779
 
      // Compute value(s).
 
776
      // Compute value(s)
780
777
      for (unsigned int r = 0; r < 21; r++)
781
778
      {
782
779
        *values += coefficients0[r]*basisvalues[r];
786
783
    case 10:
787
784
      {
788
785
        
789
 
      // Array of basisvalues.
 
786
      // Array of basisvalues
790
787
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
791
788
      
792
 
      // Declare helper variables.
 
789
      // Declare helper variables
793
790
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
794
791
      double tmp1 = (1.0 - Y)/2.0;
795
792
      double tmp2 = tmp1*tmp1;
796
793
      
797
 
      // Compute basisvalues.
 
794
      // Compute basisvalues
798
795
      basisvalues[0] = 1.0;
799
796
      basisvalues[1] = tmp0;
800
797
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
838
835
      basisvalues[16] *= std::sqrt(27.0);
839
836
      basisvalues[15] *= std::sqrt(33.0);
840
837
      
841
 
      // Table(s) of coefficients.
 
838
      // Table(s) of coefficients
842
839
      static const double coefficients0[21] = \
843
840
      {0.035074741, 0.010023442, 0.15128968, 0.0, -0.050273796, 0.093151892, 0.0, 0.0, -0.050626028, 0.064303692, 0.0, 0.0, 0.0, -0.072038426, -0.023766517, 0.0, 0.0, 0.0, 0.0, -0.065761785, -0.081359109};
844
841
      
845
 
      // Compute value(s).
 
842
      // Compute value(s)
846
843
      for (unsigned int r = 0; r < 21; r++)
847
844
      {
848
845
        *values += coefficients0[r]*basisvalues[r];
852
849
    case 11:
853
850
      {
854
851
        
855
 
      // Array of basisvalues.
 
852
      // Array of basisvalues
856
853
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
857
854
      
858
 
      // Declare helper variables.
 
855
      // Declare helper variables
859
856
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
860
857
      double tmp1 = (1.0 - Y)/2.0;
861
858
      double tmp2 = tmp1*tmp1;
862
859
      
863
 
      // Compute basisvalues.
 
860
      // Compute basisvalues
864
861
      basisvalues[0] = 1.0;
865
862
      basisvalues[1] = tmp0;
866
863
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
904
901
      basisvalues[16] *= std::sqrt(27.0);
905
902
      basisvalues[15] *= std::sqrt(33.0);
906
903
      
907
 
      // Table(s) of coefficients.
 
904
      // Table(s) of coefficients
908
905
      static const double coefficients0[21] = \
909
906
      {0.035074741, -0.13603243, -0.066964286, 0.036979725, 0.078918168, 0.06007622, -0.0077332535, -0.045750545, -0.055688631, -0.037997636, -0.064169595, -0.012576059, 0.017271653, 0.02778625, 0.020201539, 0.089945879, 0.048815465, 0.020090557, 0.00242566, -0.0056367244, -0.0054239406};
910
907
      
911
 
      // Compute value(s).
 
908
      // Compute value(s)
912
909
      for (unsigned int r = 0; r < 21; r++)
913
910
      {
914
911
        *values += coefficients0[r]*basisvalues[r];
918
915
    case 12:
919
916
      {
920
917
        
921
 
      // Array of basisvalues.
 
918
      // Array of basisvalues
922
919
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
923
920
      
924
 
      // Declare helper variables.
 
921
      // Declare helper variables
925
922
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
926
923
      double tmp1 = (1.0 - Y)/2.0;
927
924
      double tmp2 = tmp1*tmp1;
928
925
      
929
 
      // Compute basisvalues.
 
926
      // Compute basisvalues
930
927
      basisvalues[0] = 1.0;
931
928
      basisvalues[1] = tmp0;
932
929
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
970
967
      basisvalues[16] *= std::sqrt(27.0);
971
968
      basisvalues[15] *= std::sqrt(33.0);
972
969
      
973
 
      // Table(s) of coefficients.
 
970
      // Table(s) of coefficients
974
971
      static const double coefficients0[21] = \
975
972
      {0.035074741, 0.12171323, -0.032242063, -0.067167256, -0.032151846, 0.027675562, 0.13919856, 0.065357921, -0.010125206, -0.029228951, 0.04277973, -0.056592264, -0.033214717, 0.0092620833, 0.020201539, -0.17989176, -0.032543643, 0.011480318, 0.0097026399, -0.0018789081, -0.0054239406};
976
973
      
977
 
      // Compute value(s).
 
974
      // Compute value(s)
978
975
      for (unsigned int r = 0; r < 21; r++)
979
976
      {
980
977
        *values += coefficients0[r]*basisvalues[r];
984
981
    case 13:
985
982
      {
986
983
        
987
 
      // Array of basisvalues.
 
984
      // Array of basisvalues
988
985
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
989
986
      
990
 
      // Declare helper variables.
 
987
      // Declare helper variables
991
988
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
992
989
      double tmp1 = (1.0 - Y)/2.0;
993
990
      double tmp2 = tmp1*tmp1;
994
991
      
995
 
      // Compute basisvalues.
 
992
      // Compute basisvalues
996
993
      basisvalues[0] = 1.0;
997
994
      basisvalues[1] = tmp0;
998
995
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1036
1033
      basisvalues[16] *= std::sqrt(27.0);
1037
1034
      basisvalues[15] *= std::sqrt(33.0);
1038
1035
      
1039
 
      // Table(s) of coefficients.
 
1036
      // Table(s) of coefficients
1040
1037
      static const double coefficients0[21] = \
1041
1038
      {0.035074741, -0.12171323, -0.032242063, -0.067167256, 0.032151846, 0.027675562, -0.13919856, 0.065357921, 0.010125206, -0.029228951, 0.04277973, 0.056592264, -0.033214717, -0.0092620833, 0.020201539, 0.17989176, -0.032543643, -0.011480318, 0.0097026399, 0.0018789081, -0.0054239406};
1042
1039
      
1043
 
      // Compute value(s).
 
1040
      // Compute value(s)
1044
1041
      for (unsigned int r = 0; r < 21; r++)
1045
1042
      {
1046
1043
        *values += coefficients0[r]*basisvalues[r];
1050
1047
    case 14:
1051
1048
      {
1052
1049
        
1053
 
      // Array of basisvalues.
 
1050
      // Array of basisvalues
1054
1051
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1055
1052
      
1056
 
      // Declare helper variables.
 
1053
      // Declare helper variables
1057
1054
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1058
1055
      double tmp1 = (1.0 - Y)/2.0;
1059
1056
      double tmp2 = tmp1*tmp1;
1060
1057
      
1061
 
      // Compute basisvalues.
 
1058
      // Compute basisvalues
1062
1059
      basisvalues[0] = 1.0;
1063
1060
      basisvalues[1] = tmp0;
1064
1061
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1102
1099
      basisvalues[16] *= std::sqrt(27.0);
1103
1100
      basisvalues[15] *= std::sqrt(33.0);
1104
1101
      
1105
 
      // Table(s) of coefficients.
 
1102
      // Table(s) of coefficients
1106
1103
      static const double coefficients0[21] = \
1107
1104
      {0.035074741, 0.13603243, -0.066964286, 0.036979725, -0.078918168, 0.06007622, 0.0077332535, -0.045750545, 0.055688631, -0.037997636, -0.064169595, 0.012576059, 0.017271653, -0.02778625, 0.020201539, -0.089945879, 0.048815465, -0.020090557, 0.00242566, 0.0056367244, -0.0054239406};
1108
1105
      
1109
 
      // Compute value(s).
 
1106
      // Compute value(s)
1110
1107
      for (unsigned int r = 0; r < 21; r++)
1111
1108
      {
1112
1109
        *values += coefficients0[r]*basisvalues[r];
1116
1113
    case 15:
1117
1114
      {
1118
1115
        
1119
 
      // Array of basisvalues.
 
1116
      // Array of basisvalues
1120
1117
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1121
1118
      
1122
 
      // Declare helper variables.
 
1119
      // Declare helper variables
1123
1120
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1124
1121
      double tmp1 = (1.0 - Y)/2.0;
1125
1122
      double tmp2 = tmp1*tmp1;
1126
1123
      
1127
 
      // Compute basisvalues.
 
1124
      // Compute basisvalues
1128
1125
      basisvalues[0] = 1.0;
1129
1126
      basisvalues[1] = tmp0;
1130
1127
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1168
1165
      basisvalues[16] *= std::sqrt(27.0);
1169
1166
      basisvalues[15] *= std::sqrt(33.0);
1170
1167
      
1171
 
      // Table(s) of coefficients.
 
1168
      // Table(s) of coefficients
1172
1169
      static const double coefficients0[21] = \
1173
1170
      {0.28059793, -0.21478805, -0.12400794, 0.090562592, 0.1169158, 0.0, 0.15466507, 0.078429505, 0.10125206, 0.081841063, -0.1425991, 0.025152117, -0.031886128, -0.10702852, -0.095066067, 0.0, -0.13017457, -0.10332286, -0.02910792, 0.026304714, 0.032543643};
1174
1171
      
1175
 
      // Compute value(s).
 
1172
      // Compute value(s)
1176
1173
      for (unsigned int r = 0; r < 21; r++)
1177
1174
      {
1178
1175
        *values += coefficients0[r]*basisvalues[r];
1182
1179
    case 16:
1183
1180
      {
1184
1181
        
1185
 
      // Array of basisvalues.
 
1182
      // Array of basisvalues
1186
1183
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1187
1184
      
1188
 
      // Declare helper variables.
 
1185
      // Declare helper variables
1189
1186
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1190
1187
      double tmp1 = (1.0 - Y)/2.0;
1191
1188
      double tmp2 = tmp1*tmp1;
1192
1189
      
1193
 
      // Compute basisvalues.
 
1190
      // Compute basisvalues
1194
1191
      basisvalues[0] = 1.0;
1195
1192
      basisvalues[1] = tmp0;
1196
1193
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1234
1231
      basisvalues[16] *= std::sqrt(27.0);
1235
1232
      basisvalues[15] *= std::sqrt(33.0);
1236
1233
      
1237
 
      // Table(s) of coefficients.
 
1234
      // Table(s) of coefficients
1238
1235
      static const double coefficients0[21] = \
1239
1236
      {0.035074741, 0.0, -0.12400794, -0.38489102, 0.0, 0.10125206, 0.0, -0.10457267, 0.0, 0.046766322, 0.21389865, 0.0, 0.1408304, 0.0, -0.095066067, 0.0, 0.19526186, 0.0, -0.058215839, 0.0, 0.032543643};
1240
1237
      
1241
 
      // Compute value(s).
 
1238
      // Compute value(s)
1242
1239
      for (unsigned int r = 0; r < 21; r++)
1243
1240
      {
1244
1241
        *values += coefficients0[r]*basisvalues[r];
1248
1245
    case 17:
1249
1246
      {
1250
1247
        
1251
 
      // Array of basisvalues.
 
1248
      // Array of basisvalues
1252
1249
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1253
1250
      
1254
 
      // Declare helper variables.
 
1251
      // Declare helper variables
1255
1252
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1256
1253
      double tmp1 = (1.0 - Y)/2.0;
1257
1254
      double tmp2 = tmp1*tmp1;
1258
1255
      
1259
 
      // Compute basisvalues.
 
1256
      // Compute basisvalues
1260
1257
      basisvalues[0] = 1.0;
1261
1258
      basisvalues[1] = tmp0;
1262
1259
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1300
1297
      basisvalues[16] *= std::sqrt(27.0);
1301
1298
      basisvalues[15] *= std::sqrt(33.0);
1302
1299
      
1303
 
      // Table(s) of coefficients.
 
1300
      // Table(s) of coefficients
1304
1301
      static const double coefficients0[21] = \
1305
1302
      {0.28059793, 0.21478805, -0.12400794, 0.090562592, -0.1169158, 0.0, -0.15466507, 0.078429505, -0.10125206, 0.081841063, -0.1425991, -0.025152117, -0.031886128, 0.10702852, -0.095066067, 0.0, -0.13017457, 0.10332286, -0.02910792, -0.026304714, 0.032543643};
1306
1303
      
1307
 
      // Compute value(s).
 
1304
      // Compute value(s)
1308
1305
      for (unsigned int r = 0; r < 21; r++)
1309
1306
      {
1310
1307
        *values += coefficients0[r]*basisvalues[r];
1314
1311
    case 18:
1315
1312
      {
1316
1313
        
1317
 
      // Array of basisvalues.
 
1314
      // Array of basisvalues
1318
1315
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1319
1316
      
1320
 
      // Declare helper variables.
 
1317
      // Declare helper variables
1321
1318
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1322
1319
      double tmp1 = (1.0 - Y)/2.0;
1323
1320
      double tmp2 = tmp1*tmp1;
1324
1321
      
1325
 
      // Compute basisvalues.
 
1322
      // Compute basisvalues
1326
1323
      basisvalues[0] = 1.0;
1327
1324
      basisvalues[1] = tmp0;
1328
1325
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1366
1363
      basisvalues[16] *= std::sqrt(27.0);
1367
1364
      basisvalues[15] *= std::sqrt(33.0);
1368
1365
      
1369
 
      // Table(s) of coefficients.
 
1366
      // Table(s) of coefficients
1370
1367
      static const double coefficients0[21] = \
1371
1368
      {0.035074741, -0.10739402, 0.062003968, 0.011320324, -0.30690398, -0.25313014, 0.0, -0.10457267, 0.0, 0.046766322, 0.0, 0.16977679, 0.0093001208, 0.1265818, 0.17230725, 0.0, 0.0, 0.1549843, 0.10187772, -0.039457071, -0.081359109};
1372
1369
      
1373
 
      // Compute value(s).
 
1370
      // Compute value(s)
1374
1371
      for (unsigned int r = 0; r < 21; r++)
1375
1372
      {
1376
1373
        *values += coefficients0[r]*basisvalues[r];
1380
1377
    case 19:
1381
1378
      {
1382
1379
        
1383
 
      // Array of basisvalues.
 
1380
      // Array of basisvalues
1384
1381
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1385
1382
      
1386
 
      // Declare helper variables.
 
1383
      // Declare helper variables
1387
1384
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1388
1385
      double tmp1 = (1.0 - Y)/2.0;
1389
1386
      double tmp2 = tmp1*tmp1;
1390
1387
      
1391
 
      // Compute basisvalues.
 
1388
      // Compute basisvalues
1392
1389
      basisvalues[0] = 1.0;
1393
1390
      basisvalues[1] = tmp0;
1394
1391
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1432
1429
      basisvalues[16] *= std::sqrt(27.0);
1433
1430
      basisvalues[15] *= std::sqrt(33.0);
1434
1431
      
1435
 
      // Table(s) of coefficients.
 
1432
      // Table(s) of coefficients
1436
1433
      static const double coefficients0[21] = \
1437
1434
      {0.035074741, 0.10739402, 0.062003968, 0.011320324, 0.30690398, -0.25313014, 0.0, -0.10457267, 0.0, 0.046766322, 0.0, -0.16977679, 0.0093001208, -0.1265818, 0.17230725, 0.0, 0.0, -0.1549843, 0.10187772, 0.039457071, -0.081359109};
1438
1435
      
1439
 
      // Compute value(s).
 
1436
      // Compute value(s)
1440
1437
      for (unsigned int r = 0; r < 21; r++)
1441
1438
      {
1442
1439
        *values += coefficients0[r]*basisvalues[r];
1446
1443
    case 20:
1447
1444
      {
1448
1445
        
1449
 
      // Array of basisvalues.
 
1446
      // Array of basisvalues
1450
1447
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1451
1448
      
1452
 
      // Declare helper variables.
 
1449
      // Declare helper variables
1453
1450
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1454
1451
      double tmp1 = (1.0 - Y)/2.0;
1455
1452
      double tmp2 = tmp1*tmp1;
1456
1453
      
1457
 
      // Compute basisvalues.
 
1454
      // Compute basisvalues
1458
1455
      basisvalues[0] = 1.0;
1459
1456
      basisvalues[1] = tmp0;
1460
1457
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1498
1495
      basisvalues[16] *= std::sqrt(27.0);
1499
1496
      basisvalues[15] *= std::sqrt(33.0);
1500
1497
      
1501
 
      // Table(s) of coefficients.
 
1498
      // Table(s) of coefficients
1502
1499
      static const double coefficients0[21] = \
1503
1500
      {0.28059793, 0.0, 0.24801587, -0.060375061, 0.0, 0.13500274, 0.0, -0.052286337, 0.0, -0.21044845, 0.0, 0.0, -0.14880193, 0.0, -0.1425991, 0.0, 0.0, 0.0, -0.13583696, 0.0, 0.10847881};
1504
1501
      
1505
 
      // Compute value(s).
 
1502
      // Compute value(s)
1506
1503
      for (unsigned int r = 0; r < 21; r++)
1507
1504
      {
1508
1505
        *values += coefficients0[r]*basisvalues[r];
1513
1510
    
1514
1511
  }
1515
1512
 
1516
 
  /// Evaluate all basis functions at given point in cell
 
1513
  /// Evaluate all basis functions at given point x in cell
1517
1514
  virtual void evaluate_basis_all(double* values,
1518
 
                                  const double* coordinates,
1519
 
                                  const ufc::cell& c) const
 
1515
                                  const double* x,
 
1516
                                  const double* vertex_coordinates,
 
1517
                                  int cell_orientation) const
1520
1518
  {
1521
1519
    // Helper variable to hold values of a single dof.
1522
1520
    double dof_values = 0.0;
1523
1521
    
1524
 
    // Loop dofs and call evaluate_basis.
 
1522
    // Loop dofs and call evaluate_basis
1525
1523
    for (unsigned int r = 0; r < 21; r++)
1526
1524
    {
1527
 
      evaluate_basis(r, &dof_values, coordinates, c);
 
1525
      evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
1528
1526
      values[r] = dof_values;
1529
1527
    }// end loop over 'r'
1530
1528
  }
1531
1529
 
1532
 
  /// Evaluate order n derivatives of basis function i at given point in cell
1533
 
  virtual void evaluate_basis_derivatives(unsigned int i,
1534
 
                                          unsigned int n,
 
1530
  /// Evaluate order n derivatives of basis function i at given point x in cell
 
1531
  virtual void evaluate_basis_derivatives(std::size_t i,
 
1532
                                          std::size_t n,
1535
1533
                                          double* values,
1536
 
                                          const double* coordinates,
1537
 
                                          const ufc::cell& c) const
 
1534
                                          const double* x,
 
1535
                                          const double* vertex_coordinates,
 
1536
                                          int cell_orientation) const
1538
1537
  {
1539
 
    // Extract vertex coordinates
1540
 
    const double * const * x = c.coordinates;
1541
 
    
1542
 
    // Compute Jacobian of affine map from reference cell
1543
 
    const double J_00 = x[1][0] - x[0][0];
1544
 
    const double J_01 = x[2][0] - x[0][0];
1545
 
    const double J_10 = x[1][1] - x[0][1];
1546
 
    const double J_11 = x[2][1] - x[0][1];
1547
 
    
1548
 
    // Compute determinant of Jacobian
1549
 
    double detJ = J_00*J_11 - J_01*J_10;
1550
 
    
1551
 
    // Compute inverse of Jacobian
1552
 
    const double K_00 =  J_11 / detJ;
1553
 
    const double K_01 = -J_01 / detJ;
1554
 
    const double K_10 = -J_10 / detJ;
1555
 
    const double K_11 =  J_00 / detJ;
 
1538
    // Compute Jacobian
 
1539
    double J[4];
 
1540
    compute_jacobian_triangle_2d(J, vertex_coordinates);
 
1541
    
 
1542
    // Compute Jacobian inverse and determinant
 
1543
    double K[4];
 
1544
    double detJ;
 
1545
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
 
1546
    
1556
1547
    
1557
1548
    // Compute constants
1558
 
    const double C0 = x[1][0] + x[2][0];
1559
 
    const double C1 = x[1][1] + x[2][1];
 
1549
    const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
 
1550
    const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1560
1551
    
1561
1552
    // Get coordinates and map to the reference (FIAT) element
1562
 
    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1563
 
    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
 
1553
    double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
 
1554
    double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1564
1555
    
1565
1556
    // Compute number of derivatives.
1566
1557
    unsigned int num_derivatives = 1;
1597
1588
    }
1598
1589
    
1599
1590
    // Compute inverse of Jacobian
1600
 
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
 
1591
    const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
1601
1592
    
1602
1593
    // Declare transformation matrix
1603
1594
    // Declare pointer to two dimensional array and initialise
1631
1622
    case 0:
1632
1623
      {
1633
1624
        
1634
 
      // Array of basisvalues.
 
1625
      // Array of basisvalues
1635
1626
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1636
1627
      
1637
 
      // Declare helper variables.
 
1628
      // Declare helper variables
1638
1629
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1639
1630
      double tmp1 = (1.0 - Y)/2.0;
1640
1631
      double tmp2 = tmp1*tmp1;
1641
1632
      
1642
 
      // Compute basisvalues.
 
1633
      // Compute basisvalues
1643
1634
      basisvalues[0] = 1.0;
1644
1635
      basisvalues[1] = tmp0;
1645
1636
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1683
1674
      basisvalues[16] *= std::sqrt(27.0);
1684
1675
      basisvalues[15] *= std::sqrt(33.0);
1685
1676
      
1686
 
      // Table(s) of coefficients.
 
1677
      // Table(s) of coefficients
1687
1678
      static const double coefficients0[21] = \
1688
1679
      {0.015432886, -0.009450674, -0.0054563492, 0.030187531, 0.023383161, 0.013500274, -0.023199761, -0.019607376, -0.015187808, -0.0087686853, 0.021389865, 0.018864088, 0.015943064, 0.012349444, 0.007129955, -0.017989176, -0.016271822, -0.014350398, -0.0121283, -0.0093945406, -0.0054239406};
1689
1680
      
1887
1878
    case 1:
1888
1879
      {
1889
1880
        
1890
 
      // Array of basisvalues.
 
1881
      // Array of basisvalues
1891
1882
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1892
1883
      
1893
 
      // Declare helper variables.
 
1884
      // Declare helper variables
1894
1885
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1895
1886
      double tmp1 = (1.0 - Y)/2.0;
1896
1887
      double tmp2 = tmp1*tmp1;
1897
1888
      
1898
 
      // Compute basisvalues.
 
1889
      // Compute basisvalues
1899
1890
      basisvalues[0] = 1.0;
1900
1891
      basisvalues[1] = tmp0;
1901
1892
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1939
1930
      basisvalues[16] *= std::sqrt(27.0);
1940
1931
      basisvalues[15] *= std::sqrt(33.0);
1941
1932
      
1942
 
      // Table(s) of coefficients.
 
1933
      // Table(s) of coefficients
1943
1934
      static const double coefficients0[21] = \
1944
1935
      {0.015432886, 0.009450674, -0.0054563492, 0.030187531, -0.023383161, 0.013500274, 0.023199761, -0.019607376, 0.015187808, -0.0087686853, 0.021389865, -0.018864088, 0.015943064, -0.012349444, 0.007129955, 0.017989176, -0.016271822, 0.014350398, -0.0121283, 0.0093945406, -0.0054239406};
1945
1936
      
2143
2134
    case 2:
2144
2135
      {
2145
2136
        
2146
 
      // Array of basisvalues.
 
2137
      // Array of basisvalues
2147
2138
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2148
2139
      
2149
 
      // Declare helper variables.
 
2140
      // Declare helper variables
2150
2141
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2151
2142
      double tmp1 = (1.0 - Y)/2.0;
2152
2143
      double tmp2 = tmp1*tmp1;
2153
2144
      
2154
 
      // Compute basisvalues.
 
2145
      // Compute basisvalues
2155
2146
      basisvalues[0] = 1.0;
2156
2147
      basisvalues[1] = tmp0;
2157
2148
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2195
2186
      basisvalues[16] *= std::sqrt(27.0);
2196
2187
      basisvalues[15] *= std::sqrt(33.0);
2197
2188
      
2198
 
      // Table(s) of coefficients.
 
2189
      // Table(s) of coefficients
2199
2190
      static const double coefficients0[21] = \
2200
2191
      {0.015432886, 0.0, 0.010912698, 0.0, 0.0, 0.040500822, 0.0, 0.0, 0.0, 0.035074741, 0.0, 0.0, 0.0, 0.0, 0.035649775, 0.0, 0.0, 0.0, 0.0, 0.0, 0.032543643};
2201
2192
      
2399
2390
    case 3:
2400
2391
      {
2401
2392
        
2402
 
      // Array of basisvalues.
 
2393
      // Array of basisvalues
2403
2394
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2404
2395
      
2405
 
      // Declare helper variables.
 
2396
      // Declare helper variables
2406
2397
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2407
2398
      double tmp1 = (1.0 - Y)/2.0;
2408
2399
      double tmp2 = tmp1*tmp1;
2409
2400
      
2410
 
      // Compute basisvalues.
 
2401
      // Compute basisvalues
2411
2402
      basisvalues[0] = 1.0;
2412
2403
      basisvalues[1] = tmp0;
2413
2404
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2451
2442
      basisvalues[16] *= std::sqrt(27.0);
2452
2443
      basisvalues[15] *= std::sqrt(33.0);
2453
2444
      
2454
 
      // Table(s) of coefficients.
 
2445
      // Table(s) of coefficients
2455
2446
      static const double coefficients0[21] = \
2456
2447
      {0.035074741, 0.12600899, -0.084325397, 0.10188292, -0.028644372, 0.0020250411, 0.077332535, -0.026143168, 0.0, 0.0058457902, 0.035649775, 0.012576059, -0.038529072, 0.044252176, -0.029708146, 0.0, 0.032543643, -0.051661432, 0.058215839, -0.052609428, 0.032543643};
2457
2448
      
2655
2646
    case 4:
2656
2647
      {
2657
2648
        
2658
 
      // Array of basisvalues.
 
2649
      // Array of basisvalues
2659
2650
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2660
2651
      
2661
 
      // Declare helper variables.
 
2652
      // Declare helper variables
2662
2653
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2663
2654
      double tmp1 = (1.0 - Y)/2.0;
2664
2655
      double tmp2 = tmp1*tmp1;
2665
2656
      
2666
 
      // Compute basisvalues.
 
2657
      // Compute basisvalues
2667
2658
      basisvalues[0] = 1.0;
2668
2659
      basisvalues[1] = tmp0;
2669
2660
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2707
2698
      basisvalues[16] *= std::sqrt(27.0);
2708
2699
      basisvalues[15] *= std::sqrt(33.0);
2709
2700
      
2710
 
      // Table(s) of coefficients.
 
2701
      // Table(s) of coefficients
2711
2702
      static const double coefficients0[21] = \
2712
2703
      {0.035074741, -0.032934167, 0.12152778, -0.011320324, 0.075410693, -0.022275452, 0.0, 0.10457267, -0.10125206, 0.058457902, 0.0, 0.056592264, -0.0093001208, -0.040135694, 0.041591404, 0.0, 0.0, 0.051661432, -0.10187772, 0.11837121, -0.081359109};
2713
2704
      
2911
2902
    case 5:
2912
2903
      {
2913
2904
        
2914
 
      // Array of basisvalues.
 
2905
      // Array of basisvalues
2915
2906
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2916
2907
      
2917
 
      // Declare helper variables.
 
2908
      // Declare helper variables
2918
2909
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2919
2910
      double tmp1 = (1.0 - Y)/2.0;
2920
2911
      double tmp2 = tmp1*tmp1;
2921
2912
      
2922
 
      // Compute basisvalues.
 
2913
      // Compute basisvalues
2923
2914
      basisvalues[0] = 1.0;
2924
2915
      basisvalues[1] = tmp0;
2925
2916
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2963
2954
      basisvalues[16] *= std::sqrt(27.0);
2964
2955
      basisvalues[15] *= std::sqrt(33.0);
2965
2956
      
2966
 
      // Table(s) of coefficients.
 
2957
      // Table(s) of coefficients
2967
2958
      static const double coefficients0[21] = \
2968
2959
      {0.035074741, 0.088779059, -0.089285714, 0.030187531, 0.043258847, -0.059401206, 0.0, 0.026143168, 0.10125206, -0.1169158, 0.0, 0.0, 0.074400966, -0.030873611, -0.011883258, 0.0, 0.0, 0.0, 0.067918479, -0.13152357, 0.10847881};
2969
2960
      
3167
3158
    case 6:
3168
3159
      {
3169
3160
        
3170
 
      // Array of basisvalues.
 
3161
      // Array of basisvalues
3171
3162
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3172
3163
      
3173
 
      // Declare helper variables.
 
3164
      // Declare helper variables
3174
3165
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3175
3166
      double tmp1 = (1.0 - Y)/2.0;
3176
3167
      double tmp2 = tmp1*tmp1;
3177
3168
      
3178
 
      // Compute basisvalues.
 
3169
      // Compute basisvalues
3179
3170
      basisvalues[0] = 1.0;
3180
3171
      basisvalues[1] = tmp0;
3181
3172
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3219
3210
      basisvalues[16] *= std::sqrt(27.0);
3220
3211
      basisvalues[15] *= std::sqrt(33.0);
3221
3212
      
3222
 
      // Table(s) of coefficients.
 
3213
      // Table(s) of coefficients
3223
3214
      static const double coefficients0[21] = \
3224
3215
      {0.035074741, -0.010023442, 0.15128968, 0.0, 0.050273796, 0.093151892, 0.0, 0.0, 0.050626028, 0.064303692, 0.0, 0.0, 0.0, 0.072038426, -0.023766517, 0.0, 0.0, 0.0, 0.0, 0.065761785, -0.081359109};
3225
3216
      
3423
3414
    case 7:
3424
3415
      {
3425
3416
        
3426
 
      // Array of basisvalues.
 
3417
      // Array of basisvalues
3427
3418
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3428
3419
      
3429
 
      // Declare helper variables.
 
3420
      // Declare helper variables
3430
3421
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3431
3422
      double tmp1 = (1.0 - Y)/2.0;
3432
3423
      double tmp2 = tmp1*tmp1;
3433
3424
      
3434
 
      // Compute basisvalues.
 
3425
      // Compute basisvalues
3435
3426
      basisvalues[0] = 1.0;
3436
3427
      basisvalues[1] = tmp0;
3437
3428
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3475
3466
      basisvalues[16] *= std::sqrt(27.0);
3476
3467
      basisvalues[15] *= std::sqrt(33.0);
3477
3468
      
3478
 
      // Table(s) of coefficients.
 
3469
      // Table(s) of coefficients
3479
3470
      static const double coefficients0[21] = \
3480
3471
      {0.035074741, -0.12600899, -0.084325397, 0.10188292, 0.028644372, 0.0020250411, -0.077332535, -0.026143168, 0.0, 0.0058457902, 0.035649775, -0.012576059, -0.038529072, -0.044252176, -0.029708146, 0.0, 0.032543643, 0.051661432, 0.058215839, 0.052609428, 0.032543643};
3481
3472
      
3679
3670
    case 8:
3680
3671
      {
3681
3672
        
3682
 
      // Array of basisvalues.
 
3673
      // Array of basisvalues
3683
3674
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3684
3675
      
3685
 
      // Declare helper variables.
 
3676
      // Declare helper variables
3686
3677
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3687
3678
      double tmp1 = (1.0 - Y)/2.0;
3688
3679
      double tmp2 = tmp1*tmp1;
3689
3680
      
3690
 
      // Compute basisvalues.
 
3681
      // Compute basisvalues
3691
3682
      basisvalues[0] = 1.0;
3692
3683
      basisvalues[1] = tmp0;
3693
3684
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3731
3722
      basisvalues[16] *= std::sqrt(27.0);
3732
3723
      basisvalues[15] *= std::sqrt(33.0);
3733
3724
      
3734
 
      // Table(s) of coefficients.
 
3725
      // Table(s) of coefficients
3735
3726
      static const double coefficients0[21] = \
3736
3727
      {0.035074741, 0.032934167, 0.12152778, -0.011320324, -0.075410693, -0.022275452, 0.0, 0.10457267, 0.10125206, 0.058457902, 0.0, -0.056592264, -0.0093001208, 0.040135694, 0.041591404, 0.0, 0.0, -0.051661432, -0.10187772, -0.11837121, -0.081359109};
3737
3728
      
3935
3926
    case 9:
3936
3927
      {
3937
3928
        
3938
 
      // Array of basisvalues.
 
3929
      // Array of basisvalues
3939
3930
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
3940
3931
      
3941
 
      // Declare helper variables.
 
3932
      // Declare helper variables
3942
3933
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3943
3934
      double tmp1 = (1.0 - Y)/2.0;
3944
3935
      double tmp2 = tmp1*tmp1;
3945
3936
      
3946
 
      // Compute basisvalues.
 
3937
      // Compute basisvalues
3947
3938
      basisvalues[0] = 1.0;
3948
3939
      basisvalues[1] = tmp0;
3949
3940
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
3987
3978
      basisvalues[16] *= std::sqrt(27.0);
3988
3979
      basisvalues[15] *= std::sqrt(33.0);
3989
3980
      
3990
 
      // Table(s) of coefficients.
 
3981
      // Table(s) of coefficients
3991
3982
      static const double coefficients0[21] = \
3992
3983
      {0.035074741, -0.088779059, -0.089285714, 0.030187531, -0.043258847, -0.059401206, 0.0, 0.026143168, -0.10125206, -0.1169158, 0.0, 0.0, 0.074400966, 0.030873611, -0.011883258, 0.0, 0.0, 0.0, 0.067918479, 0.13152357, 0.10847881};
3993
3984
      
4191
4182
    case 10:
4192
4183
      {
4193
4184
        
4194
 
      // Array of basisvalues.
 
4185
      // Array of basisvalues
4195
4186
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
4196
4187
      
4197
 
      // Declare helper variables.
 
4188
      // Declare helper variables
4198
4189
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4199
4190
      double tmp1 = (1.0 - Y)/2.0;
4200
4191
      double tmp2 = tmp1*tmp1;
4201
4192
      
4202
 
      // Compute basisvalues.
 
4193
      // Compute basisvalues
4203
4194
      basisvalues[0] = 1.0;
4204
4195
      basisvalues[1] = tmp0;
4205
4196
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
4243
4234
      basisvalues[16] *= std::sqrt(27.0);
4244
4235
      basisvalues[15] *= std::sqrt(33.0);
4245
4236
      
4246
 
      // Table(s) of coefficients.
 
4237
      // Table(s) of coefficients
4247
4238
      static const double coefficients0[21] = \
4248
4239
      {0.035074741, 0.010023442, 0.15128968, 0.0, -0.050273796, 0.093151892, 0.0, 0.0, -0.050626028, 0.064303692, 0.0, 0.0, 0.0, -0.072038426, -0.023766517, 0.0, 0.0, 0.0, 0.0, -0.065761785, -0.081359109};
4249
4240
      
4447
4438
    case 11:
4448
4439
      {
4449
4440
        
4450
 
      // Array of basisvalues.
 
4441
      // Array of basisvalues
4451
4442
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
4452
4443
      
4453
 
      // Declare helper variables.
 
4444
      // Declare helper variables
4454
4445
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4455
4446
      double tmp1 = (1.0 - Y)/2.0;
4456
4447
      double tmp2 = tmp1*tmp1;
4457
4448
      
4458
 
      // Compute basisvalues.
 
4449
      // Compute basisvalues
4459
4450
      basisvalues[0] = 1.0;
4460
4451
      basisvalues[1] = tmp0;
4461
4452
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
4499
4490
      basisvalues[16] *= std::sqrt(27.0);
4500
4491
      basisvalues[15] *= std::sqrt(33.0);
4501
4492
      
4502
 
      // Table(s) of coefficients.
 
4493
      // Table(s) of coefficients
4503
4494
      static const double coefficients0[21] = \
4504
4495
      {0.035074741, -0.13603243, -0.066964286, 0.036979725, 0.078918168, 0.06007622, -0.0077332535, -0.045750545, -0.055688631, -0.037997636, -0.064169595, -0.012576059, 0.017271653, 0.02778625, 0.020201539, 0.089945879, 0.048815465, 0.020090557, 0.00242566, -0.0056367244, -0.0054239406};
4505
4496
      
4703
4694
    case 12:
4704
4695
      {
4705
4696
        
4706
 
      // Array of basisvalues.
 
4697
      // Array of basisvalues
4707
4698
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
4708
4699
      
4709
 
      // Declare helper variables.
 
4700
      // Declare helper variables
4710
4701
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4711
4702
      double tmp1 = (1.0 - Y)/2.0;
4712
4703
      double tmp2 = tmp1*tmp1;
4713
4704
      
4714
 
      // Compute basisvalues.
 
4705
      // Compute basisvalues
4715
4706
      basisvalues[0] = 1.0;
4716
4707
      basisvalues[1] = tmp0;
4717
4708
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
4755
4746
      basisvalues[16] *= std::sqrt(27.0);
4756
4747
      basisvalues[15] *= std::sqrt(33.0);
4757
4748
      
4758
 
      // Table(s) of coefficients.
 
4749
      // Table(s) of coefficients
4759
4750
      static const double coefficients0[21] = \
4760
4751
      {0.035074741, 0.12171323, -0.032242063, -0.067167256, -0.032151846, 0.027675562, 0.13919856, 0.065357921, -0.010125206, -0.029228951, 0.04277973, -0.056592264, -0.033214717, 0.0092620833, 0.020201539, -0.17989176, -0.032543643, 0.011480318, 0.0097026399, -0.0018789081, -0.0054239406};
4761
4752
      
4959
4950
    case 13:
4960
4951
      {
4961
4952
        
4962
 
      // Array of basisvalues.
 
4953
      // Array of basisvalues
4963
4954
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
4964
4955
      
4965
 
      // Declare helper variables.
 
4956
      // Declare helper variables
4966
4957
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
4967
4958
      double tmp1 = (1.0 - Y)/2.0;
4968
4959
      double tmp2 = tmp1*tmp1;
4969
4960
      
4970
 
      // Compute basisvalues.
 
4961
      // Compute basisvalues
4971
4962
      basisvalues[0] = 1.0;
4972
4963
      basisvalues[1] = tmp0;
4973
4964
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5011
5002
      basisvalues[16] *= std::sqrt(27.0);
5012
5003
      basisvalues[15] *= std::sqrt(33.0);
5013
5004
      
5014
 
      // Table(s) of coefficients.
 
5005
      // Table(s) of coefficients
5015
5006
      static const double coefficients0[21] = \
5016
5007
      {0.035074741, -0.12171323, -0.032242063, -0.067167256, 0.032151846, 0.027675562, -0.13919856, 0.065357921, 0.010125206, -0.029228951, 0.04277973, 0.056592264, -0.033214717, -0.0092620833, 0.020201539, 0.17989176, -0.032543643, -0.011480318, 0.0097026399, 0.0018789081, -0.0054239406};
5017
5008
      
5215
5206
    case 14:
5216
5207
      {
5217
5208
        
5218
 
      // Array of basisvalues.
 
5209
      // Array of basisvalues
5219
5210
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5220
5211
      
5221
 
      // Declare helper variables.
 
5212
      // Declare helper variables
5222
5213
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5223
5214
      double tmp1 = (1.0 - Y)/2.0;
5224
5215
      double tmp2 = tmp1*tmp1;
5225
5216
      
5226
 
      // Compute basisvalues.
 
5217
      // Compute basisvalues
5227
5218
      basisvalues[0] = 1.0;
5228
5219
      basisvalues[1] = tmp0;
5229
5220
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5267
5258
      basisvalues[16] *= std::sqrt(27.0);
5268
5259
      basisvalues[15] *= std::sqrt(33.0);
5269
5260
      
5270
 
      // Table(s) of coefficients.
 
5261
      // Table(s) of coefficients
5271
5262
      static const double coefficients0[21] = \
5272
5263
      {0.035074741, 0.13603243, -0.066964286, 0.036979725, -0.078918168, 0.06007622, 0.0077332535, -0.045750545, 0.055688631, -0.037997636, -0.064169595, 0.012576059, 0.017271653, -0.02778625, 0.020201539, -0.089945879, 0.048815465, -0.020090557, 0.00242566, 0.0056367244, -0.0054239406};
5273
5264
      
5471
5462
    case 15:
5472
5463
      {
5473
5464
        
5474
 
      // Array of basisvalues.
 
5465
      // Array of basisvalues
5475
5466
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5476
5467
      
5477
 
      // Declare helper variables.
 
5468
      // Declare helper variables
5478
5469
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5479
5470
      double tmp1 = (1.0 - Y)/2.0;
5480
5471
      double tmp2 = tmp1*tmp1;
5481
5472
      
5482
 
      // Compute basisvalues.
 
5473
      // Compute basisvalues
5483
5474
      basisvalues[0] = 1.0;
5484
5475
      basisvalues[1] = tmp0;
5485
5476
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5523
5514
      basisvalues[16] *= std::sqrt(27.0);
5524
5515
      basisvalues[15] *= std::sqrt(33.0);
5525
5516
      
5526
 
      // Table(s) of coefficients.
 
5517
      // Table(s) of coefficients
5527
5518
      static const double coefficients0[21] = \
5528
5519
      {0.28059793, -0.21478805, -0.12400794, 0.090562592, 0.1169158, 0.0, 0.15466507, 0.078429505, 0.10125206, 0.081841063, -0.1425991, 0.025152117, -0.031886128, -0.10702852, -0.095066067, 0.0, -0.13017457, -0.10332286, -0.02910792, 0.026304714, 0.032543643};
5529
5520
      
5727
5718
    case 16:
5728
5719
      {
5729
5720
        
5730
 
      // Array of basisvalues.
 
5721
      // Array of basisvalues
5731
5722
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5732
5723
      
5733
 
      // Declare helper variables.
 
5724
      // Declare helper variables
5734
5725
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5735
5726
      double tmp1 = (1.0 - Y)/2.0;
5736
5727
      double tmp2 = tmp1*tmp1;
5737
5728
      
5738
 
      // Compute basisvalues.
 
5729
      // Compute basisvalues
5739
5730
      basisvalues[0] = 1.0;
5740
5731
      basisvalues[1] = tmp0;
5741
5732
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
5779
5770
      basisvalues[16] *= std::sqrt(27.0);
5780
5771
      basisvalues[15] *= std::sqrt(33.0);
5781
5772
      
5782
 
      // Table(s) of coefficients.
 
5773
      // Table(s) of coefficients
5783
5774
      static const double coefficients0[21] = \
5784
5775
      {0.035074741, 0.0, -0.12400794, -0.38489102, 0.0, 0.10125206, 0.0, -0.10457267, 0.0, 0.046766322, 0.21389865, 0.0, 0.1408304, 0.0, -0.095066067, 0.0, 0.19526186, 0.0, -0.058215839, 0.0, 0.032543643};
5785
5776
      
5983
5974
    case 17:
5984
5975
      {
5985
5976
        
5986
 
      // Array of basisvalues.
 
5977
      // Array of basisvalues
5987
5978
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
5988
5979
      
5989
 
      // Declare helper variables.
 
5980
      // Declare helper variables
5990
5981
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
5991
5982
      double tmp1 = (1.0 - Y)/2.0;
5992
5983
      double tmp2 = tmp1*tmp1;
5993
5984
      
5994
 
      // Compute basisvalues.
 
5985
      // Compute basisvalues
5995
5986
      basisvalues[0] = 1.0;
5996
5987
      basisvalues[1] = tmp0;
5997
5988
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6035
6026
      basisvalues[16] *= std::sqrt(27.0);
6036
6027
      basisvalues[15] *= std::sqrt(33.0);
6037
6028
      
6038
 
      // Table(s) of coefficients.
 
6029
      // Table(s) of coefficients
6039
6030
      static const double coefficients0[21] = \
6040
6031
      {0.28059793, 0.21478805, -0.12400794, 0.090562592, -0.1169158, 0.0, -0.15466507, 0.078429505, -0.10125206, 0.081841063, -0.1425991, -0.025152117, -0.031886128, 0.10702852, -0.095066067, 0.0, -0.13017457, 0.10332286, -0.02910792, -0.026304714, 0.032543643};
6041
6032
      
6239
6230
    case 18:
6240
6231
      {
6241
6232
        
6242
 
      // Array of basisvalues.
 
6233
      // Array of basisvalues
6243
6234
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6244
6235
      
6245
 
      // Declare helper variables.
 
6236
      // Declare helper variables
6246
6237
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6247
6238
      double tmp1 = (1.0 - Y)/2.0;
6248
6239
      double tmp2 = tmp1*tmp1;
6249
6240
      
6250
 
      // Compute basisvalues.
 
6241
      // Compute basisvalues
6251
6242
      basisvalues[0] = 1.0;
6252
6243
      basisvalues[1] = tmp0;
6253
6244
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6291
6282
      basisvalues[16] *= std::sqrt(27.0);
6292
6283
      basisvalues[15] *= std::sqrt(33.0);
6293
6284
      
6294
 
      // Table(s) of coefficients.
 
6285
      // Table(s) of coefficients
6295
6286
      static const double coefficients0[21] = \
6296
6287
      {0.035074741, -0.10739402, 0.062003968, 0.011320324, -0.30690398, -0.25313014, 0.0, -0.10457267, 0.0, 0.046766322, 0.0, 0.16977679, 0.0093001208, 0.1265818, 0.17230725, 0.0, 0.0, 0.1549843, 0.10187772, -0.039457071, -0.081359109};
6297
6288
      
6495
6486
    case 19:
6496
6487
      {
6497
6488
        
6498
 
      // Array of basisvalues.
 
6489
      // Array of basisvalues
6499
6490
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6500
6491
      
6501
 
      // Declare helper variables.
 
6492
      // Declare helper variables
6502
6493
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6503
6494
      double tmp1 = (1.0 - Y)/2.0;
6504
6495
      double tmp2 = tmp1*tmp1;
6505
6496
      
6506
 
      // Compute basisvalues.
 
6497
      // Compute basisvalues
6507
6498
      basisvalues[0] = 1.0;
6508
6499
      basisvalues[1] = tmp0;
6509
6500
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6547
6538
      basisvalues[16] *= std::sqrt(27.0);
6548
6539
      basisvalues[15] *= std::sqrt(33.0);
6549
6540
      
6550
 
      // Table(s) of coefficients.
 
6541
      // Table(s) of coefficients
6551
6542
      static const double coefficients0[21] = \
6552
6543
      {0.035074741, 0.10739402, 0.062003968, 0.011320324, 0.30690398, -0.25313014, 0.0, -0.10457267, 0.0, 0.046766322, 0.0, -0.16977679, 0.0093001208, -0.1265818, 0.17230725, 0.0, 0.0, -0.1549843, 0.10187772, 0.039457071, -0.081359109};
6553
6544
      
6751
6742
    case 20:
6752
6743
      {
6753
6744
        
6754
 
      // Array of basisvalues.
 
6745
      // Array of basisvalues
6755
6746
      double basisvalues[21] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
6756
6747
      
6757
 
      // Declare helper variables.
 
6748
      // Declare helper variables
6758
6749
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
6759
6750
      double tmp1 = (1.0 - Y)/2.0;
6760
6751
      double tmp2 = tmp1*tmp1;
6761
6752
      
6762
 
      // Compute basisvalues.
 
6753
      // Compute basisvalues
6763
6754
      basisvalues[0] = 1.0;
6764
6755
      basisvalues[1] = tmp0;
6765
6756
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
6803
6794
      basisvalues[16] *= std::sqrt(27.0);
6804
6795
      basisvalues[15] *= std::sqrt(33.0);
6805
6796
      
6806
 
      // Table(s) of coefficients.
 
6797
      // Table(s) of coefficients
6807
6798
      static const double coefficients0[21] = \
6808
6799
      {0.28059793, 0.0, 0.24801587, -0.060375061, 0.0, 0.13500274, 0.0, -0.052286337, 0.0, -0.21044845, 0.0, 0.0, -0.14880193, 0.0, -0.1425991, 0.0, 0.0, 0.0, -0.13583696, 0.0, 0.10847881};
6809
6800
      
7008
6999
    
7009
7000
  }
7010
7001
 
7011
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
7012
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
7002
  /// Evaluate order n derivatives of all basis functions at given point x in cell
 
7003
  virtual void evaluate_basis_derivatives_all(std::size_t n,
7013
7004
                                              double* values,
7014
 
                                              const double* coordinates,
7015
 
                                              const ufc::cell& c) const
 
7005
                                              const double* x,
 
7006
                                              const double* vertex_coordinates,
 
7007
                                              int cell_orientation) const
7016
7008
  {
7017
7009
    // Compute number of derivatives.
7018
7010
    unsigned int num_derivatives = 1;
7031
7023
    // Loop dofs and call evaluate_basis_derivatives.
7032
7024
    for (unsigned int r = 0; r < 21; r++)
7033
7025
    {
7034
 
      evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
 
7026
      evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
7035
7027
      for (unsigned int s = 0; s < num_derivatives; s++)
7036
7028
      {
7037
7029
        values[r*num_derivatives + s] = dof_values[s];
7043
7035
  }
7044
7036
 
7045
7037
  /// Evaluate linear functional for dof i on the function f
7046
 
  virtual double evaluate_dof(unsigned int i,
 
7038
  virtual double evaluate_dof(std::size_t i,
7047
7039
                              const ufc::function& f,
 
7040
                              const double* vertex_coordinates,
 
7041
                              int cell_orientation,
7048
7042
                              const ufc::cell& c) const
7049
7043
  {
7050
 
    // Declare variables for result of evaluation.
 
7044
    // Declare variables for result of evaluation
7051
7045
    double vals[1];
7052
7046
    
7053
 
    // Declare variable for physical coordinates.
 
7047
    // Declare variable for physical coordinates
7054
7048
    double y[2];
7055
 
    const double * const * x = c.coordinates;
7056
7049
    switch (i)
7057
7050
    {
7058
7051
    case 0:
7059
7052
      {
7060
 
        y[0] = x[0][0];
7061
 
      y[1] = x[0][1];
 
7053
        y[0] = vertex_coordinates[0];
 
7054
      y[1] = vertex_coordinates[1];
7062
7055
      f.evaluate(vals, y, c);
7063
7056
      return vals[0];
7064
7057
        break;
7065
7058
      }
7066
7059
    case 1:
7067
7060
      {
7068
 
        y[0] = x[1][0];
7069
 
      y[1] = x[1][1];
 
7061
        y[0] = vertex_coordinates[2];
 
7062
      y[1] = vertex_coordinates[3];
7070
7063
      f.evaluate(vals, y, c);
7071
7064
      return vals[0];
7072
7065
        break;
7073
7066
      }
7074
7067
    case 2:
7075
7068
      {
7076
 
        y[0] = x[2][0];
7077
 
      y[1] = x[2][1];
 
7069
        y[0] = vertex_coordinates[4];
 
7070
      y[1] = vertex_coordinates[5];
7078
7071
      f.evaluate(vals, y, c);
7079
7072
      return vals[0];
7080
7073
        break;
7081
7074
      }
7082
7075
    case 3:
7083
7076
      {
7084
 
        y[0] = 0.8*x[1][0] + 0.2*x[2][0];
7085
 
      y[1] = 0.8*x[1][1] + 0.2*x[2][1];
 
7077
        y[0] = 0.8*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7078
      y[1] = 0.8*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7086
7079
      f.evaluate(vals, y, c);
7087
7080
      return vals[0];
7088
7081
        break;
7089
7082
      }
7090
7083
    case 4:
7091
7084
      {
7092
 
        y[0] = 0.6*x[1][0] + 0.4*x[2][0];
7093
 
      y[1] = 0.6*x[1][1] + 0.4*x[2][1];
 
7085
        y[0] = 0.6*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7086
      y[1] = 0.6*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
7094
7087
      f.evaluate(vals, y, c);
7095
7088
      return vals[0];
7096
7089
        break;
7097
7090
      }
7098
7091
    case 5:
7099
7092
      {
7100
 
        y[0] = 0.4*x[1][0] + 0.6*x[2][0];
7101
 
      y[1] = 0.4*x[1][1] + 0.6*x[2][1];
 
7093
        y[0] = 0.4*vertex_coordinates[2] + 0.6*vertex_coordinates[4];
 
7094
      y[1] = 0.4*vertex_coordinates[3] + 0.6*vertex_coordinates[5];
7102
7095
      f.evaluate(vals, y, c);
7103
7096
      return vals[0];
7104
7097
        break;
7105
7098
      }
7106
7099
    case 6:
7107
7100
      {
7108
 
        y[0] = 0.2*x[1][0] + 0.8*x[2][0];
7109
 
      y[1] = 0.2*x[1][1] + 0.8*x[2][1];
 
7101
        y[0] = 0.2*vertex_coordinates[2] + 0.8*vertex_coordinates[4];
 
7102
      y[1] = 0.2*vertex_coordinates[3] + 0.8*vertex_coordinates[5];
7110
7103
      f.evaluate(vals, y, c);
7111
7104
      return vals[0];
7112
7105
        break;
7113
7106
      }
7114
7107
    case 7:
7115
7108
      {
7116
 
        y[0] = 0.8*x[0][0] + 0.2*x[2][0];
7117
 
      y[1] = 0.8*x[0][1] + 0.2*x[2][1];
 
7109
        y[0] = 0.8*vertex_coordinates[0] + 0.2*vertex_coordinates[4];
 
7110
      y[1] = 0.8*vertex_coordinates[1] + 0.2*vertex_coordinates[5];
7118
7111
      f.evaluate(vals, y, c);
7119
7112
      return vals[0];
7120
7113
        break;
7121
7114
      }
7122
7115
    case 8:
7123
7116
      {
7124
 
        y[0] = 0.6*x[0][0] + 0.4*x[2][0];
7125
 
      y[1] = 0.6*x[0][1] + 0.4*x[2][1];
 
7117
        y[0] = 0.6*vertex_coordinates[0] + 0.4*vertex_coordinates[4];
 
7118
      y[1] = 0.6*vertex_coordinates[1] + 0.4*vertex_coordinates[5];
7126
7119
      f.evaluate(vals, y, c);
7127
7120
      return vals[0];
7128
7121
        break;
7129
7122
      }
7130
7123
    case 9:
7131
7124
      {
7132
 
        y[0] = 0.4*x[0][0] + 0.6*x[2][0];
7133
 
      y[1] = 0.4*x[0][1] + 0.6*x[2][1];
 
7125
        y[0] = 0.4*vertex_coordinates[0] + 0.6*vertex_coordinates[4];
 
7126
      y[1] = 0.4*vertex_coordinates[1] + 0.6*vertex_coordinates[5];
7134
7127
      f.evaluate(vals, y, c);
7135
7128
      return vals[0];
7136
7129
        break;
7137
7130
      }
7138
7131
    case 10:
7139
7132
      {
7140
 
        y[0] = 0.2*x[0][0] + 0.8*x[2][0];
7141
 
      y[1] = 0.2*x[0][1] + 0.8*x[2][1];
 
7133
        y[0] = 0.2*vertex_coordinates[0] + 0.8*vertex_coordinates[4];
 
7134
      y[1] = 0.2*vertex_coordinates[1] + 0.8*vertex_coordinates[5];
7142
7135
      f.evaluate(vals, y, c);
7143
7136
      return vals[0];
7144
7137
        break;
7145
7138
      }
7146
7139
    case 11:
7147
7140
      {
7148
 
        y[0] = 0.8*x[0][0] + 0.2*x[1][0];
7149
 
      y[1] = 0.8*x[0][1] + 0.2*x[1][1];
 
7141
        y[0] = 0.8*vertex_coordinates[0] + 0.2*vertex_coordinates[2];
 
7142
      y[1] = 0.8*vertex_coordinates[1] + 0.2*vertex_coordinates[3];
7150
7143
      f.evaluate(vals, y, c);
7151
7144
      return vals[0];
7152
7145
        break;
7153
7146
      }
7154
7147
    case 12:
7155
7148
      {
7156
 
        y[0] = 0.6*x[0][0] + 0.4*x[1][0];
7157
 
      y[1] = 0.6*x[0][1] + 0.4*x[1][1];
 
7149
        y[0] = 0.6*vertex_coordinates[0] + 0.4*vertex_coordinates[2];
 
7150
      y[1] = 0.6*vertex_coordinates[1] + 0.4*vertex_coordinates[3];
7158
7151
      f.evaluate(vals, y, c);
7159
7152
      return vals[0];
7160
7153
        break;
7161
7154
      }
7162
7155
    case 13:
7163
7156
      {
7164
 
        y[0] = 0.4*x[0][0] + 0.6*x[1][0];
7165
 
      y[1] = 0.4*x[0][1] + 0.6*x[1][1];
 
7157
        y[0] = 0.4*vertex_coordinates[0] + 0.6*vertex_coordinates[2];
 
7158
      y[1] = 0.4*vertex_coordinates[1] + 0.6*vertex_coordinates[3];
7166
7159
      f.evaluate(vals, y, c);
7167
7160
      return vals[0];
7168
7161
        break;
7169
7162
      }
7170
7163
    case 14:
7171
7164
      {
7172
 
        y[0] = 0.2*x[0][0] + 0.8*x[1][0];
7173
 
      y[1] = 0.2*x[0][1] + 0.8*x[1][1];
 
7165
        y[0] = 0.2*vertex_coordinates[0] + 0.8*vertex_coordinates[2];
 
7166
      y[1] = 0.2*vertex_coordinates[1] + 0.8*vertex_coordinates[3];
7174
7167
      f.evaluate(vals, y, c);
7175
7168
      return vals[0];
7176
7169
        break;
7177
7170
      }
7178
7171
    case 15:
7179
7172
      {
7180
 
        y[0] = 0.6*x[0][0] + 0.2*x[1][0] + 0.2*x[2][0];
7181
 
      y[1] = 0.6*x[0][1] + 0.2*x[1][1] + 0.2*x[2][1];
 
7173
        y[0] = 0.6*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7174
      y[1] = 0.6*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7182
7175
      f.evaluate(vals, y, c);
7183
7176
      return vals[0];
7184
7177
        break;
7185
7178
      }
7186
7179
    case 16:
7187
7180
      {
7188
 
        y[0] = 0.4*x[0][0] + 0.4*x[1][0] + 0.2*x[2][0];
7189
 
      y[1] = 0.4*x[0][1] + 0.4*x[1][1] + 0.2*x[2][1];
 
7181
        y[0] = 0.4*vertex_coordinates[0] + 0.4*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7182
      y[1] = 0.4*vertex_coordinates[1] + 0.4*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7190
7183
      f.evaluate(vals, y, c);
7191
7184
      return vals[0];
7192
7185
        break;
7193
7186
      }
7194
7187
    case 17:
7195
7188
      {
7196
 
        y[0] = 0.2*x[0][0] + 0.6*x[1][0] + 0.2*x[2][0];
7197
 
      y[1] = 0.2*x[0][1] + 0.6*x[1][1] + 0.2*x[2][1];
 
7189
        y[0] = 0.2*vertex_coordinates[0] + 0.6*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7190
      y[1] = 0.2*vertex_coordinates[1] + 0.6*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7198
7191
      f.evaluate(vals, y, c);
7199
7192
      return vals[0];
7200
7193
        break;
7201
7194
      }
7202
7195
    case 18:
7203
7196
      {
7204
 
        y[0] = 0.4*x[0][0] + 0.2*x[1][0] + 0.4*x[2][0];
7205
 
      y[1] = 0.4*x[0][1] + 0.2*x[1][1] + 0.4*x[2][1];
 
7197
        y[0] = 0.4*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7198
      y[1] = 0.4*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
7206
7199
      f.evaluate(vals, y, c);
7207
7200
      return vals[0];
7208
7201
        break;
7209
7202
      }
7210
7203
    case 19:
7211
7204
      {
7212
 
        y[0] = 0.2*x[0][0] + 0.4*x[1][0] + 0.4*x[2][0];
7213
 
      y[1] = 0.2*x[0][1] + 0.4*x[1][1] + 0.4*x[2][1];
 
7205
        y[0] = 0.2*vertex_coordinates[0] + 0.4*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7206
      y[1] = 0.2*vertex_coordinates[1] + 0.4*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
7214
7207
      f.evaluate(vals, y, c);
7215
7208
      return vals[0];
7216
7209
        break;
7217
7210
      }
7218
7211
    case 20:
7219
7212
      {
7220
 
        y[0] = 0.2*x[0][0] + 0.2*x[1][0] + 0.6*x[2][0];
7221
 
      y[1] = 0.2*x[0][1] + 0.2*x[1][1] + 0.6*x[2][1];
 
7213
        y[0] = 0.2*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.6*vertex_coordinates[4];
 
7214
      y[1] = 0.2*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.6*vertex_coordinates[5];
7222
7215
      f.evaluate(vals, y, c);
7223
7216
      return vals[0];
7224
7217
        break;
7231
7224
  /// Evaluate linear functionals for all dofs on the function f
7232
7225
  virtual void evaluate_dofs(double* values,
7233
7226
                             const ufc::function& f,
 
7227
                             const double* vertex_coordinates,
 
7228
                             int cell_orientation,
7234
7229
                             const ufc::cell& c) const
7235
7230
  {
7236
 
    // Declare variables for result of evaluation.
 
7231
    // Declare variables for result of evaluation
7237
7232
    double vals[1];
7238
7233
    
7239
 
    // Declare variable for physical coordinates.
 
7234
    // Declare variable for physical coordinates
7240
7235
    double y[2];
7241
 
    const double * const * x = c.coordinates;
7242
 
    y[0] = x[0][0];
7243
 
    y[1] = x[0][1];
 
7236
    y[0] = vertex_coordinates[0];
 
7237
    y[1] = vertex_coordinates[1];
7244
7238
    f.evaluate(vals, y, c);
7245
7239
    values[0] = vals[0];
7246
 
    y[0] = x[1][0];
7247
 
    y[1] = x[1][1];
 
7240
    y[0] = vertex_coordinates[2];
 
7241
    y[1] = vertex_coordinates[3];
7248
7242
    f.evaluate(vals, y, c);
7249
7243
    values[1] = vals[0];
7250
 
    y[0] = x[2][0];
7251
 
    y[1] = x[2][1];
 
7244
    y[0] = vertex_coordinates[4];
 
7245
    y[1] = vertex_coordinates[5];
7252
7246
    f.evaluate(vals, y, c);
7253
7247
    values[2] = vals[0];
7254
 
    y[0] = 0.8*x[1][0] + 0.2*x[2][0];
7255
 
    y[1] = 0.8*x[1][1] + 0.2*x[2][1];
 
7248
    y[0] = 0.8*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7249
    y[1] = 0.8*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7256
7250
    f.evaluate(vals, y, c);
7257
7251
    values[3] = vals[0];
7258
 
    y[0] = 0.6*x[1][0] + 0.4*x[2][0];
7259
 
    y[1] = 0.6*x[1][1] + 0.4*x[2][1];
 
7252
    y[0] = 0.6*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7253
    y[1] = 0.6*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
7260
7254
    f.evaluate(vals, y, c);
7261
7255
    values[4] = vals[0];
7262
 
    y[0] = 0.4*x[1][0] + 0.6*x[2][0];
7263
 
    y[1] = 0.4*x[1][1] + 0.6*x[2][1];
 
7256
    y[0] = 0.4*vertex_coordinates[2] + 0.6*vertex_coordinates[4];
 
7257
    y[1] = 0.4*vertex_coordinates[3] + 0.6*vertex_coordinates[5];
7264
7258
    f.evaluate(vals, y, c);
7265
7259
    values[5] = vals[0];
7266
 
    y[0] = 0.2*x[1][0] + 0.8*x[2][0];
7267
 
    y[1] = 0.2*x[1][1] + 0.8*x[2][1];
 
7260
    y[0] = 0.2*vertex_coordinates[2] + 0.8*vertex_coordinates[4];
 
7261
    y[1] = 0.2*vertex_coordinates[3] + 0.8*vertex_coordinates[5];
7268
7262
    f.evaluate(vals, y, c);
7269
7263
    values[6] = vals[0];
7270
 
    y[0] = 0.8*x[0][0] + 0.2*x[2][0];
7271
 
    y[1] = 0.8*x[0][1] + 0.2*x[2][1];
 
7264
    y[0] = 0.8*vertex_coordinates[0] + 0.2*vertex_coordinates[4];
 
7265
    y[1] = 0.8*vertex_coordinates[1] + 0.2*vertex_coordinates[5];
7272
7266
    f.evaluate(vals, y, c);
7273
7267
    values[7] = vals[0];
7274
 
    y[0] = 0.6*x[0][0] + 0.4*x[2][0];
7275
 
    y[1] = 0.6*x[0][1] + 0.4*x[2][1];
 
7268
    y[0] = 0.6*vertex_coordinates[0] + 0.4*vertex_coordinates[4];
 
7269
    y[1] = 0.6*vertex_coordinates[1] + 0.4*vertex_coordinates[5];
7276
7270
    f.evaluate(vals, y, c);
7277
7271
    values[8] = vals[0];
7278
 
    y[0] = 0.4*x[0][0] + 0.6*x[2][0];
7279
 
    y[1] = 0.4*x[0][1] + 0.6*x[2][1];
 
7272
    y[0] = 0.4*vertex_coordinates[0] + 0.6*vertex_coordinates[4];
 
7273
    y[1] = 0.4*vertex_coordinates[1] + 0.6*vertex_coordinates[5];
7280
7274
    f.evaluate(vals, y, c);
7281
7275
    values[9] = vals[0];
7282
 
    y[0] = 0.2*x[0][0] + 0.8*x[2][0];
7283
 
    y[1] = 0.2*x[0][1] + 0.8*x[2][1];
 
7276
    y[0] = 0.2*vertex_coordinates[0] + 0.8*vertex_coordinates[4];
 
7277
    y[1] = 0.2*vertex_coordinates[1] + 0.8*vertex_coordinates[5];
7284
7278
    f.evaluate(vals, y, c);
7285
7279
    values[10] = vals[0];
7286
 
    y[0] = 0.8*x[0][0] + 0.2*x[1][0];
7287
 
    y[1] = 0.8*x[0][1] + 0.2*x[1][1];
 
7280
    y[0] = 0.8*vertex_coordinates[0] + 0.2*vertex_coordinates[2];
 
7281
    y[1] = 0.8*vertex_coordinates[1] + 0.2*vertex_coordinates[3];
7288
7282
    f.evaluate(vals, y, c);
7289
7283
    values[11] = vals[0];
7290
 
    y[0] = 0.6*x[0][0] + 0.4*x[1][0];
7291
 
    y[1] = 0.6*x[0][1] + 0.4*x[1][1];
 
7284
    y[0] = 0.6*vertex_coordinates[0] + 0.4*vertex_coordinates[2];
 
7285
    y[1] = 0.6*vertex_coordinates[1] + 0.4*vertex_coordinates[3];
7292
7286
    f.evaluate(vals, y, c);
7293
7287
    values[12] = vals[0];
7294
 
    y[0] = 0.4*x[0][0] + 0.6*x[1][0];
7295
 
    y[1] = 0.4*x[0][1] + 0.6*x[1][1];
 
7288
    y[0] = 0.4*vertex_coordinates[0] + 0.6*vertex_coordinates[2];
 
7289
    y[1] = 0.4*vertex_coordinates[1] + 0.6*vertex_coordinates[3];
7296
7290
    f.evaluate(vals, y, c);
7297
7291
    values[13] = vals[0];
7298
 
    y[0] = 0.2*x[0][0] + 0.8*x[1][0];
7299
 
    y[1] = 0.2*x[0][1] + 0.8*x[1][1];
 
7292
    y[0] = 0.2*vertex_coordinates[0] + 0.8*vertex_coordinates[2];
 
7293
    y[1] = 0.2*vertex_coordinates[1] + 0.8*vertex_coordinates[3];
7300
7294
    f.evaluate(vals, y, c);
7301
7295
    values[14] = vals[0];
7302
 
    y[0] = 0.6*x[0][0] + 0.2*x[1][0] + 0.2*x[2][0];
7303
 
    y[1] = 0.6*x[0][1] + 0.2*x[1][1] + 0.2*x[2][1];
 
7296
    y[0] = 0.6*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7297
    y[1] = 0.6*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7304
7298
    f.evaluate(vals, y, c);
7305
7299
    values[15] = vals[0];
7306
 
    y[0] = 0.4*x[0][0] + 0.4*x[1][0] + 0.2*x[2][0];
7307
 
    y[1] = 0.4*x[0][1] + 0.4*x[1][1] + 0.2*x[2][1];
 
7300
    y[0] = 0.4*vertex_coordinates[0] + 0.4*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7301
    y[1] = 0.4*vertex_coordinates[1] + 0.4*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7308
7302
    f.evaluate(vals, y, c);
7309
7303
    values[16] = vals[0];
7310
 
    y[0] = 0.2*x[0][0] + 0.6*x[1][0] + 0.2*x[2][0];
7311
 
    y[1] = 0.2*x[0][1] + 0.6*x[1][1] + 0.2*x[2][1];
 
7304
    y[0] = 0.2*vertex_coordinates[0] + 0.6*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7305
    y[1] = 0.2*vertex_coordinates[1] + 0.6*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
7312
7306
    f.evaluate(vals, y, c);
7313
7307
    values[17] = vals[0];
7314
 
    y[0] = 0.4*x[0][0] + 0.2*x[1][0] + 0.4*x[2][0];
7315
 
    y[1] = 0.4*x[0][1] + 0.2*x[1][1] + 0.4*x[2][1];
 
7308
    y[0] = 0.4*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7309
    y[1] = 0.4*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
7316
7310
    f.evaluate(vals, y, c);
7317
7311
    values[18] = vals[0];
7318
 
    y[0] = 0.2*x[0][0] + 0.4*x[1][0] + 0.4*x[2][0];
7319
 
    y[1] = 0.2*x[0][1] + 0.4*x[1][1] + 0.4*x[2][1];
 
7312
    y[0] = 0.2*vertex_coordinates[0] + 0.4*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7313
    y[1] = 0.2*vertex_coordinates[1] + 0.4*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
7320
7314
    f.evaluate(vals, y, c);
7321
7315
    values[19] = vals[0];
7322
 
    y[0] = 0.2*x[0][0] + 0.2*x[1][0] + 0.6*x[2][0];
7323
 
    y[1] = 0.2*x[0][1] + 0.2*x[1][1] + 0.6*x[2][1];
 
7316
    y[0] = 0.2*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.6*vertex_coordinates[4];
 
7317
    y[1] = 0.2*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.6*vertex_coordinates[5];
7324
7318
    f.evaluate(vals, y, c);
7325
7319
    values[20] = vals[0];
7326
7320
  }
7328
7322
  /// Interpolate vertex values from dof values
7329
7323
  virtual void interpolate_vertex_values(double* vertex_values,
7330
7324
                                         const double* dof_values,
 
7325
                                         const double* vertex_coordinates,
 
7326
                                         int cell_orientation,
7331
7327
                                         const ufc::cell& c) const
7332
7328
  {
7333
7329
    // Evaluate function and change variables
7341
7337
                                       const double* xhat,
7342
7338
                                       const ufc::cell& c) const
7343
7339
  {
7344
 
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
7340
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
7345
7341
  }
7346
7342
 
7347
7343
  /// Map from coordinate x in cell to coordinate xhat in reference cell
7349
7345
                                     const double* x,
7350
7346
                                     const ufc::cell& c) const
7351
7347
  {
7352
 
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
7348
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
7353
7349
  }
7354
7350
 
7355
7351
  /// Return the number of sub elements (for a mixed element)
7356
 
  virtual unsigned int num_sub_elements() const
 
7352
  virtual std::size_t num_sub_elements() const
7357
7353
  {
7358
7354
    return 0;
7359
7355
  }
7360
7356
 
7361
7357
  /// Create a new finite element for sub element i (for a mixed element)
7362
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
7358
  virtual ufc::finite_element* create_sub_element(std::size_t i) const
7363
7359
  {
7364
7360
    return 0;
7365
7361
  }
7377
7373
 
7378
7374
class p5tri_dofmap_0: public ufc::dofmap
7379
7375
{
7380
 
private:
7381
 
 
7382
 
  unsigned int _global_dimension;
7383
7376
public:
7384
7377
 
7385
7378
  /// Constructor
7386
7379
  p5tri_dofmap_0() : ufc::dofmap()
7387
7380
  {
7388
 
    _global_dimension = 0;
 
7381
    // Do nothing
7389
7382
  }
7390
7383
 
7391
7384
  /// Destructor
7397
7390
  /// Return a string identifying the dofmap
7398
7391
  virtual const char* signature() const
7399
7392
  {
7400
 
    return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', Space(2)), 5, None)";
 
7393
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 5, None)";
7401
7394
  }
7402
7395
 
7403
7396
  /// Return true iff mesh entities of topological dimension d are needed
7404
 
  virtual bool needs_mesh_entities(unsigned int d) const
 
7397
  virtual bool needs_mesh_entities(std::size_t d) const
7405
7398
  {
7406
7399
    switch (d)
7407
7400
    {
7425
7418
    return false;
7426
7419
  }
7427
7420
 
7428
 
  /// Initialize dofmap for mesh (return true iff init_cell() is needed)
7429
 
  virtual bool init_mesh(const ufc::mesh& m)
7430
 
  {
7431
 
    _global_dimension = m.num_entities[0] + 4*m.num_entities[1] + 6*m.num_entities[2];
7432
 
    return false;
7433
 
  }
7434
 
 
7435
 
  /// Initialize dofmap for given cell
7436
 
  virtual void init_cell(const ufc::mesh& m,
7437
 
                         const ufc::cell& c)
7438
 
  {
7439
 
    // Do nothing
7440
 
  }
7441
 
 
7442
 
  /// Finish initialization of dofmap for cells
7443
 
  virtual void init_cell_finalize()
7444
 
  {
7445
 
    // Do nothing
7446
 
  }
7447
 
 
7448
7421
  /// Return the topological dimension of the associated cell shape
7449
 
  virtual unsigned int topological_dimension() const
 
7422
  virtual std::size_t topological_dimension() const
7450
7423
  {
7451
7424
    return 2;
7452
7425
  }
7453
7426
 
7454
7427
  /// Return the geometric dimension of the associated cell shape
7455
 
  virtual unsigned int geometric_dimension() const
 
7428
  virtual std::size_t geometric_dimension() const
7456
7429
  {
7457
7430
    return 2;
7458
7431
  }
7459
7432
 
7460
7433
  /// Return the dimension of the global finite element function space
7461
 
  virtual unsigned int global_dimension() const
 
7434
  virtual std::size_t global_dimension(const std::vector<std::size_t>&
 
7435
                                       num_global_entities) const
7462
7436
  {
7463
 
    return _global_dimension;
 
7437
    return num_global_entities[0] + 4*num_global_entities[1] + 6*num_global_entities[2];
7464
7438
  }
7465
7439
 
7466
7440
  /// Return the dimension of the local finite element function space for a cell
7467
 
  virtual unsigned int local_dimension(const ufc::cell& c) const
 
7441
  virtual std::size_t local_dimension(const ufc::cell& c) const
7468
7442
  {
7469
7443
    return 21;
7470
7444
  }
7471
7445
 
7472
7446
  /// Return the maximum dimension of the local finite element function space
7473
 
  virtual unsigned int max_local_dimension() const
 
7447
  virtual std::size_t max_local_dimension() const
7474
7448
  {
7475
7449
    return 21;
7476
7450
  }
7477
7451
 
7478
7452
  /// Return the number of dofs on each cell facet
7479
 
  virtual unsigned int num_facet_dofs() const
 
7453
  virtual std::size_t num_facet_dofs() const
7480
7454
  {
7481
7455
    return 6;
7482
7456
  }
7483
7457
 
7484
7458
  /// Return the number of dofs associated with each cell entity of dimension d
7485
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
7459
  virtual std::size_t num_entity_dofs(std::size_t d) const
7486
7460
  {
7487
7461
    switch (d)
7488
7462
    {
7507
7481
  }
7508
7482
 
7509
7483
  /// Tabulate the local-to-global mapping of dofs on a cell
7510
 
  virtual void tabulate_dofs(unsigned int* dofs,
7511
 
                             const ufc::mesh& m,
 
7484
  virtual void tabulate_dofs(std::size_t* dofs,
 
7485
                             const std::vector<std::size_t>& num_global_entities,
7512
7486
                             const ufc::cell& c) const
7513
7487
  {
7514
7488
    unsigned int offset = 0;
7515
7489
    dofs[0] = offset + c.entity_indices[0][0];
7516
7490
    dofs[1] = offset + c.entity_indices[0][1];
7517
7491
    dofs[2] = offset + c.entity_indices[0][2];
7518
 
    offset += m.num_entities[0];
 
7492
    offset += num_global_entities[0];
7519
7493
    dofs[3] = offset + 4*c.entity_indices[1][0];
7520
7494
    dofs[4] = offset + 4*c.entity_indices[1][0] + 1;
7521
7495
    dofs[5] = offset + 4*c.entity_indices[1][0] + 2;
7528
7502
    dofs[12] = offset + 4*c.entity_indices[1][2] + 1;
7529
7503
    dofs[13] = offset + 4*c.entity_indices[1][2] + 2;
7530
7504
    dofs[14] = offset + 4*c.entity_indices[1][2] + 3;
7531
 
    offset += 4*m.num_entities[1];
 
7505
    offset += 4*num_global_entities[1];
7532
7506
    dofs[15] = offset + 6*c.entity_indices[2][0];
7533
7507
    dofs[16] = offset + 6*c.entity_indices[2][0] + 1;
7534
7508
    dofs[17] = offset + 6*c.entity_indices[2][0] + 2;
7535
7509
    dofs[18] = offset + 6*c.entity_indices[2][0] + 3;
7536
7510
    dofs[19] = offset + 6*c.entity_indices[2][0] + 4;
7537
7511
    dofs[20] = offset + 6*c.entity_indices[2][0] + 5;
7538
 
    offset += 6*m.num_entities[2];
 
7512
    offset += 6*num_global_entities[2];
7539
7513
  }
7540
7514
 
7541
7515
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
7542
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
7543
 
                                   unsigned int facet) const
 
7516
  virtual void tabulate_facet_dofs(std::size_t* dofs,
 
7517
                                   std::size_t facet) const
7544
7518
  {
7545
7519
    switch (facet)
7546
7520
    {
7579
7553
  }
7580
7554
 
7581
7555
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
7582
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
7583
 
                                    unsigned int d, unsigned int i) const
 
7556
  virtual void tabulate_entity_dofs(std::size_t* dofs,
 
7557
                                    std::size_t d, std::size_t i) const
7584
7558
  {
7585
7559
    if (d > 2)
7586
7560
    {
7674
7648
  }
7675
7649
 
7676
7650
  /// Tabulate the coordinates of all dofs on a cell
7677
 
  virtual void tabulate_coordinates(double** coordinates,
7678
 
                                    const ufc::cell& c) const
 
7651
  virtual void tabulate_coordinates(double** dof_coordinates,
 
7652
                                    const double* vertex_coordinates) const
7679
7653
  {
7680
 
    const double * const * x = c.coordinates;
7681
 
    
7682
 
    coordinates[0][0] = x[0][0];
7683
 
    coordinates[0][1] = x[0][1];
7684
 
    coordinates[1][0] = x[1][0];
7685
 
    coordinates[1][1] = x[1][1];
7686
 
    coordinates[2][0] = x[2][0];
7687
 
    coordinates[2][1] = x[2][1];
7688
 
    coordinates[3][0] = 0.8*x[1][0] + 0.2*x[2][0];
7689
 
    coordinates[3][1] = 0.8*x[1][1] + 0.2*x[2][1];
7690
 
    coordinates[4][0] = 0.6*x[1][0] + 0.4*x[2][0];
7691
 
    coordinates[4][1] = 0.6*x[1][1] + 0.4*x[2][1];
7692
 
    coordinates[5][0] = 0.4*x[1][0] + 0.6*x[2][0];
7693
 
    coordinates[5][1] = 0.4*x[1][1] + 0.6*x[2][1];
7694
 
    coordinates[6][0] = 0.2*x[1][0] + 0.8*x[2][0];
7695
 
    coordinates[6][1] = 0.2*x[1][1] + 0.8*x[2][1];
7696
 
    coordinates[7][0] = 0.8*x[0][0] + 0.2*x[2][0];
7697
 
    coordinates[7][1] = 0.8*x[0][1] + 0.2*x[2][1];
7698
 
    coordinates[8][0] = 0.6*x[0][0] + 0.4*x[2][0];
7699
 
    coordinates[8][1] = 0.6*x[0][1] + 0.4*x[2][1];
7700
 
    coordinates[9][0] = 0.4*x[0][0] + 0.6*x[2][0];
7701
 
    coordinates[9][1] = 0.4*x[0][1] + 0.6*x[2][1];
7702
 
    coordinates[10][0] = 0.2*x[0][0] + 0.8*x[2][0];
7703
 
    coordinates[10][1] = 0.2*x[0][1] + 0.8*x[2][1];
7704
 
    coordinates[11][0] = 0.8*x[0][0] + 0.2*x[1][0];
7705
 
    coordinates[11][1] = 0.8*x[0][1] + 0.2*x[1][1];
7706
 
    coordinates[12][0] = 0.6*x[0][0] + 0.4*x[1][0];
7707
 
    coordinates[12][1] = 0.6*x[0][1] + 0.4*x[1][1];
7708
 
    coordinates[13][0] = 0.4*x[0][0] + 0.6*x[1][0];
7709
 
    coordinates[13][1] = 0.4*x[0][1] + 0.6*x[1][1];
7710
 
    coordinates[14][0] = 0.2*x[0][0] + 0.8*x[1][0];
7711
 
    coordinates[14][1] = 0.2*x[0][1] + 0.8*x[1][1];
7712
 
    coordinates[15][0] = 0.6*x[0][0] + 0.2*x[1][0] + 0.2*x[2][0];
7713
 
    coordinates[15][1] = 0.6*x[0][1] + 0.2*x[1][1] + 0.2*x[2][1];
7714
 
    coordinates[16][0] = 0.4*x[0][0] + 0.4*x[1][0] + 0.2*x[2][0];
7715
 
    coordinates[16][1] = 0.4*x[0][1] + 0.4*x[1][1] + 0.2*x[2][1];
7716
 
    coordinates[17][0] = 0.2*x[0][0] + 0.6*x[1][0] + 0.2*x[2][0];
7717
 
    coordinates[17][1] = 0.2*x[0][1] + 0.6*x[1][1] + 0.2*x[2][1];
7718
 
    coordinates[18][0] = 0.4*x[0][0] + 0.2*x[1][0] + 0.4*x[2][0];
7719
 
    coordinates[18][1] = 0.4*x[0][1] + 0.2*x[1][1] + 0.4*x[2][1];
7720
 
    coordinates[19][0] = 0.2*x[0][0] + 0.4*x[1][0] + 0.4*x[2][0];
7721
 
    coordinates[19][1] = 0.2*x[0][1] + 0.4*x[1][1] + 0.4*x[2][1];
7722
 
    coordinates[20][0] = 0.2*x[0][0] + 0.2*x[1][0] + 0.6*x[2][0];
7723
 
    coordinates[20][1] = 0.2*x[0][1] + 0.2*x[1][1] + 0.6*x[2][1];
 
7654
    dof_coordinates[0][0] = vertex_coordinates[0];
 
7655
    dof_coordinates[0][1] = vertex_coordinates[1];
 
7656
    dof_coordinates[1][0] = vertex_coordinates[2];
 
7657
    dof_coordinates[1][1] = vertex_coordinates[3];
 
7658
    dof_coordinates[2][0] = vertex_coordinates[4];
 
7659
    dof_coordinates[2][1] = vertex_coordinates[5];
 
7660
    dof_coordinates[3][0] = 0.8*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7661
    dof_coordinates[3][1] = 0.8*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
 
7662
    dof_coordinates[4][0] = 0.6*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7663
    dof_coordinates[4][1] = 0.6*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
 
7664
    dof_coordinates[5][0] = 0.4*vertex_coordinates[2] + 0.6*vertex_coordinates[4];
 
7665
    dof_coordinates[5][1] = 0.4*vertex_coordinates[3] + 0.6*vertex_coordinates[5];
 
7666
    dof_coordinates[6][0] = 0.2*vertex_coordinates[2] + 0.8*vertex_coordinates[4];
 
7667
    dof_coordinates[6][1] = 0.2*vertex_coordinates[3] + 0.8*vertex_coordinates[5];
 
7668
    dof_coordinates[7][0] = 0.8*vertex_coordinates[0] + 0.2*vertex_coordinates[4];
 
7669
    dof_coordinates[7][1] = 0.8*vertex_coordinates[1] + 0.2*vertex_coordinates[5];
 
7670
    dof_coordinates[8][0] = 0.6*vertex_coordinates[0] + 0.4*vertex_coordinates[4];
 
7671
    dof_coordinates[8][1] = 0.6*vertex_coordinates[1] + 0.4*vertex_coordinates[5];
 
7672
    dof_coordinates[9][0] = 0.4*vertex_coordinates[0] + 0.6*vertex_coordinates[4];
 
7673
    dof_coordinates[9][1] = 0.4*vertex_coordinates[1] + 0.6*vertex_coordinates[5];
 
7674
    dof_coordinates[10][0] = 0.2*vertex_coordinates[0] + 0.8*vertex_coordinates[4];
 
7675
    dof_coordinates[10][1] = 0.2*vertex_coordinates[1] + 0.8*vertex_coordinates[5];
 
7676
    dof_coordinates[11][0] = 0.8*vertex_coordinates[0] + 0.2*vertex_coordinates[2];
 
7677
    dof_coordinates[11][1] = 0.8*vertex_coordinates[1] + 0.2*vertex_coordinates[3];
 
7678
    dof_coordinates[12][0] = 0.6*vertex_coordinates[0] + 0.4*vertex_coordinates[2];
 
7679
    dof_coordinates[12][1] = 0.6*vertex_coordinates[1] + 0.4*vertex_coordinates[3];
 
7680
    dof_coordinates[13][0] = 0.4*vertex_coordinates[0] + 0.6*vertex_coordinates[2];
 
7681
    dof_coordinates[13][1] = 0.4*vertex_coordinates[1] + 0.6*vertex_coordinates[3];
 
7682
    dof_coordinates[14][0] = 0.2*vertex_coordinates[0] + 0.8*vertex_coordinates[2];
 
7683
    dof_coordinates[14][1] = 0.2*vertex_coordinates[1] + 0.8*vertex_coordinates[3];
 
7684
    dof_coordinates[15][0] = 0.6*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7685
    dof_coordinates[15][1] = 0.6*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
 
7686
    dof_coordinates[16][0] = 0.4*vertex_coordinates[0] + 0.4*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7687
    dof_coordinates[16][1] = 0.4*vertex_coordinates[1] + 0.4*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
 
7688
    dof_coordinates[17][0] = 0.2*vertex_coordinates[0] + 0.6*vertex_coordinates[2] + 0.2*vertex_coordinates[4];
 
7689
    dof_coordinates[17][1] = 0.2*vertex_coordinates[1] + 0.6*vertex_coordinates[3] + 0.2*vertex_coordinates[5];
 
7690
    dof_coordinates[18][0] = 0.4*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7691
    dof_coordinates[18][1] = 0.4*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
 
7692
    dof_coordinates[19][0] = 0.2*vertex_coordinates[0] + 0.4*vertex_coordinates[2] + 0.4*vertex_coordinates[4];
 
7693
    dof_coordinates[19][1] = 0.2*vertex_coordinates[1] + 0.4*vertex_coordinates[3] + 0.4*vertex_coordinates[5];
 
7694
    dof_coordinates[20][0] = 0.2*vertex_coordinates[0] + 0.2*vertex_coordinates[2] + 0.6*vertex_coordinates[4];
 
7695
    dof_coordinates[20][1] = 0.2*vertex_coordinates[1] + 0.2*vertex_coordinates[3] + 0.6*vertex_coordinates[5];
7724
7696
  }
7725
7697
 
7726
7698
  /// Return the number of sub dofmaps (for a mixed element)
7727
 
  virtual unsigned int num_sub_dofmaps() const
 
7699
  virtual std::size_t num_sub_dofmaps() const
7728
7700
  {
7729
7701
    return 0;
7730
7702
  }
7731
7703
 
7732
7704
  /// Create a new dofmap for sub dofmap i (for a mixed element)
7733
 
  virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
 
7705
  virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
7734
7706
  {
7735
7707
    return 0;
7736
7708
  }