← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to demo/pde/lift-drag/cpp/Lift.h

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

Show diffs side-by-side

added added

removed removed

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