~ubuntu-branches/ubuntu/vivid/ffc/vivid

« back to all changes in this revision

Viewing changes to test/regression/references/r_auto/SubDomain.h

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2014-01-10 13:56:45 UTC
  • mfrom: (1.1.14)
  • Revision ID: package-import@ubuntu.com-20140110135645-4ozcd71y1oggj44z
Tags: 1.3.0-1
* New upstream release.
* debian/watch: Update URL for move to Bitbucket.
* debian/docs: README -> README.rst and remove TODO.
* debian/control:
  - Add python-numpy to Build-Depends.
  - Replace python-all with python-all-dev in Build-Depends.
  - Add ${shlibs:Depends} to Depends.
  - Change to Architecture: any.
  - Bump Standards-Version to 3.9.5 (no changes needed).
* debian/rules: Call dh_numpy in override_dh_python2.

Show diffs side-by-side

added added

removed removed

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