← 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/Conditional.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/Conditional.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)), 2, None)";
56
 
  }
57
 
 
58
 
  /// Return the cell shape
59
 
  virtual ufc::shape cell_shape() const
60
 
  {
61
 
    return ufc::triangle;
62
 
  }
63
 
 
64
 
  /// Return the topological dimension of the cell shape
65
 
  virtual unsigned int topological_dimension() const
66
 
  {
67
 
    return 2;
68
 
  }
69
 
 
70
 
  /// Return the geometric dimension of the cell shape
71
 
  virtual unsigned int geometric_dimension() const
72
 
  {
73
 
    return 2;
74
 
  }
75
 
 
76
 
  /// Return the dimension of the finite element function space
77
 
  virtual unsigned int space_dimension() const
 
55
    return "FiniteElement('Real', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
 
56
  }
 
57
 
 
58
  /// Return the cell shape
 
59
  virtual ufc::shape cell_shape() const
 
60
  {
 
61
    return ufc::triangle;
 
62
  }
 
63
 
 
64
  /// Return the topological dimension of the cell shape
 
65
  virtual std::size_t topological_dimension() const
 
66
  {
 
67
    return 2;
 
68
  }
 
69
 
 
70
  /// Return the geometric dimension of the cell shape
 
71
  virtual std::size_t geometric_dimension() const
 
72
  {
 
73
    return 2;
 
74
  }
 
75
 
 
76
  /// Return the dimension of the finite element function space
 
77
  virtual std::size_t space_dimension() const
 
78
  {
 
79
    return 1;
 
80
  }
 
81
 
 
82
  /// Return the rank of the value space
 
83
  virtual std::size_t value_rank() const
 
84
  {
 
85
    return 0;
 
86
  }
 
87
 
 
88
  /// Return the dimension of the value space for axis i
 
89
  virtual std::size_t value_dimension(std::size_t i) const
 
90
  {
 
91
    return 1;
 
92
  }
 
93
 
 
94
  /// Evaluate basis function i at given point x in cell
 
95
  virtual void evaluate_basis(std::size_t i,
 
96
                              double* values,
 
97
                              const double* x,
 
98
                              const double* vertex_coordinates,
 
99
                              int cell_orientation) const
 
100
  {
 
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
    
 
110
    
 
111
    // Compute constants
 
112
    
 
113
    // Get coordinates and map to the reference (FIAT) element
 
114
    
 
115
    // Reset values
 
116
    *values = 0.0;
 
117
    
 
118
    // Array of basisvalues
 
119
    double basisvalues[1] = {0.0};
 
120
    
 
121
    // Declare helper variables
 
122
    
 
123
    // Compute basisvalues
 
124
    basisvalues[0] = 1.0;
 
125
    
 
126
    // Table(s) of coefficients
 
127
    static const double coefficients0[1] = \
 
128
    {1.0};
 
129
    
 
130
    // Compute value(s)
 
131
    for (unsigned int r = 0; r < 1; r++)
 
132
    {
 
133
      *values += coefficients0[r]*basisvalues[r];
 
134
    }// end loop over 'r'
 
135
  }
 
136
 
 
137
  /// Evaluate all basis functions at given point x in cell
 
138
  virtual void evaluate_basis_all(double* values,
 
139
                                  const double* x,
 
140
                                  const double* vertex_coordinates,
 
141
                                  int cell_orientation) const
 
142
  {
 
143
    // Element is constant, calling evaluate_basis.
 
144
    evaluate_basis(0, values, x, vertex_coordinates, cell_orientation);
 
145
  }
 
146
 
 
147
  /// Evaluate order n derivatives of basis function i at given point x in cell
 
148
  virtual void evaluate_basis_derivatives(std::size_t i,
 
149
                                          std::size_t n,
 
150
                                          double* values,
 
151
                                          const double* x,
 
152
                                          const double* vertex_coordinates,
 
153
                                          int cell_orientation) const
 
154
  {
 
155
    // Compute Jacobian
 
156
    double J[4];
 
157
    compute_jacobian_triangle_2d(J, vertex_coordinates);
 
158
    
 
159
    // Compute Jacobian inverse and determinant
 
160
    double K[4];
 
161
    double detJ;
 
162
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
 
163
    
 
164
    
 
165
    // Compute constants
 
166
    
 
167
    // Get coordinates and map to the reference (FIAT) element
 
168
    
 
169
    // Compute number of derivatives.
 
170
    unsigned int num_derivatives = 1;
 
171
    for (unsigned int r = 0; r < n; r++)
 
172
    {
 
173
      num_derivatives *= 2;
 
174
    }// end loop over 'r'
 
175
    
 
176
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
177
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
178
    for (unsigned int row = 0; row < num_derivatives; row++)
 
179
    {
 
180
      combinations[row] = new unsigned int [n];
 
181
      for (unsigned int col = 0; col < n; col++)
 
182
        combinations[row][col] = 0;
 
183
    }
 
184
    
 
185
    // Generate combinations of derivatives
 
186
    for (unsigned int row = 1; row < num_derivatives; row++)
 
187
    {
 
188
      for (unsigned int num = 0; num < row; num++)
 
189
      {
 
190
        for (unsigned int col = n-1; col+1 > 0; col--)
 
191
        {
 
192
          if (combinations[row][col] + 1 > 1)
 
193
            combinations[row][col] = 0;
 
194
          else
 
195
          {
 
196
            combinations[row][col] += 1;
 
197
            break;
 
198
          }
 
199
        }
 
200
      }
 
201
    }
 
202
    
 
203
    // Compute inverse of Jacobian
 
204
    const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
 
205
    
 
206
    // Declare transformation matrix
 
207
    // Declare pointer to two dimensional array and initialise
 
208
    double **transform = new double *[num_derivatives];
 
209
    
 
210
    for (unsigned int j = 0; j < num_derivatives; j++)
 
211
    {
 
212
      transform[j] = new double [num_derivatives];
 
213
      for (unsigned int k = 0; k < num_derivatives; k++)
 
214
        transform[j][k] = 1;
 
215
    }
 
216
    
 
217
    // Construct transformation matrix
 
218
    for (unsigned int row = 0; row < num_derivatives; row++)
 
219
    {
 
220
      for (unsigned int col = 0; col < num_derivatives; col++)
 
221
      {
 
222
        for (unsigned int k = 0; k < n; k++)
 
223
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
224
      }
 
225
    }
 
226
    
 
227
    // Reset values. Assuming that values is always an array.
 
228
    for (unsigned int r = 0; r < num_derivatives; r++)
 
229
    {
 
230
      values[r] = 0.0;
 
231
    }// end loop over 'r'
 
232
    
 
233
    
 
234
    // Array of basisvalues
 
235
    double basisvalues[1] = {0.0};
 
236
    
 
237
    // Declare helper variables
 
238
    
 
239
    // Compute basisvalues
 
240
    basisvalues[0] = 1.0;
 
241
    
 
242
    // Table(s) of coefficients
 
243
    static const double coefficients0[1] = \
 
244
    {1.0};
 
245
    
 
246
    // Tables of derivatives of the polynomial base (transpose).
 
247
    static const double dmats0[1][1] = \
 
248
    {{0.0}};
 
249
    
 
250
    static const double dmats1[1][1] = \
 
251
    {{0.0}};
 
252
    
 
253
    // Compute reference derivatives.
 
254
    // Declare pointer to array of derivatives on FIAT element.
 
255
    double *derivatives = new double[num_derivatives];
 
256
    for (unsigned int r = 0; r < num_derivatives; r++)
 
257
    {
 
258
      derivatives[r] = 0.0;
 
259
    }// end loop over 'r'
 
260
    
 
261
    // Declare derivative matrix (of polynomial basis).
 
262
    double dmats[1][1] = \
 
263
    {{1.0}};
 
264
    
 
265
    // Declare (auxiliary) derivative matrix (of polynomial basis).
 
266
    double dmats_old[1][1] = \
 
267
    {{1.0}};
 
268
    
 
269
    // Loop possible derivatives.
 
270
    for (unsigned int r = 0; r < num_derivatives; r++)
 
271
    {
 
272
      // Resetting dmats values to compute next derivative.
 
273
      for (unsigned int t = 0; t < 1; t++)
 
274
      {
 
275
        for (unsigned int u = 0; u < 1; u++)
 
276
        {
 
277
          dmats[t][u] = 0.0;
 
278
          if (t == u)
 
279
          {
 
280
          dmats[t][u] = 1.0;
 
281
          }
 
282
          
 
283
        }// end loop over 'u'
 
284
      }// end loop over 't'
 
285
      
 
286
      // Looping derivative order to generate dmats.
 
287
      for (unsigned int s = 0; s < n; s++)
 
288
      {
 
289
        // Updating dmats_old with new values and resetting dmats.
 
290
        for (unsigned int t = 0; t < 1; t++)
 
291
        {
 
292
          for (unsigned int u = 0; u < 1; u++)
 
293
          {
 
294
            dmats_old[t][u] = dmats[t][u];
 
295
            dmats[t][u] = 0.0;
 
296
          }// end loop over 'u'
 
297
        }// end loop over 't'
 
298
        
 
299
        // Update dmats using an inner product.
 
300
        if (combinations[r][s] == 0)
 
301
        {
 
302
        for (unsigned int t = 0; t < 1; t++)
 
303
        {
 
304
          for (unsigned int u = 0; u < 1; u++)
 
305
          {
 
306
            for (unsigned int tu = 0; tu < 1; tu++)
 
307
            {
 
308
              dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
309
            }// end loop over 'tu'
 
310
          }// end loop over 'u'
 
311
        }// end loop over 't'
 
312
        }
 
313
        
 
314
        if (combinations[r][s] == 1)
 
315
        {
 
316
        for (unsigned int t = 0; t < 1; t++)
 
317
        {
 
318
          for (unsigned int u = 0; u < 1; u++)
 
319
          {
 
320
            for (unsigned int tu = 0; tu < 1; tu++)
 
321
            {
 
322
              dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
323
            }// end loop over 'tu'
 
324
          }// end loop over 'u'
 
325
        }// end loop over 't'
 
326
        }
 
327
        
 
328
      }// end loop over 's'
 
329
      for (unsigned int s = 0; s < 1; s++)
 
330
      {
 
331
        for (unsigned int t = 0; t < 1; t++)
 
332
        {
 
333
          derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
334
        }// end loop over 't'
 
335
      }// end loop over 's'
 
336
    }// end loop over 'r'
 
337
    
 
338
    // Transform derivatives back to physical element
 
339
    for (unsigned int r = 0; r < num_derivatives; r++)
 
340
    {
 
341
      for (unsigned int s = 0; s < num_derivatives; s++)
 
342
      {
 
343
        values[r] += transform[r][s]*derivatives[s];
 
344
      }// end loop over 's'
 
345
    }// end loop over 'r'
 
346
    
 
347
    // Delete pointer to array of derivatives on FIAT element
 
348
    delete [] derivatives;
 
349
    
 
350
    // Delete pointer to array of combinations of derivatives and transform
 
351
    for (unsigned int r = 0; r < num_derivatives; r++)
 
352
    {
 
353
      delete [] combinations[r];
 
354
    }// end loop over 'r'
 
355
    delete [] combinations;
 
356
    for (unsigned int r = 0; r < num_derivatives; r++)
 
357
    {
 
358
      delete [] transform[r];
 
359
    }// end loop over 'r'
 
360
    delete [] transform;
 
361
  }
 
362
 
 
363
  /// Evaluate order n derivatives of all basis functions at given point x in cell
 
364
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
365
                                              double* values,
 
366
                                              const double* x,
 
367
                                              const double* vertex_coordinates,
 
368
                                              int cell_orientation) const
 
369
  {
 
370
    // Element is constant, calling evaluate_basis_derivatives.
 
371
    evaluate_basis_derivatives(0, n, values, x, vertex_coordinates, cell_orientation);
 
372
  }
 
373
 
 
374
  /// Evaluate linear functional for dof i on the function f
 
375
  virtual double evaluate_dof(std::size_t i,
 
376
                              const ufc::function& f,
 
377
                              const double* vertex_coordinates,
 
378
                              int cell_orientation,
 
379
                              const ufc::cell& c) const
 
380
  {
 
381
    // Declare variables for result of evaluation
 
382
    double vals[1];
 
383
    
 
384
    // Declare variable for physical coordinates
 
385
    double y[2];
 
386
    switch (i)
 
387
    {
 
388
    case 0:
 
389
      {
 
390
        y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
 
391
      y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
 
392
      f.evaluate(vals, y, c);
 
393
      return vals[0];
 
394
        break;
 
395
      }
 
396
    }
 
397
    
 
398
    return 0.0;
 
399
  }
 
400
 
 
401
  /// Evaluate linear functionals for all dofs on the function f
 
402
  virtual void evaluate_dofs(double* values,
 
403
                             const ufc::function& f,
 
404
                             const double* vertex_coordinates,
 
405
                             int cell_orientation,
 
406
                             const ufc::cell& c) const
 
407
  {
 
408
    // Declare variables for result of evaluation
 
409
    double vals[1];
 
410
    
 
411
    // Declare variable for physical coordinates
 
412
    double y[2];
 
413
    y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
 
414
    y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
 
415
    f.evaluate(vals, y, c);
 
416
    values[0] = vals[0];
 
417
  }
 
418
 
 
419
  /// Interpolate vertex values from dof values
 
420
  virtual void interpolate_vertex_values(double* vertex_values,
 
421
                                         const double* dof_values,
 
422
                                         const double* vertex_coordinates,
 
423
                                         int cell_orientation,
 
424
                                         const ufc::cell& c) const
 
425
  {
 
426
    // Evaluate function and change variables
 
427
    vertex_values[0] = dof_values[0];
 
428
    vertex_values[1] = dof_values[0];
 
429
    vertex_values[2] = dof_values[0];
 
430
  }
 
431
 
 
432
  /// Map coordinate xhat from reference cell to coordinate x in cell
 
433
  virtual void map_from_reference_cell(double* x,
 
434
                                       const double* xhat,
 
435
                                       const ufc::cell& c) const
 
436
  {
 
437
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
 
438
  }
 
439
 
 
440
  /// Map from coordinate x in cell to coordinate xhat in reference cell
 
441
  virtual void map_to_reference_cell(double* xhat,
 
442
                                     const double* x,
 
443
                                     const ufc::cell& c) const
 
444
  {
 
445
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
 
446
  }
 
447
 
 
448
  /// Return the number of sub elements (for a mixed element)
 
449
  virtual std::size_t num_sub_elements() const
 
450
  {
 
451
    return 0;
 
452
  }
 
453
 
 
454
  /// Create a new finite element for sub element i (for a mixed element)
 
455
  virtual ufc::finite_element* create_sub_element(std::size_t i) const
 
456
  {
 
457
    return 0;
 
458
  }
 
459
 
 
460
  /// Create a new class instance
 
461
  virtual ufc::finite_element* create() const
 
462
  {
 
463
    return new conditional_finite_element_0();
 
464
  }
 
465
 
 
466
};
 
467
 
 
468
/// This class defines the interface for a finite element.
 
469
 
 
470
class conditional_finite_element_1: public ufc::finite_element
 
471
{
 
472
public:
 
473
 
 
474
  /// Constructor
 
475
  conditional_finite_element_1() : ufc::finite_element()
 
476
  {
 
477
    // Do nothing
 
478
  }
 
479
 
 
480
  /// Destructor
 
481
  virtual ~conditional_finite_element_1()
 
482
  {
 
483
    // Do nothing
 
484
  }
 
485
 
 
486
  /// Return a string identifying the finite element
 
487
  virtual const char* signature() const
 
488
  {
 
489
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
 
490
  }
 
491
 
 
492
  /// Return the cell shape
 
493
  virtual ufc::shape cell_shape() const
 
494
  {
 
495
    return ufc::triangle;
 
496
  }
 
497
 
 
498
  /// Return the topological dimension of the cell shape
 
499
  virtual std::size_t topological_dimension() const
 
500
  {
 
501
    return 2;
 
502
  }
 
503
 
 
504
  /// Return the geometric dimension of the cell shape
 
505
  virtual std::size_t geometric_dimension() const
 
506
  {
 
507
    return 2;
 
508
  }
 
509
 
 
510
  /// Return the dimension of the finite element function space
 
511
  virtual std::size_t space_dimension() const
78
512
  {
79
513
    return 6;
80
514
  }
81
515
 
82
516
  /// Return the rank of the value space
83
 
  virtual unsigned int value_rank() const
 
517
  virtual std::size_t value_rank() const
84
518
  {
85
519
    return 0;
86
520
  }
87
521
 
88
522
  /// Return the dimension of the value space for axis i
89
 
  virtual unsigned int value_dimension(unsigned int i) const
 
523
  virtual std::size_t value_dimension(std::size_t i) const
90
524
  {
91
525
    return 1;
92
526
  }
93
527
 
94
 
  /// Evaluate basis function i at given point in cell
95
 
  virtual void evaluate_basis(unsigned int i,
 
528
  /// Evaluate basis function i at given point x in cell
 
529
  virtual void evaluate_basis(std::size_t i,
96
530
                              double* values,
97
 
                              const double* coordinates,
98
 
                              const ufc::cell& c) const
 
531
                              const double* x,
 
532
                              const double* vertex_coordinates,
 
533
                              int cell_orientation) const
99
534
  {
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
 
535
    // Compute Jacobian
 
536
    double J[4];
 
537
    compute_jacobian_triangle_2d(J, vertex_coordinates);
 
538
    
 
539
    // Compute Jacobian inverse and determinant
 
540
    double K[4];
 
541
    double detJ;
 
542
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
 
543
    
113
544
    
114
545
    // Compute constants
115
 
    const double C0 = x[1][0] + x[2][0];
116
 
    const double C1 = x[1][1] + x[2][1];
 
546
    const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
 
547
    const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
117
548
    
118
549
    // 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;
 
550
    double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
 
551
    double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
121
552
    
122
 
    // Reset values.
 
553
    // Reset values
123
554
    *values = 0.0;
124
555
    switch (i)
125
556
    {
126
557
    case 0:
127
558
      {
128
559
        
129
 
      // Array of basisvalues.
 
560
      // Array of basisvalues
130
561
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
131
562
      
132
 
      // Declare helper variables.
 
563
      // Declare helper variables
133
564
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
134
565
      double tmp1 = (1.0 - Y)/2.0;
135
566
      double tmp2 = tmp1*tmp1;
136
567
      
137
 
      // Compute basisvalues.
 
568
      // Compute basisvalues
138
569
      basisvalues[0] = 1.0;
139
570
      basisvalues[1] = tmp0;
140
571
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
148
579
      basisvalues[4] *= std::sqrt(4.5);
149
580
      basisvalues[3] *= std::sqrt(7.5);
150
581
      
151
 
      // Table(s) of coefficients.
 
582
      // Table(s) of coefficients
152
583
      static const double coefficients0[6] = \
153
584
      {0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
154
585
      
155
 
      // Compute value(s).
 
586
      // Compute value(s)
156
587
      for (unsigned int r = 0; r < 6; r++)
157
588
      {
158
589
        *values += coefficients0[r]*basisvalues[r];
162
593
    case 1:
163
594
      {
164
595
        
165
 
      // Array of basisvalues.
 
596
      // Array of basisvalues
166
597
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
167
598
      
168
 
      // Declare helper variables.
 
599
      // Declare helper variables
169
600
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
170
601
      double tmp1 = (1.0 - Y)/2.0;
171
602
      double tmp2 = tmp1*tmp1;
172
603
      
173
 
      // Compute basisvalues.
 
604
      // Compute basisvalues
174
605
      basisvalues[0] = 1.0;
175
606
      basisvalues[1] = tmp0;
176
607
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
184
615
      basisvalues[4] *= std::sqrt(4.5);
185
616
      basisvalues[3] *= std::sqrt(7.5);
186
617
      
187
 
      // Table(s) of coefficients.
 
618
      // Table(s) of coefficients
188
619
      static const double coefficients0[6] = \
189
620
      {0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
190
621
      
191
 
      // Compute value(s).
 
622
      // Compute value(s)
192
623
      for (unsigned int r = 0; r < 6; r++)
193
624
      {
194
625
        *values += coefficients0[r]*basisvalues[r];
198
629
    case 2:
199
630
      {
200
631
        
201
 
      // Array of basisvalues.
 
632
      // Array of basisvalues
202
633
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
203
634
      
204
 
      // Declare helper variables.
 
635
      // Declare helper variables
205
636
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
206
637
      double tmp1 = (1.0 - Y)/2.0;
207
638
      double tmp2 = tmp1*tmp1;
208
639
      
209
 
      // Compute basisvalues.
 
640
      // Compute basisvalues
210
641
      basisvalues[0] = 1.0;
211
642
      basisvalues[1] = tmp0;
212
643
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
220
651
      basisvalues[4] *= std::sqrt(4.5);
221
652
      basisvalues[3] *= std::sqrt(7.5);
222
653
      
223
 
      // Table(s) of coefficients.
 
654
      // Table(s) of coefficients
224
655
      static const double coefficients0[6] = \
225
656
      {0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
226
657
      
227
 
      // Compute value(s).
 
658
      // Compute value(s)
228
659
      for (unsigned int r = 0; r < 6; r++)
229
660
      {
230
661
        *values += coefficients0[r]*basisvalues[r];
234
665
    case 3:
235
666
      {
236
667
        
237
 
      // Array of basisvalues.
 
668
      // Array of basisvalues
238
669
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
239
670
      
240
 
      // Declare helper variables.
 
671
      // Declare helper variables
241
672
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
242
673
      double tmp1 = (1.0 - Y)/2.0;
243
674
      double tmp2 = tmp1*tmp1;
244
675
      
245
 
      // Compute basisvalues.
 
676
      // Compute basisvalues
246
677
      basisvalues[0] = 1.0;
247
678
      basisvalues[1] = tmp0;
248
679
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
256
687
      basisvalues[4] *= std::sqrt(4.5);
257
688
      basisvalues[3] *= std::sqrt(7.5);
258
689
      
259
 
      // Table(s) of coefficients.
 
690
      // Table(s) of coefficients
260
691
      static const double coefficients0[6] = \
261
692
      {0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
262
693
      
263
 
      // Compute value(s).
 
694
      // Compute value(s)
264
695
      for (unsigned int r = 0; r < 6; r++)
265
696
      {
266
697
        *values += coefficients0[r]*basisvalues[r];
270
701
    case 4:
271
702
      {
272
703
        
273
 
      // Array of basisvalues.
 
704
      // Array of basisvalues
274
705
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
275
706
      
276
 
      // Declare helper variables.
 
707
      // Declare helper variables
277
708
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
278
709
      double tmp1 = (1.0 - Y)/2.0;
279
710
      double tmp2 = tmp1*tmp1;
280
711
      
281
 
      // Compute basisvalues.
 
712
      // Compute basisvalues
282
713
      basisvalues[0] = 1.0;
283
714
      basisvalues[1] = tmp0;
284
715
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
292
723
      basisvalues[4] *= std::sqrt(4.5);
293
724
      basisvalues[3] *= std::sqrt(7.5);
294
725
      
295
 
      // Table(s) of coefficients.
 
726
      // Table(s) of coefficients
296
727
      static const double coefficients0[6] = \
297
728
      {0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
298
729
      
299
 
      // Compute value(s).
 
730
      // Compute value(s)
300
731
      for (unsigned int r = 0; r < 6; r++)
301
732
      {
302
733
        *values += coefficients0[r]*basisvalues[r];
306
737
    case 5:
307
738
      {
308
739
        
309
 
      // Array of basisvalues.
 
740
      // Array of basisvalues
310
741
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
311
742
      
312
 
      // Declare helper variables.
 
743
      // Declare helper variables
313
744
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
314
745
      double tmp1 = (1.0 - Y)/2.0;
315
746
      double tmp2 = tmp1*tmp1;
316
747
      
317
 
      // Compute basisvalues.
 
748
      // Compute basisvalues
318
749
      basisvalues[0] = 1.0;
319
750
      basisvalues[1] = tmp0;
320
751
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
328
759
      basisvalues[4] *= std::sqrt(4.5);
329
760
      basisvalues[3] *= std::sqrt(7.5);
330
761
      
331
 
      // Table(s) of coefficients.
 
762
      // Table(s) of coefficients
332
763
      static const double coefficients0[6] = \
333
764
      {0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
334
765
      
335
 
      // Compute value(s).
 
766
      // Compute value(s)
336
767
      for (unsigned int r = 0; r < 6; r++)
337
768
      {
338
769
        *values += coefficients0[r]*basisvalues[r];
343
774
    
344
775
  }
345
776
 
346
 
  /// Evaluate all basis functions at given point in cell
 
777
  /// Evaluate all basis functions at given point x in cell
347
778
  virtual void evaluate_basis_all(double* values,
348
 
                                  const double* coordinates,
349
 
                                  const ufc::cell& c) const
 
779
                                  const double* x,
 
780
                                  const double* vertex_coordinates,
 
781
                                  int cell_orientation) const
350
782
  {
351
783
    // Helper variable to hold values of a single dof.
352
784
    double dof_values = 0.0;
353
785
    
354
 
    // Loop dofs and call evaluate_basis.
 
786
    // Loop dofs and call evaluate_basis
355
787
    for (unsigned int r = 0; r < 6; r++)
356
788
    {
357
 
      evaluate_basis(r, &dof_values, coordinates, c);
 
789
      evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
358
790
      values[r] = dof_values;
359
791
    }// end loop over 'r'
360
792
  }
361
793
 
362
 
  /// Evaluate order n derivatives of basis function i at given point in cell
363
 
  virtual void evaluate_basis_derivatives(unsigned int i,
364
 
                                          unsigned int n,
 
794
  /// Evaluate order n derivatives of basis function i at given point x in cell
 
795
  virtual void evaluate_basis_derivatives(std::size_t i,
 
796
                                          std::size_t n,
365
797
                                          double* values,
366
 
                                          const double* coordinates,
367
 
                                          const ufc::cell& c) const
 
798
                                          const double* x,
 
799
                                          const double* vertex_coordinates,
 
800
                                          int cell_orientation) const
368
801
  {
369
 
    // Extract vertex coordinates
370
 
    const double * const * x = c.coordinates;
371
 
    
372
 
    // Compute Jacobian of affine map from reference cell
373
 
    const double J_00 = x[1][0] - x[0][0];
374
 
    const double J_01 = x[2][0] - x[0][0];
375
 
    const double J_10 = x[1][1] - x[0][1];
376
 
    const double J_11 = x[2][1] - x[0][1];
377
 
    
378
 
    // Compute determinant of Jacobian
379
 
    double detJ = J_00*J_11 - J_01*J_10;
380
 
    
381
 
    // Compute inverse of Jacobian
382
 
    const double K_00 =  J_11 / detJ;
383
 
    const double K_01 = -J_01 / detJ;
384
 
    const double K_10 = -J_10 / detJ;
385
 
    const double K_11 =  J_00 / detJ;
 
802
    // Compute Jacobian
 
803
    double J[4];
 
804
    compute_jacobian_triangle_2d(J, vertex_coordinates);
 
805
    
 
806
    // Compute Jacobian inverse and determinant
 
807
    double K[4];
 
808
    double detJ;
 
809
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
 
810
    
386
811
    
387
812
    // Compute constants
388
 
    const double C0 = x[1][0] + x[2][0];
389
 
    const double C1 = x[1][1] + x[2][1];
 
813
    const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
 
814
    const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
390
815
    
391
816
    // Get coordinates and map to the reference (FIAT) element
392
 
    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
393
 
    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
 
817
    double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
 
818
    double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
394
819
    
395
820
    // Compute number of derivatives.
396
821
    unsigned int num_derivatives = 1;
427
852
    }
428
853
    
429
854
    // Compute inverse of Jacobian
430
 
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
 
855
    const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
431
856
    
432
857
    // Declare transformation matrix
433
858
    // Declare pointer to two dimensional array and initialise
461
886
    case 0:
462
887
      {
463
888
        
464
 
      // Array of basisvalues.
 
889
      // Array of basisvalues
465
890
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
466
891
      
467
 
      // Declare helper variables.
 
892
      // Declare helper variables
468
893
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
469
894
      double tmp1 = (1.0 - Y)/2.0;
470
895
      double tmp2 = tmp1*tmp1;
471
896
      
472
 
      // Compute basisvalues.
 
897
      // Compute basisvalues
473
898
      basisvalues[0] = 1.0;
474
899
      basisvalues[1] = tmp0;
475
900
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
483
908
      basisvalues[4] *= std::sqrt(4.5);
484
909
      basisvalues[3] *= std::sqrt(7.5);
485
910
      
486
 
      // Table(s) of coefficients.
 
911
      // Table(s) of coefficients
487
912
      static const double coefficients0[6] = \
488
913
      {0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
489
914
      
627
1052
    case 1:
628
1053
      {
629
1054
        
630
 
      // Array of basisvalues.
 
1055
      // Array of basisvalues
631
1056
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
632
1057
      
633
 
      // Declare helper variables.
 
1058
      // Declare helper variables
634
1059
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
635
1060
      double tmp1 = (1.0 - Y)/2.0;
636
1061
      double tmp2 = tmp1*tmp1;
637
1062
      
638
 
      // Compute basisvalues.
 
1063
      // Compute basisvalues
639
1064
      basisvalues[0] = 1.0;
640
1065
      basisvalues[1] = tmp0;
641
1066
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
649
1074
      basisvalues[4] *= std::sqrt(4.5);
650
1075
      basisvalues[3] *= std::sqrt(7.5);
651
1076
      
652
 
      // Table(s) of coefficients.
 
1077
      // Table(s) of coefficients
653
1078
      static const double coefficients0[6] = \
654
1079
      {0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
655
1080
      
793
1218
    case 2:
794
1219
      {
795
1220
        
796
 
      // Array of basisvalues.
 
1221
      // Array of basisvalues
797
1222
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
798
1223
      
799
 
      // Declare helper variables.
 
1224
      // Declare helper variables
800
1225
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
801
1226
      double tmp1 = (1.0 - Y)/2.0;
802
1227
      double tmp2 = tmp1*tmp1;
803
1228
      
804
 
      // Compute basisvalues.
 
1229
      // Compute basisvalues
805
1230
      basisvalues[0] = 1.0;
806
1231
      basisvalues[1] = tmp0;
807
1232
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
815
1240
      basisvalues[4] *= std::sqrt(4.5);
816
1241
      basisvalues[3] *= std::sqrt(7.5);
817
1242
      
818
 
      // Table(s) of coefficients.
 
1243
      // Table(s) of coefficients
819
1244
      static const double coefficients0[6] = \
820
1245
      {0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
821
1246
      
959
1384
    case 3:
960
1385
      {
961
1386
        
962
 
      // Array of basisvalues.
 
1387
      // Array of basisvalues
963
1388
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
964
1389
      
965
 
      // Declare helper variables.
 
1390
      // Declare helper variables
966
1391
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
967
1392
      double tmp1 = (1.0 - Y)/2.0;
968
1393
      double tmp2 = tmp1*tmp1;
969
1394
      
970
 
      // Compute basisvalues.
 
1395
      // Compute basisvalues
971
1396
      basisvalues[0] = 1.0;
972
1397
      basisvalues[1] = tmp0;
973
1398
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
981
1406
      basisvalues[4] *= std::sqrt(4.5);
982
1407
      basisvalues[3] *= std::sqrt(7.5);
983
1408
      
984
 
      // Table(s) of coefficients.
 
1409
      // Table(s) of coefficients
985
1410
      static const double coefficients0[6] = \
986
1411
      {0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
987
1412
      
1125
1550
    case 4:
1126
1551
      {
1127
1552
        
1128
 
      // Array of basisvalues.
 
1553
      // Array of basisvalues
1129
1554
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1130
1555
      
1131
 
      // Declare helper variables.
 
1556
      // Declare helper variables
1132
1557
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1133
1558
      double tmp1 = (1.0 - Y)/2.0;
1134
1559
      double tmp2 = tmp1*tmp1;
1135
1560
      
1136
 
      // Compute basisvalues.
 
1561
      // Compute basisvalues
1137
1562
      basisvalues[0] = 1.0;
1138
1563
      basisvalues[1] = tmp0;
1139
1564
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1147
1572
      basisvalues[4] *= std::sqrt(4.5);
1148
1573
      basisvalues[3] *= std::sqrt(7.5);
1149
1574
      
1150
 
      // Table(s) of coefficients.
 
1575
      // Table(s) of coefficients
1151
1576
      static const double coefficients0[6] = \
1152
1577
      {0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
1153
1578
      
1291
1716
    case 5:
1292
1717
      {
1293
1718
        
1294
 
      // Array of basisvalues.
 
1719
      // Array of basisvalues
1295
1720
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1296
1721
      
1297
 
      // Declare helper variables.
 
1722
      // Declare helper variables
1298
1723
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1299
1724
      double tmp1 = (1.0 - Y)/2.0;
1300
1725
      double tmp2 = tmp1*tmp1;
1301
1726
      
1302
 
      // Compute basisvalues.
 
1727
      // Compute basisvalues
1303
1728
      basisvalues[0] = 1.0;
1304
1729
      basisvalues[1] = tmp0;
1305
1730
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1313
1738
      basisvalues[4] *= std::sqrt(4.5);
1314
1739
      basisvalues[3] *= std::sqrt(7.5);
1315
1740
      
1316
 
      // Table(s) of coefficients.
 
1741
      // Table(s) of coefficients
1317
1742
      static const double coefficients0[6] = \
1318
1743
      {0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
1319
1744
      
1458
1883
    
1459
1884
  }
1460
1885
 
1461
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
1462
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
1886
  /// Evaluate order n derivatives of all basis functions at given point x in cell
 
1887
  virtual void evaluate_basis_derivatives_all(std::size_t n,
1463
1888
                                              double* values,
1464
 
                                              const double* coordinates,
1465
 
                                              const ufc::cell& c) const
 
1889
                                              const double* x,
 
1890
                                              const double* vertex_coordinates,
 
1891
                                              int cell_orientation) const
1466
1892
  {
1467
1893
    // Compute number of derivatives.
1468
1894
    unsigned int num_derivatives = 1;
1481
1907
    // Loop dofs and call evaluate_basis_derivatives.
1482
1908
    for (unsigned int r = 0; r < 6; r++)
1483
1909
    {
1484
 
      evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
 
1910
      evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1485
1911
      for (unsigned int s = 0; s < num_derivatives; s++)
1486
1912
      {
1487
1913
        values[r*num_derivatives + s] = dof_values[s];
1493
1919
  }
1494
1920
 
1495
1921
  /// Evaluate linear functional for dof i on the function f
1496
 
  virtual double evaluate_dof(unsigned int i,
 
1922
  virtual double evaluate_dof(std::size_t i,
1497
1923
                              const ufc::function& f,
 
1924
                              const double* vertex_coordinates,
 
1925
                              int cell_orientation,
1498
1926
                              const ufc::cell& c) const
1499
1927
  {
1500
 
    // Declare variables for result of evaluation.
 
1928
    // Declare variables for result of evaluation
1501
1929
    double vals[1];
1502
1930
    
1503
 
    // Declare variable for physical coordinates.
 
1931
    // Declare variable for physical coordinates
1504
1932
    double y[2];
1505
 
    const double * const * x = c.coordinates;
1506
1933
    switch (i)
1507
1934
    {
1508
1935
    case 0:
1509
1936
      {
1510
 
        y[0] = x[0][0];
1511
 
      y[1] = x[0][1];
 
1937
        y[0] = vertex_coordinates[0];
 
1938
      y[1] = vertex_coordinates[1];
1512
1939
      f.evaluate(vals, y, c);
1513
1940
      return vals[0];
1514
1941
        break;
1515
1942
      }
1516
1943
    case 1:
1517
1944
      {
1518
 
        y[0] = x[1][0];
1519
 
      y[1] = x[1][1];
 
1945
        y[0] = vertex_coordinates[2];
 
1946
      y[1] = vertex_coordinates[3];
1520
1947
      f.evaluate(vals, y, c);
1521
1948
      return vals[0];
1522
1949
        break;
1523
1950
      }
1524
1951
    case 2:
1525
1952
      {
1526
 
        y[0] = x[2][0];
1527
 
      y[1] = x[2][1];
 
1953
        y[0] = vertex_coordinates[4];
 
1954
      y[1] = vertex_coordinates[5];
1528
1955
      f.evaluate(vals, y, c);
1529
1956
      return vals[0];
1530
1957
        break;
1531
1958
      }
1532
1959
    case 3:
1533
1960
      {
1534
 
        y[0] = 0.5*x[1][0] + 0.5*x[2][0];
1535
 
      y[1] = 0.5*x[1][1] + 0.5*x[2][1];
 
1961
        y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
 
1962
      y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1536
1963
      f.evaluate(vals, y, c);
1537
1964
      return vals[0];
1538
1965
        break;
1539
1966
      }
1540
1967
    case 4:
1541
1968
      {
1542
 
        y[0] = 0.5*x[0][0] + 0.5*x[2][0];
1543
 
      y[1] = 0.5*x[0][1] + 0.5*x[2][1];
 
1969
        y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
 
1970
      y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1544
1971
      f.evaluate(vals, y, c);
1545
1972
      return vals[0];
1546
1973
        break;
1547
1974
      }
1548
1975
    case 5:
1549
1976
      {
1550
 
        y[0] = 0.5*x[0][0] + 0.5*x[1][0];
1551
 
      y[1] = 0.5*x[0][1] + 0.5*x[1][1];
 
1977
        y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
 
1978
      y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1552
1979
      f.evaluate(vals, y, c);
1553
1980
      return vals[0];
1554
1981
        break;
1561
1988
  /// Evaluate linear functionals for all dofs on the function f
1562
1989
  virtual void evaluate_dofs(double* values,
1563
1990
                             const ufc::function& f,
 
1991
                             const double* vertex_coordinates,
 
1992
                             int cell_orientation,
1564
1993
                             const ufc::cell& c) const
1565
1994
  {
1566
 
    // Declare variables for result of evaluation.
 
1995
    // Declare variables for result of evaluation
1567
1996
    double vals[1];
1568
1997
    
1569
 
    // Declare variable for physical coordinates.
 
1998
    // Declare variable for physical coordinates
1570
1999
    double y[2];
1571
 
    const double * const * x = c.coordinates;
1572
 
    y[0] = x[0][0];
1573
 
    y[1] = x[0][1];
 
2000
    y[0] = vertex_coordinates[0];
 
2001
    y[1] = vertex_coordinates[1];
1574
2002
    f.evaluate(vals, y, c);
1575
2003
    values[0] = vals[0];
1576
 
    y[0] = x[1][0];
1577
 
    y[1] = x[1][1];
 
2004
    y[0] = vertex_coordinates[2];
 
2005
    y[1] = vertex_coordinates[3];
1578
2006
    f.evaluate(vals, y, c);
1579
2007
    values[1] = vals[0];
1580
 
    y[0] = x[2][0];
1581
 
    y[1] = x[2][1];
 
2008
    y[0] = vertex_coordinates[4];
 
2009
    y[1] = vertex_coordinates[5];
1582
2010
    f.evaluate(vals, y, c);
1583
2011
    values[2] = vals[0];
1584
 
    y[0] = 0.5*x[1][0] + 0.5*x[2][0];
1585
 
    y[1] = 0.5*x[1][1] + 0.5*x[2][1];
 
2012
    y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
 
2013
    y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1586
2014
    f.evaluate(vals, y, c);
1587
2015
    values[3] = vals[0];
1588
 
    y[0] = 0.5*x[0][0] + 0.5*x[2][0];
1589
 
    y[1] = 0.5*x[0][1] + 0.5*x[2][1];
 
2016
    y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
 
2017
    y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1590
2018
    f.evaluate(vals, y, c);
1591
2019
    values[4] = vals[0];
1592
 
    y[0] = 0.5*x[0][0] + 0.5*x[1][0];
1593
 
    y[1] = 0.5*x[0][1] + 0.5*x[1][1];
 
2020
    y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
 
2021
    y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1594
2022
    f.evaluate(vals, y, c);
1595
2023
    values[5] = vals[0];
1596
2024
  }
1598
2026
  /// Interpolate vertex values from dof values
1599
2027
  virtual void interpolate_vertex_values(double* vertex_values,
1600
2028
                                         const double* dof_values,
 
2029
                                         const double* vertex_coordinates,
 
2030
                                         int cell_orientation,
1601
2031
                                         const ufc::cell& c) const
1602
2032
  {
1603
2033
    // Evaluate function and change variables
1611
2041
                                       const double* xhat,
1612
2042
                                       const ufc::cell& c) const
1613
2043
  {
1614
 
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
2044
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1615
2045
  }
1616
2046
 
1617
2047
  /// Map from coordinate x in cell to coordinate xhat in reference cell
1619
2049
                                     const double* x,
1620
2050
                                     const ufc::cell& c) const
1621
2051
  {
1622
 
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
 
2052
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1623
2053
  }
1624
2054
 
1625
2055
  /// Return the number of sub elements (for a mixed element)
1626
 
  virtual unsigned int num_sub_elements() const
 
2056
  virtual std::size_t num_sub_elements() const
1627
2057
  {
1628
2058
    return 0;
1629
2059
  }
1630
2060
 
1631
2061
  /// Create a new finite element for sub element i (for a mixed element)
1632
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
2062
  virtual ufc::finite_element* create_sub_element(std::size_t i) const
1633
2063
  {
1634
2064
    return 0;
1635
2065
  }
1637
2067
  /// Create a new class instance
1638
2068
  virtual ufc::finite_element* create() const
1639
2069
  {
1640
 
    return new conditional_finite_element_0();
 
2070
    return new conditional_finite_element_1();
1641
2071
  }
1642
2072
 
1643
2073
};
1647
2077
 
1648
2078
class conditional_dofmap_0: public ufc::dofmap
1649
2079
{
1650
 
private:
1651
 
 
1652
 
  unsigned int _global_dimension;
1653
2080
public:
1654
2081
 
1655
2082
  /// Constructor
1656
2083
  conditional_dofmap_0() : ufc::dofmap()
1657
2084
  {
1658
 
    _global_dimension = 0;
 
2085
    // Do nothing
1659
2086
  }
1660
2087
 
1661
2088
  /// Destructor
1667
2094
  /// Return a string identifying the dofmap
1668
2095
  virtual const char* signature() const
1669
2096
  {
1670
 
    return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', Space(2)), 2, None)";
1671
 
  }
1672
 
 
1673
 
  /// Return true iff mesh entities of topological dimension d are needed
1674
 
  virtual bool needs_mesh_entities(unsigned int d) const
1675
 
  {
1676
 
    switch (d)
1677
 
    {
1678
 
    case 0:
1679
 
      {
1680
 
        return true;
1681
 
        break;
1682
 
      }
1683
 
    case 1:
1684
 
      {
1685
 
        return true;
1686
 
        break;
1687
 
      }
1688
 
    case 2:
1689
 
      {
1690
 
        return false;
1691
 
        break;
1692
 
      }
1693
 
    }
1694
 
    
1695
 
    return false;
1696
 
  }
1697
 
 
1698
 
  /// Initialize dofmap for mesh (return true iff init_cell() is needed)
1699
 
  virtual bool init_mesh(const ufc::mesh& m)
1700
 
  {
1701
 
    _global_dimension = m.num_entities[0] + m.num_entities[1];
1702
 
    return false;
1703
 
  }
1704
 
 
1705
 
  /// Initialize dofmap for given cell
1706
 
  virtual void init_cell(const ufc::mesh& m,
1707
 
                         const ufc::cell& c)
1708
 
  {
1709
 
    // Do nothing
1710
 
  }
1711
 
 
1712
 
  /// Finish initialization of dofmap for cells
1713
 
  virtual void init_cell_finalize()
1714
 
  {
1715
 
    // Do nothing
1716
 
  }
1717
 
 
1718
 
  /// Return the topological dimension of the associated cell shape
1719
 
  virtual unsigned int topological_dimension() const
1720
 
  {
1721
 
    return 2;
1722
 
  }
1723
 
 
1724
 
  /// Return the geometric dimension of the associated cell shape
1725
 
  virtual unsigned int geometric_dimension() const
1726
 
  {
1727
 
    return 2;
1728
 
  }
1729
 
 
1730
 
  /// Return the dimension of the global finite element function space
1731
 
  virtual unsigned int global_dimension() const
1732
 
  {
1733
 
    return _global_dimension;
1734
 
  }
1735
 
 
1736
 
  /// Return the dimension of the local finite element function space for a cell
1737
 
  virtual unsigned int local_dimension(const ufc::cell& c) const
1738
 
  {
1739
 
    return 6;
1740
 
  }
1741
 
 
1742
 
  /// Return the maximum dimension of the local finite element function space
1743
 
  virtual unsigned int max_local_dimension() const
1744
 
  {
1745
 
    return 6;
1746
 
  }
1747
 
 
1748
 
  /// Return the number of dofs on each cell facet
1749
 
  virtual unsigned int num_facet_dofs() const
 
2097
    return "FFC dofmap for FiniteElement('Real', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
 
2098
  }
 
2099
 
 
2100
  /// Return true iff mesh entities of topological dimension d are needed
 
2101
  virtual bool needs_mesh_entities(std::size_t d) const
 
2102
  {
 
2103
    switch (d)
 
2104
    {
 
2105
    case 0:
 
2106
      {
 
2107
        return false;
 
2108
        break;
 
2109
      }
 
2110
    case 1:
 
2111
      {
 
2112
        return false;
 
2113
        break;
 
2114
      }
 
2115
    case 2:
 
2116
      {
 
2117
        return false;
 
2118
        break;
 
2119
      }
 
2120
    }
 
2121
    
 
2122
    return false;
 
2123
  }
 
2124
 
 
2125
  /// Return the topological dimension of the associated cell shape
 
2126
  virtual std::size_t topological_dimension() const
 
2127
  {
 
2128
    return 2;
 
2129
  }
 
2130
 
 
2131
  /// Return the geometric dimension of the associated cell shape
 
2132
  virtual std::size_t geometric_dimension() const
 
2133
  {
 
2134
    return 2;
 
2135
  }
 
2136
 
 
2137
  /// Return the dimension of the global finite element function space
 
2138
  virtual std::size_t global_dimension(const std::vector<std::size_t>&
 
2139
                                       num_global_entities) const
 
2140
  {
 
2141
    return 1;
 
2142
  }
 
2143
 
 
2144
  /// Return the dimension of the local finite element function space for a cell
 
2145
  virtual std::size_t local_dimension(const ufc::cell& c) const
 
2146
  {
 
2147
    return 1;
 
2148
  }
 
2149
 
 
2150
  /// Return the maximum dimension of the local finite element function space
 
2151
  virtual std::size_t max_local_dimension() const
 
2152
  {
 
2153
    return 1;
 
2154
  }
 
2155
 
 
2156
  /// Return the number of dofs on each cell facet
 
2157
  virtual std::size_t num_facet_dofs() const
 
2158
  {
 
2159
    return 0;
 
2160
  }
 
2161
 
 
2162
  /// Return the number of dofs associated with each cell entity of dimension d
 
2163
  virtual std::size_t num_entity_dofs(std::size_t d) const
 
2164
  {
 
2165
    switch (d)
 
2166
    {
 
2167
    case 0:
 
2168
      {
 
2169
        return 0;
 
2170
        break;
 
2171
      }
 
2172
    case 1:
 
2173
      {
 
2174
        return 0;
 
2175
        break;
 
2176
      }
 
2177
    case 2:
 
2178
      {
 
2179
        return 1;
 
2180
        break;
 
2181
      }
 
2182
    }
 
2183
    
 
2184
    return 0;
 
2185
  }
 
2186
 
 
2187
  /// Tabulate the local-to-global mapping of dofs on a cell
 
2188
  virtual void tabulate_dofs(std::size_t* dofs,
 
2189
                             const std::vector<std::size_t>& num_global_entities,
 
2190
                             const ufc::cell& c) const
 
2191
  {
 
2192
    dofs[0] = 0;
 
2193
  }
 
2194
 
 
2195
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
2196
  virtual void tabulate_facet_dofs(std::size_t* dofs,
 
2197
                                   std::size_t facet) const
 
2198
  {
 
2199
    switch (facet)
 
2200
    {
 
2201
    case 0:
 
2202
      {
 
2203
        
 
2204
        break;
 
2205
      }
 
2206
    case 1:
 
2207
      {
 
2208
        
 
2209
        break;
 
2210
      }
 
2211
    case 2:
 
2212
      {
 
2213
        
 
2214
        break;
 
2215
      }
 
2216
    }
 
2217
    
 
2218
  }
 
2219
 
 
2220
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
2221
  virtual void tabulate_entity_dofs(std::size_t* dofs,
 
2222
                                    std::size_t d, std::size_t i) const
 
2223
  {
 
2224
    if (d > 2)
 
2225
    {
 
2226
    std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
 
2227
    }
 
2228
    
 
2229
    switch (d)
 
2230
    {
 
2231
    case 0:
 
2232
      {
 
2233
        
 
2234
        break;
 
2235
      }
 
2236
    case 1:
 
2237
      {
 
2238
        
 
2239
        break;
 
2240
      }
 
2241
    case 2:
 
2242
      {
 
2243
        if (i > 0)
 
2244
      {
 
2245
      std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
 
2246
      }
 
2247
      
 
2248
      dofs[0] = 0;
 
2249
        break;
 
2250
      }
 
2251
    }
 
2252
    
 
2253
  }
 
2254
 
 
2255
  /// Tabulate the coordinates of all dofs on a cell
 
2256
  virtual void tabulate_coordinates(double** dof_coordinates,
 
2257
                                    const double* vertex_coordinates) const
 
2258
  {
 
2259
    dof_coordinates[0][0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
 
2260
    dof_coordinates[0][1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
 
2261
  }
 
2262
 
 
2263
  /// Return the number of sub dofmaps (for a mixed element)
 
2264
  virtual std::size_t num_sub_dofmaps() const
 
2265
  {
 
2266
    return 0;
 
2267
  }
 
2268
 
 
2269
  /// Create a new dofmap for sub dofmap i (for a mixed element)
 
2270
  virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
 
2271
  {
 
2272
    return 0;
 
2273
  }
 
2274
 
 
2275
  /// Create a new class instance
 
2276
  virtual ufc::dofmap* create() const
 
2277
  {
 
2278
    return new conditional_dofmap_0();
 
2279
  }
 
2280
 
 
2281
};
 
2282
 
 
2283
/// This class defines the interface for a local-to-global mapping of
 
2284
/// degrees of freedom (dofs).
 
2285
 
 
2286
class conditional_dofmap_1: public ufc::dofmap
 
2287
{
 
2288
public:
 
2289
 
 
2290
  /// Constructor
 
2291
  conditional_dofmap_1() : ufc::dofmap()
 
2292
  {
 
2293
    // Do nothing
 
2294
  }
 
2295
 
 
2296
  /// Destructor
 
2297
  virtual ~conditional_dofmap_1()
 
2298
  {
 
2299
    // Do nothing
 
2300
  }
 
2301
 
 
2302
  /// Return a string identifying the dofmap
 
2303
  virtual const char* signature() const
 
2304
  {
 
2305
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
 
2306
  }
 
2307
 
 
2308
  /// Return true iff mesh entities of topological dimension d are needed
 
2309
  virtual bool needs_mesh_entities(std::size_t d) const
 
2310
  {
 
2311
    switch (d)
 
2312
    {
 
2313
    case 0:
 
2314
      {
 
2315
        return true;
 
2316
        break;
 
2317
      }
 
2318
    case 1:
 
2319
      {
 
2320
        return true;
 
2321
        break;
 
2322
      }
 
2323
    case 2:
 
2324
      {
 
2325
        return false;
 
2326
        break;
 
2327
      }
 
2328
    }
 
2329
    
 
2330
    return false;
 
2331
  }
 
2332
 
 
2333
  /// Return the topological dimension of the associated cell shape
 
2334
  virtual std::size_t topological_dimension() const
 
2335
  {
 
2336
    return 2;
 
2337
  }
 
2338
 
 
2339
  /// Return the geometric dimension of the associated cell shape
 
2340
  virtual std::size_t geometric_dimension() const
 
2341
  {
 
2342
    return 2;
 
2343
  }
 
2344
 
 
2345
  /// Return the dimension of the global finite element function space
 
2346
  virtual std::size_t global_dimension(const std::vector<std::size_t>&
 
2347
                                       num_global_entities) const
 
2348
  {
 
2349
    return num_global_entities[0] + num_global_entities[1];
 
2350
  }
 
2351
 
 
2352
  /// Return the dimension of the local finite element function space for a cell
 
2353
  virtual std::size_t local_dimension(const ufc::cell& c) const
 
2354
  {
 
2355
    return 6;
 
2356
  }
 
2357
 
 
2358
  /// Return the maximum dimension of the local finite element function space
 
2359
  virtual std::size_t max_local_dimension() const
 
2360
  {
 
2361
    return 6;
 
2362
  }
 
2363
 
 
2364
  /// Return the number of dofs on each cell facet
 
2365
  virtual std::size_t num_facet_dofs() const
1750
2366
  {
1751
2367
    return 3;
1752
2368
  }
1753
2369
 
1754
2370
  /// Return the number of dofs associated with each cell entity of dimension d
1755
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
2371
  virtual std::size_t num_entity_dofs(std::size_t d) const
1756
2372
  {
1757
2373
    switch (d)
1758
2374
    {
1777
2393
  }
1778
2394
 
1779
2395
  /// Tabulate the local-to-global mapping of dofs on a cell
1780
 
  virtual void tabulate_dofs(unsigned int* dofs,
1781
 
                             const ufc::mesh& m,
 
2396
  virtual void tabulate_dofs(std::size_t* dofs,
 
2397
                             const std::vector<std::size_t>& num_global_entities,
1782
2398
                             const ufc::cell& c) const
1783
2399
  {
1784
2400
    unsigned int offset = 0;
1785
2401
    dofs[0] = offset + c.entity_indices[0][0];
1786
2402
    dofs[1] = offset + c.entity_indices[0][1];
1787
2403
    dofs[2] = offset + c.entity_indices[0][2];
1788
 
    offset += m.num_entities[0];
 
2404
    offset += num_global_entities[0];
1789
2405
    dofs[3] = offset + c.entity_indices[1][0];
1790
2406
    dofs[4] = offset + c.entity_indices[1][1];
1791
2407
    dofs[5] = offset + c.entity_indices[1][2];
1792
 
    offset += m.num_entities[1];
 
2408
    offset += num_global_entities[1];
1793
2409
  }
1794
2410
 
1795
2411
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
1796
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
1797
 
                                   unsigned int facet) const
 
2412
  virtual void tabulate_facet_dofs(std::size_t* dofs,
 
2413
                                   std::size_t facet) const
1798
2414
  {
1799
2415
    switch (facet)
1800
2416
    {
1824
2440
  }
1825
2441
 
1826
2442
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
1827
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
1828
 
                                    unsigned int d, unsigned int i) const
 
2443
  virtual void tabulate_entity_dofs(std::size_t* dofs,
 
2444
                                    std::size_t d, std::size_t i) const
1829
2445
  {
1830
2446
    if (d > 2)
1831
2447
    {
1900
2516
  }
1901
2517
 
1902
2518
  /// Tabulate the coordinates of all dofs on a cell
1903
 
  virtual void tabulate_coordinates(double** coordinates,
1904
 
                                    const ufc::cell& c) const
 
2519
  virtual void tabulate_coordinates(double** dof_coordinates,
 
2520
                                    const double* vertex_coordinates) const
1905
2521
  {
1906
 
    const double * const * x = c.coordinates;
1907
 
    
1908
 
    coordinates[0][0] = x[0][0];
1909
 
    coordinates[0][1] = x[0][1];
1910
 
    coordinates[1][0] = x[1][0];
1911
 
    coordinates[1][1] = x[1][1];
1912
 
    coordinates[2][0] = x[2][0];
1913
 
    coordinates[2][1] = x[2][1];
1914
 
    coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1915
 
    coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1916
 
    coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1917
 
    coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1918
 
    coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1919
 
    coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
 
2522
    dof_coordinates[0][0] = vertex_coordinates[0];
 
2523
    dof_coordinates[0][1] = vertex_coordinates[1];
 
2524
    dof_coordinates[1][0] = vertex_coordinates[2];
 
2525
    dof_coordinates[1][1] = vertex_coordinates[3];
 
2526
    dof_coordinates[2][0] = vertex_coordinates[4];
 
2527
    dof_coordinates[2][1] = vertex_coordinates[5];
 
2528
    dof_coordinates[3][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
 
2529
    dof_coordinates[3][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
 
2530
    dof_coordinates[4][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
 
2531
    dof_coordinates[4][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
 
2532
    dof_coordinates[5][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
 
2533
    dof_coordinates[5][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1920
2534
  }
1921
2535
 
1922
2536
  /// Return the number of sub dofmaps (for a mixed element)
1923
 
  virtual unsigned int num_sub_dofmaps() const
 
2537
  virtual std::size_t num_sub_dofmaps() const
1924
2538
  {
1925
2539
    return 0;
1926
2540
  }
1927
2541
 
1928
2542
  /// Create a new dofmap for sub dofmap i (for a mixed element)
1929
 
  virtual ufc::dofmap* create_sub_dofmap(unsigned int i) const
 
2543
  virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
1930
2544
  {
1931
2545
    return 0;
1932
2546
  }
1934
2548
  /// Create a new class instance
1935
2549
  virtual ufc::dofmap* create() const
1936
2550
  {
1937
 
    return new conditional_dofmap_0();
 
2551
    return new conditional_dofmap_1();
1938
2552
  }
1939
2553
 
1940
2554
};
1943
2557
/// tensor corresponding to the local contribution to a form from
1944
2558
/// the integral over a cell.
1945
2559
 
1946
 
class conditional_cell_integral_0_0: public ufc::cell_integral
 
2560
class conditional_cell_integral_0_otherwise: public ufc::cell_integral
1947
2561
{
1948
2562
public:
1949
2563
 
1950
2564
  /// Constructor
1951
 
  conditional_cell_integral_0_0() : ufc::cell_integral()
 
2565
  conditional_cell_integral_0_otherwise() : ufc::cell_integral()
1952
2566
  {
1953
2567
    // Do nothing
1954
2568
  }
1955
2569
 
1956
2570
  /// Destructor
1957
 
  virtual ~conditional_cell_integral_0_0()
 
2571
  virtual ~conditional_cell_integral_0_otherwise()
1958
2572
  {
1959
2573
    // Do nothing
1960
2574
  }
1962
2576
  /// Tabulate the tensor for the contribution from a local cell
1963
2577
  virtual void tabulate_tensor(double* A,
1964
2578
                               const double * const * w,
1965
 
                               const ufc::cell& c) const
 
2579
                               const double* vertex_coordinates,
 
2580
                               int cell_orientation) const
1966
2581
  {
1967
 
    // Extract vertex coordinates
1968
 
    const double * const * x = c.coordinates;
1969
 
    
1970
 
    // Compute Jacobian of affine map from reference cell
1971
 
    const double J_00 = x[1][0] - x[0][0];
1972
 
    const double J_01 = x[2][0] - x[0][0];
1973
 
    const double J_10 = x[1][1] - x[0][1];
1974
 
    const double J_11 = x[2][1] - x[0][1];
1975
 
    
1976
 
    // Compute determinant of Jacobian
1977
 
    double detJ = J_00*J_11 - J_01*J_10;
1978
 
    
1979
 
    // Compute inverse of Jacobian
 
2582
    // Compute Jacobian
 
2583
    double J[4];
 
2584
    compute_jacobian_triangle_2d(J, vertex_coordinates);
 
2585
    
 
2586
    // Compute Jacobian inverse and determinant
 
2587
    double K[4];
 
2588
    double detJ;
 
2589
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
1980
2590
    
1981
2591
    // Set scale factor
1982
2592
    const double det = std::abs(detJ);
1983
2593
    
1984
 
    // Cell Volume.
1985
 
    
1986
 
    // Compute circumradius, assuming triangle is embedded in 2D.
1987
 
    
1988
 
    
1989
 
    // Facet Area.
 
2594
    // Cell volume
 
2595
    
 
2596
    // Compute circumradius of triangle in 2D
 
2597
    
 
2598
    
 
2599
    // Facet area
1990
2600
    
1991
2601
    // Array of quadrature weights.
1992
 
    static const double W6[6] = {0.054975872, 0.054975872, 0.054975872, 0.11169079, 0.11169079, 0.11169079};
1993
 
    // Quadrature points on the UFC reference element: (0.81684757, 0.091576214), (0.091576214, 0.81684757), (0.091576214, 0.091576214), (0.10810302, 0.44594849), (0.44594849, 0.10810302), (0.44594849, 0.44594849)
 
2602
    static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
 
2603
    // Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
1994
2604
    
1995
2605
    // Value of basis functions at quadrature points.
1996
 
    static const double FE0[6][6] = \
1997
 
    {{-0.074803808, 0.51763234, -0.074803808, 0.29921523, 0.033544812, 0.29921523},
1998
 
    {-0.074803808, -0.074803808, 0.51763234, 0.29921523, 0.29921523, 0.033544812},
1999
 
    {0.51763234, -0.074803808, -0.074803808, 0.033544812, 0.29921523, 0.29921523},
2000
 
    {-0.048208378, -0.084730493, -0.048208378, 0.19283351, 0.79548023, 0.19283351},
2001
 
    {-0.048208378, -0.048208378, -0.084730493, 0.19283351, 0.19283351, 0.79548023},
2002
 
    {-0.084730493, -0.048208378, -0.048208378, 0.79548023, 0.19283351, 0.19283351}};
 
2606
    static const double FE0[3][6] = \
 
2607
    {{0.22222222, -0.11111111, -0.11111111, 0.11111111, 0.44444444, 0.44444444},
 
2608
    {-0.11111111, -0.11111111, 0.22222222, 0.44444444, 0.44444444, 0.11111111},
 
2609
    {-0.11111111, 0.22222222, -0.11111111, 0.44444444, 0.11111111, 0.44444444}};
2003
2610
    
2004
 
    static const double FEA6_f0[6][3] = \
2005
 
    {{0.091576214, 0.81684757, 0.091576214},
2006
 
    {0.091576214, 0.091576214, 0.81684757},
2007
 
    {0.81684757, 0.091576214, 0.091576214},
2008
 
    {0.44594849, 0.10810302, 0.44594849},
2009
 
    {0.44594849, 0.44594849, 0.10810302},
2010
 
    {0.10810302, 0.44594849, 0.44594849}};
 
2611
    static const double FEA3_f0[3][3] = \
 
2612
    {{0.66666667, 0.16666667, 0.16666667},
 
2613
    {0.16666667, 0.16666667, 0.66666667},
 
2614
    {0.16666667, 0.66666667, 0.16666667}};
2011
2615
    
2012
2616
    // Reset values in the element tensor.
2013
2617
    for (unsigned int r = 0; r < 6; r++)
2021
2625
    // Loop quadrature points for integral.
2022
2626
    
2023
2627
    // Declare array to hold physical coordinate of quadrature point.
2024
 
    double X6[2];
2025
 
    // Number of operations to compute element tensor for following IP loop = 318
2026
 
    for (unsigned int ip = 0; ip < 6; ip++)
 
2628
    double X3[2];
 
2629
    // Number of operations to compute element tensor for following IP loop = 168
 
2630
    for (unsigned int ip = 0; ip < 3; ip++)
2027
2631
    {
2028
2632
      
2029
2633
      // Compute physical coordinate of quadrature point, operations: 10.
2030
 
      X6[0] = FEA6_f0[ip][0]*x[0][0] + FEA6_f0[ip][1]*x[1][0] + FEA6_f0[ip][2]*x[2][0];
2031
 
      X6[1] = FEA6_f0[ip][0]*x[0][1] + FEA6_f0[ip][1]*x[1][1] + FEA6_f0[ip][2]*x[2][1];
2032
 
      double C[3];
 
2634
      X3[0] = FEA3_f0[ip][0]*vertex_coordinates[0] +                  FEA3_f0[ip][1]*vertex_coordinates[2] + FEA3_f0[ip][2]*vertex_coordinates[4];
 
2635
      X3[1] = FEA3_f0[ip][0]*vertex_coordinates[1] +                  FEA3_f0[ip][1]*vertex_coordinates[3] + FEA3_f0[ip][2]*vertex_coordinates[5];
 
2636
      double C[4];
 
2637
      // Compute conditional, operations: 2.
 
2638
      C[0] = ((1.0 > 0.0)) ? w[0][0] : (1.0 + w[0][0]);
2033
2639
      // Compute conditional, operations: 12.
2034
 
      C[0] = ((0.5 + X6[1] - X6[0]) >= 0.0 && X6[0] >= 0.55 && X6[0] <= 0.95 && !(X6[1] < 0.05 || X6[1] > 0.45)) ? -1.0 : 0.0;
2035
 
      // Compute conditional, operations: 8.
2036
 
      C[1] = (((X6[0]-0.33)*(X6[0]-0.33) + (X6[1]-0.67)*(X6[1]-0.67)) <= 0.015) ? -1.0 : 5.0;
2037
 
      // Compute conditional, operations: 8.
2038
 
      C[2] = (((X6[0]-0.33)*(X6[0]-0.33) + (X6[1]-0.67)*(X6[1]-0.67)) <= 0.025) ? C[1] : 0.0;
 
2640
      C[1] = ((((0.5 + X3[1] - X3[0]) >= 0.0) && (((X3[0] >= 0.55) && (X3[0] <= 0.95)) && !(((X3[1] < 0.05) || (X3[1] > 0.45)))))) ? -1.0 : 0.0;
 
2641
      // Compute conditional, operations: 8.
 
2642
      C[2] = ((((X3[0]-0.33)*(X3[0]-0.33) + (X3[1]-0.67)*(X3[1]-0.67)) <= 0.015)) ? -1.0 : 5.0;
 
2643
      // Compute conditional, operations: 8.
 
2644
      C[3] = ((((X3[0]-0.33)*(X3[0]-0.33) + (X3[1]-0.67)*(X3[1]-0.67)) <= 0.025)) ? C[2] : 0.0;
2039
2645
      
2040
 
      // Number of operations to compute ip constants: 3
 
2646
      // Number of operations to compute ip constants: 4
2041
2647
      double I[1];
2042
 
      // Number of operations: 3
2043
 
      I[0] = W6[ip]*det*(C[0] + C[2]);
 
2648
      // Number of operations: 4
 
2649
      I[0] = W3[ip]*det*(C[0] + C[1] + C[3]);
2044
2650
      
2045
2651
      
2046
2652
      // Number of operations for primary indices: 12
2052
2658
    }// end loop over 'ip'
2053
2659
  }
2054
2660
 
2055
 
  /// Tabulate the tensor for the contribution from a local cell
2056
 
  /// using the specified reference cell quadrature points/weights
2057
 
  virtual void tabulate_tensor(double* A,
2058
 
                               const double * const * w,
2059
 
                               const ufc::cell& c,
2060
 
                               unsigned int num_quadrature_points,
2061
 
                               const double * const * quadrature_points,
2062
 
                               const double* quadrature_weights) const
2063
 
  {
2064
 
    std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
2065
 
  }
2066
 
 
2067
2661
};
2068
2662
 
2069
2663
/// This class defines the interface for the assembly of the global
2100
2694
  /// Return a string identifying the form
2101
2695
  virtual const char* signature() const
2102
2696
  {
2103
 
    return "Form([Integral(Product(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)), 2, None), 0), Sum(Conditional(AndCondition(GE(Sum(FloatValue(0.55, (), (), {}), Sum(FloatValue(-0.05, (), (), {}), Sum(Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(1),), {})), Product(IntValue(-1, (), (), {}), Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(0),), {})))))), Zero((), (), {})), AndCondition(AndCondition(GE(Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(0),), {})), FloatValue(0.55, (), (), {})), LE(Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(0),), {})), FloatValue(0.95, (), (), {}))), NotCondition(OrCondition(LT(Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(1),), {})), FloatValue(0.05, (), (), {})), GT(Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(1),), {})), FloatValue(0.45, (), (), {})))))), FloatValue(-1, (), (), {}), Zero((), (), {})), Conditional(LE(Sum(Power(Sum(FloatValue(-0.33, (), (), {}), Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(0),), {}))), IntValue(2, (), (), {})), Power(Sum(FloatValue(-0.67, (), (), {}), Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(1),), {}))), IntValue(2, (), (), {}))), FloatValue(0.025, (), (), {})), Conditional(LE(Sum(Power(Sum(FloatValue(-0.33, (), (), {}), Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(0),), {}))), IntValue(2, (), (), {})), Power(Sum(FloatValue(-0.67, (), (), {}), Indexed(SpatialCoordinate(Cell('triangle', Space(2))), MultiIndex((FixedIndex(1),), {}))), IntValue(2, (), (), {}))), FloatValue(0.015, (), (), {})), FloatValue(-1, (), (), {}), FloatValue(5, (), (), {})), Zero((), (), {})))), Measure('cell', 0, None))])";
 
2697
    return "d58fed1bc260f7f4defd7a416ccd3bfaa17fe02fb2f4fe3afd9221e7afab62ff17ceb9ccd1ed8d48f90cf118e951b44c22a33afcde5104ad00188a048194ce7e";
2104
2698
  }
2105
2699
 
2106
2700
  /// Return the rank of the global tensor (r)
2107
 
  virtual unsigned int rank() const
 
2701
  virtual std::size_t rank() const
2108
2702
  {
2109
2703
    return 1;
2110
2704
  }
2111
2705
 
2112
2706
  /// Return the number of coefficients (n)
2113
 
  virtual unsigned int num_coefficients() const
 
2707
  virtual std::size_t num_coefficients() const
2114
2708
  {
2115
 
    return 0;
 
2709
    return 1;
2116
2710
  }
2117
2711
 
2118
2712
  /// Return the number of cell domains
2119
 
  virtual unsigned int num_cell_domains() const
 
2713
  virtual std::size_t num_cell_domains() const
2120
2714
  {
2121
 
    return 1;
 
2715
    return 0;
2122
2716
  }
2123
2717
 
2124
2718
  /// Return the number of exterior facet domains
2125
 
  virtual unsigned int num_exterior_facet_domains() const
 
2719
  virtual std::size_t num_exterior_facet_domains() const
2126
2720
  {
2127
2721
    return 0;
2128
2722
  }
2129
2723
 
2130
2724
  /// Return the number of interior facet domains
2131
 
  virtual unsigned int num_interior_facet_domains() const
2132
 
  {
2133
 
    return 0;
 
2725
  virtual std::size_t num_interior_facet_domains() const
 
2726
  {
 
2727
    return 0;
 
2728
  }
 
2729
 
 
2730
  /// Return the number of point domains
 
2731
  virtual std::size_t num_point_domains() const
 
2732
  {
 
2733
    return 0;
 
2734
  }
 
2735
 
 
2736
  /// Return whether the form has any cell integrals
 
2737
  virtual bool has_cell_integrals() const
 
2738
  {
 
2739
    return true;
 
2740
  }
 
2741
 
 
2742
  /// Return whether the form has any exterior facet integrals
 
2743
  virtual bool has_exterior_facet_integrals() const
 
2744
  {
 
2745
    return false;
 
2746
  }
 
2747
 
 
2748
  /// Return whether the form has any interior facet integrals
 
2749
  virtual bool has_interior_facet_integrals() const
 
2750
  {
 
2751
    return false;
 
2752
  }
 
2753
 
 
2754
  /// Return whether the form has any point integrals
 
2755
  virtual bool has_point_integrals() const
 
2756
  {
 
2757
    return false;
2134
2758
  }
2135
2759
 
2136
2760
  /// Create a new finite element for argument function i
2137
 
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
2761
  virtual ufc::finite_element* create_finite_element(std::size_t i) const
2138
2762
  {
2139
2763
    switch (i)
2140
2764
    {
2141
2765
    case 0:
2142
2766
      {
 
2767
        return new conditional_finite_element_1();
 
2768
        break;
 
2769
      }
 
2770
    case 1:
 
2771
      {
2143
2772
        return new conditional_finite_element_0();
2144
2773
        break;
2145
2774
      }
2149
2778
  }
2150
2779
 
2151
2780
  /// Create a new dofmap for argument function i
2152
 
  virtual ufc::dofmap* create_dofmap(unsigned int i) const
 
2781
  virtual ufc::dofmap* create_dofmap(std::size_t i) const
2153
2782
  {
2154
2783
    switch (i)
2155
2784
    {
2156
2785
    case 0:
2157
2786
      {
 
2787
        return new conditional_dofmap_1();
 
2788
        break;
 
2789
      }
 
2790
    case 1:
 
2791
      {
2158
2792
        return new conditional_dofmap_0();
2159
2793
        break;
2160
2794
      }
2164
2798
  }
2165
2799
 
2166
2800
  /// Create a new cell integral on sub domain i
2167
 
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
2801
  virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
2168
2802
  {
2169
 
    switch (i)
2170
 
    {
2171
 
    case 0:
2172
 
      {
2173
 
        return new conditional_cell_integral_0_0();
2174
 
        break;
2175
 
      }
2176
 
    }
2177
 
    
2178
2803
    return 0;
2179
2804
  }
2180
2805
 
2181
2806
  /// Create a new exterior facet integral on sub domain i
2182
 
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
2807
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
2183
2808
  {
2184
2809
    return 0;
2185
2810
  }
2186
2811
 
2187
2812
  /// Create a new interior facet integral on sub domain i
2188
 
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
2813
  virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
 
2814
  {
 
2815
    return 0;
 
2816
  }
 
2817
 
 
2818
  /// Create a new point integral on sub domain i
 
2819
  virtual ufc::point_integral* create_point_integral(std::size_t i) const
 
2820
  {
 
2821
    return 0;
 
2822
  }
 
2823
 
 
2824
  /// Create a new cell integral on everywhere else
 
2825
  virtual ufc::cell_integral* create_default_cell_integral() const
 
2826
  {
 
2827
    return new conditional_cell_integral_0_otherwise();
 
2828
  }
 
2829
 
 
2830
  /// Create a new exterior facet integral on everywhere else
 
2831
  virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
 
2832
  {
 
2833
    return 0;
 
2834
  }
 
2835
 
 
2836
  /// Create a new interior facet integral on everywhere else
 
2837
  virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
 
2838
  {
 
2839
    return 0;
 
2840
  }
 
2841
 
 
2842
  /// Create a new point integral on everywhere else
 
2843
  virtual ufc::point_integral* create_default_point_integral() const
2189
2844
  {
2190
2845
    return 0;
2191
2846
  }