~corrado-maurini/dolfin/tao

« back to all changes in this revision

Viewing changes to demo/undocumented/restriction/cpp/Poisson.h

  • Committer: corrado maurini
  • Date: 2012-12-18 12:16:08 UTC
  • mfrom: (6685.78.207 trunk)
  • Revision ID: corrado.maurini@upmc.fr-20121218121608-nk82ly9jgsld9u84
updating with trunk, fix uint in TAO solver and hacking the check for tao FindTAO.cmake

Show diffs side-by-side

added added

removed removed

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