← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to demo/pde/nonlinear-poisson/cpp/NonlinearPoisson.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/nonlinear-poisson/cpp/NonlinearPoisson.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 __NONLINEARPOISSON_H
 
8
#define __NONLINEARPOISSON_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_NonlinearPoissonBilinearForm_finite_element_0: public ufc::finite_element
 
18
{
 
19
public:
 
20
 
 
21
  /// Constructor
 
22
  UFC_NonlinearPoissonBilinearForm_finite_element_0() : ufc::finite_element()
 
23
  {
 
24
    // Do nothing
 
25
  }
 
26
 
 
27
  /// Destructor
 
28
  virtual ~UFC_NonlinearPoissonBilinearForm_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_NonlinearPoissonBilinearForm_finite_element_0();
 
432
  }
 
433
 
 
434
};
 
435
 
 
436
/// This class defines the interface for a finite element.
 
437
 
 
438
class UFC_NonlinearPoissonBilinearForm_finite_element_1: public ufc::finite_element
 
439
{
 
440
public:
 
441
 
 
442
  /// Constructor
 
443
  UFC_NonlinearPoissonBilinearForm_finite_element_1() : ufc::finite_element()
 
444
  {
 
445
    // Do nothing
 
446
  }
 
447
 
 
448
  /// Destructor
 
449
  virtual ~UFC_NonlinearPoissonBilinearForm_finite_element_1()
 
450
  {
 
451
    // Do nothing
 
452
  }
 
453
 
 
454
  /// Return a string identifying the finite element
 
455
  virtual const char* signature() const
 
456
  {
 
457
    return "Lagrange finite element of degree 1 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 3;
 
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
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
528
    
 
529
    // Compute psitilde_a
 
530
    const double psitilde_a_0 = 1;
 
531
    const double psitilde_a_1 = x;
 
532
    
 
533
    // Compute psitilde_bs
 
534
    const double psitilde_bs_0_0 = 1;
 
535
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
536
    const double psitilde_bs_1_0 = 1;
 
537
    
 
538
    // Compute basisvalues
 
539
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
540
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
541
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
542
    
 
543
    // Table(s) of coefficients
 
544
    const static double coefficients0[3][3] = \
 
545
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
546
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
547
    {0.471404520791032, 0, 0.333333333333333}};
 
548
    
 
549
    // Extract relevant coefficients
 
550
    const double coeff0_0 = coefficients0[dof][0];
 
551
    const double coeff0_1 = coefficients0[dof][1];
 
552
    const double coeff0_2 = coefficients0[dof][2];
 
553
    
 
554
    // Compute value(s)
 
555
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
556
  }
 
557
 
 
558
  /// Evaluate all basis functions at given point in cell
 
559
  virtual void evaluate_basis_all(double* values,
 
560
                                  const double* coordinates,
 
561
                                  const ufc::cell& c) const
 
562
  {
 
563
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
564
  }
 
565
 
 
566
  /// Evaluate order n derivatives of basis function i at given point in cell
 
567
  virtual void evaluate_basis_derivatives(unsigned int i,
 
568
                                          unsigned int n,
 
569
                                          double* values,
 
570
                                          const double* coordinates,
 
571
                                          const ufc::cell& c) const
 
572
  {
 
573
    // Extract vertex coordinates
 
574
    const double * const * element_coordinates = c.coordinates;
 
575
    
 
576
    // Compute Jacobian of affine map from reference cell
 
577
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
578
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
579
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
580
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
581
      
 
582
    // Compute determinant of Jacobian
 
583
    const double detJ = J_00*J_11 - J_01*J_10;
 
584
    
 
585
    // Compute inverse of Jacobian
 
586
    
 
587
    // Get coordinates and map to the reference (UFC) element
 
588
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
589
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
590
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
591
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
592
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
593
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
594
    
 
595
    // Map coordinates to the reference square
 
596
    if (std::abs(y - 1.0) < 1e-14)
 
597
      x = -1.0;
 
598
    else
 
599
      x = 2.0 *x/(1.0 - y) - 1.0;
 
600
    y = 2.0*y - 1.0;
 
601
    
 
602
    // Compute number of derivatives
 
603
    unsigned int num_derivatives = 1;
 
604
    
 
605
    for (unsigned int j = 0; j < n; j++)
 
606
      num_derivatives *= 2;
 
607
    
 
608
    
 
609
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
610
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
611
        
 
612
    for (unsigned int j = 0; j < num_derivatives; j++)
 
613
    {
 
614
      combinations[j] = new unsigned int [n];
 
615
      for (unsigned int k = 0; k < n; k++)
 
616
        combinations[j][k] = 0;
 
617
    }
 
618
        
 
619
    // Generate combinations of derivatives
 
620
    for (unsigned int row = 1; row < num_derivatives; row++)
 
621
    {
 
622
      for (unsigned int num = 0; num < row; num++)
 
623
      {
 
624
        for (unsigned int col = n-1; col+1 > 0; col--)
 
625
        {
 
626
          if (combinations[row][col] + 1 > 1)
 
627
            combinations[row][col] = 0;
 
628
          else
 
629
          {
 
630
            combinations[row][col] += 1;
 
631
            break;
 
632
          }
 
633
        }
 
634
      }
 
635
    }
 
636
    
 
637
    // Compute inverse of Jacobian
 
638
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
639
    
 
640
    // Declare transformation matrix
 
641
    // Declare pointer to two dimensional array and initialise
 
642
    double **transform = new double *[num_derivatives];
 
643
        
 
644
    for (unsigned int j = 0; j < num_derivatives; j++)
 
645
    {
 
646
      transform[j] = new double [num_derivatives];
 
647
      for (unsigned int k = 0; k < num_derivatives; k++)
 
648
        transform[j][k] = 1;
 
649
    }
 
650
    
 
651
    // Construct transformation matrix
 
652
    for (unsigned int row = 0; row < num_derivatives; row++)
 
653
    {
 
654
      for (unsigned int col = 0; col < num_derivatives; col++)
 
655
      {
 
656
        for (unsigned int k = 0; k < n; k++)
 
657
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
658
      }
 
659
    }
 
660
    
 
661
    // Reset values
 
662
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
663
      values[j] = 0;
 
664
    
 
665
    // Map degree of freedom to element degree of freedom
 
666
    const unsigned int dof = i;
 
667
    
 
668
    // Generate scalings
 
669
    const double scalings_y_0 = 1;
 
670
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
671
    
 
672
    // Compute psitilde_a
 
673
    const double psitilde_a_0 = 1;
 
674
    const double psitilde_a_1 = x;
 
675
    
 
676
    // Compute psitilde_bs
 
677
    const double psitilde_bs_0_0 = 1;
 
678
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
679
    const double psitilde_bs_1_0 = 1;
 
680
    
 
681
    // Compute basisvalues
 
682
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
683
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
684
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
685
    
 
686
    // Table(s) of coefficients
 
687
    const static double coefficients0[3][3] = \
 
688
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
689
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
690
    {0.471404520791032, 0, 0.333333333333333}};
 
691
    
 
692
    // Interesting (new) part
 
693
    // Tables of derivatives of the polynomial base (transpose)
 
694
    const static double dmats0[3][3] = \
 
695
    {{0, 0, 0},
 
696
    {4.89897948556636, 0, 0},
 
697
    {0, 0, 0}};
 
698
    
 
699
    const static double dmats1[3][3] = \
 
700
    {{0, 0, 0},
 
701
    {2.44948974278318, 0, 0},
 
702
    {4.24264068711928, 0, 0}};
 
703
    
 
704
    // Compute reference derivatives
 
705
    // Declare pointer to array of derivatives on FIAT element
 
706
    double *derivatives = new double [num_derivatives];
 
707
    
 
708
    // Declare coefficients
 
709
    double coeff0_0 = 0;
 
710
    double coeff0_1 = 0;
 
711
    double coeff0_2 = 0;
 
712
    
 
713
    // Declare new coefficients
 
714
    double new_coeff0_0 = 0;
 
715
    double new_coeff0_1 = 0;
 
716
    double new_coeff0_2 = 0;
 
717
    
 
718
    // Loop possible derivatives
 
719
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
720
    {
 
721
      // Get values from coefficients array
 
722
      new_coeff0_0 = coefficients0[dof][0];
 
723
      new_coeff0_1 = coefficients0[dof][1];
 
724
      new_coeff0_2 = coefficients0[dof][2];
 
725
    
 
726
      // Loop derivative order
 
727
      for (unsigned int j = 0; j < n; j++)
 
728
      {
 
729
        // Update old coefficients
 
730
        coeff0_0 = new_coeff0_0;
 
731
        coeff0_1 = new_coeff0_1;
 
732
        coeff0_2 = new_coeff0_2;
 
733
    
 
734
        if(combinations[deriv_num][j] == 0)
 
735
        {
 
736
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
737
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
738
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
739
        }
 
740
        if(combinations[deriv_num][j] == 1)
 
741
        {
 
742
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
743
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
744
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
745
        }
 
746
    
 
747
      }
 
748
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
749
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
750
    }
 
751
    
 
752
    // Transform derivatives back to physical element
 
753
    for (unsigned int row = 0; row < num_derivatives; row++)
 
754
    {
 
755
      for (unsigned int col = 0; col < num_derivatives; col++)
 
756
      {
 
757
        values[row] += transform[row][col]*derivatives[col];
 
758
      }
 
759
    }
 
760
    // Delete pointer to array of derivatives on FIAT element
 
761
    delete [] derivatives;
 
762
    
 
763
    // Delete pointer to array of combinations of derivatives and transform
 
764
    for (unsigned int row = 0; row < num_derivatives; row++)
 
765
    {
 
766
      delete [] combinations[row];
 
767
      delete [] transform[row];
 
768
    }
 
769
    
 
770
    delete [] combinations;
 
771
    delete [] transform;
 
772
  }
 
773
 
 
774
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
775
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
776
                                              double* values,
 
777
                                              const double* coordinates,
 
778
                                              const ufc::cell& c) const
 
779
  {
 
780
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
781
  }
 
782
 
 
783
  /// Evaluate linear functional for dof i on the function f
 
784
  virtual double evaluate_dof(unsigned int i,
 
785
                              const ufc::function& f,
 
786
                              const ufc::cell& c) const
 
787
  {
 
788
    // The reference points, direction and weights:
 
789
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
790
    const static double W[3][1] = {{1}, {1}, {1}};
 
791
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
792
    
 
793
    const double * const * x = c.coordinates;
 
794
    double result = 0.0;
 
795
    // Iterate over the points:
 
796
    // Evaluate basis functions for affine mapping
 
797
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
798
    const double w1 = X[i][0][0];
 
799
    const double w2 = X[i][0][1];
 
800
    
 
801
    // Compute affine mapping y = F(X)
 
802
    double y[2];
 
803
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
804
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
805
    
 
806
    // Evaluate function at physical points
 
807
    double values[1];
 
808
    f.evaluate(values, y, c);
 
809
    
 
810
    // Map function values using appropriate mapping
 
811
    // Affine map: Do nothing
 
812
    
 
813
    // Note that we do not map the weights (yet).
 
814
    
 
815
    // Take directional components
 
816
    for(int k = 0; k < 1; k++)
 
817
      result += values[k]*D[i][0][k];
 
818
    // Multiply by weights 
 
819
    result *= W[i][0];
 
820
    
 
821
    return result;
 
822
  }
 
823
 
 
824
  /// Evaluate linear functionals for all dofs on the function f
 
825
  virtual void evaluate_dofs(double* values,
 
826
                             const ufc::function& f,
 
827
                             const ufc::cell& c) const
 
828
  {
 
829
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
830
  }
 
831
 
 
832
  /// Interpolate vertex values from dof values
 
833
  virtual void interpolate_vertex_values(double* vertex_values,
 
834
                                         const double* dof_values,
 
835
                                         const ufc::cell& c) const
 
836
  {
 
837
    // Evaluate at vertices and use affine mapping
 
838
    vertex_values[0] = dof_values[0];
 
839
    vertex_values[1] = dof_values[1];
 
840
    vertex_values[2] = dof_values[2];
 
841
  }
 
842
 
 
843
  /// Return the number of sub elements (for a mixed element)
 
844
  virtual unsigned int num_sub_elements() const
 
845
  {
 
846
    return 1;
 
847
  }
 
848
 
 
849
  /// Create a new finite element for sub element i (for a mixed element)
 
850
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
851
  {
 
852
    return new UFC_NonlinearPoissonBilinearForm_finite_element_1();
 
853
  }
 
854
 
 
855
};
 
856
 
 
857
/// This class defines the interface for a finite element.
 
858
 
 
859
class UFC_NonlinearPoissonBilinearForm_finite_element_2: public ufc::finite_element
 
860
{
 
861
public:
 
862
 
 
863
  /// Constructor
 
864
  UFC_NonlinearPoissonBilinearForm_finite_element_2() : ufc::finite_element()
 
865
  {
 
866
    // Do nothing
 
867
  }
 
868
 
 
869
  /// Destructor
 
870
  virtual ~UFC_NonlinearPoissonBilinearForm_finite_element_2()
 
871
  {
 
872
    // Do nothing
 
873
  }
 
874
 
 
875
  /// Return a string identifying the finite element
 
876
  virtual const char* signature() const
 
877
  {
 
878
    return "Lagrange finite element of degree 1 on a triangle";
 
879
  }
 
880
 
 
881
  /// Return the cell shape
 
882
  virtual ufc::shape cell_shape() const
 
883
  {
 
884
    return ufc::triangle;
 
885
  }
 
886
 
 
887
  /// Return the dimension of the finite element function space
 
888
  virtual unsigned int space_dimension() const
 
889
  {
 
890
    return 3;
 
891
  }
 
892
 
 
893
  /// Return the rank of the value space
 
894
  virtual unsigned int value_rank() const
 
895
  {
 
896
    return 0;
 
897
  }
 
898
 
 
899
  /// Return the dimension of the value space for axis i
 
900
  virtual unsigned int value_dimension(unsigned int i) const
 
901
  {
 
902
    return 1;
 
903
  }
 
904
 
 
905
  /// Evaluate basis function i at given point in cell
 
906
  virtual void evaluate_basis(unsigned int i,
 
907
                              double* values,
 
908
                              const double* coordinates,
 
909
                              const ufc::cell& c) const
 
910
  {
 
911
    // Extract vertex coordinates
 
912
    const double * const * element_coordinates = c.coordinates;
 
913
    
 
914
    // Compute Jacobian of affine map from reference cell
 
915
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
916
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
917
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
918
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
919
      
 
920
    // Compute determinant of Jacobian
 
921
    const double detJ = J_00*J_11 - J_01*J_10;
 
922
    
 
923
    // Compute inverse of Jacobian
 
924
    
 
925
    // Get coordinates and map to the reference (UFC) element
 
926
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
927
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
928
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
929
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
930
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
931
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
932
    
 
933
    // Map coordinates to the reference square
 
934
    if (std::abs(y - 1.0) < 1e-14)
 
935
      x = -1.0;
 
936
    else
 
937
      x = 2.0 *x/(1.0 - y) - 1.0;
 
938
    y = 2.0*y - 1.0;
 
939
    
 
940
    // Reset values
 
941
    *values = 0;
 
942
    
 
943
    // Map degree of freedom to element degree of freedom
 
944
    const unsigned int dof = i;
 
945
    
 
946
    // Generate scalings
 
947
    const double scalings_y_0 = 1;
 
948
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
949
    
 
950
    // Compute psitilde_a
 
951
    const double psitilde_a_0 = 1;
 
952
    const double psitilde_a_1 = x;
 
953
    
 
954
    // Compute psitilde_bs
 
955
    const double psitilde_bs_0_0 = 1;
 
956
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
957
    const double psitilde_bs_1_0 = 1;
 
958
    
 
959
    // Compute basisvalues
 
960
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
961
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
962
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
963
    
 
964
    // Table(s) of coefficients
 
965
    const static double coefficients0[3][3] = \
 
966
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
967
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
968
    {0.471404520791032, 0, 0.333333333333333}};
 
969
    
 
970
    // Extract relevant coefficients
 
971
    const double coeff0_0 = coefficients0[dof][0];
 
972
    const double coeff0_1 = coefficients0[dof][1];
 
973
    const double coeff0_2 = coefficients0[dof][2];
 
974
    
 
975
    // Compute value(s)
 
976
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
977
  }
 
978
 
 
979
  /// Evaluate all basis functions at given point in cell
 
980
  virtual void evaluate_basis_all(double* values,
 
981
                                  const double* coordinates,
 
982
                                  const ufc::cell& c) const
 
983
  {
 
984
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
985
  }
 
986
 
 
987
  /// Evaluate order n derivatives of basis function i at given point in cell
 
988
  virtual void evaluate_basis_derivatives(unsigned int i,
 
989
                                          unsigned int n,
 
990
                                          double* values,
 
991
                                          const double* coordinates,
 
992
                                          const ufc::cell& c) const
 
993
  {
 
994
    // Extract vertex coordinates
 
995
    const double * const * element_coordinates = c.coordinates;
 
996
    
 
997
    // Compute Jacobian of affine map from reference cell
 
998
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
999
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
1000
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
1001
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
1002
      
 
1003
    // Compute determinant of Jacobian
 
1004
    const double detJ = J_00*J_11 - J_01*J_10;
 
1005
    
 
1006
    // Compute inverse of Jacobian
 
1007
    
 
1008
    // Get coordinates and map to the reference (UFC) element
 
1009
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
1010
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
1011
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
1012
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
1013
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
1014
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
1015
    
 
1016
    // Map coordinates to the reference square
 
1017
    if (std::abs(y - 1.0) < 1e-14)
 
1018
      x = -1.0;
 
1019
    else
 
1020
      x = 2.0 *x/(1.0 - y) - 1.0;
 
1021
    y = 2.0*y - 1.0;
 
1022
    
 
1023
    // Compute number of derivatives
 
1024
    unsigned int num_derivatives = 1;
 
1025
    
 
1026
    for (unsigned int j = 0; j < n; j++)
 
1027
      num_derivatives *= 2;
 
1028
    
 
1029
    
 
1030
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
1031
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
1032
        
 
1033
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1034
    {
 
1035
      combinations[j] = new unsigned int [n];
 
1036
      for (unsigned int k = 0; k < n; k++)
 
1037
        combinations[j][k] = 0;
 
1038
    }
 
1039
        
 
1040
    // Generate combinations of derivatives
 
1041
    for (unsigned int row = 1; row < num_derivatives; row++)
 
1042
    {
 
1043
      for (unsigned int num = 0; num < row; num++)
 
1044
      {
 
1045
        for (unsigned int col = n-1; col+1 > 0; col--)
 
1046
        {
 
1047
          if (combinations[row][col] + 1 > 1)
 
1048
            combinations[row][col] = 0;
 
1049
          else
 
1050
          {
 
1051
            combinations[row][col] += 1;
 
1052
            break;
 
1053
          }
 
1054
        }
 
1055
      }
 
1056
    }
 
1057
    
 
1058
    // Compute inverse of Jacobian
 
1059
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
1060
    
 
1061
    // Declare transformation matrix
 
1062
    // Declare pointer to two dimensional array and initialise
 
1063
    double **transform = new double *[num_derivatives];
 
1064
        
 
1065
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1066
    {
 
1067
      transform[j] = new double [num_derivatives];
 
1068
      for (unsigned int k = 0; k < num_derivatives; k++)
 
1069
        transform[j][k] = 1;
 
1070
    }
 
1071
    
 
1072
    // Construct transformation matrix
 
1073
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1074
    {
 
1075
      for (unsigned int col = 0; col < num_derivatives; col++)
 
1076
      {
 
1077
        for (unsigned int k = 0; k < n; k++)
 
1078
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
1079
      }
 
1080
    }
 
1081
    
 
1082
    // Reset values
 
1083
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
1084
      values[j] = 0;
 
1085
    
 
1086
    // Map degree of freedom to element degree of freedom
 
1087
    const unsigned int dof = i;
 
1088
    
 
1089
    // Generate scalings
 
1090
    const double scalings_y_0 = 1;
 
1091
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
1092
    
 
1093
    // Compute psitilde_a
 
1094
    const double psitilde_a_0 = 1;
 
1095
    const double psitilde_a_1 = x;
 
1096
    
 
1097
    // Compute psitilde_bs
 
1098
    const double psitilde_bs_0_0 = 1;
 
1099
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
1100
    const double psitilde_bs_1_0 = 1;
 
1101
    
 
1102
    // Compute basisvalues
 
1103
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
1104
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
1105
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
1106
    
 
1107
    // Table(s) of coefficients
 
1108
    const static double coefficients0[3][3] = \
 
1109
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
1110
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
1111
    {0.471404520791032, 0, 0.333333333333333}};
 
1112
    
 
1113
    // Interesting (new) part
 
1114
    // Tables of derivatives of the polynomial base (transpose)
 
1115
    const static double dmats0[3][3] = \
 
1116
    {{0, 0, 0},
 
1117
    {4.89897948556636, 0, 0},
 
1118
    {0, 0, 0}};
 
1119
    
 
1120
    const static double dmats1[3][3] = \
 
1121
    {{0, 0, 0},
 
1122
    {2.44948974278318, 0, 0},
 
1123
    {4.24264068711928, 0, 0}};
 
1124
    
 
1125
    // Compute reference derivatives
 
1126
    // Declare pointer to array of derivatives on FIAT element
 
1127
    double *derivatives = new double [num_derivatives];
 
1128
    
 
1129
    // Declare coefficients
 
1130
    double coeff0_0 = 0;
 
1131
    double coeff0_1 = 0;
 
1132
    double coeff0_2 = 0;
 
1133
    
 
1134
    // Declare new coefficients
 
1135
    double new_coeff0_0 = 0;
 
1136
    double new_coeff0_1 = 0;
 
1137
    double new_coeff0_2 = 0;
 
1138
    
 
1139
    // Loop possible derivatives
 
1140
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
1141
    {
 
1142
      // Get values from coefficients array
 
1143
      new_coeff0_0 = coefficients0[dof][0];
 
1144
      new_coeff0_1 = coefficients0[dof][1];
 
1145
      new_coeff0_2 = coefficients0[dof][2];
 
1146
    
 
1147
      // Loop derivative order
 
1148
      for (unsigned int j = 0; j < n; j++)
 
1149
      {
 
1150
        // Update old coefficients
 
1151
        coeff0_0 = new_coeff0_0;
 
1152
        coeff0_1 = new_coeff0_1;
 
1153
        coeff0_2 = new_coeff0_2;
 
1154
    
 
1155
        if(combinations[deriv_num][j] == 0)
 
1156
        {
 
1157
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
1158
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
1159
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
1160
        }
 
1161
        if(combinations[deriv_num][j] == 1)
 
1162
        {
 
1163
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
1164
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
1165
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
1166
        }
 
1167
    
 
1168
      }
 
1169
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
1170
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
1171
    }
 
1172
    
 
1173
    // Transform derivatives back to physical element
 
1174
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1175
    {
 
1176
      for (unsigned int col = 0; col < num_derivatives; col++)
 
1177
      {
 
1178
        values[row] += transform[row][col]*derivatives[col];
 
1179
      }
 
1180
    }
 
1181
    // Delete pointer to array of derivatives on FIAT element
 
1182
    delete [] derivatives;
 
1183
    
 
1184
    // Delete pointer to array of combinations of derivatives and transform
 
1185
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1186
    {
 
1187
      delete [] combinations[row];
 
1188
      delete [] transform[row];
 
1189
    }
 
1190
    
 
1191
    delete [] combinations;
 
1192
    delete [] transform;
 
1193
  }
 
1194
 
 
1195
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
1196
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
1197
                                              double* values,
 
1198
                                              const double* coordinates,
 
1199
                                              const ufc::cell& c) const
 
1200
  {
 
1201
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
1202
  }
 
1203
 
 
1204
  /// Evaluate linear functional for dof i on the function f
 
1205
  virtual double evaluate_dof(unsigned int i,
 
1206
                              const ufc::function& f,
 
1207
                              const ufc::cell& c) const
 
1208
  {
 
1209
    // The reference points, direction and weights:
 
1210
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
1211
    const static double W[3][1] = {{1}, {1}, {1}};
 
1212
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
1213
    
 
1214
    const double * const * x = c.coordinates;
 
1215
    double result = 0.0;
 
1216
    // Iterate over the points:
 
1217
    // Evaluate basis functions for affine mapping
 
1218
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
1219
    const double w1 = X[i][0][0];
 
1220
    const double w2 = X[i][0][1];
 
1221
    
 
1222
    // Compute affine mapping y = F(X)
 
1223
    double y[2];
 
1224
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
1225
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
1226
    
 
1227
    // Evaluate function at physical points
 
1228
    double values[1];
 
1229
    f.evaluate(values, y, c);
 
1230
    
 
1231
    // Map function values using appropriate mapping
 
1232
    // Affine map: Do nothing
 
1233
    
 
1234
    // Note that we do not map the weights (yet).
 
1235
    
 
1236
    // Take directional components
 
1237
    for(int k = 0; k < 1; k++)
 
1238
      result += values[k]*D[i][0][k];
 
1239
    // Multiply by weights 
 
1240
    result *= W[i][0];
 
1241
    
 
1242
    return result;
 
1243
  }
 
1244
 
 
1245
  /// Evaluate linear functionals for all dofs on the function f
 
1246
  virtual void evaluate_dofs(double* values,
 
1247
                             const ufc::function& f,
 
1248
                             const ufc::cell& c) const
 
1249
  {
 
1250
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1251
  }
 
1252
 
 
1253
  /// Interpolate vertex values from dof values
 
1254
  virtual void interpolate_vertex_values(double* vertex_values,
 
1255
                                         const double* dof_values,
 
1256
                                         const ufc::cell& c) const
 
1257
  {
 
1258
    // Evaluate at vertices and use affine mapping
 
1259
    vertex_values[0] = dof_values[0];
 
1260
    vertex_values[1] = dof_values[1];
 
1261
    vertex_values[2] = dof_values[2];
 
1262
  }
 
1263
 
 
1264
  /// Return the number of sub elements (for a mixed element)
 
1265
  virtual unsigned int num_sub_elements() const
 
1266
  {
 
1267
    return 1;
 
1268
  }
 
1269
 
 
1270
  /// Create a new finite element for sub element i (for a mixed element)
 
1271
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
1272
  {
 
1273
    return new UFC_NonlinearPoissonBilinearForm_finite_element_2();
 
1274
  }
 
1275
 
 
1276
};
 
1277
 
 
1278
/// This class defines the interface for a local-to-global mapping of
 
1279
/// degrees of freedom (dofs).
 
1280
 
 
1281
class UFC_NonlinearPoissonBilinearForm_dof_map_0: public ufc::dof_map
 
1282
{
 
1283
private:
 
1284
 
 
1285
  unsigned int __global_dimension;
 
1286
 
 
1287
public:
 
1288
 
 
1289
  /// Constructor
 
1290
  UFC_NonlinearPoissonBilinearForm_dof_map_0() : ufc::dof_map()
 
1291
  {
 
1292
    __global_dimension = 0;
 
1293
  }
 
1294
 
 
1295
  /// Destructor
 
1296
  virtual ~UFC_NonlinearPoissonBilinearForm_dof_map_0()
 
1297
  {
 
1298
    // Do nothing
 
1299
  }
 
1300
 
 
1301
  /// Return a string identifying the dof map
 
1302
  virtual const char* signature() const
 
1303
  {
 
1304
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
 
1305
  }
 
1306
 
 
1307
  /// Return true iff mesh entities of topological dimension d are needed
 
1308
  virtual bool needs_mesh_entities(unsigned int d) const
 
1309
  {
 
1310
    switch ( d )
 
1311
    {
 
1312
    case 0:
 
1313
      return true;
 
1314
      break;
 
1315
    case 1:
 
1316
      return false;
 
1317
      break;
 
1318
    case 2:
 
1319
      return false;
 
1320
      break;
 
1321
    }
 
1322
    return false;
 
1323
  }
 
1324
 
 
1325
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1326
  virtual bool init_mesh(const ufc::mesh& m)
 
1327
  {
 
1328
    __global_dimension = m.num_entities[0];
 
1329
    return false;
 
1330
  }
 
1331
 
 
1332
  /// Initialize dof map for given cell
 
1333
  virtual void init_cell(const ufc::mesh& m,
 
1334
                         const ufc::cell& c)
 
1335
  {
 
1336
    // Do nothing
 
1337
  }
 
1338
 
 
1339
  /// Finish initialization of dof map for cells
 
1340
  virtual void init_cell_finalize()
 
1341
  {
 
1342
    // Do nothing
 
1343
  }
 
1344
 
 
1345
  /// Return the dimension of the global finite element function space
 
1346
  virtual unsigned int global_dimension() const
 
1347
  {
 
1348
    return __global_dimension;
 
1349
  }
 
1350
 
 
1351
  /// Return the dimension of the local finite element function space
 
1352
  virtual unsigned int local_dimension() const
 
1353
  {
 
1354
    return 3;
 
1355
  }
 
1356
 
 
1357
  // Return the geometric dimension of the coordinates this dof map provides
 
1358
  virtual unsigned int geometric_dimension() const
 
1359
  {
 
1360
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1361
  }
 
1362
 
 
1363
  /// Return the number of dofs on each cell facet
 
1364
  virtual unsigned int num_facet_dofs() const
 
1365
  {
 
1366
    return 2;
 
1367
  }
 
1368
 
 
1369
  /// Return the number of dofs associated with each cell entity of dimension d
 
1370
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1371
  {
 
1372
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1373
  }
 
1374
 
 
1375
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1376
  virtual void tabulate_dofs(unsigned int* dofs,
 
1377
                             const ufc::mesh& m,
 
1378
                             const ufc::cell& c) const
 
1379
  {
 
1380
    dofs[0] = c.entity_indices[0][0];
 
1381
    dofs[1] = c.entity_indices[0][1];
 
1382
    dofs[2] = c.entity_indices[0][2];
 
1383
  }
 
1384
 
 
1385
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1386
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1387
                                   unsigned int facet) const
 
1388
  {
 
1389
    switch ( facet )
 
1390
    {
 
1391
    case 0:
 
1392
      dofs[0] = 1;
 
1393
      dofs[1] = 2;
 
1394
      break;
 
1395
    case 1:
 
1396
      dofs[0] = 0;
 
1397
      dofs[1] = 2;
 
1398
      break;
 
1399
    case 2:
 
1400
      dofs[0] = 0;
 
1401
      dofs[1] = 1;
 
1402
      break;
 
1403
    }
 
1404
  }
 
1405
 
 
1406
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1407
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1408
                                    unsigned int d, unsigned int i) const
 
1409
  {
 
1410
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1411
  }
 
1412
 
 
1413
  /// Tabulate the coordinates of all dofs on a cell
 
1414
  virtual void tabulate_coordinates(double** coordinates,
 
1415
                                    const ufc::cell& c) const
 
1416
  {
 
1417
    const double * const * x = c.coordinates;
 
1418
    coordinates[0][0] = x[0][0];
 
1419
    coordinates[0][1] = x[0][1];
 
1420
    coordinates[1][0] = x[1][0];
 
1421
    coordinates[1][1] = x[1][1];
 
1422
    coordinates[2][0] = x[2][0];
 
1423
    coordinates[2][1] = x[2][1];
 
1424
  }
 
1425
 
 
1426
  /// Return the number of sub dof maps (for a mixed element)
 
1427
  virtual unsigned int num_sub_dof_maps() const
 
1428
  {
 
1429
    return 1;
 
1430
  }
 
1431
 
 
1432
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1433
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1434
  {
 
1435
    return new UFC_NonlinearPoissonBilinearForm_dof_map_0();
 
1436
  }
 
1437
 
 
1438
};
 
1439
 
 
1440
/// This class defines the interface for a local-to-global mapping of
 
1441
/// degrees of freedom (dofs).
 
1442
 
 
1443
class UFC_NonlinearPoissonBilinearForm_dof_map_1: public ufc::dof_map
 
1444
{
 
1445
private:
 
1446
 
 
1447
  unsigned int __global_dimension;
 
1448
 
 
1449
public:
 
1450
 
 
1451
  /// Constructor
 
1452
  UFC_NonlinearPoissonBilinearForm_dof_map_1() : ufc::dof_map()
 
1453
  {
 
1454
    __global_dimension = 0;
 
1455
  }
 
1456
 
 
1457
  /// Destructor
 
1458
  virtual ~UFC_NonlinearPoissonBilinearForm_dof_map_1()
 
1459
  {
 
1460
    // Do nothing
 
1461
  }
 
1462
 
 
1463
  /// Return a string identifying the dof map
 
1464
  virtual const char* signature() const
 
1465
  {
 
1466
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
 
1467
  }
 
1468
 
 
1469
  /// Return true iff mesh entities of topological dimension d are needed
 
1470
  virtual bool needs_mesh_entities(unsigned int d) const
 
1471
  {
 
1472
    switch ( d )
 
1473
    {
 
1474
    case 0:
 
1475
      return true;
 
1476
      break;
 
1477
    case 1:
 
1478
      return false;
 
1479
      break;
 
1480
    case 2:
 
1481
      return false;
 
1482
      break;
 
1483
    }
 
1484
    return false;
 
1485
  }
 
1486
 
 
1487
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1488
  virtual bool init_mesh(const ufc::mesh& m)
 
1489
  {
 
1490
    __global_dimension = m.num_entities[0];
 
1491
    return false;
 
1492
  }
 
1493
 
 
1494
  /// Initialize dof map for given cell
 
1495
  virtual void init_cell(const ufc::mesh& m,
 
1496
                         const ufc::cell& c)
 
1497
  {
 
1498
    // Do nothing
 
1499
  }
 
1500
 
 
1501
  /// Finish initialization of dof map for cells
 
1502
  virtual void init_cell_finalize()
 
1503
  {
 
1504
    // Do nothing
 
1505
  }
 
1506
 
 
1507
  /// Return the dimension of the global finite element function space
 
1508
  virtual unsigned int global_dimension() const
 
1509
  {
 
1510
    return __global_dimension;
 
1511
  }
 
1512
 
 
1513
  /// Return the dimension of the local finite element function space
 
1514
  virtual unsigned int local_dimension() const
 
1515
  {
 
1516
    return 3;
 
1517
  }
 
1518
 
 
1519
  // Return the geometric dimension of the coordinates this dof map provides
 
1520
  virtual unsigned int geometric_dimension() const
 
1521
  {
 
1522
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1523
  }
 
1524
 
 
1525
  /// Return the number of dofs on each cell facet
 
1526
  virtual unsigned int num_facet_dofs() const
 
1527
  {
 
1528
    return 2;
 
1529
  }
 
1530
 
 
1531
  /// Return the number of dofs associated with each cell entity of dimension d
 
1532
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1533
  {
 
1534
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1535
  }
 
1536
 
 
1537
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1538
  virtual void tabulate_dofs(unsigned int* dofs,
 
1539
                             const ufc::mesh& m,
 
1540
                             const ufc::cell& c) const
 
1541
  {
 
1542
    dofs[0] = c.entity_indices[0][0];
 
1543
    dofs[1] = c.entity_indices[0][1];
 
1544
    dofs[2] = c.entity_indices[0][2];
 
1545
  }
 
1546
 
 
1547
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1548
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1549
                                   unsigned int facet) const
 
1550
  {
 
1551
    switch ( facet )
 
1552
    {
 
1553
    case 0:
 
1554
      dofs[0] = 1;
 
1555
      dofs[1] = 2;
 
1556
      break;
 
1557
    case 1:
 
1558
      dofs[0] = 0;
 
1559
      dofs[1] = 2;
 
1560
      break;
 
1561
    case 2:
 
1562
      dofs[0] = 0;
 
1563
      dofs[1] = 1;
 
1564
      break;
 
1565
    }
 
1566
  }
 
1567
 
 
1568
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1569
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1570
                                    unsigned int d, unsigned int i) const
 
1571
  {
 
1572
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1573
  }
 
1574
 
 
1575
  /// Tabulate the coordinates of all dofs on a cell
 
1576
  virtual void tabulate_coordinates(double** coordinates,
 
1577
                                    const ufc::cell& c) const
 
1578
  {
 
1579
    const double * const * x = c.coordinates;
 
1580
    coordinates[0][0] = x[0][0];
 
1581
    coordinates[0][1] = x[0][1];
 
1582
    coordinates[1][0] = x[1][0];
 
1583
    coordinates[1][1] = x[1][1];
 
1584
    coordinates[2][0] = x[2][0];
 
1585
    coordinates[2][1] = x[2][1];
 
1586
  }
 
1587
 
 
1588
  /// Return the number of sub dof maps (for a mixed element)
 
1589
  virtual unsigned int num_sub_dof_maps() const
 
1590
  {
 
1591
    return 1;
 
1592
  }
 
1593
 
 
1594
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1595
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1596
  {
 
1597
    return new UFC_NonlinearPoissonBilinearForm_dof_map_1();
 
1598
  }
 
1599
 
 
1600
};
 
1601
 
 
1602
/// This class defines the interface for a local-to-global mapping of
 
1603
/// degrees of freedom (dofs).
 
1604
 
 
1605
class UFC_NonlinearPoissonBilinearForm_dof_map_2: public ufc::dof_map
 
1606
{
 
1607
private:
 
1608
 
 
1609
  unsigned int __global_dimension;
 
1610
 
 
1611
public:
 
1612
 
 
1613
  /// Constructor
 
1614
  UFC_NonlinearPoissonBilinearForm_dof_map_2() : ufc::dof_map()
 
1615
  {
 
1616
    __global_dimension = 0;
 
1617
  }
 
1618
 
 
1619
  /// Destructor
 
1620
  virtual ~UFC_NonlinearPoissonBilinearForm_dof_map_2()
 
1621
  {
 
1622
    // Do nothing
 
1623
  }
 
1624
 
 
1625
  /// Return a string identifying the dof map
 
1626
  virtual const char* signature() const
 
1627
  {
 
1628
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
 
1629
  }
 
1630
 
 
1631
  /// Return true iff mesh entities of topological dimension d are needed
 
1632
  virtual bool needs_mesh_entities(unsigned int d) const
 
1633
  {
 
1634
    switch ( d )
 
1635
    {
 
1636
    case 0:
 
1637
      return true;
 
1638
      break;
 
1639
    case 1:
 
1640
      return false;
 
1641
      break;
 
1642
    case 2:
 
1643
      return false;
 
1644
      break;
 
1645
    }
 
1646
    return false;
 
1647
  }
 
1648
 
 
1649
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1650
  virtual bool init_mesh(const ufc::mesh& m)
 
1651
  {
 
1652
    __global_dimension = m.num_entities[0];
 
1653
    return false;
 
1654
  }
 
1655
 
 
1656
  /// Initialize dof map for given cell
 
1657
  virtual void init_cell(const ufc::mesh& m,
 
1658
                         const ufc::cell& c)
 
1659
  {
 
1660
    // Do nothing
 
1661
  }
 
1662
 
 
1663
  /// Finish initialization of dof map for cells
 
1664
  virtual void init_cell_finalize()
 
1665
  {
 
1666
    // Do nothing
 
1667
  }
 
1668
 
 
1669
  /// Return the dimension of the global finite element function space
 
1670
  virtual unsigned int global_dimension() const
 
1671
  {
 
1672
    return __global_dimension;
 
1673
  }
 
1674
 
 
1675
  /// Return the dimension of the local finite element function space
 
1676
  virtual unsigned int local_dimension() const
 
1677
  {
 
1678
    return 3;
 
1679
  }
 
1680
 
 
1681
  // Return the geometric dimension of the coordinates this dof map provides
 
1682
  virtual unsigned int geometric_dimension() const
 
1683
  {
 
1684
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1685
  }
 
1686
 
 
1687
  /// Return the number of dofs on each cell facet
 
1688
  virtual unsigned int num_facet_dofs() const
 
1689
  {
 
1690
    return 2;
 
1691
  }
 
1692
 
 
1693
  /// Return the number of dofs associated with each cell entity of dimension d
 
1694
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1695
  {
 
1696
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1697
  }
 
1698
 
 
1699
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1700
  virtual void tabulate_dofs(unsigned int* dofs,
 
1701
                             const ufc::mesh& m,
 
1702
                             const ufc::cell& c) const
 
1703
  {
 
1704
    dofs[0] = c.entity_indices[0][0];
 
1705
    dofs[1] = c.entity_indices[0][1];
 
1706
    dofs[2] = c.entity_indices[0][2];
 
1707
  }
 
1708
 
 
1709
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1710
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1711
                                   unsigned int facet) const
 
1712
  {
 
1713
    switch ( facet )
 
1714
    {
 
1715
    case 0:
 
1716
      dofs[0] = 1;
 
1717
      dofs[1] = 2;
 
1718
      break;
 
1719
    case 1:
 
1720
      dofs[0] = 0;
 
1721
      dofs[1] = 2;
 
1722
      break;
 
1723
    case 2:
 
1724
      dofs[0] = 0;
 
1725
      dofs[1] = 1;
 
1726
      break;
 
1727
    }
 
1728
  }
 
1729
 
 
1730
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1731
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1732
                                    unsigned int d, unsigned int i) const
 
1733
  {
 
1734
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1735
  }
 
1736
 
 
1737
  /// Tabulate the coordinates of all dofs on a cell
 
1738
  virtual void tabulate_coordinates(double** coordinates,
 
1739
                                    const ufc::cell& c) const
 
1740
  {
 
1741
    const double * const * x = c.coordinates;
 
1742
    coordinates[0][0] = x[0][0];
 
1743
    coordinates[0][1] = x[0][1];
 
1744
    coordinates[1][0] = x[1][0];
 
1745
    coordinates[1][1] = x[1][1];
 
1746
    coordinates[2][0] = x[2][0];
 
1747
    coordinates[2][1] = x[2][1];
 
1748
  }
 
1749
 
 
1750
  /// Return the number of sub dof maps (for a mixed element)
 
1751
  virtual unsigned int num_sub_dof_maps() const
 
1752
  {
 
1753
    return 1;
 
1754
  }
 
1755
 
 
1756
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1757
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1758
  {
 
1759
    return new UFC_NonlinearPoissonBilinearForm_dof_map_2();
 
1760
  }
 
1761
 
 
1762
};
 
1763
 
 
1764
/// This class defines the interface for the tabulation of the cell
 
1765
/// tensor corresponding to the local contribution to a form from
 
1766
/// the integral over a cell.
 
1767
 
 
1768
class UFC_NonlinearPoissonBilinearForm_cell_integral_0: public ufc::cell_integral
 
1769
{
 
1770
public:
 
1771
 
 
1772
  /// Constructor
 
1773
  UFC_NonlinearPoissonBilinearForm_cell_integral_0() : ufc::cell_integral()
 
1774
  {
 
1775
    // Do nothing
 
1776
  }
 
1777
 
 
1778
  /// Destructor
 
1779
  virtual ~UFC_NonlinearPoissonBilinearForm_cell_integral_0()
 
1780
  {
 
1781
    // Do nothing
 
1782
  }
 
1783
 
 
1784
  /// Tabulate the tensor for the contribution from a local cell
 
1785
  virtual void tabulate_tensor(double* A,
 
1786
                               const double * const * w,
 
1787
                               const ufc::cell& c) const
 
1788
  {
 
1789
    // Extract vertex coordinates
 
1790
    const double * const * x = c.coordinates;
 
1791
    
 
1792
    // Compute Jacobian of affine map from reference cell
 
1793
    const double J_00 = x[1][0] - x[0][0];
 
1794
    const double J_01 = x[2][0] - x[0][0];
 
1795
    const double J_10 = x[1][1] - x[0][1];
 
1796
    const double J_11 = x[2][1] - x[0][1];
 
1797
      
 
1798
    // Compute determinant of Jacobian
 
1799
    double detJ = J_00*J_11 - J_01*J_10;
 
1800
      
 
1801
    // Compute inverse of Jacobian
 
1802
    const double Jinv_00 =  J_11 / detJ;
 
1803
    const double Jinv_01 = -J_01 / detJ;
 
1804
    const double Jinv_10 = -J_10 / detJ;
 
1805
    const double Jinv_11 =  J_00 / detJ;
 
1806
    
 
1807
    // Set scale factor
 
1808
    const double det = std::abs(detJ);
 
1809
    
 
1810
    // Compute coefficients
 
1811
    const double c0_0_0_0 = w[0][0];
 
1812
    const double c0_0_0_1 = w[0][1];
 
1813
    const double c0_0_0_2 = w[0][2];
 
1814
    const double c0_0_1_0 = w[0][0];
 
1815
    const double c0_0_1_1 = w[0][1];
 
1816
    const double c0_0_1_2 = w[0][2];
 
1817
    const double c0_2_0_0 = w[0][0];
 
1818
    const double c0_2_0_1 = w[0][1];
 
1819
    const double c0_2_0_2 = w[0][2];
 
1820
    const double c0_2_1_0 = w[0][0];
 
1821
    const double c0_2_1_1 = w[0][1];
 
1822
    const double c0_2_1_2 = w[0][2];
 
1823
    
 
1824
    // Compute geometry tensors
 
1825
    const double G0_0_0_0_0 = det*c0_0_0_0*c0_0_1_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1826
    const double G0_0_0_0_1 = det*c0_0_0_0*c0_0_1_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1827
    const double G0_0_0_1_0 = det*c0_0_0_0*c0_0_1_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1828
    const double G0_0_0_1_1 = det*c0_0_0_0*c0_0_1_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1829
    const double G0_0_1_0_0 = det*c0_0_0_0*c0_0_1_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1830
    const double G0_0_1_0_1 = det*c0_0_0_0*c0_0_1_1*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1831
    const double G0_0_1_1_0 = det*c0_0_0_0*c0_0_1_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1832
    const double G0_0_1_1_1 = det*c0_0_0_0*c0_0_1_1*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1833
    const double G0_0_2_0_0 = det*c0_0_0_0*c0_0_1_2*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1834
    const double G0_0_2_0_1 = det*c0_0_0_0*c0_0_1_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1835
    const double G0_0_2_1_0 = det*c0_0_0_0*c0_0_1_2*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1836
    const double G0_0_2_1_1 = det*c0_0_0_0*c0_0_1_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1837
    const double G0_1_0_0_0 = det*c0_0_0_1*c0_0_1_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1838
    const double G0_1_0_0_1 = det*c0_0_0_1*c0_0_1_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1839
    const double G0_1_0_1_0 = det*c0_0_0_1*c0_0_1_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1840
    const double G0_1_0_1_1 = det*c0_0_0_1*c0_0_1_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1841
    const double G0_1_1_0_0 = det*c0_0_0_1*c0_0_1_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1842
    const double G0_1_1_0_1 = det*c0_0_0_1*c0_0_1_1*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1843
    const double G0_1_1_1_0 = det*c0_0_0_1*c0_0_1_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1844
    const double G0_1_1_1_1 = det*c0_0_0_1*c0_0_1_1*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1845
    const double G0_1_2_0_0 = det*c0_0_0_1*c0_0_1_2*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1846
    const double G0_1_2_0_1 = det*c0_0_0_1*c0_0_1_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1847
    const double G0_1_2_1_0 = det*c0_0_0_1*c0_0_1_2*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1848
    const double G0_1_2_1_1 = det*c0_0_0_1*c0_0_1_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1849
    const double G0_2_0_0_0 = det*c0_0_0_2*c0_0_1_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1850
    const double G0_2_0_0_1 = det*c0_0_0_2*c0_0_1_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1851
    const double G0_2_0_1_0 = det*c0_0_0_2*c0_0_1_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1852
    const double G0_2_0_1_1 = det*c0_0_0_2*c0_0_1_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1853
    const double G0_2_1_0_0 = det*c0_0_0_2*c0_0_1_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1854
    const double G0_2_1_0_1 = det*c0_0_0_2*c0_0_1_1*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1855
    const double G0_2_1_1_0 = det*c0_0_0_2*c0_0_1_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1856
    const double G0_2_1_1_1 = det*c0_0_0_2*c0_0_1_1*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1857
    const double G0_2_2_0_0 = det*c0_0_0_2*c0_0_1_2*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1858
    const double G0_2_2_0_1 = det*c0_0_0_2*c0_0_1_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1859
    const double G0_2_2_1_0 = det*c0_0_0_2*c0_0_1_2*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1860
    const double G0_2_2_1_1 = det*c0_0_0_2*c0_0_1_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1861
    const double G1_0_0 = det*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1862
    const double G1_0_1 = det*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1863
    const double G1_1_0 = det*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1864
    const double G1_1_1 = det*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1865
    const double G2_0_0_0_0 = det*c0_2_0_0*c0_2_1_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1866
    const double G2_0_0_0_1 = det*c0_2_0_0*c0_2_1_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1867
    const double G2_0_0_1_0 = det*c0_2_0_0*c0_2_1_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1868
    const double G2_0_0_2_1 = det*c0_2_0_0*c0_2_1_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1869
    const double G2_0_1_0_0 = det*c0_2_0_0*c0_2_1_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1870
    const double G2_0_1_0_1 = det*c0_2_0_0*c0_2_1_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1871
    const double G2_0_1_1_0 = det*c0_2_0_0*c0_2_1_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1872
    const double G2_0_1_2_1 = det*c0_2_0_0*c0_2_1_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1873
    const double G2_1_0_0_0 = det*c0_2_0_1*c0_2_1_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1874
    const double G2_1_0_0_1 = det*c0_2_0_1*c0_2_1_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1875
    const double G2_1_0_1_0 = det*c0_2_0_1*c0_2_1_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1876
    const double G2_1_0_2_1 = det*c0_2_0_1*c0_2_1_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1877
    const double G2_1_1_0_0 = det*c0_2_0_1*c0_2_1_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1878
    const double G2_1_1_0_1 = det*c0_2_0_1*c0_2_1_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1879
    const double G2_1_1_1_0 = det*c0_2_0_1*c0_2_1_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1880
    const double G2_1_1_2_1 = det*c0_2_0_1*c0_2_1_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1881
    const double G2_2_0_0_0 = det*c0_2_0_2*c0_2_1_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1882
    const double G2_2_0_0_1 = det*c0_2_0_2*c0_2_1_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1883
    const double G2_2_0_1_0 = det*c0_2_0_2*c0_2_1_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
1884
    const double G2_2_0_2_1 = det*c0_2_0_2*c0_2_1_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
1885
    const double G2_2_1_0_0 = det*c0_2_0_2*c0_2_1_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1886
    const double G2_2_1_0_1 = det*c0_2_0_2*c0_2_1_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1887
    const double G2_2_1_1_0 = det*c0_2_0_2*c0_2_1_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
1888
    const double G2_2_1_2_1 = det*c0_2_0_2*c0_2_1_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
1889
    
 
1890
    // Compute element tensor
 
1891
    A[0] = 0.0833333333333332*G0_0_0_0_0 + 0.0833333333333332*G0_0_0_0_1 + 0.0833333333333332*G0_0_0_1_0 + 0.0833333333333332*G0_0_0_1_1 + 0.0416666666666666*G0_0_1_0_0 + 0.0416666666666666*G0_0_1_0_1 + 0.0416666666666666*G0_0_1_1_0 + 0.0416666666666666*G0_0_1_1_1 + 0.0416666666666666*G0_0_2_0_0 + 0.0416666666666666*G0_0_2_0_1 + 0.0416666666666666*G0_0_2_1_0 + 0.0416666666666666*G0_0_2_1_1 + 0.0416666666666666*G0_1_0_0_0 + 0.0416666666666666*G0_1_0_0_1 + 0.0416666666666666*G0_1_0_1_0 + 0.0416666666666666*G0_1_0_1_1 + 0.0833333333333332*G0_1_1_0_0 + 0.0833333333333332*G0_1_1_0_1 + 0.0833333333333332*G0_1_1_1_0 + 0.0833333333333332*G0_1_1_1_1 + 0.0416666666666666*G0_1_2_0_0 + 0.0416666666666666*G0_1_2_0_1 + 0.0416666666666666*G0_1_2_1_0 + 0.0416666666666666*G0_1_2_1_1 + 0.0416666666666666*G0_2_0_0_0 + 0.0416666666666666*G0_2_0_0_1 + 0.0416666666666666*G0_2_0_1_0 + 0.0416666666666666*G0_2_0_1_1 + 0.0416666666666666*G0_2_1_0_0 + 0.0416666666666666*G0_2_1_0_1 + 0.0416666666666666*G0_2_1_1_0 + 0.0416666666666666*G0_2_1_1_1 + 0.0833333333333332*G0_2_2_0_0 + 0.0833333333333332*G0_2_2_0_1 + 0.0833333333333332*G0_2_2_1_0 + 0.0833333333333332*G0_2_2_1_1 + 0.5*G1_0_0 + 0.5*G1_0_1 + 0.5*G1_1_0 + 0.5*G1_1_1 + 0.166666666666666*G2_0_0_0_0 + 0.166666666666666*G2_0_0_0_1 - 0.166666666666666*G2_0_0_1_0 - 0.166666666666666*G2_0_0_2_1 + 0.166666666666666*G2_0_1_0_0 + 0.166666666666666*G2_0_1_0_1 - 0.166666666666666*G2_0_1_1_0 - 0.166666666666666*G2_0_1_2_1 + 0.0833333333333332*G2_1_0_0_0 + 0.0833333333333332*G2_1_0_0_1 - 0.0833333333333332*G2_1_0_1_0 - 0.0833333333333332*G2_1_0_2_1 + 0.0833333333333332*G2_1_1_0_0 + 0.0833333333333332*G2_1_1_0_1 - 0.0833333333333332*G2_1_1_1_0 - 0.0833333333333332*G2_1_1_2_1 + 0.0833333333333332*G2_2_0_0_0 + 0.0833333333333332*G2_2_0_0_1 - 0.0833333333333332*G2_2_0_1_0 - 0.0833333333333332*G2_2_0_2_1 + 0.0833333333333332*G2_2_1_0_0 + 0.0833333333333332*G2_2_1_0_1 - 0.0833333333333332*G2_2_1_1_0 - 0.0833333333333332*G2_2_1_2_1;
 
1892
    A[1] = -0.0833333333333332*G0_0_0_0_0 - 0.0833333333333332*G0_0_0_1_0 - 0.0416666666666666*G0_0_1_0_0 - 0.0416666666666666*G0_0_1_1_0 - 0.0416666666666666*G0_0_2_0_0 - 0.0416666666666666*G0_0_2_1_0 - 0.0416666666666666*G0_1_0_0_0 - 0.0416666666666666*G0_1_0_1_0 - 0.0833333333333332*G0_1_1_0_0 - 0.0833333333333332*G0_1_1_1_0 - 0.0416666666666666*G0_1_2_0_0 - 0.0416666666666666*G0_1_2_1_0 - 0.0416666666666666*G0_2_0_0_0 - 0.0416666666666666*G0_2_0_1_0 - 0.0416666666666666*G0_2_1_0_0 - 0.0416666666666666*G0_2_1_1_0 - 0.0833333333333332*G0_2_2_0_0 - 0.0833333333333332*G0_2_2_1_0 - 0.5*G1_0_0 - 0.5*G1_1_0 + 0.0833333333333332*G2_0_0_0_0 + 0.0833333333333332*G2_0_0_0_1 - 0.0833333333333332*G2_0_0_1_0 - 0.0833333333333332*G2_0_0_2_1 + 0.0833333333333332*G2_0_1_0_0 + 0.0833333333333332*G2_0_1_0_1 - 0.0833333333333332*G2_0_1_1_0 - 0.0833333333333332*G2_0_1_2_1 + 0.166666666666666*G2_1_0_0_0 + 0.166666666666666*G2_1_0_0_1 - 0.166666666666666*G2_1_0_1_0 - 0.166666666666666*G2_1_0_2_1 + 0.166666666666666*G2_1_1_0_0 + 0.166666666666666*G2_1_1_0_1 - 0.166666666666666*G2_1_1_1_0 - 0.166666666666666*G2_1_1_2_1 + 0.0833333333333332*G2_2_0_0_0 + 0.0833333333333332*G2_2_0_0_1 - 0.0833333333333332*G2_2_0_1_0 - 0.0833333333333332*G2_2_0_2_1 + 0.0833333333333332*G2_2_1_0_0 + 0.0833333333333332*G2_2_1_0_1 - 0.0833333333333332*G2_2_1_1_0 - 0.0833333333333332*G2_2_1_2_1;
 
1893
    A[2] = -0.0833333333333332*G0_0_0_0_1 - 0.0833333333333332*G0_0_0_1_1 - 0.0416666666666666*G0_0_1_0_1 - 0.0416666666666666*G0_0_1_1_1 - 0.0416666666666666*G0_0_2_0_1 - 0.0416666666666666*G0_0_2_1_1 - 0.0416666666666666*G0_1_0_0_1 - 0.0416666666666666*G0_1_0_1_1 - 0.0833333333333332*G0_1_1_0_1 - 0.0833333333333332*G0_1_1_1_1 - 0.0416666666666666*G0_1_2_0_1 - 0.0416666666666666*G0_1_2_1_1 - 0.0416666666666666*G0_2_0_0_1 - 0.0416666666666666*G0_2_0_1_1 - 0.0416666666666666*G0_2_1_0_1 - 0.0416666666666666*G0_2_1_1_1 - 0.0833333333333332*G0_2_2_0_1 - 0.0833333333333332*G0_2_2_1_1 - 0.5*G1_0_1 - 0.5*G1_1_1 + 0.0833333333333332*G2_0_0_0_0 + 0.0833333333333332*G2_0_0_0_1 - 0.0833333333333332*G2_0_0_1_0 - 0.0833333333333332*G2_0_0_2_1 + 0.0833333333333332*G2_0_1_0_0 + 0.0833333333333332*G2_0_1_0_1 - 0.0833333333333332*G2_0_1_1_0 - 0.0833333333333332*G2_0_1_2_1 + 0.0833333333333332*G2_1_0_0_0 + 0.0833333333333332*G2_1_0_0_1 - 0.0833333333333332*G2_1_0_1_0 - 0.0833333333333332*G2_1_0_2_1 + 0.0833333333333332*G2_1_1_0_0 + 0.0833333333333332*G2_1_1_0_1 - 0.0833333333333332*G2_1_1_1_0 - 0.0833333333333332*G2_1_1_2_1 + 0.166666666666666*G2_2_0_0_0 + 0.166666666666666*G2_2_0_0_1 - 0.166666666666666*G2_2_0_1_0 - 0.166666666666666*G2_2_0_2_1 + 0.166666666666666*G2_2_1_0_0 + 0.166666666666666*G2_2_1_0_1 - 0.166666666666666*G2_2_1_1_0 - 0.166666666666666*G2_2_1_2_1;
 
1894
    A[3] = -0.0833333333333332*G0_0_0_0_0 - 0.0833333333333332*G0_0_0_0_1 - 0.0416666666666666*G0_0_1_0_0 - 0.0416666666666666*G0_0_1_0_1 - 0.0416666666666666*G0_0_2_0_0 - 0.0416666666666666*G0_0_2_0_1 - 0.0416666666666666*G0_1_0_0_0 - 0.0416666666666666*G0_1_0_0_1 - 0.0833333333333332*G0_1_1_0_0 - 0.0833333333333332*G0_1_1_0_1 - 0.0416666666666666*G0_1_2_0_0 - 0.0416666666666666*G0_1_2_0_1 - 0.0416666666666666*G0_2_0_0_0 - 0.0416666666666666*G0_2_0_0_1 - 0.0416666666666666*G0_2_1_0_0 - 0.0416666666666666*G0_2_1_0_1 - 0.0833333333333332*G0_2_2_0_0 - 0.0833333333333332*G0_2_2_0_1 - 0.5*G1_0_0 - 0.5*G1_0_1 - 0.166666666666666*G2_0_0_0_0 - 0.166666666666666*G2_0_0_0_1 + 0.166666666666666*G2_0_0_1_0 + 0.166666666666666*G2_0_0_2_1 - 0.0833333333333332*G2_1_0_0_0 - 0.0833333333333332*G2_1_0_0_1 + 0.0833333333333332*G2_1_0_1_0 + 0.0833333333333332*G2_1_0_2_1 - 0.0833333333333332*G2_2_0_0_0 - 0.0833333333333332*G2_2_0_0_1 + 0.0833333333333332*G2_2_0_1_0 + 0.0833333333333332*G2_2_0_2_1;
 
1895
    A[4] = 0.0833333333333332*G0_0_0_0_0 + 0.0416666666666666*G0_0_1_0_0 + 0.0416666666666666*G0_0_2_0_0 + 0.0416666666666666*G0_1_0_0_0 + 0.0833333333333332*G0_1_1_0_0 + 0.0416666666666666*G0_1_2_0_0 + 0.0416666666666666*G0_2_0_0_0 + 0.0416666666666666*G0_2_1_0_0 + 0.0833333333333332*G0_2_2_0_0 + 0.5*G1_0_0 - 0.0833333333333332*G2_0_0_0_0 - 0.0833333333333332*G2_0_0_0_1 + 0.0833333333333332*G2_0_0_1_0 + 0.0833333333333332*G2_0_0_2_1 - 0.166666666666666*G2_1_0_0_0 - 0.166666666666666*G2_1_0_0_1 + 0.166666666666666*G2_1_0_1_0 + 0.166666666666666*G2_1_0_2_1 - 0.0833333333333332*G2_2_0_0_0 - 0.0833333333333332*G2_2_0_0_1 + 0.0833333333333332*G2_2_0_1_0 + 0.0833333333333332*G2_2_0_2_1;
 
1896
    A[5] = 0.0833333333333332*G0_0_0_0_1 + 0.0416666666666666*G0_0_1_0_1 + 0.0416666666666666*G0_0_2_0_1 + 0.0416666666666666*G0_1_0_0_1 + 0.0833333333333332*G0_1_1_0_1 + 0.0416666666666666*G0_1_2_0_1 + 0.0416666666666666*G0_2_0_0_1 + 0.0416666666666666*G0_2_1_0_1 + 0.0833333333333332*G0_2_2_0_1 + 0.5*G1_0_1 - 0.0833333333333332*G2_0_0_0_0 - 0.0833333333333332*G2_0_0_0_1 + 0.0833333333333332*G2_0_0_1_0 + 0.0833333333333332*G2_0_0_2_1 - 0.0833333333333332*G2_1_0_0_0 - 0.0833333333333332*G2_1_0_0_1 + 0.0833333333333332*G2_1_0_1_0 + 0.0833333333333332*G2_1_0_2_1 - 0.166666666666666*G2_2_0_0_0 - 0.166666666666666*G2_2_0_0_1 + 0.166666666666666*G2_2_0_1_0 + 0.166666666666666*G2_2_0_2_1;
 
1897
    A[6] = -0.0833333333333332*G0_0_0_1_0 - 0.0833333333333332*G0_0_0_1_1 - 0.0416666666666666*G0_0_1_1_0 - 0.0416666666666666*G0_0_1_1_1 - 0.0416666666666666*G0_0_2_1_0 - 0.0416666666666666*G0_0_2_1_1 - 0.0416666666666666*G0_1_0_1_0 - 0.0416666666666666*G0_1_0_1_1 - 0.0833333333333332*G0_1_1_1_0 - 0.0833333333333332*G0_1_1_1_1 - 0.0416666666666666*G0_1_2_1_0 - 0.0416666666666666*G0_1_2_1_1 - 0.0416666666666666*G0_2_0_1_0 - 0.0416666666666666*G0_2_0_1_1 - 0.0416666666666666*G0_2_1_1_0 - 0.0416666666666666*G0_2_1_1_1 - 0.0833333333333332*G0_2_2_1_0 - 0.0833333333333332*G0_2_2_1_1 - 0.5*G1_1_0 - 0.5*G1_1_1 - 0.166666666666666*G2_0_1_0_0 - 0.166666666666666*G2_0_1_0_1 + 0.166666666666666*G2_0_1_1_0 + 0.166666666666666*G2_0_1_2_1 - 0.0833333333333332*G2_1_1_0_0 - 0.0833333333333332*G2_1_1_0_1 + 0.0833333333333332*G2_1_1_1_0 + 0.0833333333333332*G2_1_1_2_1 - 0.0833333333333332*G2_2_1_0_0 - 0.0833333333333332*G2_2_1_0_1 + 0.0833333333333332*G2_2_1_1_0 + 0.0833333333333332*G2_2_1_2_1;
 
1898
    A[7] = 0.0833333333333332*G0_0_0_1_0 + 0.0416666666666666*G0_0_1_1_0 + 0.0416666666666666*G0_0_2_1_0 + 0.0416666666666666*G0_1_0_1_0 + 0.0833333333333332*G0_1_1_1_0 + 0.0416666666666666*G0_1_2_1_0 + 0.0416666666666666*G0_2_0_1_0 + 0.0416666666666666*G0_2_1_1_0 + 0.0833333333333332*G0_2_2_1_0 + 0.5*G1_1_0 - 0.0833333333333332*G2_0_1_0_0 - 0.0833333333333332*G2_0_1_0_1 + 0.0833333333333332*G2_0_1_1_0 + 0.0833333333333332*G2_0_1_2_1 - 0.166666666666666*G2_1_1_0_0 - 0.166666666666666*G2_1_1_0_1 + 0.166666666666666*G2_1_1_1_0 + 0.166666666666666*G2_1_1_2_1 - 0.0833333333333332*G2_2_1_0_0 - 0.0833333333333332*G2_2_1_0_1 + 0.0833333333333332*G2_2_1_1_0 + 0.0833333333333332*G2_2_1_2_1;
 
1899
    A[8] = 0.0833333333333332*G0_0_0_1_1 + 0.0416666666666666*G0_0_1_1_1 + 0.0416666666666666*G0_0_2_1_1 + 0.0416666666666666*G0_1_0_1_1 + 0.0833333333333332*G0_1_1_1_1 + 0.0416666666666666*G0_1_2_1_1 + 0.0416666666666666*G0_2_0_1_1 + 0.0416666666666666*G0_2_1_1_1 + 0.0833333333333332*G0_2_2_1_1 + 0.5*G1_1_1 - 0.0833333333333332*G2_0_1_0_0 - 0.0833333333333332*G2_0_1_0_1 + 0.0833333333333332*G2_0_1_1_0 + 0.0833333333333332*G2_0_1_2_1 - 0.0833333333333332*G2_1_1_0_0 - 0.0833333333333332*G2_1_1_0_1 + 0.0833333333333332*G2_1_1_1_0 + 0.0833333333333332*G2_1_1_2_1 - 0.166666666666666*G2_2_1_0_0 - 0.166666666666666*G2_2_1_0_1 + 0.166666666666666*G2_2_1_1_0 + 0.166666666666666*G2_2_1_2_1;
 
1900
  }
 
1901
 
 
1902
};
 
1903
 
 
1904
/// This class defines the interface for the assembly of the global
 
1905
/// tensor corresponding to a form with r + n arguments, that is, a
 
1906
/// mapping
 
1907
///
 
1908
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
1909
///
 
1910
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
1911
/// global tensor A is defined by
 
1912
///
 
1913
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
1914
///
 
1915
/// where each argument Vj represents the application to the
 
1916
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
1917
/// fixed functions (coefficients).
 
1918
 
 
1919
class UFC_NonlinearPoissonBilinearForm: public ufc::form
 
1920
{
 
1921
public:
 
1922
 
 
1923
  /// Constructor
 
1924
  UFC_NonlinearPoissonBilinearForm() : ufc::form()
 
1925
  {
 
1926
    // Do nothing
 
1927
  }
 
1928
 
 
1929
  /// Destructor
 
1930
  virtual ~UFC_NonlinearPoissonBilinearForm()
 
1931
  {
 
1932
    // Do nothing
 
1933
  }
 
1934
 
 
1935
  /// Return a string identifying the form
 
1936
  virtual const char* signature() const
 
1937
  {
 
1938
    return "w0_a0w0_a1(dXa2/dxb0)(dXa3/dxb0) | ((d/dXa2)vi0)*va0*va1*((d/dXa3)vi1)*dX(0) + (dXa0/dxb0)(dXa1/dxb0) | ((d/dXa0)vi0)*((d/dXa1)vi1)*dX(0) + 2.0w0_a0w0_a2(dXa1/dxb0)(dXa3/dxb0) | ((d/dXa1)vi0)*va0*vi1*((d/dXa3)va2)*dX(0)";
 
1939
  }
 
1940
 
 
1941
  /// Return the rank of the global tensor (r)
 
1942
  virtual unsigned int rank() const
 
1943
  {
 
1944
    return 2;
 
1945
  }
 
1946
 
 
1947
  /// Return the number of coefficients (n)
 
1948
  virtual unsigned int num_coefficients() const
 
1949
  {
 
1950
    return 1;
 
1951
  }
 
1952
 
 
1953
  /// Return the number of cell integrals
 
1954
  virtual unsigned int num_cell_integrals() const
 
1955
  {
 
1956
    return 1;
 
1957
  }
 
1958
  
 
1959
  /// Return the number of exterior facet integrals
 
1960
  virtual unsigned int num_exterior_facet_integrals() const
 
1961
  {
 
1962
    return 0;
 
1963
  }
 
1964
  
 
1965
  /// Return the number of interior facet integrals
 
1966
  virtual unsigned int num_interior_facet_integrals() const
 
1967
  {
 
1968
    return 0;
 
1969
  }
 
1970
    
 
1971
  /// Create a new finite element for argument function i
 
1972
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
1973
  {
 
1974
    switch ( i )
 
1975
    {
 
1976
    case 0:
 
1977
      return new UFC_NonlinearPoissonBilinearForm_finite_element_0();
 
1978
      break;
 
1979
    case 1:
 
1980
      return new UFC_NonlinearPoissonBilinearForm_finite_element_1();
 
1981
      break;
 
1982
    case 2:
 
1983
      return new UFC_NonlinearPoissonBilinearForm_finite_element_2();
 
1984
      break;
 
1985
    }
 
1986
    return 0;
 
1987
  }
 
1988
  
 
1989
  /// Create a new dof map for argument function i
 
1990
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
1991
  {
 
1992
    switch ( i )
 
1993
    {
 
1994
    case 0:
 
1995
      return new UFC_NonlinearPoissonBilinearForm_dof_map_0();
 
1996
      break;
 
1997
    case 1:
 
1998
      return new UFC_NonlinearPoissonBilinearForm_dof_map_1();
 
1999
      break;
 
2000
    case 2:
 
2001
      return new UFC_NonlinearPoissonBilinearForm_dof_map_2();
 
2002
      break;
 
2003
    }
 
2004
    return 0;
 
2005
  }
 
2006
 
 
2007
  /// Create a new cell integral on sub domain i
 
2008
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
2009
  {
 
2010
    return new UFC_NonlinearPoissonBilinearForm_cell_integral_0();
 
2011
  }
 
2012
 
 
2013
  /// Create a new exterior facet integral on sub domain i
 
2014
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
2015
  {
 
2016
    return 0;
 
2017
  }
 
2018
 
 
2019
  /// Create a new interior facet integral on sub domain i
 
2020
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
2021
  {
 
2022
    return 0;
 
2023
  }
 
2024
 
 
2025
};
 
2026
 
 
2027
/// This class defines the interface for a finite element.
 
2028
 
 
2029
class UFC_NonlinearPoissonLinearForm_finite_element_0: public ufc::finite_element
 
2030
{
 
2031
public:
 
2032
 
 
2033
  /// Constructor
 
2034
  UFC_NonlinearPoissonLinearForm_finite_element_0() : ufc::finite_element()
 
2035
  {
 
2036
    // Do nothing
 
2037
  }
 
2038
 
 
2039
  /// Destructor
 
2040
  virtual ~UFC_NonlinearPoissonLinearForm_finite_element_0()
 
2041
  {
 
2042
    // Do nothing
 
2043
  }
 
2044
 
 
2045
  /// Return a string identifying the finite element
 
2046
  virtual const char* signature() const
 
2047
  {
 
2048
    return "Lagrange finite element of degree 1 on a triangle";
 
2049
  }
 
2050
 
 
2051
  /// Return the cell shape
 
2052
  virtual ufc::shape cell_shape() const
 
2053
  {
 
2054
    return ufc::triangle;
 
2055
  }
 
2056
 
 
2057
  /// Return the dimension of the finite element function space
 
2058
  virtual unsigned int space_dimension() const
 
2059
  {
 
2060
    return 3;
 
2061
  }
 
2062
 
 
2063
  /// Return the rank of the value space
 
2064
  virtual unsigned int value_rank() const
 
2065
  {
 
2066
    return 0;
 
2067
  }
 
2068
 
 
2069
  /// Return the dimension of the value space for axis i
 
2070
  virtual unsigned int value_dimension(unsigned int i) const
 
2071
  {
 
2072
    return 1;
 
2073
  }
 
2074
 
 
2075
  /// Evaluate basis function i at given point in cell
 
2076
  virtual void evaluate_basis(unsigned int i,
 
2077
                              double* values,
 
2078
                              const double* coordinates,
 
2079
                              const ufc::cell& c) const
 
2080
  {
 
2081
    // Extract vertex coordinates
 
2082
    const double * const * element_coordinates = c.coordinates;
 
2083
    
 
2084
    // Compute Jacobian of affine map from reference cell
 
2085
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
2086
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
2087
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
2088
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
2089
      
 
2090
    // Compute determinant of Jacobian
 
2091
    const double detJ = J_00*J_11 - J_01*J_10;
 
2092
    
 
2093
    // Compute inverse of Jacobian
 
2094
    
 
2095
    // Get coordinates and map to the reference (UFC) element
 
2096
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
2097
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
2098
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
2099
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
2100
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
2101
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
2102
    
 
2103
    // Map coordinates to the reference square
 
2104
    if (std::abs(y - 1.0) < 1e-14)
 
2105
      x = -1.0;
 
2106
    else
 
2107
      x = 2.0 *x/(1.0 - y) - 1.0;
 
2108
    y = 2.0*y - 1.0;
 
2109
    
 
2110
    // Reset values
 
2111
    *values = 0;
 
2112
    
 
2113
    // Map degree of freedom to element degree of freedom
 
2114
    const unsigned int dof = i;
 
2115
    
 
2116
    // Generate scalings
 
2117
    const double scalings_y_0 = 1;
 
2118
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2119
    
 
2120
    // Compute psitilde_a
 
2121
    const double psitilde_a_0 = 1;
 
2122
    const double psitilde_a_1 = x;
 
2123
    
 
2124
    // Compute psitilde_bs
 
2125
    const double psitilde_bs_0_0 = 1;
 
2126
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2127
    const double psitilde_bs_1_0 = 1;
 
2128
    
 
2129
    // Compute basisvalues
 
2130
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
2131
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
2132
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
2133
    
 
2134
    // Table(s) of coefficients
 
2135
    const static double coefficients0[3][3] = \
 
2136
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
2137
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
2138
    {0.471404520791032, 0, 0.333333333333333}};
 
2139
    
 
2140
    // Extract relevant coefficients
 
2141
    const double coeff0_0 = coefficients0[dof][0];
 
2142
    const double coeff0_1 = coefficients0[dof][1];
 
2143
    const double coeff0_2 = coefficients0[dof][2];
 
2144
    
 
2145
    // Compute value(s)
 
2146
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
2147
  }
 
2148
 
 
2149
  /// Evaluate all basis functions at given point in cell
 
2150
  virtual void evaluate_basis_all(double* values,
 
2151
                                  const double* coordinates,
 
2152
                                  const ufc::cell& c) const
 
2153
  {
 
2154
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
2155
  }
 
2156
 
 
2157
  /// Evaluate order n derivatives of basis function i at given point in cell
 
2158
  virtual void evaluate_basis_derivatives(unsigned int i,
 
2159
                                          unsigned int n,
 
2160
                                          double* values,
 
2161
                                          const double* coordinates,
 
2162
                                          const ufc::cell& c) const
 
2163
  {
 
2164
    // Extract vertex coordinates
 
2165
    const double * const * element_coordinates = c.coordinates;
 
2166
    
 
2167
    // Compute Jacobian of affine map from reference cell
 
2168
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
2169
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
2170
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
2171
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
2172
      
 
2173
    // Compute determinant of Jacobian
 
2174
    const double detJ = J_00*J_11 - J_01*J_10;
 
2175
    
 
2176
    // Compute inverse of Jacobian
 
2177
    
 
2178
    // Get coordinates and map to the reference (UFC) element
 
2179
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
2180
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
2181
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
2182
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
2183
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
2184
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
2185
    
 
2186
    // Map coordinates to the reference square
 
2187
    if (std::abs(y - 1.0) < 1e-14)
 
2188
      x = -1.0;
 
2189
    else
 
2190
      x = 2.0 *x/(1.0 - y) - 1.0;
 
2191
    y = 2.0*y - 1.0;
 
2192
    
 
2193
    // Compute number of derivatives
 
2194
    unsigned int num_derivatives = 1;
 
2195
    
 
2196
    for (unsigned int j = 0; j < n; j++)
 
2197
      num_derivatives *= 2;
 
2198
    
 
2199
    
 
2200
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
2201
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
2202
        
 
2203
    for (unsigned int j = 0; j < num_derivatives; j++)
 
2204
    {
 
2205
      combinations[j] = new unsigned int [n];
 
2206
      for (unsigned int k = 0; k < n; k++)
 
2207
        combinations[j][k] = 0;
 
2208
    }
 
2209
        
 
2210
    // Generate combinations of derivatives
 
2211
    for (unsigned int row = 1; row < num_derivatives; row++)
 
2212
    {
 
2213
      for (unsigned int num = 0; num < row; num++)
 
2214
      {
 
2215
        for (unsigned int col = n-1; col+1 > 0; col--)
 
2216
        {
 
2217
          if (combinations[row][col] + 1 > 1)
 
2218
            combinations[row][col] = 0;
 
2219
          else
 
2220
          {
 
2221
            combinations[row][col] += 1;
 
2222
            break;
 
2223
          }
 
2224
        }
 
2225
      }
 
2226
    }
 
2227
    
 
2228
    // Compute inverse of Jacobian
 
2229
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
2230
    
 
2231
    // Declare transformation matrix
 
2232
    // Declare pointer to two dimensional array and initialise
 
2233
    double **transform = new double *[num_derivatives];
 
2234
        
 
2235
    for (unsigned int j = 0; j < num_derivatives; j++)
 
2236
    {
 
2237
      transform[j] = new double [num_derivatives];
 
2238
      for (unsigned int k = 0; k < num_derivatives; k++)
 
2239
        transform[j][k] = 1;
 
2240
    }
 
2241
    
 
2242
    // Construct transformation matrix
 
2243
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2244
    {
 
2245
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2246
      {
 
2247
        for (unsigned int k = 0; k < n; k++)
 
2248
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
2249
      }
 
2250
    }
 
2251
    
 
2252
    // Reset values
 
2253
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
2254
      values[j] = 0;
 
2255
    
 
2256
    // Map degree of freedom to element degree of freedom
 
2257
    const unsigned int dof = i;
 
2258
    
 
2259
    // Generate scalings
 
2260
    const double scalings_y_0 = 1;
 
2261
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2262
    
 
2263
    // Compute psitilde_a
 
2264
    const double psitilde_a_0 = 1;
 
2265
    const double psitilde_a_1 = x;
 
2266
    
 
2267
    // Compute psitilde_bs
 
2268
    const double psitilde_bs_0_0 = 1;
 
2269
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2270
    const double psitilde_bs_1_0 = 1;
 
2271
    
 
2272
    // Compute basisvalues
 
2273
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
2274
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
2275
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
2276
    
 
2277
    // Table(s) of coefficients
 
2278
    const static double coefficients0[3][3] = \
 
2279
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
2280
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
2281
    {0.471404520791032, 0, 0.333333333333333}};
 
2282
    
 
2283
    // Interesting (new) part
 
2284
    // Tables of derivatives of the polynomial base (transpose)
 
2285
    const static double dmats0[3][3] = \
 
2286
    {{0, 0, 0},
 
2287
    {4.89897948556636, 0, 0},
 
2288
    {0, 0, 0}};
 
2289
    
 
2290
    const static double dmats1[3][3] = \
 
2291
    {{0, 0, 0},
 
2292
    {2.44948974278318, 0, 0},
 
2293
    {4.24264068711928, 0, 0}};
 
2294
    
 
2295
    // Compute reference derivatives
 
2296
    // Declare pointer to array of derivatives on FIAT element
 
2297
    double *derivatives = new double [num_derivatives];
 
2298
    
 
2299
    // Declare coefficients
 
2300
    double coeff0_0 = 0;
 
2301
    double coeff0_1 = 0;
 
2302
    double coeff0_2 = 0;
 
2303
    
 
2304
    // Declare new coefficients
 
2305
    double new_coeff0_0 = 0;
 
2306
    double new_coeff0_1 = 0;
 
2307
    double new_coeff0_2 = 0;
 
2308
    
 
2309
    // Loop possible derivatives
 
2310
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
2311
    {
 
2312
      // Get values from coefficients array
 
2313
      new_coeff0_0 = coefficients0[dof][0];
 
2314
      new_coeff0_1 = coefficients0[dof][1];
 
2315
      new_coeff0_2 = coefficients0[dof][2];
 
2316
    
 
2317
      // Loop derivative order
 
2318
      for (unsigned int j = 0; j < n; j++)
 
2319
      {
 
2320
        // Update old coefficients
 
2321
        coeff0_0 = new_coeff0_0;
 
2322
        coeff0_1 = new_coeff0_1;
 
2323
        coeff0_2 = new_coeff0_2;
 
2324
    
 
2325
        if(combinations[deriv_num][j] == 0)
 
2326
        {
 
2327
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
2328
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
2329
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
2330
        }
 
2331
        if(combinations[deriv_num][j] == 1)
 
2332
        {
 
2333
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
2334
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
2335
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
2336
        }
 
2337
    
 
2338
      }
 
2339
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
2340
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
2341
    }
 
2342
    
 
2343
    // Transform derivatives back to physical element
 
2344
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2345
    {
 
2346
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2347
      {
 
2348
        values[row] += transform[row][col]*derivatives[col];
 
2349
      }
 
2350
    }
 
2351
    // Delete pointer to array of derivatives on FIAT element
 
2352
    delete [] derivatives;
 
2353
    
 
2354
    // Delete pointer to array of combinations of derivatives and transform
 
2355
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2356
    {
 
2357
      delete [] combinations[row];
 
2358
      delete [] transform[row];
 
2359
    }
 
2360
    
 
2361
    delete [] combinations;
 
2362
    delete [] transform;
 
2363
  }
 
2364
 
 
2365
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
2366
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
2367
                                              double* values,
 
2368
                                              const double* coordinates,
 
2369
                                              const ufc::cell& c) const
 
2370
  {
 
2371
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
2372
  }
 
2373
 
 
2374
  /// Evaluate linear functional for dof i on the function f
 
2375
  virtual double evaluate_dof(unsigned int i,
 
2376
                              const ufc::function& f,
 
2377
                              const ufc::cell& c) const
 
2378
  {
 
2379
    // The reference points, direction and weights:
 
2380
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
2381
    const static double W[3][1] = {{1}, {1}, {1}};
 
2382
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
2383
    
 
2384
    const double * const * x = c.coordinates;
 
2385
    double result = 0.0;
 
2386
    // Iterate over the points:
 
2387
    // Evaluate basis functions for affine mapping
 
2388
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
2389
    const double w1 = X[i][0][0];
 
2390
    const double w2 = X[i][0][1];
 
2391
    
 
2392
    // Compute affine mapping y = F(X)
 
2393
    double y[2];
 
2394
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
2395
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
2396
    
 
2397
    // Evaluate function at physical points
 
2398
    double values[1];
 
2399
    f.evaluate(values, y, c);
 
2400
    
 
2401
    // Map function values using appropriate mapping
 
2402
    // Affine map: Do nothing
 
2403
    
 
2404
    // Note that we do not map the weights (yet).
 
2405
    
 
2406
    // Take directional components
 
2407
    for(int k = 0; k < 1; k++)
 
2408
      result += values[k]*D[i][0][k];
 
2409
    // Multiply by weights 
 
2410
    result *= W[i][0];
 
2411
    
 
2412
    return result;
 
2413
  }
 
2414
 
 
2415
  /// Evaluate linear functionals for all dofs on the function f
 
2416
  virtual void evaluate_dofs(double* values,
 
2417
                             const ufc::function& f,
 
2418
                             const ufc::cell& c) const
 
2419
  {
 
2420
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2421
  }
 
2422
 
 
2423
  /// Interpolate vertex values from dof values
 
2424
  virtual void interpolate_vertex_values(double* vertex_values,
 
2425
                                         const double* dof_values,
 
2426
                                         const ufc::cell& c) const
 
2427
  {
 
2428
    // Evaluate at vertices and use affine mapping
 
2429
    vertex_values[0] = dof_values[0];
 
2430
    vertex_values[1] = dof_values[1];
 
2431
    vertex_values[2] = dof_values[2];
 
2432
  }
 
2433
 
 
2434
  /// Return the number of sub elements (for a mixed element)
 
2435
  virtual unsigned int num_sub_elements() const
 
2436
  {
 
2437
    return 1;
 
2438
  }
 
2439
 
 
2440
  /// Create a new finite element for sub element i (for a mixed element)
 
2441
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
2442
  {
 
2443
    return new UFC_NonlinearPoissonLinearForm_finite_element_0();
 
2444
  }
 
2445
 
 
2446
};
 
2447
 
 
2448
/// This class defines the interface for a finite element.
 
2449
 
 
2450
class UFC_NonlinearPoissonLinearForm_finite_element_1: public ufc::finite_element
 
2451
{
 
2452
public:
 
2453
 
 
2454
  /// Constructor
 
2455
  UFC_NonlinearPoissonLinearForm_finite_element_1() : ufc::finite_element()
 
2456
  {
 
2457
    // Do nothing
 
2458
  }
 
2459
 
 
2460
  /// Destructor
 
2461
  virtual ~UFC_NonlinearPoissonLinearForm_finite_element_1()
 
2462
  {
 
2463
    // Do nothing
 
2464
  }
 
2465
 
 
2466
  /// Return a string identifying the finite element
 
2467
  virtual const char* signature() const
 
2468
  {
 
2469
    return "Lagrange finite element of degree 1 on a triangle";
 
2470
  }
 
2471
 
 
2472
  /// Return the cell shape
 
2473
  virtual ufc::shape cell_shape() const
 
2474
  {
 
2475
    return ufc::triangle;
 
2476
  }
 
2477
 
 
2478
  /// Return the dimension of the finite element function space
 
2479
  virtual unsigned int space_dimension() const
 
2480
  {
 
2481
    return 3;
 
2482
  }
 
2483
 
 
2484
  /// Return the rank of the value space
 
2485
  virtual unsigned int value_rank() const
 
2486
  {
 
2487
    return 0;
 
2488
  }
 
2489
 
 
2490
  /// Return the dimension of the value space for axis i
 
2491
  virtual unsigned int value_dimension(unsigned int i) const
 
2492
  {
 
2493
    return 1;
 
2494
  }
 
2495
 
 
2496
  /// Evaluate basis function i at given point in cell
 
2497
  virtual void evaluate_basis(unsigned int i,
 
2498
                              double* values,
 
2499
                              const double* coordinates,
 
2500
                              const ufc::cell& c) const
 
2501
  {
 
2502
    // Extract vertex coordinates
 
2503
    const double * const * element_coordinates = c.coordinates;
 
2504
    
 
2505
    // Compute Jacobian of affine map from reference cell
 
2506
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
2507
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
2508
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
2509
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
2510
      
 
2511
    // Compute determinant of Jacobian
 
2512
    const double detJ = J_00*J_11 - J_01*J_10;
 
2513
    
 
2514
    // Compute inverse of Jacobian
 
2515
    
 
2516
    // Get coordinates and map to the reference (UFC) element
 
2517
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
2518
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
2519
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
2520
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
2521
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
2522
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
2523
    
 
2524
    // Map coordinates to the reference square
 
2525
    if (std::abs(y - 1.0) < 1e-14)
 
2526
      x = -1.0;
 
2527
    else
 
2528
      x = 2.0 *x/(1.0 - y) - 1.0;
 
2529
    y = 2.0*y - 1.0;
 
2530
    
 
2531
    // Reset values
 
2532
    *values = 0;
 
2533
    
 
2534
    // Map degree of freedom to element degree of freedom
 
2535
    const unsigned int dof = i;
 
2536
    
 
2537
    // Generate scalings
 
2538
    const double scalings_y_0 = 1;
 
2539
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2540
    
 
2541
    // Compute psitilde_a
 
2542
    const double psitilde_a_0 = 1;
 
2543
    const double psitilde_a_1 = x;
 
2544
    
 
2545
    // Compute psitilde_bs
 
2546
    const double psitilde_bs_0_0 = 1;
 
2547
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2548
    const double psitilde_bs_1_0 = 1;
 
2549
    
 
2550
    // Compute basisvalues
 
2551
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
2552
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
2553
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
2554
    
 
2555
    // Table(s) of coefficients
 
2556
    const static double coefficients0[3][3] = \
 
2557
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
2558
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
2559
    {0.471404520791032, 0, 0.333333333333333}};
 
2560
    
 
2561
    // Extract relevant coefficients
 
2562
    const double coeff0_0 = coefficients0[dof][0];
 
2563
    const double coeff0_1 = coefficients0[dof][1];
 
2564
    const double coeff0_2 = coefficients0[dof][2];
 
2565
    
 
2566
    // Compute value(s)
 
2567
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
2568
  }
 
2569
 
 
2570
  /// Evaluate all basis functions at given point in cell
 
2571
  virtual void evaluate_basis_all(double* values,
 
2572
                                  const double* coordinates,
 
2573
                                  const ufc::cell& c) const
 
2574
  {
 
2575
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
2576
  }
 
2577
 
 
2578
  /// Evaluate order n derivatives of basis function i at given point in cell
 
2579
  virtual void evaluate_basis_derivatives(unsigned int i,
 
2580
                                          unsigned int n,
 
2581
                                          double* values,
 
2582
                                          const double* coordinates,
 
2583
                                          const ufc::cell& c) const
 
2584
  {
 
2585
    // Extract vertex coordinates
 
2586
    const double * const * element_coordinates = c.coordinates;
 
2587
    
 
2588
    // Compute Jacobian of affine map from reference cell
 
2589
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
2590
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
2591
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
2592
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
2593
      
 
2594
    // Compute determinant of Jacobian
 
2595
    const double detJ = J_00*J_11 - J_01*J_10;
 
2596
    
 
2597
    // Compute inverse of Jacobian
 
2598
    
 
2599
    // Get coordinates and map to the reference (UFC) element
 
2600
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
2601
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
2602
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
2603
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
2604
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
2605
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
2606
    
 
2607
    // Map coordinates to the reference square
 
2608
    if (std::abs(y - 1.0) < 1e-14)
 
2609
      x = -1.0;
 
2610
    else
 
2611
      x = 2.0 *x/(1.0 - y) - 1.0;
 
2612
    y = 2.0*y - 1.0;
 
2613
    
 
2614
    // Compute number of derivatives
 
2615
    unsigned int num_derivatives = 1;
 
2616
    
 
2617
    for (unsigned int j = 0; j < n; j++)
 
2618
      num_derivatives *= 2;
 
2619
    
 
2620
    
 
2621
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
2622
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
2623
        
 
2624
    for (unsigned int j = 0; j < num_derivatives; j++)
 
2625
    {
 
2626
      combinations[j] = new unsigned int [n];
 
2627
      for (unsigned int k = 0; k < n; k++)
 
2628
        combinations[j][k] = 0;
 
2629
    }
 
2630
        
 
2631
    // Generate combinations of derivatives
 
2632
    for (unsigned int row = 1; row < num_derivatives; row++)
 
2633
    {
 
2634
      for (unsigned int num = 0; num < row; num++)
 
2635
      {
 
2636
        for (unsigned int col = n-1; col+1 > 0; col--)
 
2637
        {
 
2638
          if (combinations[row][col] + 1 > 1)
 
2639
            combinations[row][col] = 0;
 
2640
          else
 
2641
          {
 
2642
            combinations[row][col] += 1;
 
2643
            break;
 
2644
          }
 
2645
        }
 
2646
      }
 
2647
    }
 
2648
    
 
2649
    // Compute inverse of Jacobian
 
2650
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
2651
    
 
2652
    // Declare transformation matrix
 
2653
    // Declare pointer to two dimensional array and initialise
 
2654
    double **transform = new double *[num_derivatives];
 
2655
        
 
2656
    for (unsigned int j = 0; j < num_derivatives; j++)
 
2657
    {
 
2658
      transform[j] = new double [num_derivatives];
 
2659
      for (unsigned int k = 0; k < num_derivatives; k++)
 
2660
        transform[j][k] = 1;
 
2661
    }
 
2662
    
 
2663
    // Construct transformation matrix
 
2664
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2665
    {
 
2666
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2667
      {
 
2668
        for (unsigned int k = 0; k < n; k++)
 
2669
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
2670
      }
 
2671
    }
 
2672
    
 
2673
    // Reset values
 
2674
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
2675
      values[j] = 0;
 
2676
    
 
2677
    // Map degree of freedom to element degree of freedom
 
2678
    const unsigned int dof = i;
 
2679
    
 
2680
    // Generate scalings
 
2681
    const double scalings_y_0 = 1;
 
2682
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2683
    
 
2684
    // Compute psitilde_a
 
2685
    const double psitilde_a_0 = 1;
 
2686
    const double psitilde_a_1 = x;
 
2687
    
 
2688
    // Compute psitilde_bs
 
2689
    const double psitilde_bs_0_0 = 1;
 
2690
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2691
    const double psitilde_bs_1_0 = 1;
 
2692
    
 
2693
    // Compute basisvalues
 
2694
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
2695
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
2696
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
2697
    
 
2698
    // Table(s) of coefficients
 
2699
    const static double coefficients0[3][3] = \
 
2700
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
2701
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
2702
    {0.471404520791032, 0, 0.333333333333333}};
 
2703
    
 
2704
    // Interesting (new) part
 
2705
    // Tables of derivatives of the polynomial base (transpose)
 
2706
    const static double dmats0[3][3] = \
 
2707
    {{0, 0, 0},
 
2708
    {4.89897948556636, 0, 0},
 
2709
    {0, 0, 0}};
 
2710
    
 
2711
    const static double dmats1[3][3] = \
 
2712
    {{0, 0, 0},
 
2713
    {2.44948974278318, 0, 0},
 
2714
    {4.24264068711928, 0, 0}};
 
2715
    
 
2716
    // Compute reference derivatives
 
2717
    // Declare pointer to array of derivatives on FIAT element
 
2718
    double *derivatives = new double [num_derivatives];
 
2719
    
 
2720
    // Declare coefficients
 
2721
    double coeff0_0 = 0;
 
2722
    double coeff0_1 = 0;
 
2723
    double coeff0_2 = 0;
 
2724
    
 
2725
    // Declare new coefficients
 
2726
    double new_coeff0_0 = 0;
 
2727
    double new_coeff0_1 = 0;
 
2728
    double new_coeff0_2 = 0;
 
2729
    
 
2730
    // Loop possible derivatives
 
2731
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
2732
    {
 
2733
      // Get values from coefficients array
 
2734
      new_coeff0_0 = coefficients0[dof][0];
 
2735
      new_coeff0_1 = coefficients0[dof][1];
 
2736
      new_coeff0_2 = coefficients0[dof][2];
 
2737
    
 
2738
      // Loop derivative order
 
2739
      for (unsigned int j = 0; j < n; j++)
 
2740
      {
 
2741
        // Update old coefficients
 
2742
        coeff0_0 = new_coeff0_0;
 
2743
        coeff0_1 = new_coeff0_1;
 
2744
        coeff0_2 = new_coeff0_2;
 
2745
    
 
2746
        if(combinations[deriv_num][j] == 0)
 
2747
        {
 
2748
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
2749
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
2750
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
2751
        }
 
2752
        if(combinations[deriv_num][j] == 1)
 
2753
        {
 
2754
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
2755
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
2756
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
2757
        }
 
2758
    
 
2759
      }
 
2760
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
2761
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
2762
    }
 
2763
    
 
2764
    // Transform derivatives back to physical element
 
2765
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2766
    {
 
2767
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2768
      {
 
2769
        values[row] += transform[row][col]*derivatives[col];
 
2770
      }
 
2771
    }
 
2772
    // Delete pointer to array of derivatives on FIAT element
 
2773
    delete [] derivatives;
 
2774
    
 
2775
    // Delete pointer to array of combinations of derivatives and transform
 
2776
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2777
    {
 
2778
      delete [] combinations[row];
 
2779
      delete [] transform[row];
 
2780
    }
 
2781
    
 
2782
    delete [] combinations;
 
2783
    delete [] transform;
 
2784
  }
 
2785
 
 
2786
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
2787
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
2788
                                              double* values,
 
2789
                                              const double* coordinates,
 
2790
                                              const ufc::cell& c) const
 
2791
  {
 
2792
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
2793
  }
 
2794
 
 
2795
  /// Evaluate linear functional for dof i on the function f
 
2796
  virtual double evaluate_dof(unsigned int i,
 
2797
                              const ufc::function& f,
 
2798
                              const ufc::cell& c) const
 
2799
  {
 
2800
    // The reference points, direction and weights:
 
2801
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
2802
    const static double W[3][1] = {{1}, {1}, {1}};
 
2803
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
2804
    
 
2805
    const double * const * x = c.coordinates;
 
2806
    double result = 0.0;
 
2807
    // Iterate over the points:
 
2808
    // Evaluate basis functions for affine mapping
 
2809
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
2810
    const double w1 = X[i][0][0];
 
2811
    const double w2 = X[i][0][1];
 
2812
    
 
2813
    // Compute affine mapping y = F(X)
 
2814
    double y[2];
 
2815
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
2816
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
2817
    
 
2818
    // Evaluate function at physical points
 
2819
    double values[1];
 
2820
    f.evaluate(values, y, c);
 
2821
    
 
2822
    // Map function values using appropriate mapping
 
2823
    // Affine map: Do nothing
 
2824
    
 
2825
    // Note that we do not map the weights (yet).
 
2826
    
 
2827
    // Take directional components
 
2828
    for(int k = 0; k < 1; k++)
 
2829
      result += values[k]*D[i][0][k];
 
2830
    // Multiply by weights 
 
2831
    result *= W[i][0];
 
2832
    
 
2833
    return result;
 
2834
  }
 
2835
 
 
2836
  /// Evaluate linear functionals for all dofs on the function f
 
2837
  virtual void evaluate_dofs(double* values,
 
2838
                             const ufc::function& f,
 
2839
                             const ufc::cell& c) const
 
2840
  {
 
2841
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2842
  }
 
2843
 
 
2844
  /// Interpolate vertex values from dof values
 
2845
  virtual void interpolate_vertex_values(double* vertex_values,
 
2846
                                         const double* dof_values,
 
2847
                                         const ufc::cell& c) const
 
2848
  {
 
2849
    // Evaluate at vertices and use affine mapping
 
2850
    vertex_values[0] = dof_values[0];
 
2851
    vertex_values[1] = dof_values[1];
 
2852
    vertex_values[2] = dof_values[2];
 
2853
  }
 
2854
 
 
2855
  /// Return the number of sub elements (for a mixed element)
 
2856
  virtual unsigned int num_sub_elements() const
 
2857
  {
 
2858
    return 1;
 
2859
  }
 
2860
 
 
2861
  /// Create a new finite element for sub element i (for a mixed element)
 
2862
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
2863
  {
 
2864
    return new UFC_NonlinearPoissonLinearForm_finite_element_1();
 
2865
  }
 
2866
 
 
2867
};
 
2868
 
 
2869
/// This class defines the interface for a finite element.
 
2870
 
 
2871
class UFC_NonlinearPoissonLinearForm_finite_element_2: public ufc::finite_element
 
2872
{
 
2873
public:
 
2874
 
 
2875
  /// Constructor
 
2876
  UFC_NonlinearPoissonLinearForm_finite_element_2() : ufc::finite_element()
 
2877
  {
 
2878
    // Do nothing
 
2879
  }
 
2880
 
 
2881
  /// Destructor
 
2882
  virtual ~UFC_NonlinearPoissonLinearForm_finite_element_2()
 
2883
  {
 
2884
    // Do nothing
 
2885
  }
 
2886
 
 
2887
  /// Return a string identifying the finite element
 
2888
  virtual const char* signature() const
 
2889
  {
 
2890
    return "Lagrange finite element of degree 1 on a triangle";
 
2891
  }
 
2892
 
 
2893
  /// Return the cell shape
 
2894
  virtual ufc::shape cell_shape() const
 
2895
  {
 
2896
    return ufc::triangle;
 
2897
  }
 
2898
 
 
2899
  /// Return the dimension of the finite element function space
 
2900
  virtual unsigned int space_dimension() const
 
2901
  {
 
2902
    return 3;
 
2903
  }
 
2904
 
 
2905
  /// Return the rank of the value space
 
2906
  virtual unsigned int value_rank() const
 
2907
  {
 
2908
    return 0;
 
2909
  }
 
2910
 
 
2911
  /// Return the dimension of the value space for axis i
 
2912
  virtual unsigned int value_dimension(unsigned int i) const
 
2913
  {
 
2914
    return 1;
 
2915
  }
 
2916
 
 
2917
  /// Evaluate basis function i at given point in cell
 
2918
  virtual void evaluate_basis(unsigned int i,
 
2919
                              double* values,
 
2920
                              const double* coordinates,
 
2921
                              const ufc::cell& c) const
 
2922
  {
 
2923
    // Extract vertex coordinates
 
2924
    const double * const * element_coordinates = c.coordinates;
 
2925
    
 
2926
    // Compute Jacobian of affine map from reference cell
 
2927
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
2928
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
2929
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
2930
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
2931
      
 
2932
    // Compute determinant of Jacobian
 
2933
    const double detJ = J_00*J_11 - J_01*J_10;
 
2934
    
 
2935
    // Compute inverse of Jacobian
 
2936
    
 
2937
    // Get coordinates and map to the reference (UFC) element
 
2938
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
2939
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
2940
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
2941
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
2942
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
2943
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
2944
    
 
2945
    // Map coordinates to the reference square
 
2946
    if (std::abs(y - 1.0) < 1e-14)
 
2947
      x = -1.0;
 
2948
    else
 
2949
      x = 2.0 *x/(1.0 - y) - 1.0;
 
2950
    y = 2.0*y - 1.0;
 
2951
    
 
2952
    // Reset values
 
2953
    *values = 0;
 
2954
    
 
2955
    // Map degree of freedom to element degree of freedom
 
2956
    const unsigned int dof = i;
 
2957
    
 
2958
    // Generate scalings
 
2959
    const double scalings_y_0 = 1;
 
2960
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2961
    
 
2962
    // Compute psitilde_a
 
2963
    const double psitilde_a_0 = 1;
 
2964
    const double psitilde_a_1 = x;
 
2965
    
 
2966
    // Compute psitilde_bs
 
2967
    const double psitilde_bs_0_0 = 1;
 
2968
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2969
    const double psitilde_bs_1_0 = 1;
 
2970
    
 
2971
    // Compute basisvalues
 
2972
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
2973
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
2974
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
2975
    
 
2976
    // Table(s) of coefficients
 
2977
    const static double coefficients0[3][3] = \
 
2978
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
2979
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
2980
    {0.471404520791032, 0, 0.333333333333333}};
 
2981
    
 
2982
    // Extract relevant coefficients
 
2983
    const double coeff0_0 = coefficients0[dof][0];
 
2984
    const double coeff0_1 = coefficients0[dof][1];
 
2985
    const double coeff0_2 = coefficients0[dof][2];
 
2986
    
 
2987
    // Compute value(s)
 
2988
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
 
2989
  }
 
2990
 
 
2991
  /// Evaluate all basis functions at given point in cell
 
2992
  virtual void evaluate_basis_all(double* values,
 
2993
                                  const double* coordinates,
 
2994
                                  const ufc::cell& c) const
 
2995
  {
 
2996
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
2997
  }
 
2998
 
 
2999
  /// Evaluate order n derivatives of basis function i at given point in cell
 
3000
  virtual void evaluate_basis_derivatives(unsigned int i,
 
3001
                                          unsigned int n,
 
3002
                                          double* values,
 
3003
                                          const double* coordinates,
 
3004
                                          const ufc::cell& c) const
 
3005
  {
 
3006
    // Extract vertex coordinates
 
3007
    const double * const * element_coordinates = c.coordinates;
 
3008
    
 
3009
    // Compute Jacobian of affine map from reference cell
 
3010
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
3011
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
3012
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
3013
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
3014
      
 
3015
    // Compute determinant of Jacobian
 
3016
    const double detJ = J_00*J_11 - J_01*J_10;
 
3017
    
 
3018
    // Compute inverse of Jacobian
 
3019
    
 
3020
    // Get coordinates and map to the reference (UFC) element
 
3021
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
 
3022
                element_coordinates[0][0]*element_coordinates[2][1] +\
 
3023
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
 
3024
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
 
3025
                element_coordinates[1][0]*element_coordinates[0][1] -\
 
3026
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
 
3027
    
 
3028
    // Map coordinates to the reference square
 
3029
    if (std::abs(y - 1.0) < 1e-14)
 
3030
      x = -1.0;
 
3031
    else
 
3032
      x = 2.0 *x/(1.0 - y) - 1.0;
 
3033
    y = 2.0*y - 1.0;
 
3034
    
 
3035
    // Compute number of derivatives
 
3036
    unsigned int num_derivatives = 1;
 
3037
    
 
3038
    for (unsigned int j = 0; j < n; j++)
 
3039
      num_derivatives *= 2;
 
3040
    
 
3041
    
 
3042
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
3043
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
3044
        
 
3045
    for (unsigned int j = 0; j < num_derivatives; j++)
 
3046
    {
 
3047
      combinations[j] = new unsigned int [n];
 
3048
      for (unsigned int k = 0; k < n; k++)
 
3049
        combinations[j][k] = 0;
 
3050
    }
 
3051
        
 
3052
    // Generate combinations of derivatives
 
3053
    for (unsigned int row = 1; row < num_derivatives; row++)
 
3054
    {
 
3055
      for (unsigned int num = 0; num < row; num++)
 
3056
      {
 
3057
        for (unsigned int col = n-1; col+1 > 0; col--)
 
3058
        {
 
3059
          if (combinations[row][col] + 1 > 1)
 
3060
            combinations[row][col] = 0;
 
3061
          else
 
3062
          {
 
3063
            combinations[row][col] += 1;
 
3064
            break;
 
3065
          }
 
3066
        }
 
3067
      }
 
3068
    }
 
3069
    
 
3070
    // Compute inverse of Jacobian
 
3071
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
 
3072
    
 
3073
    // Declare transformation matrix
 
3074
    // Declare pointer to two dimensional array and initialise
 
3075
    double **transform = new double *[num_derivatives];
 
3076
        
 
3077
    for (unsigned int j = 0; j < num_derivatives; j++)
 
3078
    {
 
3079
      transform[j] = new double [num_derivatives];
 
3080
      for (unsigned int k = 0; k < num_derivatives; k++)
 
3081
        transform[j][k] = 1;
 
3082
    }
 
3083
    
 
3084
    // Construct transformation matrix
 
3085
    for (unsigned int row = 0; row < num_derivatives; row++)
 
3086
    {
 
3087
      for (unsigned int col = 0; col < num_derivatives; col++)
 
3088
      {
 
3089
        for (unsigned int k = 0; k < n; k++)
 
3090
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
3091
      }
 
3092
    }
 
3093
    
 
3094
    // Reset values
 
3095
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
3096
      values[j] = 0;
 
3097
    
 
3098
    // Map degree of freedom to element degree of freedom
 
3099
    const unsigned int dof = i;
 
3100
    
 
3101
    // Generate scalings
 
3102
    const double scalings_y_0 = 1;
 
3103
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
3104
    
 
3105
    // Compute psitilde_a
 
3106
    const double psitilde_a_0 = 1;
 
3107
    const double psitilde_a_1 = x;
 
3108
    
 
3109
    // Compute psitilde_bs
 
3110
    const double psitilde_bs_0_0 = 1;
 
3111
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
3112
    const double psitilde_bs_1_0 = 1;
 
3113
    
 
3114
    // Compute basisvalues
 
3115
    const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
 
3116
    const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
 
3117
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
 
3118
    
 
3119
    // Table(s) of coefficients
 
3120
    const static double coefficients0[3][3] = \
 
3121
    {{0.471404520791032, -0.288675134594813, -0.166666666666667},
 
3122
    {0.471404520791032, 0.288675134594813, -0.166666666666667},
 
3123
    {0.471404520791032, 0, 0.333333333333333}};
 
3124
    
 
3125
    // Interesting (new) part
 
3126
    // Tables of derivatives of the polynomial base (transpose)
 
3127
    const static double dmats0[3][3] = \
 
3128
    {{0, 0, 0},
 
3129
    {4.89897948556636, 0, 0},
 
3130
    {0, 0, 0}};
 
3131
    
 
3132
    const static double dmats1[3][3] = \
 
3133
    {{0, 0, 0},
 
3134
    {2.44948974278318, 0, 0},
 
3135
    {4.24264068711928, 0, 0}};
 
3136
    
 
3137
    // Compute reference derivatives
 
3138
    // Declare pointer to array of derivatives on FIAT element
 
3139
    double *derivatives = new double [num_derivatives];
 
3140
    
 
3141
    // Declare coefficients
 
3142
    double coeff0_0 = 0;
 
3143
    double coeff0_1 = 0;
 
3144
    double coeff0_2 = 0;
 
3145
    
 
3146
    // Declare new coefficients
 
3147
    double new_coeff0_0 = 0;
 
3148
    double new_coeff0_1 = 0;
 
3149
    double new_coeff0_2 = 0;
 
3150
    
 
3151
    // Loop possible derivatives
 
3152
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
3153
    {
 
3154
      // Get values from coefficients array
 
3155
      new_coeff0_0 = coefficients0[dof][0];
 
3156
      new_coeff0_1 = coefficients0[dof][1];
 
3157
      new_coeff0_2 = coefficients0[dof][2];
 
3158
    
 
3159
      // Loop derivative order
 
3160
      for (unsigned int j = 0; j < n; j++)
 
3161
      {
 
3162
        // Update old coefficients
 
3163
        coeff0_0 = new_coeff0_0;
 
3164
        coeff0_1 = new_coeff0_1;
 
3165
        coeff0_2 = new_coeff0_2;
 
3166
    
 
3167
        if(combinations[deriv_num][j] == 0)
 
3168
        {
 
3169
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
 
3170
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
 
3171
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
 
3172
        }
 
3173
        if(combinations[deriv_num][j] == 1)
 
3174
        {
 
3175
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
 
3176
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
 
3177
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
 
3178
        }
 
3179
    
 
3180
      }
 
3181
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
3182
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
 
3183
    }
 
3184
    
 
3185
    // Transform derivatives back to physical element
 
3186
    for (unsigned int row = 0; row < num_derivatives; row++)
 
3187
    {
 
3188
      for (unsigned int col = 0; col < num_derivatives; col++)
 
3189
      {
 
3190
        values[row] += transform[row][col]*derivatives[col];
 
3191
      }
 
3192
    }
 
3193
    // Delete pointer to array of derivatives on FIAT element
 
3194
    delete [] derivatives;
 
3195
    
 
3196
    // Delete pointer to array of combinations of derivatives and transform
 
3197
    for (unsigned int row = 0; row < num_derivatives; row++)
 
3198
    {
 
3199
      delete [] combinations[row];
 
3200
      delete [] transform[row];
 
3201
    }
 
3202
    
 
3203
    delete [] combinations;
 
3204
    delete [] transform;
 
3205
  }
 
3206
 
 
3207
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
3208
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
3209
                                              double* values,
 
3210
                                              const double* coordinates,
 
3211
                                              const ufc::cell& c) const
 
3212
  {
 
3213
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
3214
  }
 
3215
 
 
3216
  /// Evaluate linear functional for dof i on the function f
 
3217
  virtual double evaluate_dof(unsigned int i,
 
3218
                              const ufc::function& f,
 
3219
                              const ufc::cell& c) const
 
3220
  {
 
3221
    // The reference points, direction and weights:
 
3222
    const static double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
 
3223
    const static double W[3][1] = {{1}, {1}, {1}};
 
3224
    const static double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
 
3225
    
 
3226
    const double * const * x = c.coordinates;
 
3227
    double result = 0.0;
 
3228
    // Iterate over the points:
 
3229
    // Evaluate basis functions for affine mapping
 
3230
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
 
3231
    const double w1 = X[i][0][0];
 
3232
    const double w2 = X[i][0][1];
 
3233
    
 
3234
    // Compute affine mapping y = F(X)
 
3235
    double y[2];
 
3236
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
 
3237
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
 
3238
    
 
3239
    // Evaluate function at physical points
 
3240
    double values[1];
 
3241
    f.evaluate(values, y, c);
 
3242
    
 
3243
    // Map function values using appropriate mapping
 
3244
    // Affine map: Do nothing
 
3245
    
 
3246
    // Note that we do not map the weights (yet).
 
3247
    
 
3248
    // Take directional components
 
3249
    for(int k = 0; k < 1; k++)
 
3250
      result += values[k]*D[i][0][k];
 
3251
    // Multiply by weights 
 
3252
    result *= W[i][0];
 
3253
    
 
3254
    return result;
 
3255
  }
 
3256
 
 
3257
  /// Evaluate linear functionals for all dofs on the function f
 
3258
  virtual void evaluate_dofs(double* values,
 
3259
                             const ufc::function& f,
 
3260
                             const ufc::cell& c) const
 
3261
  {
 
3262
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3263
  }
 
3264
 
 
3265
  /// Interpolate vertex values from dof values
 
3266
  virtual void interpolate_vertex_values(double* vertex_values,
 
3267
                                         const double* dof_values,
 
3268
                                         const ufc::cell& c) const
 
3269
  {
 
3270
    // Evaluate at vertices and use affine mapping
 
3271
    vertex_values[0] = dof_values[0];
 
3272
    vertex_values[1] = dof_values[1];
 
3273
    vertex_values[2] = dof_values[2];
 
3274
  }
 
3275
 
 
3276
  /// Return the number of sub elements (for a mixed element)
 
3277
  virtual unsigned int num_sub_elements() const
 
3278
  {
 
3279
    return 1;
 
3280
  }
 
3281
 
 
3282
  /// Create a new finite element for sub element i (for a mixed element)
 
3283
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
3284
  {
 
3285
    return new UFC_NonlinearPoissonLinearForm_finite_element_2();
 
3286
  }
 
3287
 
 
3288
};
 
3289
 
 
3290
/// This class defines the interface for a local-to-global mapping of
 
3291
/// degrees of freedom (dofs).
 
3292
 
 
3293
class UFC_NonlinearPoissonLinearForm_dof_map_0: public ufc::dof_map
 
3294
{
 
3295
private:
 
3296
 
 
3297
  unsigned int __global_dimension;
 
3298
 
 
3299
public:
 
3300
 
 
3301
  /// Constructor
 
3302
  UFC_NonlinearPoissonLinearForm_dof_map_0() : ufc::dof_map()
 
3303
  {
 
3304
    __global_dimension = 0;
 
3305
  }
 
3306
 
 
3307
  /// Destructor
 
3308
  virtual ~UFC_NonlinearPoissonLinearForm_dof_map_0()
 
3309
  {
 
3310
    // Do nothing
 
3311
  }
 
3312
 
 
3313
  /// Return a string identifying the dof map
 
3314
  virtual const char* signature() const
 
3315
  {
 
3316
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
 
3317
  }
 
3318
 
 
3319
  /// Return true iff mesh entities of topological dimension d are needed
 
3320
  virtual bool needs_mesh_entities(unsigned int d) const
 
3321
  {
 
3322
    switch ( d )
 
3323
    {
 
3324
    case 0:
 
3325
      return true;
 
3326
      break;
 
3327
    case 1:
 
3328
      return false;
 
3329
      break;
 
3330
    case 2:
 
3331
      return false;
 
3332
      break;
 
3333
    }
 
3334
    return false;
 
3335
  }
 
3336
 
 
3337
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
3338
  virtual bool init_mesh(const ufc::mesh& m)
 
3339
  {
 
3340
    __global_dimension = m.num_entities[0];
 
3341
    return false;
 
3342
  }
 
3343
 
 
3344
  /// Initialize dof map for given cell
 
3345
  virtual void init_cell(const ufc::mesh& m,
 
3346
                         const ufc::cell& c)
 
3347
  {
 
3348
    // Do nothing
 
3349
  }
 
3350
 
 
3351
  /// Finish initialization of dof map for cells
 
3352
  virtual void init_cell_finalize()
 
3353
  {
 
3354
    // Do nothing
 
3355
  }
 
3356
 
 
3357
  /// Return the dimension of the global finite element function space
 
3358
  virtual unsigned int global_dimension() const
 
3359
  {
 
3360
    return __global_dimension;
 
3361
  }
 
3362
 
 
3363
  /// Return the dimension of the local finite element function space
 
3364
  virtual unsigned int local_dimension() const
 
3365
  {
 
3366
    return 3;
 
3367
  }
 
3368
 
 
3369
  // Return the geometric dimension of the coordinates this dof map provides
 
3370
  virtual unsigned int geometric_dimension() const
 
3371
  {
 
3372
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3373
  }
 
3374
 
 
3375
  /// Return the number of dofs on each cell facet
 
3376
  virtual unsigned int num_facet_dofs() const
 
3377
  {
 
3378
    return 2;
 
3379
  }
 
3380
 
 
3381
  /// Return the number of dofs associated with each cell entity of dimension d
 
3382
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
3383
  {
 
3384
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3385
  }
 
3386
 
 
3387
  /// Tabulate the local-to-global mapping of dofs on a cell
 
3388
  virtual void tabulate_dofs(unsigned int* dofs,
 
3389
                             const ufc::mesh& m,
 
3390
                             const ufc::cell& c) const
 
3391
  {
 
3392
    dofs[0] = c.entity_indices[0][0];
 
3393
    dofs[1] = c.entity_indices[0][1];
 
3394
    dofs[2] = c.entity_indices[0][2];
 
3395
  }
 
3396
 
 
3397
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
3398
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
3399
                                   unsigned int facet) const
 
3400
  {
 
3401
    switch ( facet )
 
3402
    {
 
3403
    case 0:
 
3404
      dofs[0] = 1;
 
3405
      dofs[1] = 2;
 
3406
      break;
 
3407
    case 1:
 
3408
      dofs[0] = 0;
 
3409
      dofs[1] = 2;
 
3410
      break;
 
3411
    case 2:
 
3412
      dofs[0] = 0;
 
3413
      dofs[1] = 1;
 
3414
      break;
 
3415
    }
 
3416
  }
 
3417
 
 
3418
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
3419
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
3420
                                    unsigned int d, unsigned int i) const
 
3421
  {
 
3422
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3423
  }
 
3424
 
 
3425
  /// Tabulate the coordinates of all dofs on a cell
 
3426
  virtual void tabulate_coordinates(double** coordinates,
 
3427
                                    const ufc::cell& c) const
 
3428
  {
 
3429
    const double * const * x = c.coordinates;
 
3430
    coordinates[0][0] = x[0][0];
 
3431
    coordinates[0][1] = x[0][1];
 
3432
    coordinates[1][0] = x[1][0];
 
3433
    coordinates[1][1] = x[1][1];
 
3434
    coordinates[2][0] = x[2][0];
 
3435
    coordinates[2][1] = x[2][1];
 
3436
  }
 
3437
 
 
3438
  /// Return the number of sub dof maps (for a mixed element)
 
3439
  virtual unsigned int num_sub_dof_maps() const
 
3440
  {
 
3441
    return 1;
 
3442
  }
 
3443
 
 
3444
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
3445
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
3446
  {
 
3447
    return new UFC_NonlinearPoissonLinearForm_dof_map_0();
 
3448
  }
 
3449
 
 
3450
};
 
3451
 
 
3452
/// This class defines the interface for a local-to-global mapping of
 
3453
/// degrees of freedom (dofs).
 
3454
 
 
3455
class UFC_NonlinearPoissonLinearForm_dof_map_1: public ufc::dof_map
 
3456
{
 
3457
private:
 
3458
 
 
3459
  unsigned int __global_dimension;
 
3460
 
 
3461
public:
 
3462
 
 
3463
  /// Constructor
 
3464
  UFC_NonlinearPoissonLinearForm_dof_map_1() : ufc::dof_map()
 
3465
  {
 
3466
    __global_dimension = 0;
 
3467
  }
 
3468
 
 
3469
  /// Destructor
 
3470
  virtual ~UFC_NonlinearPoissonLinearForm_dof_map_1()
 
3471
  {
 
3472
    // Do nothing
 
3473
  }
 
3474
 
 
3475
  /// Return a string identifying the dof map
 
3476
  virtual const char* signature() const
 
3477
  {
 
3478
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
 
3479
  }
 
3480
 
 
3481
  /// Return true iff mesh entities of topological dimension d are needed
 
3482
  virtual bool needs_mesh_entities(unsigned int d) const
 
3483
  {
 
3484
    switch ( d )
 
3485
    {
 
3486
    case 0:
 
3487
      return true;
 
3488
      break;
 
3489
    case 1:
 
3490
      return false;
 
3491
      break;
 
3492
    case 2:
 
3493
      return false;
 
3494
      break;
 
3495
    }
 
3496
    return false;
 
3497
  }
 
3498
 
 
3499
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
3500
  virtual bool init_mesh(const ufc::mesh& m)
 
3501
  {
 
3502
    __global_dimension = m.num_entities[0];
 
3503
    return false;
 
3504
  }
 
3505
 
 
3506
  /// Initialize dof map for given cell
 
3507
  virtual void init_cell(const ufc::mesh& m,
 
3508
                         const ufc::cell& c)
 
3509
  {
 
3510
    // Do nothing
 
3511
  }
 
3512
 
 
3513
  /// Finish initialization of dof map for cells
 
3514
  virtual void init_cell_finalize()
 
3515
  {
 
3516
    // Do nothing
 
3517
  }
 
3518
 
 
3519
  /// Return the dimension of the global finite element function space
 
3520
  virtual unsigned int global_dimension() const
 
3521
  {
 
3522
    return __global_dimension;
 
3523
  }
 
3524
 
 
3525
  /// Return the dimension of the local finite element function space
 
3526
  virtual unsigned int local_dimension() const
 
3527
  {
 
3528
    return 3;
 
3529
  }
 
3530
 
 
3531
  // Return the geometric dimension of the coordinates this dof map provides
 
3532
  virtual unsigned int geometric_dimension() const
 
3533
  {
 
3534
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3535
  }
 
3536
 
 
3537
  /// Return the number of dofs on each cell facet
 
3538
  virtual unsigned int num_facet_dofs() const
 
3539
  {
 
3540
    return 2;
 
3541
  }
 
3542
 
 
3543
  /// Return the number of dofs associated with each cell entity of dimension d
 
3544
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
3545
  {
 
3546
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3547
  }
 
3548
 
 
3549
  /// Tabulate the local-to-global mapping of dofs on a cell
 
3550
  virtual void tabulate_dofs(unsigned int* dofs,
 
3551
                             const ufc::mesh& m,
 
3552
                             const ufc::cell& c) const
 
3553
  {
 
3554
    dofs[0] = c.entity_indices[0][0];
 
3555
    dofs[1] = c.entity_indices[0][1];
 
3556
    dofs[2] = c.entity_indices[0][2];
 
3557
  }
 
3558
 
 
3559
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
3560
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
3561
                                   unsigned int facet) const
 
3562
  {
 
3563
    switch ( facet )
 
3564
    {
 
3565
    case 0:
 
3566
      dofs[0] = 1;
 
3567
      dofs[1] = 2;
 
3568
      break;
 
3569
    case 1:
 
3570
      dofs[0] = 0;
 
3571
      dofs[1] = 2;
 
3572
      break;
 
3573
    case 2:
 
3574
      dofs[0] = 0;
 
3575
      dofs[1] = 1;
 
3576
      break;
 
3577
    }
 
3578
  }
 
3579
 
 
3580
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
3581
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
3582
                                    unsigned int d, unsigned int i) const
 
3583
  {
 
3584
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3585
  }
 
3586
 
 
3587
  /// Tabulate the coordinates of all dofs on a cell
 
3588
  virtual void tabulate_coordinates(double** coordinates,
 
3589
                                    const ufc::cell& c) const
 
3590
  {
 
3591
    const double * const * x = c.coordinates;
 
3592
    coordinates[0][0] = x[0][0];
 
3593
    coordinates[0][1] = x[0][1];
 
3594
    coordinates[1][0] = x[1][0];
 
3595
    coordinates[1][1] = x[1][1];
 
3596
    coordinates[2][0] = x[2][0];
 
3597
    coordinates[2][1] = x[2][1];
 
3598
  }
 
3599
 
 
3600
  /// Return the number of sub dof maps (for a mixed element)
 
3601
  virtual unsigned int num_sub_dof_maps() const
 
3602
  {
 
3603
    return 1;
 
3604
  }
 
3605
 
 
3606
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
3607
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
3608
  {
 
3609
    return new UFC_NonlinearPoissonLinearForm_dof_map_1();
 
3610
  }
 
3611
 
 
3612
};
 
3613
 
 
3614
/// This class defines the interface for a local-to-global mapping of
 
3615
/// degrees of freedom (dofs).
 
3616
 
 
3617
class UFC_NonlinearPoissonLinearForm_dof_map_2: public ufc::dof_map
 
3618
{
 
3619
private:
 
3620
 
 
3621
  unsigned int __global_dimension;
 
3622
 
 
3623
public:
 
3624
 
 
3625
  /// Constructor
 
3626
  UFC_NonlinearPoissonLinearForm_dof_map_2() : ufc::dof_map()
 
3627
  {
 
3628
    __global_dimension = 0;
 
3629
  }
 
3630
 
 
3631
  /// Destructor
 
3632
  virtual ~UFC_NonlinearPoissonLinearForm_dof_map_2()
 
3633
  {
 
3634
    // Do nothing
 
3635
  }
 
3636
 
 
3637
  /// Return a string identifying the dof map
 
3638
  virtual const char* signature() const
 
3639
  {
 
3640
    return "FFC dof map for Lagrange finite element of degree 1 on a triangle";
 
3641
  }
 
3642
 
 
3643
  /// Return true iff mesh entities of topological dimension d are needed
 
3644
  virtual bool needs_mesh_entities(unsigned int d) const
 
3645
  {
 
3646
    switch ( d )
 
3647
    {
 
3648
    case 0:
 
3649
      return true;
 
3650
      break;
 
3651
    case 1:
 
3652
      return false;
 
3653
      break;
 
3654
    case 2:
 
3655
      return false;
 
3656
      break;
 
3657
    }
 
3658
    return false;
 
3659
  }
 
3660
 
 
3661
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
3662
  virtual bool init_mesh(const ufc::mesh& m)
 
3663
  {
 
3664
    __global_dimension = m.num_entities[0];
 
3665
    return false;
 
3666
  }
 
3667
 
 
3668
  /// Initialize dof map for given cell
 
3669
  virtual void init_cell(const ufc::mesh& m,
 
3670
                         const ufc::cell& c)
 
3671
  {
 
3672
    // Do nothing
 
3673
  }
 
3674
 
 
3675
  /// Finish initialization of dof map for cells
 
3676
  virtual void init_cell_finalize()
 
3677
  {
 
3678
    // Do nothing
 
3679
  }
 
3680
 
 
3681
  /// Return the dimension of the global finite element function space
 
3682
  virtual unsigned int global_dimension() const
 
3683
  {
 
3684
    return __global_dimension;
 
3685
  }
 
3686
 
 
3687
  /// Return the dimension of the local finite element function space
 
3688
  virtual unsigned int local_dimension() const
 
3689
  {
 
3690
    return 3;
 
3691
  }
 
3692
 
 
3693
  // Return the geometric dimension of the coordinates this dof map provides
 
3694
  virtual unsigned int geometric_dimension() const
 
3695
  {
 
3696
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3697
  }
 
3698
 
 
3699
  /// Return the number of dofs on each cell facet
 
3700
  virtual unsigned int num_facet_dofs() const
 
3701
  {
 
3702
    return 2;
 
3703
  }
 
3704
 
 
3705
  /// Return the number of dofs associated with each cell entity of dimension d
 
3706
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
3707
  {
 
3708
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3709
  }
 
3710
 
 
3711
  /// Tabulate the local-to-global mapping of dofs on a cell
 
3712
  virtual void tabulate_dofs(unsigned int* dofs,
 
3713
                             const ufc::mesh& m,
 
3714
                             const ufc::cell& c) const
 
3715
  {
 
3716
    dofs[0] = c.entity_indices[0][0];
 
3717
    dofs[1] = c.entity_indices[0][1];
 
3718
    dofs[2] = c.entity_indices[0][2];
 
3719
  }
 
3720
 
 
3721
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
3722
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
3723
                                   unsigned int facet) const
 
3724
  {
 
3725
    switch ( facet )
 
3726
    {
 
3727
    case 0:
 
3728
      dofs[0] = 1;
 
3729
      dofs[1] = 2;
 
3730
      break;
 
3731
    case 1:
 
3732
      dofs[0] = 0;
 
3733
      dofs[1] = 2;
 
3734
      break;
 
3735
    case 2:
 
3736
      dofs[0] = 0;
 
3737
      dofs[1] = 1;
 
3738
      break;
 
3739
    }
 
3740
  }
 
3741
 
 
3742
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
3743
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
3744
                                    unsigned int d, unsigned int i) const
 
3745
  {
 
3746
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3747
  }
 
3748
 
 
3749
  /// Tabulate the coordinates of all dofs on a cell
 
3750
  virtual void tabulate_coordinates(double** coordinates,
 
3751
                                    const ufc::cell& c) const
 
3752
  {
 
3753
    const double * const * x = c.coordinates;
 
3754
    coordinates[0][0] = x[0][0];
 
3755
    coordinates[0][1] = x[0][1];
 
3756
    coordinates[1][0] = x[1][0];
 
3757
    coordinates[1][1] = x[1][1];
 
3758
    coordinates[2][0] = x[2][0];
 
3759
    coordinates[2][1] = x[2][1];
 
3760
  }
 
3761
 
 
3762
  /// Return the number of sub dof maps (for a mixed element)
 
3763
  virtual unsigned int num_sub_dof_maps() const
 
3764
  {
 
3765
    return 1;
 
3766
  }
 
3767
 
 
3768
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
3769
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
3770
  {
 
3771
    return new UFC_NonlinearPoissonLinearForm_dof_map_2();
 
3772
  }
 
3773
 
 
3774
};
 
3775
 
 
3776
/// This class defines the interface for the tabulation of the cell
 
3777
/// tensor corresponding to the local contribution to a form from
 
3778
/// the integral over a cell.
 
3779
 
 
3780
class UFC_NonlinearPoissonLinearForm_cell_integral_0: public ufc::cell_integral
 
3781
{
 
3782
public:
 
3783
 
 
3784
  /// Constructor
 
3785
  UFC_NonlinearPoissonLinearForm_cell_integral_0() : ufc::cell_integral()
 
3786
  {
 
3787
    // Do nothing
 
3788
  }
 
3789
 
 
3790
  /// Destructor
 
3791
  virtual ~UFC_NonlinearPoissonLinearForm_cell_integral_0()
 
3792
  {
 
3793
    // Do nothing
 
3794
  }
 
3795
 
 
3796
  /// Tabulate the tensor for the contribution from a local cell
 
3797
  virtual void tabulate_tensor(double* A,
 
3798
                               const double * const * w,
 
3799
                               const ufc::cell& c) const
 
3800
  {
 
3801
    // Extract vertex coordinates
 
3802
    const double * const * x = c.coordinates;
 
3803
    
 
3804
    // Compute Jacobian of affine map from reference cell
 
3805
    const double J_00 = x[1][0] - x[0][0];
 
3806
    const double J_01 = x[2][0] - x[0][0];
 
3807
    const double J_10 = x[1][1] - x[0][1];
 
3808
    const double J_11 = x[2][1] - x[0][1];
 
3809
      
 
3810
    // Compute determinant of Jacobian
 
3811
    double detJ = J_00*J_11 - J_01*J_10;
 
3812
      
 
3813
    // Compute inverse of Jacobian
 
3814
    const double Jinv_00 =  J_11 / detJ;
 
3815
    const double Jinv_01 = -J_01 / detJ;
 
3816
    const double Jinv_10 = -J_10 / detJ;
 
3817
    const double Jinv_11 =  J_00 / detJ;
 
3818
    
 
3819
    // Set scale factor
 
3820
    const double det = std::abs(detJ);
 
3821
    
 
3822
    // Compute coefficients
 
3823
    const double c1_0_0_0 = w[1][0];
 
3824
    const double c1_0_0_1 = w[1][1];
 
3825
    const double c1_0_0_2 = w[1][2];
 
3826
    const double c0_1_0_0 = w[0][0];
 
3827
    const double c0_1_0_1 = w[0][1];
 
3828
    const double c0_1_0_2 = w[0][2];
 
3829
    const double c0_1_1_0 = w[0][0];
 
3830
    const double c0_1_1_1 = w[0][1];
 
3831
    const double c0_1_1_2 = w[0][2];
 
3832
    const double c0_1_2_0 = w[0][0];
 
3833
    const double c0_1_2_1 = w[0][1];
 
3834
    const double c0_1_2_2 = w[0][2];
 
3835
    const double c0_2_0_0 = w[0][0];
 
3836
    const double c0_2_0_1 = w[0][1];
 
3837
    const double c0_2_0_2 = w[0][2];
 
3838
    
 
3839
    // Compute geometry tensors
 
3840
    const double G0_0 = det*c1_0_0_0;
 
3841
    const double G0_1 = det*c1_0_0_1;
 
3842
    const double G0_2 = det*c1_0_0_2;
 
3843
    const double G1_0_0_0_0_0 = det*c0_1_0_0*c0_1_1_0*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3844
    const double G1_0_0_0_0_1 = det*c0_1_0_0*c0_1_1_0*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3845
    const double G1_0_0_0_1_0 = det*c0_1_0_0*c0_1_1_0*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3846
    const double G1_0_0_0_2_1 = det*c0_1_0_0*c0_1_1_0*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3847
    const double G1_0_0_1_0_0 = det*c0_1_0_0*c0_1_1_0*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3848
    const double G1_0_0_1_0_1 = det*c0_1_0_0*c0_1_1_0*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3849
    const double G1_0_0_1_1_0 = det*c0_1_0_0*c0_1_1_0*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3850
    const double G1_0_0_1_2_1 = det*c0_1_0_0*c0_1_1_0*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3851
    const double G1_0_1_0_0_0 = det*c0_1_0_0*c0_1_1_1*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3852
    const double G1_0_1_0_0_1 = det*c0_1_0_0*c0_1_1_1*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3853
    const double G1_0_1_0_1_0 = det*c0_1_0_0*c0_1_1_1*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3854
    const double G1_0_1_0_2_1 = det*c0_1_0_0*c0_1_1_1*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3855
    const double G1_0_1_1_0_0 = det*c0_1_0_0*c0_1_1_1*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3856
    const double G1_0_1_1_0_1 = det*c0_1_0_0*c0_1_1_1*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3857
    const double G1_0_1_1_1_0 = det*c0_1_0_0*c0_1_1_1*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3858
    const double G1_0_1_1_2_1 = det*c0_1_0_0*c0_1_1_1*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3859
    const double G1_0_2_0_0_0 = det*c0_1_0_0*c0_1_1_2*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3860
    const double G1_0_2_0_0_1 = det*c0_1_0_0*c0_1_1_2*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3861
    const double G1_0_2_0_1_0 = det*c0_1_0_0*c0_1_1_2*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3862
    const double G1_0_2_0_2_1 = det*c0_1_0_0*c0_1_1_2*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3863
    const double G1_0_2_1_0_0 = det*c0_1_0_0*c0_1_1_2*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3864
    const double G1_0_2_1_0_1 = det*c0_1_0_0*c0_1_1_2*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3865
    const double G1_0_2_1_1_0 = det*c0_1_0_0*c0_1_1_2*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3866
    const double G1_0_2_1_2_1 = det*c0_1_0_0*c0_1_1_2*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3867
    const double G1_1_0_0_0_0 = det*c0_1_0_1*c0_1_1_0*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3868
    const double G1_1_0_0_0_1 = det*c0_1_0_1*c0_1_1_0*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3869
    const double G1_1_0_0_1_0 = det*c0_1_0_1*c0_1_1_0*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3870
    const double G1_1_0_0_2_1 = det*c0_1_0_1*c0_1_1_0*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3871
    const double G1_1_0_1_0_0 = det*c0_1_0_1*c0_1_1_0*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3872
    const double G1_1_0_1_0_1 = det*c0_1_0_1*c0_1_1_0*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3873
    const double G1_1_0_1_1_0 = det*c0_1_0_1*c0_1_1_0*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3874
    const double G1_1_0_1_2_1 = det*c0_1_0_1*c0_1_1_0*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3875
    const double G1_1_1_0_0_0 = det*c0_1_0_1*c0_1_1_1*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3876
    const double G1_1_1_0_0_1 = det*c0_1_0_1*c0_1_1_1*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3877
    const double G1_1_1_0_1_0 = det*c0_1_0_1*c0_1_1_1*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3878
    const double G1_1_1_0_2_1 = det*c0_1_0_1*c0_1_1_1*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3879
    const double G1_1_1_1_0_0 = det*c0_1_0_1*c0_1_1_1*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3880
    const double G1_1_1_1_0_1 = det*c0_1_0_1*c0_1_1_1*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3881
    const double G1_1_1_1_1_0 = det*c0_1_0_1*c0_1_1_1*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3882
    const double G1_1_1_1_2_1 = det*c0_1_0_1*c0_1_1_1*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3883
    const double G1_1_2_0_0_0 = det*c0_1_0_1*c0_1_1_2*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3884
    const double G1_1_2_0_0_1 = det*c0_1_0_1*c0_1_1_2*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3885
    const double G1_1_2_0_1_0 = det*c0_1_0_1*c0_1_1_2*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3886
    const double G1_1_2_0_2_1 = det*c0_1_0_1*c0_1_1_2*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3887
    const double G1_1_2_1_0_0 = det*c0_1_0_1*c0_1_1_2*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3888
    const double G1_1_2_1_0_1 = det*c0_1_0_1*c0_1_1_2*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3889
    const double G1_1_2_1_1_0 = det*c0_1_0_1*c0_1_1_2*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3890
    const double G1_1_2_1_2_1 = det*c0_1_0_1*c0_1_1_2*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3891
    const double G1_2_0_0_0_0 = det*c0_1_0_2*c0_1_1_0*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3892
    const double G1_2_0_0_0_1 = det*c0_1_0_2*c0_1_1_0*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3893
    const double G1_2_0_0_1_0 = det*c0_1_0_2*c0_1_1_0*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3894
    const double G1_2_0_0_2_1 = det*c0_1_0_2*c0_1_1_0*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3895
    const double G1_2_0_1_0_0 = det*c0_1_0_2*c0_1_1_0*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3896
    const double G1_2_0_1_0_1 = det*c0_1_0_2*c0_1_1_0*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3897
    const double G1_2_0_1_1_0 = det*c0_1_0_2*c0_1_1_0*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3898
    const double G1_2_0_1_2_1 = det*c0_1_0_2*c0_1_1_0*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3899
    const double G1_2_1_0_0_0 = det*c0_1_0_2*c0_1_1_1*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3900
    const double G1_2_1_0_0_1 = det*c0_1_0_2*c0_1_1_1*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3901
    const double G1_2_1_0_1_0 = det*c0_1_0_2*c0_1_1_1*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3902
    const double G1_2_1_0_2_1 = det*c0_1_0_2*c0_1_1_1*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3903
    const double G1_2_1_1_0_0 = det*c0_1_0_2*c0_1_1_1*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3904
    const double G1_2_1_1_0_1 = det*c0_1_0_2*c0_1_1_1*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3905
    const double G1_2_1_1_1_0 = det*c0_1_0_2*c0_1_1_1*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3906
    const double G1_2_1_1_2_1 = det*c0_1_0_2*c0_1_1_1*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3907
    const double G1_2_2_0_0_0 = det*c0_1_0_2*c0_1_1_2*c0_1_2_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3908
    const double G1_2_2_0_0_1 = det*c0_1_0_2*c0_1_1_2*c0_1_2_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3909
    const double G1_2_2_0_1_0 = det*c0_1_0_2*c0_1_1_2*c0_1_2_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3910
    const double G1_2_2_0_2_1 = det*c0_1_0_2*c0_1_1_2*c0_1_2_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3911
    const double G1_2_2_1_0_0 = det*c0_1_0_2*c0_1_1_2*c0_1_2_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3912
    const double G1_2_2_1_0_1 = det*c0_1_0_2*c0_1_1_2*c0_1_2_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3913
    const double G1_2_2_1_1_0 = det*c0_1_0_2*c0_1_1_2*c0_1_2_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3914
    const double G1_2_2_1_2_1 = det*c0_1_0_2*c0_1_1_2*c0_1_2_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3915
    const double G2_0_0_0 = det*c0_2_0_0*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3916
    const double G2_0_0_1 = det*c0_2_0_0*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3917
    const double G2_0_1_0 = det*c0_2_0_1*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);
 
3918
    const double G2_0_2_1 = det*c0_2_0_2*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);
 
3919
    const double G2_1_0_0 = det*c0_2_0_0*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3920
    const double G2_1_0_1 = det*c0_2_0_0*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3921
    const double G2_1_1_0 = det*c0_2_0_1*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);
 
3922
    const double G2_1_2_1 = det*c0_2_0_2*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);
 
3923
    
 
3924
    // Compute element tensor
 
3925
    A[0] = 0.0833333333333332*G0_0 + 0.0416666666666666*G0_1 + 0.0416666666666666*G0_2 - 0.0833333333333332*G1_0_0_0_0_0 - 0.0833333333333332*G1_0_0_0_0_1 + 0.0833333333333332*G1_0_0_0_1_0 + 0.0833333333333332*G1_0_0_0_2_1 - 0.0833333333333332*G1_0_0_1_0_0 - 0.0833333333333332*G1_0_0_1_0_1 + 0.0833333333333332*G1_0_0_1_1_0 + 0.0833333333333332*G1_0_0_1_2_1 - 0.0416666666666666*G1_0_1_0_0_0 - 0.0416666666666666*G1_0_1_0_0_1 + 0.0416666666666666*G1_0_1_0_1_0 + 0.0416666666666666*G1_0_1_0_2_1 - 0.0416666666666666*G1_0_1_1_0_0 - 0.0416666666666666*G1_0_1_1_0_1 + 0.0416666666666666*G1_0_1_1_1_0 + 0.0416666666666666*G1_0_1_1_2_1 - 0.0416666666666666*G1_0_2_0_0_0 - 0.0416666666666666*G1_0_2_0_0_1 + 0.0416666666666666*G1_0_2_0_1_0 + 0.0416666666666666*G1_0_2_0_2_1 - 0.0416666666666666*G1_0_2_1_0_0 - 0.0416666666666666*G1_0_2_1_0_1 + 0.0416666666666666*G1_0_2_1_1_0 + 0.0416666666666666*G1_0_2_1_2_1 - 0.0416666666666666*G1_1_0_0_0_0 - 0.0416666666666666*G1_1_0_0_0_1 + 0.0416666666666666*G1_1_0_0_1_0 + 0.0416666666666666*G1_1_0_0_2_1 - 0.0416666666666666*G1_1_0_1_0_0 - 0.0416666666666666*G1_1_0_1_0_1 + 0.0416666666666666*G1_1_0_1_1_0 + 0.0416666666666666*G1_1_0_1_2_1 - 0.0833333333333332*G1_1_1_0_0_0 - 0.0833333333333332*G1_1_1_0_0_1 + 0.0833333333333332*G1_1_1_0_1_0 + 0.0833333333333332*G1_1_1_0_2_1 - 0.0833333333333332*G1_1_1_1_0_0 - 0.0833333333333332*G1_1_1_1_0_1 + 0.0833333333333332*G1_1_1_1_1_0 + 0.0833333333333332*G1_1_1_1_2_1 - 0.0416666666666666*G1_1_2_0_0_0 - 0.0416666666666666*G1_1_2_0_0_1 + 0.0416666666666666*G1_1_2_0_1_0 + 0.0416666666666666*G1_1_2_0_2_1 - 0.0416666666666666*G1_1_2_1_0_0 - 0.0416666666666666*G1_1_2_1_0_1 + 0.0416666666666666*G1_1_2_1_1_0 + 0.0416666666666666*G1_1_2_1_2_1 - 0.0416666666666666*G1_2_0_0_0_0 - 0.0416666666666666*G1_2_0_0_0_1 + 0.0416666666666666*G1_2_0_0_1_0 + 0.0416666666666666*G1_2_0_0_2_1 - 0.0416666666666666*G1_2_0_1_0_0 - 0.0416666666666666*G1_2_0_1_0_1 + 0.0416666666666666*G1_2_0_1_1_0 + 0.0416666666666666*G1_2_0_1_2_1 - 0.0416666666666666*G1_2_1_0_0_0 - 0.0416666666666666*G1_2_1_0_0_1 + 0.0416666666666666*G1_2_1_0_1_0 + 0.0416666666666666*G1_2_1_0_2_1 - 0.0416666666666666*G1_2_1_1_0_0 - 0.0416666666666666*G1_2_1_1_0_1 + 0.0416666666666666*G1_2_1_1_1_0 + 0.0416666666666666*G1_2_1_1_2_1 - 0.0833333333333332*G1_2_2_0_0_0 - 0.0833333333333332*G1_2_2_0_0_1 + 0.0833333333333332*G1_2_2_0_1_0 + 0.0833333333333332*G1_2_2_0_2_1 - 0.0833333333333332*G1_2_2_1_0_0 - 0.0833333333333332*G1_2_2_1_0_1 + 0.0833333333333332*G1_2_2_1_1_0 + 0.0833333333333332*G1_2_2_1_2_1 - 0.5*G2_0_0_0 - 0.5*G2_0_0_1 + 0.5*G2_0_1_0 + 0.5*G2_0_2_1 - 0.5*G2_1_0_0 - 0.5*G2_1_0_1 + 0.5*G2_1_1_0 + 0.5*G2_1_2_1;
 
3926
    A[1] = 0.0416666666666666*G0_0 + 0.0833333333333332*G0_1 + 0.0416666666666666*G0_2 + 0.0833333333333332*G1_0_0_0_0_0 + 0.0833333333333332*G1_0_0_0_0_1 - 0.0833333333333332*G1_0_0_0_1_0 - 0.0833333333333332*G1_0_0_0_2_1 + 0.0416666666666666*G1_0_1_0_0_0 + 0.0416666666666666*G1_0_1_0_0_1 - 0.0416666666666666*G1_0_1_0_1_0 - 0.0416666666666666*G1_0_1_0_2_1 + 0.0416666666666666*G1_0_2_0_0_0 + 0.0416666666666666*G1_0_2_0_0_1 - 0.0416666666666666*G1_0_2_0_1_0 - 0.0416666666666666*G1_0_2_0_2_1 + 0.0416666666666666*G1_1_0_0_0_0 + 0.0416666666666666*G1_1_0_0_0_1 - 0.0416666666666666*G1_1_0_0_1_0 - 0.0416666666666666*G1_1_0_0_2_1 + 0.0833333333333332*G1_1_1_0_0_0 + 0.0833333333333332*G1_1_1_0_0_1 - 0.0833333333333332*G1_1_1_0_1_0 - 0.0833333333333332*G1_1_1_0_2_1 + 0.0416666666666666*G1_1_2_0_0_0 + 0.0416666666666666*G1_1_2_0_0_1 - 0.0416666666666666*G1_1_2_0_1_0 - 0.0416666666666666*G1_1_2_0_2_1 + 0.0416666666666666*G1_2_0_0_0_0 + 0.0416666666666666*G1_2_0_0_0_1 - 0.0416666666666666*G1_2_0_0_1_0 - 0.0416666666666666*G1_2_0_0_2_1 + 0.0416666666666666*G1_2_1_0_0_0 + 0.0416666666666666*G1_2_1_0_0_1 - 0.0416666666666666*G1_2_1_0_1_0 - 0.0416666666666666*G1_2_1_0_2_1 + 0.0833333333333332*G1_2_2_0_0_0 + 0.0833333333333332*G1_2_2_0_0_1 - 0.0833333333333332*G1_2_2_0_1_0 - 0.0833333333333332*G1_2_2_0_2_1 + 0.5*G2_0_0_0 + 0.5*G2_0_0_1 - 0.5*G2_0_1_0 - 0.5*G2_0_2_1;
 
3927
    A[2] = 0.0416666666666666*G0_0 + 0.0416666666666666*G0_1 + 0.0833333333333332*G0_2 + 0.0833333333333332*G1_0_0_1_0_0 + 0.0833333333333332*G1_0_0_1_0_1 - 0.0833333333333332*G1_0_0_1_1_0 - 0.0833333333333332*G1_0_0_1_2_1 + 0.0416666666666666*G1_0_1_1_0_0 + 0.0416666666666666*G1_0_1_1_0_1 - 0.0416666666666666*G1_0_1_1_1_0 - 0.0416666666666666*G1_0_1_1_2_1 + 0.0416666666666666*G1_0_2_1_0_0 + 0.0416666666666666*G1_0_2_1_0_1 - 0.0416666666666666*G1_0_2_1_1_0 - 0.0416666666666666*G1_0_2_1_2_1 + 0.0416666666666666*G1_1_0_1_0_0 + 0.0416666666666666*G1_1_0_1_0_1 - 0.0416666666666666*G1_1_0_1_1_0 - 0.0416666666666666*G1_1_0_1_2_1 + 0.0833333333333332*G1_1_1_1_0_0 + 0.0833333333333332*G1_1_1_1_0_1 - 0.0833333333333332*G1_1_1_1_1_0 - 0.0833333333333332*G1_1_1_1_2_1 + 0.0416666666666666*G1_1_2_1_0_0 + 0.0416666666666666*G1_1_2_1_0_1 - 0.0416666666666666*G1_1_2_1_1_0 - 0.0416666666666666*G1_1_2_1_2_1 + 0.0416666666666666*G1_2_0_1_0_0 + 0.0416666666666666*G1_2_0_1_0_1 - 0.0416666666666666*G1_2_0_1_1_0 - 0.0416666666666666*G1_2_0_1_2_1 + 0.0416666666666666*G1_2_1_1_0_0 + 0.0416666666666666*G1_2_1_1_0_1 - 0.0416666666666666*G1_2_1_1_1_0 - 0.0416666666666666*G1_2_1_1_2_1 + 0.0833333333333332*G1_2_2_1_0_0 + 0.0833333333333332*G1_2_2_1_0_1 - 0.0833333333333332*G1_2_2_1_1_0 - 0.0833333333333332*G1_2_2_1_2_1 + 0.5*G2_1_0_0 + 0.5*G2_1_0_1 - 0.5*G2_1_1_0 - 0.5*G2_1_2_1;
 
3928
  }
 
3929
 
 
3930
};
 
3931
 
 
3932
/// This class defines the interface for the assembly of the global
 
3933
/// tensor corresponding to a form with r + n arguments, that is, a
 
3934
/// mapping
 
3935
///
 
3936
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
3937
///
 
3938
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
3939
/// global tensor A is defined by
 
3940
///
 
3941
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
3942
///
 
3943
/// where each argument Vj represents the application to the
 
3944
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
3945
/// fixed functions (coefficients).
 
3946
 
 
3947
class UFC_NonlinearPoissonLinearForm: public ufc::form
 
3948
{
 
3949
public:
 
3950
 
 
3951
  /// Constructor
 
3952
  UFC_NonlinearPoissonLinearForm() : ufc::form()
 
3953
  {
 
3954
    // Do nothing
 
3955
  }
 
3956
 
 
3957
  /// Destructor
 
3958
  virtual ~UFC_NonlinearPoissonLinearForm()
 
3959
  {
 
3960
    // Do nothing
 
3961
  }
 
3962
 
 
3963
  /// Return a string identifying the form
 
3964
  virtual const char* signature() const
 
3965
  {
 
3966
    return "w1_a0 | vi0*va0*dX(0) + -w0_a0w0_a1w0_a3(dXa2/dxb0)(dXa4/dxb0) | ((d/dXa2)vi0)*va0*va1*((d/dXa4)va3)*dX(0) + -w0_a1(dXa0/dxb0)(dXa2/dxb0) | ((d/dXa0)vi0)*((d/dXa2)va1)*dX(0)";
 
3967
  }
 
3968
 
 
3969
  /// Return the rank of the global tensor (r)
 
3970
  virtual unsigned int rank() const
 
3971
  {
 
3972
    return 1;
 
3973
  }
 
3974
 
 
3975
  /// Return the number of coefficients (n)
 
3976
  virtual unsigned int num_coefficients() const
 
3977
  {
 
3978
    return 2;
 
3979
  }
 
3980
 
 
3981
  /// Return the number of cell integrals
 
3982
  virtual unsigned int num_cell_integrals() const
 
3983
  {
 
3984
    return 1;
 
3985
  }
 
3986
  
 
3987
  /// Return the number of exterior facet integrals
 
3988
  virtual unsigned int num_exterior_facet_integrals() const
 
3989
  {
 
3990
    return 0;
 
3991
  }
 
3992
  
 
3993
  /// Return the number of interior facet integrals
 
3994
  virtual unsigned int num_interior_facet_integrals() const
 
3995
  {
 
3996
    return 0;
 
3997
  }
 
3998
    
 
3999
  /// Create a new finite element for argument function i
 
4000
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
4001
  {
 
4002
    switch ( i )
 
4003
    {
 
4004
    case 0:
 
4005
      return new UFC_NonlinearPoissonLinearForm_finite_element_0();
 
4006
      break;
 
4007
    case 1:
 
4008
      return new UFC_NonlinearPoissonLinearForm_finite_element_1();
 
4009
      break;
 
4010
    case 2:
 
4011
      return new UFC_NonlinearPoissonLinearForm_finite_element_2();
 
4012
      break;
 
4013
    }
 
4014
    return 0;
 
4015
  }
 
4016
  
 
4017
  /// Create a new dof map for argument function i
 
4018
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
4019
  {
 
4020
    switch ( i )
 
4021
    {
 
4022
    case 0:
 
4023
      return new UFC_NonlinearPoissonLinearForm_dof_map_0();
 
4024
      break;
 
4025
    case 1:
 
4026
      return new UFC_NonlinearPoissonLinearForm_dof_map_1();
 
4027
      break;
 
4028
    case 2:
 
4029
      return new UFC_NonlinearPoissonLinearForm_dof_map_2();
 
4030
      break;
 
4031
    }
 
4032
    return 0;
 
4033
  }
 
4034
 
 
4035
  /// Create a new cell integral on sub domain i
 
4036
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
4037
  {
 
4038
    return new UFC_NonlinearPoissonLinearForm_cell_integral_0();
 
4039
  }
 
4040
 
 
4041
  /// Create a new exterior facet integral on sub domain i
 
4042
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
4043
  {
 
4044
    return 0;
 
4045
  }
 
4046
 
 
4047
  /// Create a new interior facet integral on sub domain i
 
4048
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
4049
  {
 
4050
    return 0;
 
4051
  }
 
4052
 
 
4053
};
 
4054
 
 
4055
// DOLFIN wrappers
 
4056
 
 
4057
#include <dolfin/fem/Form.h>
 
4058
 
 
4059
class NonlinearPoissonBilinearForm : public dolfin::Form
 
4060
{
 
4061
public:
 
4062
 
 
4063
  NonlinearPoissonBilinearForm(dolfin::Function& w0) : dolfin::Form()
 
4064
  {
 
4065
    __coefficients.push_back(&w0);
 
4066
  }
 
4067
 
 
4068
  /// Return UFC form
 
4069
  virtual const ufc::form& form() const
 
4070
  {
 
4071
    return __form;
 
4072
  }
 
4073
  
 
4074
  /// Return array of coefficients
 
4075
  virtual const dolfin::Array<dolfin::Function*>& coefficients() const
 
4076
  {
 
4077
    return __coefficients;
 
4078
  }
 
4079
 
 
4080
private:
 
4081
 
 
4082
  // UFC form
 
4083
  UFC_NonlinearPoissonBilinearForm __form;
 
4084
 
 
4085
  /// Array of coefficients
 
4086
  dolfin::Array<dolfin::Function*> __coefficients;
 
4087
 
 
4088
};
 
4089
 
 
4090
class NonlinearPoissonLinearForm : public dolfin::Form
 
4091
{
 
4092
public:
 
4093
 
 
4094
  NonlinearPoissonLinearForm(dolfin::Function& w0, dolfin::Function& w1) : dolfin::Form()
 
4095
  {
 
4096
    __coefficients.push_back(&w0);
 
4097
    __coefficients.push_back(&w1);
 
4098
  }
 
4099
 
 
4100
  /// Return UFC form
 
4101
  virtual const ufc::form& form() const
 
4102
  {
 
4103
    return __form;
 
4104
  }
 
4105
  
 
4106
  /// Return array of coefficients
 
4107
  virtual const dolfin::Array<dolfin::Function*>& coefficients() const
 
4108
  {
 
4109
    return __coefficients;
 
4110
  }
 
4111
 
 
4112
private:
 
4113
 
 
4114
  // UFC form
 
4115
  UFC_NonlinearPoissonLinearForm __form;
 
4116
 
 
4117
  /// Array of coefficients
 
4118
  dolfin::Array<dolfin::Function*> __coefficients;
 
4119
 
 
4120
};
 
4121
 
 
4122
#endif