← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to bench/fem/passembly/Poisson3D_1.h

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
bench/fem/passembly/Poisson3D_1.h
 
1
// This code conforms with the UFC specification version 1.0
 
2
// and was automatically generated by FFC version 0.4.3.
 
3
//
 
4
// Warning: This code was generated with the option '-l dolfin'
 
5
// and contains DOLFIN-specific wrappers that depend on DOLFIN.
 
6
 
 
7
#ifndef __POISSON3D_1_H
 
8
#define __POISSON3D_1_H
 
9
 
 
10
#include <cmath>
 
11
#include <stdexcept>
 
12
#include <fstream>
 
13
#include <ufc.h>
 
14
 
 
15
/// This class defines the interface for a finite element.
 
16
 
 
17
class UFC_Poisson3D_1BilinearForm_finite_element_0: public ufc::finite_element
 
18
{
 
19
public:
 
20
 
 
21
  /// Constructor
 
22
  UFC_Poisson3D_1BilinearForm_finite_element_0() : ufc::finite_element()
 
23
  {
 
24
    // Do nothing
 
25
  }
 
26
 
 
27
  /// Destructor
 
28
  virtual ~UFC_Poisson3D_1BilinearForm_finite_element_0()
 
29
  {
 
30
    // Do nothing
 
31
  }
 
32
 
 
33
  /// Return a string identifying the finite element
 
34
  virtual const char* signature() const
 
35
  {
 
36
    return "Lagrange finite element of degree 1 on a tetrahedron";
 
37
  }
 
38
 
 
39
  /// Return the cell shape
 
40
  virtual ufc::shape cell_shape() const
 
41
  {
 
42
    return ufc::tetrahedron;
 
43
  }
 
44
 
 
45
  /// Return the dimension of the finite element function space
 
46
  virtual unsigned int space_dimension() const
 
47
  {
 
48
    return 4;
 
49
  }
 
50
 
 
51
  /// Return the rank of the value space
 
52
  virtual unsigned int value_rank() const
 
53
  {
 
54
    return 0;
 
55
  }
 
56
 
 
57
  /// Return the dimension of the value space for axis i
 
58
  virtual unsigned int value_dimension(unsigned int i) const
 
59
  {
 
60
    return 1;
 
61
  }
 
62
 
 
63
  /// Evaluate basis function i at given point in cell
 
64
  virtual void evaluate_basis(unsigned int i,
 
65
                              double* values,
 
66
                              const double* coordinates,
 
67
                              const ufc::cell& c) const
 
68
  {
 
69
    // Extract vertex coordinates
 
70
    const double * const * element_coordinates = c.coordinates;
 
71
    
 
72
    // Compute Jacobian of affine map from reference cell
 
73
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
74
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
75
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
76
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
77
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
78
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
79
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
80
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
81
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
82
      
 
83
    // Compute sub determinants
 
84
    const double d00 = J_11*J_22 - J_12*J_21;
 
85
    const double d01 = J_12*J_20 - J_10*J_22;
 
86
    const double d02 = J_10*J_21 - J_11*J_20;
 
87
    
 
88
    const double d10 = J_02*J_21 - J_01*J_22;
 
89
    const double d11 = J_00*J_22 - J_02*J_20;
 
90
    const double d12 = J_01*J_20 - J_00*J_21;
 
91
    
 
92
    const double d20 = J_01*J_12 - J_02*J_11;
 
93
    const double d21 = J_02*J_10 - J_00*J_12;
 
94
    const double d22 = J_00*J_11 - J_01*J_10;
 
95
      
 
96
    // Compute determinant of Jacobian
 
97
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
98
    
 
99
    // Compute inverse of Jacobian
 
100
    
 
101
    // Compute constants
 
102
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
103
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
104
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
105
    
 
106
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
107
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
108
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
109
    
 
110
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
111
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
112
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
113
    
 
114
    // Get coordinates and map to the UFC reference element
 
115
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
116
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
117
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
118
    
 
119
    // Map coordinates to the reference cube
 
120
    if (std::abs(y + z - 1.0) < 1e-14)
 
121
      x = 1.0;
 
122
    else
 
123
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
124
    if (std::abs(z - 1.0) < 1e-14)
 
125
      y = -1.0;
 
126
    else
 
127
      y = 2.0 * y/(1.0 - z) - 1.0;
 
128
    z = 2.0 * z - 1.0;
 
129
    
 
130
    // Reset values
 
131
    *values = 0;
 
132
    
 
133
    // Map degree of freedom to element degree of freedom
 
134
    const unsigned int dof = i;
 
135
    
 
136
    // Generate scalings
 
137
    const double scalings_y_0 = 1;
 
138
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
139
    const double scalings_z_0 = 1;
 
140
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
141
    
 
142
    // Compute psitilde_a
 
143
    const double psitilde_a_0 = 1;
 
144
    const double psitilde_a_1 = x;
 
145
    
 
146
    // Compute psitilde_bs
 
147
    const double psitilde_bs_0_0 = 1;
 
148
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
149
    const double psitilde_bs_1_0 = 1;
 
150
    
 
151
    // Compute psitilde_cs
 
152
    const double psitilde_cs_00_0 = 1;
 
153
    const double psitilde_cs_00_1 = 2*z + 1;
 
154
    const double psitilde_cs_01_0 = 1;
 
155
    const double psitilde_cs_10_0 = 1;
 
156
    
 
157
    // Compute basisvalues
 
158
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
159
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
160
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
161
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
162
    
 
163
    // Table(s) of coefficients
 
164
    const static double coefficients0[4][4] = \
 
165
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
166
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
167
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
168
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
169
    
 
170
    // Extract relevant coefficients
 
171
    const double coeff0_0 = coefficients0[dof][0];
 
172
    const double coeff0_1 = coefficients0[dof][1];
 
173
    const double coeff0_2 = coefficients0[dof][2];
 
174
    const double coeff0_3 = coefficients0[dof][3];
 
175
    
 
176
    // Compute value(s)
 
177
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3;
 
178
  }
 
179
 
 
180
  /// Evaluate all basis functions at given point in cell
 
181
  virtual void evaluate_basis_all(double* values,
 
182
                                  const double* coordinates,
 
183
                                  const ufc::cell& c) const
 
184
  {
 
185
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
186
  }
 
187
 
 
188
  /// Evaluate order n derivatives of basis function i at given point in cell
 
189
  virtual void evaluate_basis_derivatives(unsigned int i,
 
190
                                          unsigned int n,
 
191
                                          double* values,
 
192
                                          const double* coordinates,
 
193
                                          const ufc::cell& c) const
 
194
  {
 
195
    // Extract vertex coordinates
 
196
    const double * const * element_coordinates = c.coordinates;
 
197
    
 
198
    // Compute Jacobian of affine map from reference cell
 
199
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
200
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
201
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
202
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
203
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
204
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
205
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
206
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
207
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
208
      
 
209
    // Compute sub determinants
 
210
    const double d00 = J_11*J_22 - J_12*J_21;
 
211
    const double d01 = J_12*J_20 - J_10*J_22;
 
212
    const double d02 = J_10*J_21 - J_11*J_20;
 
213
    
 
214
    const double d10 = J_02*J_21 - J_01*J_22;
 
215
    const double d11 = J_00*J_22 - J_02*J_20;
 
216
    const double d12 = J_01*J_20 - J_00*J_21;
 
217
    
 
218
    const double d20 = J_01*J_12 - J_02*J_11;
 
219
    const double d21 = J_02*J_10 - J_00*J_12;
 
220
    const double d22 = J_00*J_11 - J_01*J_10;
 
221
      
 
222
    // Compute determinant of Jacobian
 
223
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
224
    
 
225
    // Compute inverse of Jacobian
 
226
    
 
227
    // Compute constants
 
228
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
229
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
230
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
231
    
 
232
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
233
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
234
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
235
    
 
236
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
237
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
238
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
239
    
 
240
    // Get coordinates and map to the UFC reference element
 
241
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
242
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
243
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
244
    
 
245
    // Map coordinates to the reference cube
 
246
    if (std::abs(y + z - 1.0) < 1e-14)
 
247
      x = 1.0;
 
248
    else
 
249
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
250
    if (std::abs(z - 1.0) < 1e-14)
 
251
      y = -1.0;
 
252
    else
 
253
      y = 2.0 * y/(1.0 - z) - 1.0;
 
254
    z = 2.0 * z - 1.0;
 
255
    
 
256
    // Compute number of derivatives
 
257
    unsigned int num_derivatives = 1;
 
258
    
 
259
    for (unsigned int j = 0; j < n; j++)
 
260
      num_derivatives *= 3;
 
261
    
 
262
    
 
263
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
264
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
265
        
 
266
    for (unsigned int j = 0; j < num_derivatives; j++)
 
267
    {
 
268
      combinations[j] = new unsigned int [n];
 
269
      for (unsigned int k = 0; k < n; k++)
 
270
        combinations[j][k] = 0;
 
271
    }
 
272
        
 
273
    // Generate combinations of derivatives
 
274
    for (unsigned int row = 1; row < num_derivatives; row++)
 
275
    {
 
276
      for (unsigned int num = 0; num < row; num++)
 
277
      {
 
278
        for (unsigned int col = n-1; col+1 > 0; col--)
 
279
        {
 
280
          if (combinations[row][col] + 1 > 2)
 
281
            combinations[row][col] = 0;
 
282
          else
 
283
          {
 
284
            combinations[row][col] += 1;
 
285
            break;
 
286
          }
 
287
        }
 
288
      }
 
289
    }
 
290
    
 
291
    // Compute inverse of Jacobian
 
292
    const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
 
293
    
 
294
    // Declare transformation matrix
 
295
    // Declare pointer to two dimensional array and initialise
 
296
    double **transform = new double *[num_derivatives];
 
297
        
 
298
    for (unsigned int j = 0; j < num_derivatives; j++)
 
299
    {
 
300
      transform[j] = new double [num_derivatives];
 
301
      for (unsigned int k = 0; k < num_derivatives; k++)
 
302
        transform[j][k] = 1;
 
303
    }
 
304
    
 
305
    // Construct transformation matrix
 
306
    for (unsigned int row = 0; row < num_derivatives; row++)
 
307
    {
 
308
      for (unsigned int col = 0; col < num_derivatives; col++)
 
309
      {
 
310
        for (unsigned int k = 0; k < n; k++)
 
311
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
312
      }
 
313
    }
 
314
    
 
315
    // Reset values
 
316
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
317
      values[j] = 0;
 
318
    
 
319
    // Map degree of freedom to element degree of freedom
 
320
    const unsigned int dof = i;
 
321
    
 
322
    // Generate scalings
 
323
    const double scalings_y_0 = 1;
 
324
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
325
    const double scalings_z_0 = 1;
 
326
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
327
    
 
328
    // Compute psitilde_a
 
329
    const double psitilde_a_0 = 1;
 
330
    const double psitilde_a_1 = x;
 
331
    
 
332
    // Compute psitilde_bs
 
333
    const double psitilde_bs_0_0 = 1;
 
334
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
335
    const double psitilde_bs_1_0 = 1;
 
336
    
 
337
    // Compute psitilde_cs
 
338
    const double psitilde_cs_00_0 = 1;
 
339
    const double psitilde_cs_00_1 = 2*z + 1;
 
340
    const double psitilde_cs_01_0 = 1;
 
341
    const double psitilde_cs_10_0 = 1;
 
342
    
 
343
    // Compute basisvalues
 
344
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
345
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
346
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
347
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
348
    
 
349
    // Table(s) of coefficients
 
350
    const static double coefficients0[4][4] = \
 
351
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
352
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
353
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
354
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
355
    
 
356
    // Interesting (new) part
 
357
    // Tables of derivatives of the polynomial base (transpose)
 
358
    const static double dmats0[4][4] = \
 
359
    {{0, 0, 0, 0},
 
360
    {6.32455532033676, 0, 0, 0},
 
361
    {0, 0, 0, 0},
 
362
    {0, 0, 0, 0}};
 
363
    
 
364
    const static double dmats1[4][4] = \
 
365
    {{0, 0, 0, 0},
 
366
    {3.16227766016838, 0, 0, 0},
 
367
    {5.47722557505166, 0, 0, 0},
 
368
    {0, 0, 0, 0}};
 
369
    
 
370
    const static double dmats2[4][4] = \
 
371
    {{0, 0, 0, 0},
 
372
    {3.16227766016838, 0, 0, 0},
 
373
    {1.82574185835055, 0, 0, 0},
 
374
    {5.16397779494322, 0, 0, 0}};
 
375
    
 
376
    // Compute reference derivatives
 
377
    // Declare pointer to array of derivatives on FIAT element
 
378
    double *derivatives = new double [num_derivatives];
 
379
    
 
380
    // Declare coefficients
 
381
    double coeff0_0 = 0;
 
382
    double coeff0_1 = 0;
 
383
    double coeff0_2 = 0;
 
384
    double coeff0_3 = 0;
 
385
    
 
386
    // Declare new coefficients
 
387
    double new_coeff0_0 = 0;
 
388
    double new_coeff0_1 = 0;
 
389
    double new_coeff0_2 = 0;
 
390
    double new_coeff0_3 = 0;
 
391
    
 
392
    // Loop possible derivatives
 
393
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
394
    {
 
395
      // Get values from coefficients array
 
396
      new_coeff0_0 = coefficients0[dof][0];
 
397
      new_coeff0_1 = coefficients0[dof][1];
 
398
      new_coeff0_2 = coefficients0[dof][2];
 
399
      new_coeff0_3 = coefficients0[dof][3];
 
400
    
 
401
      // Loop derivative order
 
402
      for (unsigned int j = 0; j < n; j++)
 
403
      {
 
404
        // Update old coefficients
 
405
        coeff0_0 = new_coeff0_0;
 
406
        coeff0_1 = new_coeff0_1;
 
407
        coeff0_2 = new_coeff0_2;
 
408
        coeff0_3 = new_coeff0_3;
 
409
    
 
410
        if(combinations[deriv_num][j] == 0)
 
411
        {
 
412
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0];
 
413
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1];
 
414
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2];
 
415
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3];
 
416
        }
 
417
        if(combinations[deriv_num][j] == 1)
 
418
        {
 
419
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0];
 
420
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1];
 
421
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2];
 
422
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3];
 
423
        }
 
424
        if(combinations[deriv_num][j] == 2)
 
425
        {
 
426
          new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0];
 
427
          new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1];
 
428
          new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2];
 
429
          new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3];
 
430
        }
 
431
    
 
432
      }
 
433
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
434
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3;
 
435
    }
 
436
    
 
437
    // Transform derivatives back to physical element
 
438
    for (unsigned int row = 0; row < num_derivatives; row++)
 
439
    {
 
440
      for (unsigned int col = 0; col < num_derivatives; col++)
 
441
      {
 
442
        values[row] += transform[row][col]*derivatives[col];
 
443
      }
 
444
    }
 
445
    // Delete pointer to array of derivatives on FIAT element
 
446
    delete [] derivatives;
 
447
    
 
448
    // Delete pointer to array of combinations of derivatives and transform
 
449
    for (unsigned int row = 0; row < num_derivatives; row++)
 
450
    {
 
451
      delete [] combinations[row];
 
452
      delete [] transform[row];
 
453
    }
 
454
    
 
455
    delete [] combinations;
 
456
    delete [] transform;
 
457
  }
 
458
 
 
459
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
460
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
461
                                              double* values,
 
462
                                              const double* coordinates,
 
463
                                              const ufc::cell& c) const
 
464
  {
 
465
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
466
  }
 
467
 
 
468
  /// Evaluate linear functional for dof i on the function f
 
469
  virtual double evaluate_dof(unsigned int i,
 
470
                              const ufc::function& f,
 
471
                              const ufc::cell& c) const
 
472
  {
 
473
    // The reference points, direction and weights:
 
474
    const static double X[4][1][3] = {{{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}};
 
475
    const static double W[4][1] = {{1}, {1}, {1}, {1}};
 
476
    const static double D[4][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}};
 
477
    
 
478
    const double * const * x = c.coordinates;
 
479
    double result = 0.0;
 
480
    // Iterate over the points:
 
481
    // Evaluate basis functions for affine mapping
 
482
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
 
483
    const double w1 = X[i][0][0];
 
484
    const double w2 = X[i][0][1];
 
485
    const double w3 = X[i][0][2];
 
486
    
 
487
    // Compute affine mapping y = F(X)
 
488
    double y[3];
 
489
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
 
490
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
 
491
    y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
 
492
    
 
493
    // Evaluate function at physical points
 
494
    double values[1];
 
495
    f.evaluate(values, y, c);
 
496
    
 
497
    // Map function values using appropriate mapping
 
498
    // Affine map: Do nothing
 
499
    
 
500
    // Note that we do not map the weights (yet).
 
501
    
 
502
    // Take directional components
 
503
    for(int k = 0; k < 1; k++)
 
504
      result += values[k]*D[i][0][k];
 
505
    // Multiply by weights 
 
506
    result *= W[i][0];
 
507
    
 
508
    return result;
 
509
  }
 
510
 
 
511
  /// Evaluate linear functionals for all dofs on the function f
 
512
  virtual void evaluate_dofs(double* values,
 
513
                             const ufc::function& f,
 
514
                             const ufc::cell& c) const
 
515
  {
 
516
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
517
  }
 
518
 
 
519
  /// Interpolate vertex values from dof values
 
520
  virtual void interpolate_vertex_values(double* vertex_values,
 
521
                                         const double* dof_values,
 
522
                                         const ufc::cell& c) const
 
523
  {
 
524
    // Evaluate at vertices and use affine mapping
 
525
    vertex_values[0] = dof_values[0];
 
526
    vertex_values[1] = dof_values[1];
 
527
    vertex_values[2] = dof_values[2];
 
528
    vertex_values[3] = dof_values[3];
 
529
  }
 
530
 
 
531
  /// Return the number of sub elements (for a mixed element)
 
532
  virtual unsigned int num_sub_elements() const
 
533
  {
 
534
    return 1;
 
535
  }
 
536
 
 
537
  /// Create a new finite element for sub element i (for a mixed element)
 
538
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
539
  {
 
540
    return new UFC_Poisson3D_1BilinearForm_finite_element_0();
 
541
  }
 
542
 
 
543
};
 
544
 
 
545
/// This class defines the interface for a finite element.
 
546
 
 
547
class UFC_Poisson3D_1BilinearForm_finite_element_1: public ufc::finite_element
 
548
{
 
549
public:
 
550
 
 
551
  /// Constructor
 
552
  UFC_Poisson3D_1BilinearForm_finite_element_1() : ufc::finite_element()
 
553
  {
 
554
    // Do nothing
 
555
  }
 
556
 
 
557
  /// Destructor
 
558
  virtual ~UFC_Poisson3D_1BilinearForm_finite_element_1()
 
559
  {
 
560
    // Do nothing
 
561
  }
 
562
 
 
563
  /// Return a string identifying the finite element
 
564
  virtual const char* signature() const
 
565
  {
 
566
    return "Lagrange finite element of degree 1 on a tetrahedron";
 
567
  }
 
568
 
 
569
  /// Return the cell shape
 
570
  virtual ufc::shape cell_shape() const
 
571
  {
 
572
    return ufc::tetrahedron;
 
573
  }
 
574
 
 
575
  /// Return the dimension of the finite element function space
 
576
  virtual unsigned int space_dimension() const
 
577
  {
 
578
    return 4;
 
579
  }
 
580
 
 
581
  /// Return the rank of the value space
 
582
  virtual unsigned int value_rank() const
 
583
  {
 
584
    return 0;
 
585
  }
 
586
 
 
587
  /// Return the dimension of the value space for axis i
 
588
  virtual unsigned int value_dimension(unsigned int i) const
 
589
  {
 
590
    return 1;
 
591
  }
 
592
 
 
593
  /// Evaluate basis function i at given point in cell
 
594
  virtual void evaluate_basis(unsigned int i,
 
595
                              double* values,
 
596
                              const double* coordinates,
 
597
                              const ufc::cell& c) const
 
598
  {
 
599
    // Extract vertex coordinates
 
600
    const double * const * element_coordinates = c.coordinates;
 
601
    
 
602
    // Compute Jacobian of affine map from reference cell
 
603
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
604
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
605
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
606
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
607
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
608
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
609
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
610
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
611
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
612
      
 
613
    // Compute sub determinants
 
614
    const double d00 = J_11*J_22 - J_12*J_21;
 
615
    const double d01 = J_12*J_20 - J_10*J_22;
 
616
    const double d02 = J_10*J_21 - J_11*J_20;
 
617
    
 
618
    const double d10 = J_02*J_21 - J_01*J_22;
 
619
    const double d11 = J_00*J_22 - J_02*J_20;
 
620
    const double d12 = J_01*J_20 - J_00*J_21;
 
621
    
 
622
    const double d20 = J_01*J_12 - J_02*J_11;
 
623
    const double d21 = J_02*J_10 - J_00*J_12;
 
624
    const double d22 = J_00*J_11 - J_01*J_10;
 
625
      
 
626
    // Compute determinant of Jacobian
 
627
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
628
    
 
629
    // Compute inverse of Jacobian
 
630
    
 
631
    // Compute constants
 
632
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
633
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
634
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
635
    
 
636
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
637
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
638
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
639
    
 
640
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
641
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
642
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
643
    
 
644
    // Get coordinates and map to the UFC reference element
 
645
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
646
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
647
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
648
    
 
649
    // Map coordinates to the reference cube
 
650
    if (std::abs(y + z - 1.0) < 1e-14)
 
651
      x = 1.0;
 
652
    else
 
653
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
654
    if (std::abs(z - 1.0) < 1e-14)
 
655
      y = -1.0;
 
656
    else
 
657
      y = 2.0 * y/(1.0 - z) - 1.0;
 
658
    z = 2.0 * z - 1.0;
 
659
    
 
660
    // Reset values
 
661
    *values = 0;
 
662
    
 
663
    // Map degree of freedom to element degree of freedom
 
664
    const unsigned int dof = i;
 
665
    
 
666
    // Generate scalings
 
667
    const double scalings_y_0 = 1;
 
668
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
669
    const double scalings_z_0 = 1;
 
670
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
671
    
 
672
    // Compute psitilde_a
 
673
    const double psitilde_a_0 = 1;
 
674
    const double psitilde_a_1 = x;
 
675
    
 
676
    // Compute psitilde_bs
 
677
    const double psitilde_bs_0_0 = 1;
 
678
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
679
    const double psitilde_bs_1_0 = 1;
 
680
    
 
681
    // Compute psitilde_cs
 
682
    const double psitilde_cs_00_0 = 1;
 
683
    const double psitilde_cs_00_1 = 2*z + 1;
 
684
    const double psitilde_cs_01_0 = 1;
 
685
    const double psitilde_cs_10_0 = 1;
 
686
    
 
687
    // Compute basisvalues
 
688
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
689
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
690
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
691
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
692
    
 
693
    // Table(s) of coefficients
 
694
    const static double coefficients0[4][4] = \
 
695
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
696
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
697
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
698
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
699
    
 
700
    // Extract relevant coefficients
 
701
    const double coeff0_0 = coefficients0[dof][0];
 
702
    const double coeff0_1 = coefficients0[dof][1];
 
703
    const double coeff0_2 = coefficients0[dof][2];
 
704
    const double coeff0_3 = coefficients0[dof][3];
 
705
    
 
706
    // Compute value(s)
 
707
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3;
 
708
  }
 
709
 
 
710
  /// Evaluate all basis functions at given point in cell
 
711
  virtual void evaluate_basis_all(double* values,
 
712
                                  const double* coordinates,
 
713
                                  const ufc::cell& c) const
 
714
  {
 
715
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
716
  }
 
717
 
 
718
  /// Evaluate order n derivatives of basis function i at given point in cell
 
719
  virtual void evaluate_basis_derivatives(unsigned int i,
 
720
                                          unsigned int n,
 
721
                                          double* values,
 
722
                                          const double* coordinates,
 
723
                                          const ufc::cell& c) const
 
724
  {
 
725
    // Extract vertex coordinates
 
726
    const double * const * element_coordinates = c.coordinates;
 
727
    
 
728
    // Compute Jacobian of affine map from reference cell
 
729
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
730
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
731
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
732
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
733
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
734
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
735
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
736
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
737
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
738
      
 
739
    // Compute sub determinants
 
740
    const double d00 = J_11*J_22 - J_12*J_21;
 
741
    const double d01 = J_12*J_20 - J_10*J_22;
 
742
    const double d02 = J_10*J_21 - J_11*J_20;
 
743
    
 
744
    const double d10 = J_02*J_21 - J_01*J_22;
 
745
    const double d11 = J_00*J_22 - J_02*J_20;
 
746
    const double d12 = J_01*J_20 - J_00*J_21;
 
747
    
 
748
    const double d20 = J_01*J_12 - J_02*J_11;
 
749
    const double d21 = J_02*J_10 - J_00*J_12;
 
750
    const double d22 = J_00*J_11 - J_01*J_10;
 
751
      
 
752
    // Compute determinant of Jacobian
 
753
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
754
    
 
755
    // Compute inverse of Jacobian
 
756
    
 
757
    // Compute constants
 
758
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
759
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
760
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
761
    
 
762
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
763
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
764
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
765
    
 
766
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
767
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
768
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
769
    
 
770
    // Get coordinates and map to the UFC reference element
 
771
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
772
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
773
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
774
    
 
775
    // Map coordinates to the reference cube
 
776
    if (std::abs(y + z - 1.0) < 1e-14)
 
777
      x = 1.0;
 
778
    else
 
779
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
780
    if (std::abs(z - 1.0) < 1e-14)
 
781
      y = -1.0;
 
782
    else
 
783
      y = 2.0 * y/(1.0 - z) - 1.0;
 
784
    z = 2.0 * z - 1.0;
 
785
    
 
786
    // Compute number of derivatives
 
787
    unsigned int num_derivatives = 1;
 
788
    
 
789
    for (unsigned int j = 0; j < n; j++)
 
790
      num_derivatives *= 3;
 
791
    
 
792
    
 
793
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
794
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
795
        
 
796
    for (unsigned int j = 0; j < num_derivatives; j++)
 
797
    {
 
798
      combinations[j] = new unsigned int [n];
 
799
      for (unsigned int k = 0; k < n; k++)
 
800
        combinations[j][k] = 0;
 
801
    }
 
802
        
 
803
    // Generate combinations of derivatives
 
804
    for (unsigned int row = 1; row < num_derivatives; row++)
 
805
    {
 
806
      for (unsigned int num = 0; num < row; num++)
 
807
      {
 
808
        for (unsigned int col = n-1; col+1 > 0; col--)
 
809
        {
 
810
          if (combinations[row][col] + 1 > 2)
 
811
            combinations[row][col] = 0;
 
812
          else
 
813
          {
 
814
            combinations[row][col] += 1;
 
815
            break;
 
816
          }
 
817
        }
 
818
      }
 
819
    }
 
820
    
 
821
    // Compute inverse of Jacobian
 
822
    const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
 
823
    
 
824
    // Declare transformation matrix
 
825
    // Declare pointer to two dimensional array and initialise
 
826
    double **transform = new double *[num_derivatives];
 
827
        
 
828
    for (unsigned int j = 0; j < num_derivatives; j++)
 
829
    {
 
830
      transform[j] = new double [num_derivatives];
 
831
      for (unsigned int k = 0; k < num_derivatives; k++)
 
832
        transform[j][k] = 1;
 
833
    }
 
834
    
 
835
    // Construct transformation matrix
 
836
    for (unsigned int row = 0; row < num_derivatives; row++)
 
837
    {
 
838
      for (unsigned int col = 0; col < num_derivatives; col++)
 
839
      {
 
840
        for (unsigned int k = 0; k < n; k++)
 
841
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
842
      }
 
843
    }
 
844
    
 
845
    // Reset values
 
846
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
847
      values[j] = 0;
 
848
    
 
849
    // Map degree of freedom to element degree of freedom
 
850
    const unsigned int dof = i;
 
851
    
 
852
    // Generate scalings
 
853
    const double scalings_y_0 = 1;
 
854
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
855
    const double scalings_z_0 = 1;
 
856
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
857
    
 
858
    // Compute psitilde_a
 
859
    const double psitilde_a_0 = 1;
 
860
    const double psitilde_a_1 = x;
 
861
    
 
862
    // Compute psitilde_bs
 
863
    const double psitilde_bs_0_0 = 1;
 
864
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
865
    const double psitilde_bs_1_0 = 1;
 
866
    
 
867
    // Compute psitilde_cs
 
868
    const double psitilde_cs_00_0 = 1;
 
869
    const double psitilde_cs_00_1 = 2*z + 1;
 
870
    const double psitilde_cs_01_0 = 1;
 
871
    const double psitilde_cs_10_0 = 1;
 
872
    
 
873
    // Compute basisvalues
 
874
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
875
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
876
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
877
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
878
    
 
879
    // Table(s) of coefficients
 
880
    const static double coefficients0[4][4] = \
 
881
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
882
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
883
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
884
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
885
    
 
886
    // Interesting (new) part
 
887
    // Tables of derivatives of the polynomial base (transpose)
 
888
    const static double dmats0[4][4] = \
 
889
    {{0, 0, 0, 0},
 
890
    {6.32455532033676, 0, 0, 0},
 
891
    {0, 0, 0, 0},
 
892
    {0, 0, 0, 0}};
 
893
    
 
894
    const static double dmats1[4][4] = \
 
895
    {{0, 0, 0, 0},
 
896
    {3.16227766016838, 0, 0, 0},
 
897
    {5.47722557505166, 0, 0, 0},
 
898
    {0, 0, 0, 0}};
 
899
    
 
900
    const static double dmats2[4][4] = \
 
901
    {{0, 0, 0, 0},
 
902
    {3.16227766016838, 0, 0, 0},
 
903
    {1.82574185835055, 0, 0, 0},
 
904
    {5.16397779494322, 0, 0, 0}};
 
905
    
 
906
    // Compute reference derivatives
 
907
    // Declare pointer to array of derivatives on FIAT element
 
908
    double *derivatives = new double [num_derivatives];
 
909
    
 
910
    // Declare coefficients
 
911
    double coeff0_0 = 0;
 
912
    double coeff0_1 = 0;
 
913
    double coeff0_2 = 0;
 
914
    double coeff0_3 = 0;
 
915
    
 
916
    // Declare new coefficients
 
917
    double new_coeff0_0 = 0;
 
918
    double new_coeff0_1 = 0;
 
919
    double new_coeff0_2 = 0;
 
920
    double new_coeff0_3 = 0;
 
921
    
 
922
    // Loop possible derivatives
 
923
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
924
    {
 
925
      // Get values from coefficients array
 
926
      new_coeff0_0 = coefficients0[dof][0];
 
927
      new_coeff0_1 = coefficients0[dof][1];
 
928
      new_coeff0_2 = coefficients0[dof][2];
 
929
      new_coeff0_3 = coefficients0[dof][3];
 
930
    
 
931
      // Loop derivative order
 
932
      for (unsigned int j = 0; j < n; j++)
 
933
      {
 
934
        // Update old coefficients
 
935
        coeff0_0 = new_coeff0_0;
 
936
        coeff0_1 = new_coeff0_1;
 
937
        coeff0_2 = new_coeff0_2;
 
938
        coeff0_3 = new_coeff0_3;
 
939
    
 
940
        if(combinations[deriv_num][j] == 0)
 
941
        {
 
942
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0];
 
943
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1];
 
944
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2];
 
945
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3];
 
946
        }
 
947
        if(combinations[deriv_num][j] == 1)
 
948
        {
 
949
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0];
 
950
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1];
 
951
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2];
 
952
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3];
 
953
        }
 
954
        if(combinations[deriv_num][j] == 2)
 
955
        {
 
956
          new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0];
 
957
          new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1];
 
958
          new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2];
 
959
          new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3];
 
960
        }
 
961
    
 
962
      }
 
963
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
964
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3;
 
965
    }
 
966
    
 
967
    // Transform derivatives back to physical element
 
968
    for (unsigned int row = 0; row < num_derivatives; row++)
 
969
    {
 
970
      for (unsigned int col = 0; col < num_derivatives; col++)
 
971
      {
 
972
        values[row] += transform[row][col]*derivatives[col];
 
973
      }
 
974
    }
 
975
    // Delete pointer to array of derivatives on FIAT element
 
976
    delete [] derivatives;
 
977
    
 
978
    // Delete pointer to array of combinations of derivatives and transform
 
979
    for (unsigned int row = 0; row < num_derivatives; row++)
 
980
    {
 
981
      delete [] combinations[row];
 
982
      delete [] transform[row];
 
983
    }
 
984
    
 
985
    delete [] combinations;
 
986
    delete [] transform;
 
987
  }
 
988
 
 
989
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
990
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
991
                                              double* values,
 
992
                                              const double* coordinates,
 
993
                                              const ufc::cell& c) const
 
994
  {
 
995
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
996
  }
 
997
 
 
998
  /// Evaluate linear functional for dof i on the function f
 
999
  virtual double evaluate_dof(unsigned int i,
 
1000
                              const ufc::function& f,
 
1001
                              const ufc::cell& c) const
 
1002
  {
 
1003
    // The reference points, direction and weights:
 
1004
    const static double X[4][1][3] = {{{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}};
 
1005
    const static double W[4][1] = {{1}, {1}, {1}, {1}};
 
1006
    const static double D[4][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}};
 
1007
    
 
1008
    const double * const * x = c.coordinates;
 
1009
    double result = 0.0;
 
1010
    // Iterate over the points:
 
1011
    // Evaluate basis functions for affine mapping
 
1012
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
 
1013
    const double w1 = X[i][0][0];
 
1014
    const double w2 = X[i][0][1];
 
1015
    const double w3 = X[i][0][2];
 
1016
    
 
1017
    // Compute affine mapping y = F(X)
 
1018
    double y[3];
 
1019
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
 
1020
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
 
1021
    y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
 
1022
    
 
1023
    // Evaluate function at physical points
 
1024
    double values[1];
 
1025
    f.evaluate(values, y, c);
 
1026
    
 
1027
    // Map function values using appropriate mapping
 
1028
    // Affine map: Do nothing
 
1029
    
 
1030
    // Note that we do not map the weights (yet).
 
1031
    
 
1032
    // Take directional components
 
1033
    for(int k = 0; k < 1; k++)
 
1034
      result += values[k]*D[i][0][k];
 
1035
    // Multiply by weights 
 
1036
    result *= W[i][0];
 
1037
    
 
1038
    return result;
 
1039
  }
 
1040
 
 
1041
  /// Evaluate linear functionals for all dofs on the function f
 
1042
  virtual void evaluate_dofs(double* values,
 
1043
                             const ufc::function& f,
 
1044
                             const ufc::cell& c) const
 
1045
  {
 
1046
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1047
  }
 
1048
 
 
1049
  /// Interpolate vertex values from dof values
 
1050
  virtual void interpolate_vertex_values(double* vertex_values,
 
1051
                                         const double* dof_values,
 
1052
                                         const ufc::cell& c) const
 
1053
  {
 
1054
    // Evaluate at vertices and use affine mapping
 
1055
    vertex_values[0] = dof_values[0];
 
1056
    vertex_values[1] = dof_values[1];
 
1057
    vertex_values[2] = dof_values[2];
 
1058
    vertex_values[3] = dof_values[3];
 
1059
  }
 
1060
 
 
1061
  /// Return the number of sub elements (for a mixed element)
 
1062
  virtual unsigned int num_sub_elements() const
 
1063
  {
 
1064
    return 1;
 
1065
  }
 
1066
 
 
1067
  /// Create a new finite element for sub element i (for a mixed element)
 
1068
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
1069
  {
 
1070
    return new UFC_Poisson3D_1BilinearForm_finite_element_1();
 
1071
  }
 
1072
 
 
1073
};
 
1074
 
 
1075
/// This class defines the interface for a local-to-global mapping of
 
1076
/// degrees of freedom (dofs).
 
1077
 
 
1078
class UFC_Poisson3D_1BilinearForm_dof_map_0: public ufc::dof_map
 
1079
{
 
1080
private:
 
1081
 
 
1082
  unsigned int __global_dimension;
 
1083
 
 
1084
public:
 
1085
 
 
1086
  /// Constructor
 
1087
  UFC_Poisson3D_1BilinearForm_dof_map_0() : ufc::dof_map()
 
1088
  {
 
1089
    __global_dimension = 0;
 
1090
  }
 
1091
 
 
1092
  /// Destructor
 
1093
  virtual ~UFC_Poisson3D_1BilinearForm_dof_map_0()
 
1094
  {
 
1095
    // Do nothing
 
1096
  }
 
1097
 
 
1098
  /// Return a string identifying the dof map
 
1099
  virtual const char* signature() const
 
1100
  {
 
1101
    return "FFC dof map for Lagrange finite element of degree 1 on a tetrahedron";
 
1102
  }
 
1103
 
 
1104
  /// Return true iff mesh entities of topological dimension d are needed
 
1105
  virtual bool needs_mesh_entities(unsigned int d) const
 
1106
  {
 
1107
    switch ( d )
 
1108
    {
 
1109
    case 0:
 
1110
      return true;
 
1111
      break;
 
1112
    case 1:
 
1113
      return false;
 
1114
      break;
 
1115
    case 2:
 
1116
      return false;
 
1117
      break;
 
1118
    case 3:
 
1119
      return false;
 
1120
      break;
 
1121
    }
 
1122
    return false;
 
1123
  }
 
1124
 
 
1125
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1126
  virtual bool init_mesh(const ufc::mesh& m)
 
1127
  {
 
1128
    __global_dimension = m.num_entities[0];
 
1129
    return false;
 
1130
  }
 
1131
 
 
1132
  /// Initialize dof map for given cell
 
1133
  virtual void init_cell(const ufc::mesh& m,
 
1134
                         const ufc::cell& c)
 
1135
  {
 
1136
    // Do nothing
 
1137
  }
 
1138
 
 
1139
  /// Finish initialization of dof map for cells
 
1140
  virtual void init_cell_finalize()
 
1141
  {
 
1142
    // Do nothing
 
1143
  }
 
1144
 
 
1145
  /// Return the dimension of the global finite element function space
 
1146
  virtual unsigned int global_dimension() const
 
1147
  {
 
1148
    return __global_dimension;
 
1149
  }
 
1150
 
 
1151
  /// Return the dimension of the local finite element function space
 
1152
  virtual unsigned int local_dimension() const
 
1153
  {
 
1154
    return 4;
 
1155
  }
 
1156
 
 
1157
  // Return the geometric dimension of the coordinates this dof map provides
 
1158
  virtual unsigned int geometric_dimension() const
 
1159
  {
 
1160
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1161
  }
 
1162
 
 
1163
  /// Return the number of dofs on each cell facet
 
1164
  virtual unsigned int num_facet_dofs() const
 
1165
  {
 
1166
    return 3;
 
1167
  }
 
1168
 
 
1169
  /// Return the number of dofs associated with each cell entity of dimension d
 
1170
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1171
  {
 
1172
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1173
  }
 
1174
 
 
1175
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1176
  virtual void tabulate_dofs(unsigned int* dofs,
 
1177
                             const ufc::mesh& m,
 
1178
                             const ufc::cell& c) const
 
1179
  {
 
1180
    dofs[0] = c.entity_indices[0][0];
 
1181
    dofs[1] = c.entity_indices[0][1];
 
1182
    dofs[2] = c.entity_indices[0][2];
 
1183
    dofs[3] = c.entity_indices[0][3];
 
1184
  }
 
1185
 
 
1186
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1187
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1188
                                   unsigned int facet) const
 
1189
  {
 
1190
    switch ( facet )
 
1191
    {
 
1192
    case 0:
 
1193
      dofs[0] = 1;
 
1194
      dofs[1] = 2;
 
1195
      dofs[2] = 3;
 
1196
      break;
 
1197
    case 1:
 
1198
      dofs[0] = 0;
 
1199
      dofs[1] = 2;
 
1200
      dofs[2] = 3;
 
1201
      break;
 
1202
    case 2:
 
1203
      dofs[0] = 0;
 
1204
      dofs[1] = 1;
 
1205
      dofs[2] = 3;
 
1206
      break;
 
1207
    case 3:
 
1208
      dofs[0] = 0;
 
1209
      dofs[1] = 1;
 
1210
      dofs[2] = 2;
 
1211
      break;
 
1212
    }
 
1213
  }
 
1214
 
 
1215
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1216
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1217
                                    unsigned int d, unsigned int i) const
 
1218
  {
 
1219
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1220
  }
 
1221
 
 
1222
  /// Tabulate the coordinates of all dofs on a cell
 
1223
  virtual void tabulate_coordinates(double** coordinates,
 
1224
                                    const ufc::cell& c) const
 
1225
  {
 
1226
    const double * const * x = c.coordinates;
 
1227
    coordinates[0][0] = x[0][0];
 
1228
    coordinates[0][1] = x[0][1];
 
1229
    coordinates[0][2] = x[0][2];
 
1230
    coordinates[1][0] = x[1][0];
 
1231
    coordinates[1][1] = x[1][1];
 
1232
    coordinates[1][2] = x[1][2];
 
1233
    coordinates[2][0] = x[2][0];
 
1234
    coordinates[2][1] = x[2][1];
 
1235
    coordinates[2][2] = x[2][2];
 
1236
    coordinates[3][0] = x[3][0];
 
1237
    coordinates[3][1] = x[3][1];
 
1238
    coordinates[3][2] = x[3][2];
 
1239
  }
 
1240
 
 
1241
  /// Return the number of sub dof maps (for a mixed element)
 
1242
  virtual unsigned int num_sub_dof_maps() const
 
1243
  {
 
1244
    return 1;
 
1245
  }
 
1246
 
 
1247
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1248
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1249
  {
 
1250
    return new UFC_Poisson3D_1BilinearForm_dof_map_0();
 
1251
  }
 
1252
 
 
1253
};
 
1254
 
 
1255
/// This class defines the interface for a local-to-global mapping of
 
1256
/// degrees of freedom (dofs).
 
1257
 
 
1258
class UFC_Poisson3D_1BilinearForm_dof_map_1: public ufc::dof_map
 
1259
{
 
1260
private:
 
1261
 
 
1262
  unsigned int __global_dimension;
 
1263
 
 
1264
public:
 
1265
 
 
1266
  /// Constructor
 
1267
  UFC_Poisson3D_1BilinearForm_dof_map_1() : ufc::dof_map()
 
1268
  {
 
1269
    __global_dimension = 0;
 
1270
  }
 
1271
 
 
1272
  /// Destructor
 
1273
  virtual ~UFC_Poisson3D_1BilinearForm_dof_map_1()
 
1274
  {
 
1275
    // Do nothing
 
1276
  }
 
1277
 
 
1278
  /// Return a string identifying the dof map
 
1279
  virtual const char* signature() const
 
1280
  {
 
1281
    return "FFC dof map for Lagrange finite element of degree 1 on a tetrahedron";
 
1282
  }
 
1283
 
 
1284
  /// Return true iff mesh entities of topological dimension d are needed
 
1285
  virtual bool needs_mesh_entities(unsigned int d) const
 
1286
  {
 
1287
    switch ( d )
 
1288
    {
 
1289
    case 0:
 
1290
      return true;
 
1291
      break;
 
1292
    case 1:
 
1293
      return false;
 
1294
      break;
 
1295
    case 2:
 
1296
      return false;
 
1297
      break;
 
1298
    case 3:
 
1299
      return false;
 
1300
      break;
 
1301
    }
 
1302
    return false;
 
1303
  }
 
1304
 
 
1305
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
1306
  virtual bool init_mesh(const ufc::mesh& m)
 
1307
  {
 
1308
    __global_dimension = m.num_entities[0];
 
1309
    return false;
 
1310
  }
 
1311
 
 
1312
  /// Initialize dof map for given cell
 
1313
  virtual void init_cell(const ufc::mesh& m,
 
1314
                         const ufc::cell& c)
 
1315
  {
 
1316
    // Do nothing
 
1317
  }
 
1318
 
 
1319
  /// Finish initialization of dof map for cells
 
1320
  virtual void init_cell_finalize()
 
1321
  {
 
1322
    // Do nothing
 
1323
  }
 
1324
 
 
1325
  /// Return the dimension of the global finite element function space
 
1326
  virtual unsigned int global_dimension() const
 
1327
  {
 
1328
    return __global_dimension;
 
1329
  }
 
1330
 
 
1331
  /// Return the dimension of the local finite element function space
 
1332
  virtual unsigned int local_dimension() const
 
1333
  {
 
1334
    return 4;
 
1335
  }
 
1336
 
 
1337
  // Return the geometric dimension of the coordinates this dof map provides
 
1338
  virtual unsigned int geometric_dimension() const
 
1339
  {
 
1340
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1341
  }
 
1342
 
 
1343
  /// Return the number of dofs on each cell facet
 
1344
  virtual unsigned int num_facet_dofs() const
 
1345
  {
 
1346
    return 3;
 
1347
  }
 
1348
 
 
1349
  /// Return the number of dofs associated with each cell entity of dimension d
 
1350
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
1351
  {
 
1352
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1353
  }
 
1354
 
 
1355
  /// Tabulate the local-to-global mapping of dofs on a cell
 
1356
  virtual void tabulate_dofs(unsigned int* dofs,
 
1357
                             const ufc::mesh& m,
 
1358
                             const ufc::cell& c) const
 
1359
  {
 
1360
    dofs[0] = c.entity_indices[0][0];
 
1361
    dofs[1] = c.entity_indices[0][1];
 
1362
    dofs[2] = c.entity_indices[0][2];
 
1363
    dofs[3] = c.entity_indices[0][3];
 
1364
  }
 
1365
 
 
1366
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
1367
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
1368
                                   unsigned int facet) const
 
1369
  {
 
1370
    switch ( facet )
 
1371
    {
 
1372
    case 0:
 
1373
      dofs[0] = 1;
 
1374
      dofs[1] = 2;
 
1375
      dofs[2] = 3;
 
1376
      break;
 
1377
    case 1:
 
1378
      dofs[0] = 0;
 
1379
      dofs[1] = 2;
 
1380
      dofs[2] = 3;
 
1381
      break;
 
1382
    case 2:
 
1383
      dofs[0] = 0;
 
1384
      dofs[1] = 1;
 
1385
      dofs[2] = 3;
 
1386
      break;
 
1387
    case 3:
 
1388
      dofs[0] = 0;
 
1389
      dofs[1] = 1;
 
1390
      dofs[2] = 2;
 
1391
      break;
 
1392
    }
 
1393
  }
 
1394
 
 
1395
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
1396
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
1397
                                    unsigned int d, unsigned int i) const
 
1398
  {
 
1399
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
1400
  }
 
1401
 
 
1402
  /// Tabulate the coordinates of all dofs on a cell
 
1403
  virtual void tabulate_coordinates(double** coordinates,
 
1404
                                    const ufc::cell& c) const
 
1405
  {
 
1406
    const double * const * x = c.coordinates;
 
1407
    coordinates[0][0] = x[0][0];
 
1408
    coordinates[0][1] = x[0][1];
 
1409
    coordinates[0][2] = x[0][2];
 
1410
    coordinates[1][0] = x[1][0];
 
1411
    coordinates[1][1] = x[1][1];
 
1412
    coordinates[1][2] = x[1][2];
 
1413
    coordinates[2][0] = x[2][0];
 
1414
    coordinates[2][1] = x[2][1];
 
1415
    coordinates[2][2] = x[2][2];
 
1416
    coordinates[3][0] = x[3][0];
 
1417
    coordinates[3][1] = x[3][1];
 
1418
    coordinates[3][2] = x[3][2];
 
1419
  }
 
1420
 
 
1421
  /// Return the number of sub dof maps (for a mixed element)
 
1422
  virtual unsigned int num_sub_dof_maps() const
 
1423
  {
 
1424
    return 1;
 
1425
  }
 
1426
 
 
1427
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
1428
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
1429
  {
 
1430
    return new UFC_Poisson3D_1BilinearForm_dof_map_1();
 
1431
  }
 
1432
 
 
1433
};
 
1434
 
 
1435
/// This class defines the interface for the tabulation of the cell
 
1436
/// tensor corresponding to the local contribution to a form from
 
1437
/// the integral over a cell.
 
1438
 
 
1439
class UFC_Poisson3D_1BilinearForm_cell_integral_0: public ufc::cell_integral
 
1440
{
 
1441
public:
 
1442
 
 
1443
  /// Constructor
 
1444
  UFC_Poisson3D_1BilinearForm_cell_integral_0() : ufc::cell_integral()
 
1445
  {
 
1446
    // Do nothing
 
1447
  }
 
1448
 
 
1449
  /// Destructor
 
1450
  virtual ~UFC_Poisson3D_1BilinearForm_cell_integral_0()
 
1451
  {
 
1452
    // Do nothing
 
1453
  }
 
1454
 
 
1455
  /// Tabulate the tensor for the contribution from a local cell
 
1456
  virtual void tabulate_tensor(double* A,
 
1457
                               const double * const * w,
 
1458
                               const ufc::cell& c) const
 
1459
  {
 
1460
    // Extract vertex coordinates
 
1461
    const double * const * x = c.coordinates;
 
1462
    
 
1463
    // Compute Jacobian of affine map from reference cell
 
1464
    const double J_00 = x[1][0] - x[0][0];
 
1465
    const double J_01 = x[2][0] - x[0][0];
 
1466
    const double J_02 = x[3][0] - x[0][0];
 
1467
    const double J_10 = x[1][1] - x[0][1];
 
1468
    const double J_11 = x[2][1] - x[0][1];
 
1469
    const double J_12 = x[3][1] - x[0][1];
 
1470
    const double J_20 = x[1][2] - x[0][2];
 
1471
    const double J_21 = x[2][2] - x[0][2];
 
1472
    const double J_22 = x[3][2] - x[0][2];
 
1473
      
 
1474
    // Compute sub determinants
 
1475
    const double d_00 = J_11*J_22 - J_12*J_21;
 
1476
    const double d_01 = J_12*J_20 - J_10*J_22;
 
1477
    const double d_02 = J_10*J_21 - J_11*J_20;
 
1478
    
 
1479
    const double d_10 = J_02*J_21 - J_01*J_22;
 
1480
    const double d_11 = J_00*J_22 - J_02*J_20;
 
1481
    const double d_12 = J_01*J_20 - J_00*J_21;
 
1482
    
 
1483
    const double d_20 = J_01*J_12 - J_02*J_11;
 
1484
    const double d_21 = J_02*J_10 - J_00*J_12;
 
1485
    const double d_22 = J_00*J_11 - J_01*J_10;
 
1486
      
 
1487
    // Compute determinant of Jacobian
 
1488
    double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
 
1489
      
 
1490
    // Compute inverse of Jacobian
 
1491
    const double Jinv_00 = d_00 / detJ;
 
1492
    const double Jinv_01 = d_10 / detJ;
 
1493
    const double Jinv_02 = d_20 / detJ;
 
1494
    const double Jinv_10 = d_01 / detJ;
 
1495
    const double Jinv_11 = d_11 / detJ;
 
1496
    const double Jinv_12 = d_21 / detJ;
 
1497
    const double Jinv_20 = d_02 / detJ;
 
1498
    const double Jinv_21 = d_12 / detJ;
 
1499
    const double Jinv_22 = d_22 / detJ;
 
1500
    
 
1501
    // Set scale factor
 
1502
    const double det = std::abs(detJ);
 
1503
    
 
1504
    // Compute geometry tensors
 
1505
    const double G0_0_0 = det*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01 + Jinv_02*Jinv_02);
 
1506
    const double G0_0_1 = det*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11 + Jinv_02*Jinv_12);
 
1507
    const double G0_0_2 = det*(Jinv_00*Jinv_20 + Jinv_01*Jinv_21 + Jinv_02*Jinv_22);
 
1508
    const double G0_1_0 = det*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01 + Jinv_12*Jinv_02);
 
1509
    const double G0_1_1 = det*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11 + Jinv_12*Jinv_12);
 
1510
    const double G0_1_2 = det*(Jinv_10*Jinv_20 + Jinv_11*Jinv_21 + Jinv_12*Jinv_22);
 
1511
    const double G0_2_0 = det*(Jinv_20*Jinv_00 + Jinv_21*Jinv_01 + Jinv_22*Jinv_02);
 
1512
    const double G0_2_1 = det*(Jinv_20*Jinv_10 + Jinv_21*Jinv_11 + Jinv_22*Jinv_12);
 
1513
    const double G0_2_2 = det*(Jinv_20*Jinv_20 + Jinv_21*Jinv_21 + Jinv_22*Jinv_22);
 
1514
    
 
1515
    // Compute element tensor
 
1516
    A[0] = 0.166666666666666*G0_0_0 + 0.166666666666666*G0_0_1 + 0.166666666666666*G0_0_2 + 0.166666666666666*G0_1_0 + 0.166666666666666*G0_1_1 + 0.166666666666666*G0_1_2 + 0.166666666666666*G0_2_0 + 0.166666666666666*G0_2_1 + 0.166666666666666*G0_2_2;
 
1517
    A[1] = -0.166666666666666*G0_0_0 - 0.166666666666666*G0_1_0 - 0.166666666666666*G0_2_0;
 
1518
    A[2] = -0.166666666666666*G0_0_1 - 0.166666666666666*G0_1_1 - 0.166666666666666*G0_2_1;
 
1519
    A[3] = -0.166666666666666*G0_0_2 - 0.166666666666666*G0_1_2 - 0.166666666666666*G0_2_2;
 
1520
    A[4] = -0.166666666666666*G0_0_0 - 0.166666666666666*G0_0_1 - 0.166666666666666*G0_0_2;
 
1521
    A[5] = 0.166666666666666*G0_0_0;
 
1522
    A[6] = 0.166666666666666*G0_0_1;
 
1523
    A[7] = 0.166666666666666*G0_0_2;
 
1524
    A[8] = -0.166666666666666*G0_1_0 - 0.166666666666666*G0_1_1 - 0.166666666666666*G0_1_2;
 
1525
    A[9] = 0.166666666666666*G0_1_0;
 
1526
    A[10] = 0.166666666666666*G0_1_1;
 
1527
    A[11] = 0.166666666666666*G0_1_2;
 
1528
    A[12] = -0.166666666666666*G0_2_0 - 0.166666666666666*G0_2_1 - 0.166666666666666*G0_2_2;
 
1529
    A[13] = 0.166666666666666*G0_2_0;
 
1530
    A[14] = 0.166666666666666*G0_2_1;
 
1531
    A[15] = 0.166666666666666*G0_2_2;
 
1532
  }
 
1533
 
 
1534
};
 
1535
 
 
1536
/// This class defines the interface for the assembly of the global
 
1537
/// tensor corresponding to a form with r + n arguments, that is, a
 
1538
/// mapping
 
1539
///
 
1540
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
1541
///
 
1542
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
1543
/// global tensor A is defined by
 
1544
///
 
1545
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
1546
///
 
1547
/// where each argument Vj represents the application to the
 
1548
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
1549
/// fixed functions (coefficients).
 
1550
 
 
1551
class UFC_Poisson3D_1BilinearForm: public ufc::form
 
1552
{
 
1553
public:
 
1554
 
 
1555
  /// Constructor
 
1556
  UFC_Poisson3D_1BilinearForm() : ufc::form()
 
1557
  {
 
1558
    // Do nothing
 
1559
  }
 
1560
 
 
1561
  /// Destructor
 
1562
  virtual ~UFC_Poisson3D_1BilinearForm()
 
1563
  {
 
1564
    // Do nothing
 
1565
  }
 
1566
 
 
1567
  /// Return a string identifying the form
 
1568
  virtual const char* signature() const
 
1569
  {
 
1570
    return "(dXa0/dxb0)(dXa1/dxb0) | ((d/dXa0)vi0)*((d/dXa1)vi1)*dX(0)";
 
1571
  }
 
1572
 
 
1573
  /// Return the rank of the global tensor (r)
 
1574
  virtual unsigned int rank() const
 
1575
  {
 
1576
    return 2;
 
1577
  }
 
1578
 
 
1579
  /// Return the number of coefficients (n)
 
1580
  virtual unsigned int num_coefficients() const
 
1581
  {
 
1582
    return 0;
 
1583
  }
 
1584
 
 
1585
  /// Return the number of cell integrals
 
1586
  virtual unsigned int num_cell_integrals() const
 
1587
  {
 
1588
    return 1;
 
1589
  }
 
1590
  
 
1591
  /// Return the number of exterior facet integrals
 
1592
  virtual unsigned int num_exterior_facet_integrals() const
 
1593
  {
 
1594
    return 0;
 
1595
  }
 
1596
  
 
1597
  /// Return the number of interior facet integrals
 
1598
  virtual unsigned int num_interior_facet_integrals() const
 
1599
  {
 
1600
    return 0;
 
1601
  }
 
1602
    
 
1603
  /// Create a new finite element for argument function i
 
1604
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
1605
  {
 
1606
    switch ( i )
 
1607
    {
 
1608
    case 0:
 
1609
      return new UFC_Poisson3D_1BilinearForm_finite_element_0();
 
1610
      break;
 
1611
    case 1:
 
1612
      return new UFC_Poisson3D_1BilinearForm_finite_element_1();
 
1613
      break;
 
1614
    }
 
1615
    return 0;
 
1616
  }
 
1617
  
 
1618
  /// Create a new dof map for argument function i
 
1619
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
1620
  {
 
1621
    switch ( i )
 
1622
    {
 
1623
    case 0:
 
1624
      return new UFC_Poisson3D_1BilinearForm_dof_map_0();
 
1625
      break;
 
1626
    case 1:
 
1627
      return new UFC_Poisson3D_1BilinearForm_dof_map_1();
 
1628
      break;
 
1629
    }
 
1630
    return 0;
 
1631
  }
 
1632
 
 
1633
  /// Create a new cell integral on sub domain i
 
1634
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
1635
  {
 
1636
    return new UFC_Poisson3D_1BilinearForm_cell_integral_0();
 
1637
  }
 
1638
 
 
1639
  /// Create a new exterior facet integral on sub domain i
 
1640
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
1641
  {
 
1642
    return 0;
 
1643
  }
 
1644
 
 
1645
  /// Create a new interior facet integral on sub domain i
 
1646
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
1647
  {
 
1648
    return 0;
 
1649
  }
 
1650
 
 
1651
};
 
1652
 
 
1653
/// This class defines the interface for a finite element.
 
1654
 
 
1655
class UFC_Poisson3D_1LinearForm_finite_element_0: public ufc::finite_element
 
1656
{
 
1657
public:
 
1658
 
 
1659
  /// Constructor
 
1660
  UFC_Poisson3D_1LinearForm_finite_element_0() : ufc::finite_element()
 
1661
  {
 
1662
    // Do nothing
 
1663
  }
 
1664
 
 
1665
  /// Destructor
 
1666
  virtual ~UFC_Poisson3D_1LinearForm_finite_element_0()
 
1667
  {
 
1668
    // Do nothing
 
1669
  }
 
1670
 
 
1671
  /// Return a string identifying the finite element
 
1672
  virtual const char* signature() const
 
1673
  {
 
1674
    return "Lagrange finite element of degree 1 on a tetrahedron";
 
1675
  }
 
1676
 
 
1677
  /// Return the cell shape
 
1678
  virtual ufc::shape cell_shape() const
 
1679
  {
 
1680
    return ufc::tetrahedron;
 
1681
  }
 
1682
 
 
1683
  /// Return the dimension of the finite element function space
 
1684
  virtual unsigned int space_dimension() const
 
1685
  {
 
1686
    return 4;
 
1687
  }
 
1688
 
 
1689
  /// Return the rank of the value space
 
1690
  virtual unsigned int value_rank() const
 
1691
  {
 
1692
    return 0;
 
1693
  }
 
1694
 
 
1695
  /// Return the dimension of the value space for axis i
 
1696
  virtual unsigned int value_dimension(unsigned int i) const
 
1697
  {
 
1698
    return 1;
 
1699
  }
 
1700
 
 
1701
  /// Evaluate basis function i at given point in cell
 
1702
  virtual void evaluate_basis(unsigned int i,
 
1703
                              double* values,
 
1704
                              const double* coordinates,
 
1705
                              const ufc::cell& c) const
 
1706
  {
 
1707
    // Extract vertex coordinates
 
1708
    const double * const * element_coordinates = c.coordinates;
 
1709
    
 
1710
    // Compute Jacobian of affine map from reference cell
 
1711
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
1712
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
1713
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
1714
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
1715
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
1716
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
1717
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
1718
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
1719
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
1720
      
 
1721
    // Compute sub determinants
 
1722
    const double d00 = J_11*J_22 - J_12*J_21;
 
1723
    const double d01 = J_12*J_20 - J_10*J_22;
 
1724
    const double d02 = J_10*J_21 - J_11*J_20;
 
1725
    
 
1726
    const double d10 = J_02*J_21 - J_01*J_22;
 
1727
    const double d11 = J_00*J_22 - J_02*J_20;
 
1728
    const double d12 = J_01*J_20 - J_00*J_21;
 
1729
    
 
1730
    const double d20 = J_01*J_12 - J_02*J_11;
 
1731
    const double d21 = J_02*J_10 - J_00*J_12;
 
1732
    const double d22 = J_00*J_11 - J_01*J_10;
 
1733
      
 
1734
    // Compute determinant of Jacobian
 
1735
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
1736
    
 
1737
    // Compute inverse of Jacobian
 
1738
    
 
1739
    // Compute constants
 
1740
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
1741
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
1742
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
1743
    
 
1744
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
1745
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
1746
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
1747
    
 
1748
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
1749
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
1750
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
1751
    
 
1752
    // Get coordinates and map to the UFC reference element
 
1753
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
1754
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
1755
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
1756
    
 
1757
    // Map coordinates to the reference cube
 
1758
    if (std::abs(y + z - 1.0) < 1e-14)
 
1759
      x = 1.0;
 
1760
    else
 
1761
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
1762
    if (std::abs(z - 1.0) < 1e-14)
 
1763
      y = -1.0;
 
1764
    else
 
1765
      y = 2.0 * y/(1.0 - z) - 1.0;
 
1766
    z = 2.0 * z - 1.0;
 
1767
    
 
1768
    // Reset values
 
1769
    *values = 0;
 
1770
    
 
1771
    // Map degree of freedom to element degree of freedom
 
1772
    const unsigned int dof = i;
 
1773
    
 
1774
    // Generate scalings
 
1775
    const double scalings_y_0 = 1;
 
1776
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
1777
    const double scalings_z_0 = 1;
 
1778
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
1779
    
 
1780
    // Compute psitilde_a
 
1781
    const double psitilde_a_0 = 1;
 
1782
    const double psitilde_a_1 = x;
 
1783
    
 
1784
    // Compute psitilde_bs
 
1785
    const double psitilde_bs_0_0 = 1;
 
1786
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
1787
    const double psitilde_bs_1_0 = 1;
 
1788
    
 
1789
    // Compute psitilde_cs
 
1790
    const double psitilde_cs_00_0 = 1;
 
1791
    const double psitilde_cs_00_1 = 2*z + 1;
 
1792
    const double psitilde_cs_01_0 = 1;
 
1793
    const double psitilde_cs_10_0 = 1;
 
1794
    
 
1795
    // Compute basisvalues
 
1796
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
1797
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
1798
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
1799
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
1800
    
 
1801
    // Table(s) of coefficients
 
1802
    const static double coefficients0[4][4] = \
 
1803
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
1804
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
1805
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
1806
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
1807
    
 
1808
    // Extract relevant coefficients
 
1809
    const double coeff0_0 = coefficients0[dof][0];
 
1810
    const double coeff0_1 = coefficients0[dof][1];
 
1811
    const double coeff0_2 = coefficients0[dof][2];
 
1812
    const double coeff0_3 = coefficients0[dof][3];
 
1813
    
 
1814
    // Compute value(s)
 
1815
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3;
 
1816
  }
 
1817
 
 
1818
  /// Evaluate all basis functions at given point in cell
 
1819
  virtual void evaluate_basis_all(double* values,
 
1820
                                  const double* coordinates,
 
1821
                                  const ufc::cell& c) const
 
1822
  {
 
1823
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
1824
  }
 
1825
 
 
1826
  /// Evaluate order n derivatives of basis function i at given point in cell
 
1827
  virtual void evaluate_basis_derivatives(unsigned int i,
 
1828
                                          unsigned int n,
 
1829
                                          double* values,
 
1830
                                          const double* coordinates,
 
1831
                                          const ufc::cell& c) const
 
1832
  {
 
1833
    // Extract vertex coordinates
 
1834
    const double * const * element_coordinates = c.coordinates;
 
1835
    
 
1836
    // Compute Jacobian of affine map from reference cell
 
1837
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
1838
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
1839
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
1840
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
1841
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
1842
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
1843
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
1844
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
1845
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
1846
      
 
1847
    // Compute sub determinants
 
1848
    const double d00 = J_11*J_22 - J_12*J_21;
 
1849
    const double d01 = J_12*J_20 - J_10*J_22;
 
1850
    const double d02 = J_10*J_21 - J_11*J_20;
 
1851
    
 
1852
    const double d10 = J_02*J_21 - J_01*J_22;
 
1853
    const double d11 = J_00*J_22 - J_02*J_20;
 
1854
    const double d12 = J_01*J_20 - J_00*J_21;
 
1855
    
 
1856
    const double d20 = J_01*J_12 - J_02*J_11;
 
1857
    const double d21 = J_02*J_10 - J_00*J_12;
 
1858
    const double d22 = J_00*J_11 - J_01*J_10;
 
1859
      
 
1860
    // Compute determinant of Jacobian
 
1861
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
1862
    
 
1863
    // Compute inverse of Jacobian
 
1864
    
 
1865
    // Compute constants
 
1866
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
1867
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
1868
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
1869
    
 
1870
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
1871
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
1872
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
1873
    
 
1874
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
1875
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
1876
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
1877
    
 
1878
    // Get coordinates and map to the UFC reference element
 
1879
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
1880
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
1881
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
1882
    
 
1883
    // Map coordinates to the reference cube
 
1884
    if (std::abs(y + z - 1.0) < 1e-14)
 
1885
      x = 1.0;
 
1886
    else
 
1887
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
1888
    if (std::abs(z - 1.0) < 1e-14)
 
1889
      y = -1.0;
 
1890
    else
 
1891
      y = 2.0 * y/(1.0 - z) - 1.0;
 
1892
    z = 2.0 * z - 1.0;
 
1893
    
 
1894
    // Compute number of derivatives
 
1895
    unsigned int num_derivatives = 1;
 
1896
    
 
1897
    for (unsigned int j = 0; j < n; j++)
 
1898
      num_derivatives *= 3;
 
1899
    
 
1900
    
 
1901
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
1902
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
1903
        
 
1904
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1905
    {
 
1906
      combinations[j] = new unsigned int [n];
 
1907
      for (unsigned int k = 0; k < n; k++)
 
1908
        combinations[j][k] = 0;
 
1909
    }
 
1910
        
 
1911
    // Generate combinations of derivatives
 
1912
    for (unsigned int row = 1; row < num_derivatives; row++)
 
1913
    {
 
1914
      for (unsigned int num = 0; num < row; num++)
 
1915
      {
 
1916
        for (unsigned int col = n-1; col+1 > 0; col--)
 
1917
        {
 
1918
          if (combinations[row][col] + 1 > 2)
 
1919
            combinations[row][col] = 0;
 
1920
          else
 
1921
          {
 
1922
            combinations[row][col] += 1;
 
1923
            break;
 
1924
          }
 
1925
        }
 
1926
      }
 
1927
    }
 
1928
    
 
1929
    // Compute inverse of Jacobian
 
1930
    const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
 
1931
    
 
1932
    // Declare transformation matrix
 
1933
    // Declare pointer to two dimensional array and initialise
 
1934
    double **transform = new double *[num_derivatives];
 
1935
        
 
1936
    for (unsigned int j = 0; j < num_derivatives; j++)
 
1937
    {
 
1938
      transform[j] = new double [num_derivatives];
 
1939
      for (unsigned int k = 0; k < num_derivatives; k++)
 
1940
        transform[j][k] = 1;
 
1941
    }
 
1942
    
 
1943
    // Construct transformation matrix
 
1944
    for (unsigned int row = 0; row < num_derivatives; row++)
 
1945
    {
 
1946
      for (unsigned int col = 0; col < num_derivatives; col++)
 
1947
      {
 
1948
        for (unsigned int k = 0; k < n; k++)
 
1949
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
1950
      }
 
1951
    }
 
1952
    
 
1953
    // Reset values
 
1954
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
1955
      values[j] = 0;
 
1956
    
 
1957
    // Map degree of freedom to element degree of freedom
 
1958
    const unsigned int dof = i;
 
1959
    
 
1960
    // Generate scalings
 
1961
    const double scalings_y_0 = 1;
 
1962
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
1963
    const double scalings_z_0 = 1;
 
1964
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
1965
    
 
1966
    // Compute psitilde_a
 
1967
    const double psitilde_a_0 = 1;
 
1968
    const double psitilde_a_1 = x;
 
1969
    
 
1970
    // Compute psitilde_bs
 
1971
    const double psitilde_bs_0_0 = 1;
 
1972
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
1973
    const double psitilde_bs_1_0 = 1;
 
1974
    
 
1975
    // Compute psitilde_cs
 
1976
    const double psitilde_cs_00_0 = 1;
 
1977
    const double psitilde_cs_00_1 = 2*z + 1;
 
1978
    const double psitilde_cs_01_0 = 1;
 
1979
    const double psitilde_cs_10_0 = 1;
 
1980
    
 
1981
    // Compute basisvalues
 
1982
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
1983
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
1984
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
1985
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
1986
    
 
1987
    // Table(s) of coefficients
 
1988
    const static double coefficients0[4][4] = \
 
1989
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
1990
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
1991
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
1992
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
1993
    
 
1994
    // Interesting (new) part
 
1995
    // Tables of derivatives of the polynomial base (transpose)
 
1996
    const static double dmats0[4][4] = \
 
1997
    {{0, 0, 0, 0},
 
1998
    {6.32455532033676, 0, 0, 0},
 
1999
    {0, 0, 0, 0},
 
2000
    {0, 0, 0, 0}};
 
2001
    
 
2002
    const static double dmats1[4][4] = \
 
2003
    {{0, 0, 0, 0},
 
2004
    {3.16227766016838, 0, 0, 0},
 
2005
    {5.47722557505166, 0, 0, 0},
 
2006
    {0, 0, 0, 0}};
 
2007
    
 
2008
    const static double dmats2[4][4] = \
 
2009
    {{0, 0, 0, 0},
 
2010
    {3.16227766016838, 0, 0, 0},
 
2011
    {1.82574185835055, 0, 0, 0},
 
2012
    {5.16397779494322, 0, 0, 0}};
 
2013
    
 
2014
    // Compute reference derivatives
 
2015
    // Declare pointer to array of derivatives on FIAT element
 
2016
    double *derivatives = new double [num_derivatives];
 
2017
    
 
2018
    // Declare coefficients
 
2019
    double coeff0_0 = 0;
 
2020
    double coeff0_1 = 0;
 
2021
    double coeff0_2 = 0;
 
2022
    double coeff0_3 = 0;
 
2023
    
 
2024
    // Declare new coefficients
 
2025
    double new_coeff0_0 = 0;
 
2026
    double new_coeff0_1 = 0;
 
2027
    double new_coeff0_2 = 0;
 
2028
    double new_coeff0_3 = 0;
 
2029
    
 
2030
    // Loop possible derivatives
 
2031
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
2032
    {
 
2033
      // Get values from coefficients array
 
2034
      new_coeff0_0 = coefficients0[dof][0];
 
2035
      new_coeff0_1 = coefficients0[dof][1];
 
2036
      new_coeff0_2 = coefficients0[dof][2];
 
2037
      new_coeff0_3 = coefficients0[dof][3];
 
2038
    
 
2039
      // Loop derivative order
 
2040
      for (unsigned int j = 0; j < n; j++)
 
2041
      {
 
2042
        // Update old coefficients
 
2043
        coeff0_0 = new_coeff0_0;
 
2044
        coeff0_1 = new_coeff0_1;
 
2045
        coeff0_2 = new_coeff0_2;
 
2046
        coeff0_3 = new_coeff0_3;
 
2047
    
 
2048
        if(combinations[deriv_num][j] == 0)
 
2049
        {
 
2050
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0];
 
2051
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1];
 
2052
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2];
 
2053
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3];
 
2054
        }
 
2055
        if(combinations[deriv_num][j] == 1)
 
2056
        {
 
2057
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0];
 
2058
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1];
 
2059
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2];
 
2060
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3];
 
2061
        }
 
2062
        if(combinations[deriv_num][j] == 2)
 
2063
        {
 
2064
          new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0];
 
2065
          new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1];
 
2066
          new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2];
 
2067
          new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3];
 
2068
        }
 
2069
    
 
2070
      }
 
2071
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
2072
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3;
 
2073
    }
 
2074
    
 
2075
    // Transform derivatives back to physical element
 
2076
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2077
    {
 
2078
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2079
      {
 
2080
        values[row] += transform[row][col]*derivatives[col];
 
2081
      }
 
2082
    }
 
2083
    // Delete pointer to array of derivatives on FIAT element
 
2084
    delete [] derivatives;
 
2085
    
 
2086
    // Delete pointer to array of combinations of derivatives and transform
 
2087
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2088
    {
 
2089
      delete [] combinations[row];
 
2090
      delete [] transform[row];
 
2091
    }
 
2092
    
 
2093
    delete [] combinations;
 
2094
    delete [] transform;
 
2095
  }
 
2096
 
 
2097
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
2098
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
2099
                                              double* values,
 
2100
                                              const double* coordinates,
 
2101
                                              const ufc::cell& c) const
 
2102
  {
 
2103
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
2104
  }
 
2105
 
 
2106
  /// Evaluate linear functional for dof i on the function f
 
2107
  virtual double evaluate_dof(unsigned int i,
 
2108
                              const ufc::function& f,
 
2109
                              const ufc::cell& c) const
 
2110
  {
 
2111
    // The reference points, direction and weights:
 
2112
    const static double X[4][1][3] = {{{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}};
 
2113
    const static double W[4][1] = {{1}, {1}, {1}, {1}};
 
2114
    const static double D[4][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}};
 
2115
    
 
2116
    const double * const * x = c.coordinates;
 
2117
    double result = 0.0;
 
2118
    // Iterate over the points:
 
2119
    // Evaluate basis functions for affine mapping
 
2120
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
 
2121
    const double w1 = X[i][0][0];
 
2122
    const double w2 = X[i][0][1];
 
2123
    const double w3 = X[i][0][2];
 
2124
    
 
2125
    // Compute affine mapping y = F(X)
 
2126
    double y[3];
 
2127
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
 
2128
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
 
2129
    y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
 
2130
    
 
2131
    // Evaluate function at physical points
 
2132
    double values[1];
 
2133
    f.evaluate(values, y, c);
 
2134
    
 
2135
    // Map function values using appropriate mapping
 
2136
    // Affine map: Do nothing
 
2137
    
 
2138
    // Note that we do not map the weights (yet).
 
2139
    
 
2140
    // Take directional components
 
2141
    for(int k = 0; k < 1; k++)
 
2142
      result += values[k]*D[i][0][k];
 
2143
    // Multiply by weights 
 
2144
    result *= W[i][0];
 
2145
    
 
2146
    return result;
 
2147
  }
 
2148
 
 
2149
  /// Evaluate linear functionals for all dofs on the function f
 
2150
  virtual void evaluate_dofs(double* values,
 
2151
                             const ufc::function& f,
 
2152
                             const ufc::cell& c) const
 
2153
  {
 
2154
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2155
  }
 
2156
 
 
2157
  /// Interpolate vertex values from dof values
 
2158
  virtual void interpolate_vertex_values(double* vertex_values,
 
2159
                                         const double* dof_values,
 
2160
                                         const ufc::cell& c) const
 
2161
  {
 
2162
    // Evaluate at vertices and use affine mapping
 
2163
    vertex_values[0] = dof_values[0];
 
2164
    vertex_values[1] = dof_values[1];
 
2165
    vertex_values[2] = dof_values[2];
 
2166
    vertex_values[3] = dof_values[3];
 
2167
  }
 
2168
 
 
2169
  /// Return the number of sub elements (for a mixed element)
 
2170
  virtual unsigned int num_sub_elements() const
 
2171
  {
 
2172
    return 1;
 
2173
  }
 
2174
 
 
2175
  /// Create a new finite element for sub element i (for a mixed element)
 
2176
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
2177
  {
 
2178
    return new UFC_Poisson3D_1LinearForm_finite_element_0();
 
2179
  }
 
2180
 
 
2181
};
 
2182
 
 
2183
/// This class defines the interface for a finite element.
 
2184
 
 
2185
class UFC_Poisson3D_1LinearForm_finite_element_1: public ufc::finite_element
 
2186
{
 
2187
public:
 
2188
 
 
2189
  /// Constructor
 
2190
  UFC_Poisson3D_1LinearForm_finite_element_1() : ufc::finite_element()
 
2191
  {
 
2192
    // Do nothing
 
2193
  }
 
2194
 
 
2195
  /// Destructor
 
2196
  virtual ~UFC_Poisson3D_1LinearForm_finite_element_1()
 
2197
  {
 
2198
    // Do nothing
 
2199
  }
 
2200
 
 
2201
  /// Return a string identifying the finite element
 
2202
  virtual const char* signature() const
 
2203
  {
 
2204
    return "Lagrange finite element of degree 1 on a tetrahedron";
 
2205
  }
 
2206
 
 
2207
  /// Return the cell shape
 
2208
  virtual ufc::shape cell_shape() const
 
2209
  {
 
2210
    return ufc::tetrahedron;
 
2211
  }
 
2212
 
 
2213
  /// Return the dimension of the finite element function space
 
2214
  virtual unsigned int space_dimension() const
 
2215
  {
 
2216
    return 4;
 
2217
  }
 
2218
 
 
2219
  /// Return the rank of the value space
 
2220
  virtual unsigned int value_rank() const
 
2221
  {
 
2222
    return 0;
 
2223
  }
 
2224
 
 
2225
  /// Return the dimension of the value space for axis i
 
2226
  virtual unsigned int value_dimension(unsigned int i) const
 
2227
  {
 
2228
    return 1;
 
2229
  }
 
2230
 
 
2231
  /// Evaluate basis function i at given point in cell
 
2232
  virtual void evaluate_basis(unsigned int i,
 
2233
                              double* values,
 
2234
                              const double* coordinates,
 
2235
                              const ufc::cell& c) const
 
2236
  {
 
2237
    // Extract vertex coordinates
 
2238
    const double * const * element_coordinates = c.coordinates;
 
2239
    
 
2240
    // Compute Jacobian of affine map from reference cell
 
2241
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
2242
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
2243
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
2244
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
2245
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
2246
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
2247
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
2248
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
2249
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
2250
      
 
2251
    // Compute sub determinants
 
2252
    const double d00 = J_11*J_22 - J_12*J_21;
 
2253
    const double d01 = J_12*J_20 - J_10*J_22;
 
2254
    const double d02 = J_10*J_21 - J_11*J_20;
 
2255
    
 
2256
    const double d10 = J_02*J_21 - J_01*J_22;
 
2257
    const double d11 = J_00*J_22 - J_02*J_20;
 
2258
    const double d12 = J_01*J_20 - J_00*J_21;
 
2259
    
 
2260
    const double d20 = J_01*J_12 - J_02*J_11;
 
2261
    const double d21 = J_02*J_10 - J_00*J_12;
 
2262
    const double d22 = J_00*J_11 - J_01*J_10;
 
2263
      
 
2264
    // Compute determinant of Jacobian
 
2265
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
2266
    
 
2267
    // Compute inverse of Jacobian
 
2268
    
 
2269
    // Compute constants
 
2270
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
2271
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
2272
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
2273
    
 
2274
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
2275
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
2276
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
2277
    
 
2278
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
2279
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
2280
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
2281
    
 
2282
    // Get coordinates and map to the UFC reference element
 
2283
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
2284
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
2285
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
2286
    
 
2287
    // Map coordinates to the reference cube
 
2288
    if (std::abs(y + z - 1.0) < 1e-14)
 
2289
      x = 1.0;
 
2290
    else
 
2291
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
2292
    if (std::abs(z - 1.0) < 1e-14)
 
2293
      y = -1.0;
 
2294
    else
 
2295
      y = 2.0 * y/(1.0 - z) - 1.0;
 
2296
    z = 2.0 * z - 1.0;
 
2297
    
 
2298
    // Reset values
 
2299
    *values = 0;
 
2300
    
 
2301
    // Map degree of freedom to element degree of freedom
 
2302
    const unsigned int dof = i;
 
2303
    
 
2304
    // Generate scalings
 
2305
    const double scalings_y_0 = 1;
 
2306
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2307
    const double scalings_z_0 = 1;
 
2308
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
2309
    
 
2310
    // Compute psitilde_a
 
2311
    const double psitilde_a_0 = 1;
 
2312
    const double psitilde_a_1 = x;
 
2313
    
 
2314
    // Compute psitilde_bs
 
2315
    const double psitilde_bs_0_0 = 1;
 
2316
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2317
    const double psitilde_bs_1_0 = 1;
 
2318
    
 
2319
    // Compute psitilde_cs
 
2320
    const double psitilde_cs_00_0 = 1;
 
2321
    const double psitilde_cs_00_1 = 2*z + 1;
 
2322
    const double psitilde_cs_01_0 = 1;
 
2323
    const double psitilde_cs_10_0 = 1;
 
2324
    
 
2325
    // Compute basisvalues
 
2326
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
2327
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
2328
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
2329
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
2330
    
 
2331
    // Table(s) of coefficients
 
2332
    const static double coefficients0[4][4] = \
 
2333
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
2334
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
2335
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
2336
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
2337
    
 
2338
    // Extract relevant coefficients
 
2339
    const double coeff0_0 = coefficients0[dof][0];
 
2340
    const double coeff0_1 = coefficients0[dof][1];
 
2341
    const double coeff0_2 = coefficients0[dof][2];
 
2342
    const double coeff0_3 = coefficients0[dof][3];
 
2343
    
 
2344
    // Compute value(s)
 
2345
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3;
 
2346
  }
 
2347
 
 
2348
  /// Evaluate all basis functions at given point in cell
 
2349
  virtual void evaluate_basis_all(double* values,
 
2350
                                  const double* coordinates,
 
2351
                                  const ufc::cell& c) const
 
2352
  {
 
2353
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
 
2354
  }
 
2355
 
 
2356
  /// Evaluate order n derivatives of basis function i at given point in cell
 
2357
  virtual void evaluate_basis_derivatives(unsigned int i,
 
2358
                                          unsigned int n,
 
2359
                                          double* values,
 
2360
                                          const double* coordinates,
 
2361
                                          const ufc::cell& c) const
 
2362
  {
 
2363
    // Extract vertex coordinates
 
2364
    const double * const * element_coordinates = c.coordinates;
 
2365
    
 
2366
    // Compute Jacobian of affine map from reference cell
 
2367
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
 
2368
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
 
2369
    const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];
 
2370
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
 
2371
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
 
2372
    const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
 
2373
    const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
 
2374
    const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
 
2375
    const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
 
2376
      
 
2377
    // Compute sub determinants
 
2378
    const double d00 = J_11*J_22 - J_12*J_21;
 
2379
    const double d01 = J_12*J_20 - J_10*J_22;
 
2380
    const double d02 = J_10*J_21 - J_11*J_20;
 
2381
    
 
2382
    const double d10 = J_02*J_21 - J_01*J_22;
 
2383
    const double d11 = J_00*J_22 - J_02*J_20;
 
2384
    const double d12 = J_01*J_20 - J_00*J_21;
 
2385
    
 
2386
    const double d20 = J_01*J_12 - J_02*J_11;
 
2387
    const double d21 = J_02*J_10 - J_00*J_12;
 
2388
    const double d22 = J_00*J_11 - J_01*J_10;
 
2389
      
 
2390
    // Compute determinant of Jacobian
 
2391
    double detJ = J_00*d00 + J_10*d10 + J_20*d20;
 
2392
    
 
2393
    // Compute inverse of Jacobian
 
2394
    
 
2395
    // Compute constants
 
2396
    const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \
 
2397
                    + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \
 
2398
                    + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);
 
2399
    
 
2400
    const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \
 
2401
                    + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \
 
2402
                    + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);
 
2403
    
 
2404
    const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \
 
2405
                    + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \
 
2406
                    + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);
 
2407
    
 
2408
    // Get coordinates and map to the UFC reference element
 
2409
    double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
 
2410
    double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
 
2411
    double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;
 
2412
    
 
2413
    // Map coordinates to the reference cube
 
2414
    if (std::abs(y + z - 1.0) < 1e-14)
 
2415
      x = 1.0;
 
2416
    else
 
2417
      x = -2.0 * x/(y + z - 1.0) - 1.0;
 
2418
    if (std::abs(z - 1.0) < 1e-14)
 
2419
      y = -1.0;
 
2420
    else
 
2421
      y = 2.0 * y/(1.0 - z) - 1.0;
 
2422
    z = 2.0 * z - 1.0;
 
2423
    
 
2424
    // Compute number of derivatives
 
2425
    unsigned int num_derivatives = 1;
 
2426
    
 
2427
    for (unsigned int j = 0; j < n; j++)
 
2428
      num_derivatives *= 3;
 
2429
    
 
2430
    
 
2431
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
 
2432
    unsigned int **combinations = new unsigned int *[num_derivatives];
 
2433
        
 
2434
    for (unsigned int j = 0; j < num_derivatives; j++)
 
2435
    {
 
2436
      combinations[j] = new unsigned int [n];
 
2437
      for (unsigned int k = 0; k < n; k++)
 
2438
        combinations[j][k] = 0;
 
2439
    }
 
2440
        
 
2441
    // Generate combinations of derivatives
 
2442
    for (unsigned int row = 1; row < num_derivatives; row++)
 
2443
    {
 
2444
      for (unsigned int num = 0; num < row; num++)
 
2445
      {
 
2446
        for (unsigned int col = n-1; col+1 > 0; col--)
 
2447
        {
 
2448
          if (combinations[row][col] + 1 > 2)
 
2449
            combinations[row][col] = 0;
 
2450
          else
 
2451
          {
 
2452
            combinations[row][col] += 1;
 
2453
            break;
 
2454
          }
 
2455
        }
 
2456
      }
 
2457
    }
 
2458
    
 
2459
    // Compute inverse of Jacobian
 
2460
    const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};
 
2461
    
 
2462
    // Declare transformation matrix
 
2463
    // Declare pointer to two dimensional array and initialise
 
2464
    double **transform = new double *[num_derivatives];
 
2465
        
 
2466
    for (unsigned int j = 0; j < num_derivatives; j++)
 
2467
    {
 
2468
      transform[j] = new double [num_derivatives];
 
2469
      for (unsigned int k = 0; k < num_derivatives; k++)
 
2470
        transform[j][k] = 1;
 
2471
    }
 
2472
    
 
2473
    // Construct transformation matrix
 
2474
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2475
    {
 
2476
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2477
      {
 
2478
        for (unsigned int k = 0; k < n; k++)
 
2479
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
 
2480
      }
 
2481
    }
 
2482
    
 
2483
    // Reset values
 
2484
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
 
2485
      values[j] = 0;
 
2486
    
 
2487
    // Map degree of freedom to element degree of freedom
 
2488
    const unsigned int dof = i;
 
2489
    
 
2490
    // Generate scalings
 
2491
    const double scalings_y_0 = 1;
 
2492
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
 
2493
    const double scalings_z_0 = 1;
 
2494
    const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);
 
2495
    
 
2496
    // Compute psitilde_a
 
2497
    const double psitilde_a_0 = 1;
 
2498
    const double psitilde_a_1 = x;
 
2499
    
 
2500
    // Compute psitilde_bs
 
2501
    const double psitilde_bs_0_0 = 1;
 
2502
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
 
2503
    const double psitilde_bs_1_0 = 1;
 
2504
    
 
2505
    // Compute psitilde_cs
 
2506
    const double psitilde_cs_00_0 = 1;
 
2507
    const double psitilde_cs_00_1 = 2*z + 1;
 
2508
    const double psitilde_cs_01_0 = 1;
 
2509
    const double psitilde_cs_10_0 = 1;
 
2510
    
 
2511
    // Compute basisvalues
 
2512
    const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;
 
2513
    const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;
 
2514
    const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;
 
2515
    const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;
 
2516
    
 
2517
    // Table(s) of coefficients
 
2518
    const static double coefficients0[4][4] = \
 
2519
    {{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},
 
2520
    {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},
 
2521
    {0.288675134594813, 0, 0.210818510677892, -0.074535599249993},
 
2522
    {0.288675134594813, 0, 0, 0.223606797749979}};
 
2523
    
 
2524
    // Interesting (new) part
 
2525
    // Tables of derivatives of the polynomial base (transpose)
 
2526
    const static double dmats0[4][4] = \
 
2527
    {{0, 0, 0, 0},
 
2528
    {6.32455532033676, 0, 0, 0},
 
2529
    {0, 0, 0, 0},
 
2530
    {0, 0, 0, 0}};
 
2531
    
 
2532
    const static double dmats1[4][4] = \
 
2533
    {{0, 0, 0, 0},
 
2534
    {3.16227766016838, 0, 0, 0},
 
2535
    {5.47722557505166, 0, 0, 0},
 
2536
    {0, 0, 0, 0}};
 
2537
    
 
2538
    const static double dmats2[4][4] = \
 
2539
    {{0, 0, 0, 0},
 
2540
    {3.16227766016838, 0, 0, 0},
 
2541
    {1.82574185835055, 0, 0, 0},
 
2542
    {5.16397779494322, 0, 0, 0}};
 
2543
    
 
2544
    // Compute reference derivatives
 
2545
    // Declare pointer to array of derivatives on FIAT element
 
2546
    double *derivatives = new double [num_derivatives];
 
2547
    
 
2548
    // Declare coefficients
 
2549
    double coeff0_0 = 0;
 
2550
    double coeff0_1 = 0;
 
2551
    double coeff0_2 = 0;
 
2552
    double coeff0_3 = 0;
 
2553
    
 
2554
    // Declare new coefficients
 
2555
    double new_coeff0_0 = 0;
 
2556
    double new_coeff0_1 = 0;
 
2557
    double new_coeff0_2 = 0;
 
2558
    double new_coeff0_3 = 0;
 
2559
    
 
2560
    // Loop possible derivatives
 
2561
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
 
2562
    {
 
2563
      // Get values from coefficients array
 
2564
      new_coeff0_0 = coefficients0[dof][0];
 
2565
      new_coeff0_1 = coefficients0[dof][1];
 
2566
      new_coeff0_2 = coefficients0[dof][2];
 
2567
      new_coeff0_3 = coefficients0[dof][3];
 
2568
    
 
2569
      // Loop derivative order
 
2570
      for (unsigned int j = 0; j < n; j++)
 
2571
      {
 
2572
        // Update old coefficients
 
2573
        coeff0_0 = new_coeff0_0;
 
2574
        coeff0_1 = new_coeff0_1;
 
2575
        coeff0_2 = new_coeff0_2;
 
2576
        coeff0_3 = new_coeff0_3;
 
2577
    
 
2578
        if(combinations[deriv_num][j] == 0)
 
2579
        {
 
2580
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0];
 
2581
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1];
 
2582
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2];
 
2583
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3];
 
2584
        }
 
2585
        if(combinations[deriv_num][j] == 1)
 
2586
        {
 
2587
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0];
 
2588
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1];
 
2589
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2];
 
2590
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3];
 
2591
        }
 
2592
        if(combinations[deriv_num][j] == 2)
 
2593
        {
 
2594
          new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0];
 
2595
          new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1];
 
2596
          new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2];
 
2597
          new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3];
 
2598
        }
 
2599
    
 
2600
      }
 
2601
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
 
2602
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3;
 
2603
    }
 
2604
    
 
2605
    // Transform derivatives back to physical element
 
2606
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2607
    {
 
2608
      for (unsigned int col = 0; col < num_derivatives; col++)
 
2609
      {
 
2610
        values[row] += transform[row][col]*derivatives[col];
 
2611
      }
 
2612
    }
 
2613
    // Delete pointer to array of derivatives on FIAT element
 
2614
    delete [] derivatives;
 
2615
    
 
2616
    // Delete pointer to array of combinations of derivatives and transform
 
2617
    for (unsigned int row = 0; row < num_derivatives; row++)
 
2618
    {
 
2619
      delete [] combinations[row];
 
2620
      delete [] transform[row];
 
2621
    }
 
2622
    
 
2623
    delete [] combinations;
 
2624
    delete [] transform;
 
2625
  }
 
2626
 
 
2627
  /// Evaluate order n derivatives of all basis functions at given point in cell
 
2628
  virtual void evaluate_basis_derivatives_all(unsigned int n,
 
2629
                                              double* values,
 
2630
                                              const double* coordinates,
 
2631
                                              const ufc::cell& c) const
 
2632
  {
 
2633
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
 
2634
  }
 
2635
 
 
2636
  /// Evaluate linear functional for dof i on the function f
 
2637
  virtual double evaluate_dof(unsigned int i,
 
2638
                              const ufc::function& f,
 
2639
                              const ufc::cell& c) const
 
2640
  {
 
2641
    // The reference points, direction and weights:
 
2642
    const static double X[4][1][3] = {{{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}};
 
2643
    const static double W[4][1] = {{1}, {1}, {1}, {1}};
 
2644
    const static double D[4][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}};
 
2645
    
 
2646
    const double * const * x = c.coordinates;
 
2647
    double result = 0.0;
 
2648
    // Iterate over the points:
 
2649
    // Evaluate basis functions for affine mapping
 
2650
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];
 
2651
    const double w1 = X[i][0][0];
 
2652
    const double w2 = X[i][0][1];
 
2653
    const double w3 = X[i][0][2];
 
2654
    
 
2655
    // Compute affine mapping y = F(X)
 
2656
    double y[3];
 
2657
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
 
2658
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
 
2659
    y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];
 
2660
    
 
2661
    // Evaluate function at physical points
 
2662
    double values[1];
 
2663
    f.evaluate(values, y, c);
 
2664
    
 
2665
    // Map function values using appropriate mapping
 
2666
    // Affine map: Do nothing
 
2667
    
 
2668
    // Note that we do not map the weights (yet).
 
2669
    
 
2670
    // Take directional components
 
2671
    for(int k = 0; k < 1; k++)
 
2672
      result += values[k]*D[i][0][k];
 
2673
    // Multiply by weights 
 
2674
    result *= W[i][0];
 
2675
    
 
2676
    return result;
 
2677
  }
 
2678
 
 
2679
  /// Evaluate linear functionals for all dofs on the function f
 
2680
  virtual void evaluate_dofs(double* values,
 
2681
                             const ufc::function& f,
 
2682
                             const ufc::cell& c) const
 
2683
  {
 
2684
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2685
  }
 
2686
 
 
2687
  /// Interpolate vertex values from dof values
 
2688
  virtual void interpolate_vertex_values(double* vertex_values,
 
2689
                                         const double* dof_values,
 
2690
                                         const ufc::cell& c) const
 
2691
  {
 
2692
    // Evaluate at vertices and use affine mapping
 
2693
    vertex_values[0] = dof_values[0];
 
2694
    vertex_values[1] = dof_values[1];
 
2695
    vertex_values[2] = dof_values[2];
 
2696
    vertex_values[3] = dof_values[3];
 
2697
  }
 
2698
 
 
2699
  /// Return the number of sub elements (for a mixed element)
 
2700
  virtual unsigned int num_sub_elements() const
 
2701
  {
 
2702
    return 1;
 
2703
  }
 
2704
 
 
2705
  /// Create a new finite element for sub element i (for a mixed element)
 
2706
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
 
2707
  {
 
2708
    return new UFC_Poisson3D_1LinearForm_finite_element_1();
 
2709
  }
 
2710
 
 
2711
};
 
2712
 
 
2713
/// This class defines the interface for a local-to-global mapping of
 
2714
/// degrees of freedom (dofs).
 
2715
 
 
2716
class UFC_Poisson3D_1LinearForm_dof_map_0: public ufc::dof_map
 
2717
{
 
2718
private:
 
2719
 
 
2720
  unsigned int __global_dimension;
 
2721
 
 
2722
public:
 
2723
 
 
2724
  /// Constructor
 
2725
  UFC_Poisson3D_1LinearForm_dof_map_0() : ufc::dof_map()
 
2726
  {
 
2727
    __global_dimension = 0;
 
2728
  }
 
2729
 
 
2730
  /// Destructor
 
2731
  virtual ~UFC_Poisson3D_1LinearForm_dof_map_0()
 
2732
  {
 
2733
    // Do nothing
 
2734
  }
 
2735
 
 
2736
  /// Return a string identifying the dof map
 
2737
  virtual const char* signature() const
 
2738
  {
 
2739
    return "FFC dof map for Lagrange finite element of degree 1 on a tetrahedron";
 
2740
  }
 
2741
 
 
2742
  /// Return true iff mesh entities of topological dimension d are needed
 
2743
  virtual bool needs_mesh_entities(unsigned int d) const
 
2744
  {
 
2745
    switch ( d )
 
2746
    {
 
2747
    case 0:
 
2748
      return true;
 
2749
      break;
 
2750
    case 1:
 
2751
      return false;
 
2752
      break;
 
2753
    case 2:
 
2754
      return false;
 
2755
      break;
 
2756
    case 3:
 
2757
      return false;
 
2758
      break;
 
2759
    }
 
2760
    return false;
 
2761
  }
 
2762
 
 
2763
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
2764
  virtual bool init_mesh(const ufc::mesh& m)
 
2765
  {
 
2766
    __global_dimension = m.num_entities[0];
 
2767
    return false;
 
2768
  }
 
2769
 
 
2770
  /// Initialize dof map for given cell
 
2771
  virtual void init_cell(const ufc::mesh& m,
 
2772
                         const ufc::cell& c)
 
2773
  {
 
2774
    // Do nothing
 
2775
  }
 
2776
 
 
2777
  /// Finish initialization of dof map for cells
 
2778
  virtual void init_cell_finalize()
 
2779
  {
 
2780
    // Do nothing
 
2781
  }
 
2782
 
 
2783
  /// Return the dimension of the global finite element function space
 
2784
  virtual unsigned int global_dimension() const
 
2785
  {
 
2786
    return __global_dimension;
 
2787
  }
 
2788
 
 
2789
  /// Return the dimension of the local finite element function space
 
2790
  virtual unsigned int local_dimension() const
 
2791
  {
 
2792
    return 4;
 
2793
  }
 
2794
 
 
2795
  // Return the geometric dimension of the coordinates this dof map provides
 
2796
  virtual unsigned int geometric_dimension() const
 
2797
  {
 
2798
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2799
  }
 
2800
 
 
2801
  /// Return the number of dofs on each cell facet
 
2802
  virtual unsigned int num_facet_dofs() const
 
2803
  {
 
2804
    return 3;
 
2805
  }
 
2806
 
 
2807
  /// Return the number of dofs associated with each cell entity of dimension d
 
2808
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
2809
  {
 
2810
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2811
  }
 
2812
 
 
2813
  /// Tabulate the local-to-global mapping of dofs on a cell
 
2814
  virtual void tabulate_dofs(unsigned int* dofs,
 
2815
                             const ufc::mesh& m,
 
2816
                             const ufc::cell& c) const
 
2817
  {
 
2818
    dofs[0] = c.entity_indices[0][0];
 
2819
    dofs[1] = c.entity_indices[0][1];
 
2820
    dofs[2] = c.entity_indices[0][2];
 
2821
    dofs[3] = c.entity_indices[0][3];
 
2822
  }
 
2823
 
 
2824
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
2825
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
2826
                                   unsigned int facet) const
 
2827
  {
 
2828
    switch ( facet )
 
2829
    {
 
2830
    case 0:
 
2831
      dofs[0] = 1;
 
2832
      dofs[1] = 2;
 
2833
      dofs[2] = 3;
 
2834
      break;
 
2835
    case 1:
 
2836
      dofs[0] = 0;
 
2837
      dofs[1] = 2;
 
2838
      dofs[2] = 3;
 
2839
      break;
 
2840
    case 2:
 
2841
      dofs[0] = 0;
 
2842
      dofs[1] = 1;
 
2843
      dofs[2] = 3;
 
2844
      break;
 
2845
    case 3:
 
2846
      dofs[0] = 0;
 
2847
      dofs[1] = 1;
 
2848
      dofs[2] = 2;
 
2849
      break;
 
2850
    }
 
2851
  }
 
2852
 
 
2853
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
2854
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
2855
                                    unsigned int d, unsigned int i) const
 
2856
  {
 
2857
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2858
  }
 
2859
 
 
2860
  /// Tabulate the coordinates of all dofs on a cell
 
2861
  virtual void tabulate_coordinates(double** coordinates,
 
2862
                                    const ufc::cell& c) const
 
2863
  {
 
2864
    const double * const * x = c.coordinates;
 
2865
    coordinates[0][0] = x[0][0];
 
2866
    coordinates[0][1] = x[0][1];
 
2867
    coordinates[0][2] = x[0][2];
 
2868
    coordinates[1][0] = x[1][0];
 
2869
    coordinates[1][1] = x[1][1];
 
2870
    coordinates[1][2] = x[1][2];
 
2871
    coordinates[2][0] = x[2][0];
 
2872
    coordinates[2][1] = x[2][1];
 
2873
    coordinates[2][2] = x[2][2];
 
2874
    coordinates[3][0] = x[3][0];
 
2875
    coordinates[3][1] = x[3][1];
 
2876
    coordinates[3][2] = x[3][2];
 
2877
  }
 
2878
 
 
2879
  /// Return the number of sub dof maps (for a mixed element)
 
2880
  virtual unsigned int num_sub_dof_maps() const
 
2881
  {
 
2882
    return 1;
 
2883
  }
 
2884
 
 
2885
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
2886
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
2887
  {
 
2888
    return new UFC_Poisson3D_1LinearForm_dof_map_0();
 
2889
  }
 
2890
 
 
2891
};
 
2892
 
 
2893
/// This class defines the interface for a local-to-global mapping of
 
2894
/// degrees of freedom (dofs).
 
2895
 
 
2896
class UFC_Poisson3D_1LinearForm_dof_map_1: public ufc::dof_map
 
2897
{
 
2898
private:
 
2899
 
 
2900
  unsigned int __global_dimension;
 
2901
 
 
2902
public:
 
2903
 
 
2904
  /// Constructor
 
2905
  UFC_Poisson3D_1LinearForm_dof_map_1() : ufc::dof_map()
 
2906
  {
 
2907
    __global_dimension = 0;
 
2908
  }
 
2909
 
 
2910
  /// Destructor
 
2911
  virtual ~UFC_Poisson3D_1LinearForm_dof_map_1()
 
2912
  {
 
2913
    // Do nothing
 
2914
  }
 
2915
 
 
2916
  /// Return a string identifying the dof map
 
2917
  virtual const char* signature() const
 
2918
  {
 
2919
    return "FFC dof map for Lagrange finite element of degree 1 on a tetrahedron";
 
2920
  }
 
2921
 
 
2922
  /// Return true iff mesh entities of topological dimension d are needed
 
2923
  virtual bool needs_mesh_entities(unsigned int d) const
 
2924
  {
 
2925
    switch ( d )
 
2926
    {
 
2927
    case 0:
 
2928
      return true;
 
2929
      break;
 
2930
    case 1:
 
2931
      return false;
 
2932
      break;
 
2933
    case 2:
 
2934
      return false;
 
2935
      break;
 
2936
    case 3:
 
2937
      return false;
 
2938
      break;
 
2939
    }
 
2940
    return false;
 
2941
  }
 
2942
 
 
2943
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
 
2944
  virtual bool init_mesh(const ufc::mesh& m)
 
2945
  {
 
2946
    __global_dimension = m.num_entities[0];
 
2947
    return false;
 
2948
  }
 
2949
 
 
2950
  /// Initialize dof map for given cell
 
2951
  virtual void init_cell(const ufc::mesh& m,
 
2952
                         const ufc::cell& c)
 
2953
  {
 
2954
    // Do nothing
 
2955
  }
 
2956
 
 
2957
  /// Finish initialization of dof map for cells
 
2958
  virtual void init_cell_finalize()
 
2959
  {
 
2960
    // Do nothing
 
2961
  }
 
2962
 
 
2963
  /// Return the dimension of the global finite element function space
 
2964
  virtual unsigned int global_dimension() const
 
2965
  {
 
2966
    return __global_dimension;
 
2967
  }
 
2968
 
 
2969
  /// Return the dimension of the local finite element function space
 
2970
  virtual unsigned int local_dimension() const
 
2971
  {
 
2972
    return 4;
 
2973
  }
 
2974
 
 
2975
  // Return the geometric dimension of the coordinates this dof map provides
 
2976
  virtual unsigned int geometric_dimension() const
 
2977
  {
 
2978
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2979
  }
 
2980
 
 
2981
  /// Return the number of dofs on each cell facet
 
2982
  virtual unsigned int num_facet_dofs() const
 
2983
  {
 
2984
    return 3;
 
2985
  }
 
2986
 
 
2987
  /// Return the number of dofs associated with each cell entity of dimension d
 
2988
  virtual unsigned int num_entity_dofs(unsigned int d) const
 
2989
  {
 
2990
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
2991
  }
 
2992
 
 
2993
  /// Tabulate the local-to-global mapping of dofs on a cell
 
2994
  virtual void tabulate_dofs(unsigned int* dofs,
 
2995
                             const ufc::mesh& m,
 
2996
                             const ufc::cell& c) const
 
2997
  {
 
2998
    dofs[0] = c.entity_indices[0][0];
 
2999
    dofs[1] = c.entity_indices[0][1];
 
3000
    dofs[2] = c.entity_indices[0][2];
 
3001
    dofs[3] = c.entity_indices[0][3];
 
3002
  }
 
3003
 
 
3004
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
 
3005
  virtual void tabulate_facet_dofs(unsigned int* dofs,
 
3006
                                   unsigned int facet) const
 
3007
  {
 
3008
    switch ( facet )
 
3009
    {
 
3010
    case 0:
 
3011
      dofs[0] = 1;
 
3012
      dofs[1] = 2;
 
3013
      dofs[2] = 3;
 
3014
      break;
 
3015
    case 1:
 
3016
      dofs[0] = 0;
 
3017
      dofs[1] = 2;
 
3018
      dofs[2] = 3;
 
3019
      break;
 
3020
    case 2:
 
3021
      dofs[0] = 0;
 
3022
      dofs[1] = 1;
 
3023
      dofs[2] = 3;
 
3024
      break;
 
3025
    case 3:
 
3026
      dofs[0] = 0;
 
3027
      dofs[1] = 1;
 
3028
      dofs[2] = 2;
 
3029
      break;
 
3030
    }
 
3031
  }
 
3032
 
 
3033
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
 
3034
  virtual void tabulate_entity_dofs(unsigned int* dofs,
 
3035
                                    unsigned int d, unsigned int i) const
 
3036
  {
 
3037
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
 
3038
  }
 
3039
 
 
3040
  /// Tabulate the coordinates of all dofs on a cell
 
3041
  virtual void tabulate_coordinates(double** coordinates,
 
3042
                                    const ufc::cell& c) const
 
3043
  {
 
3044
    const double * const * x = c.coordinates;
 
3045
    coordinates[0][0] = x[0][0];
 
3046
    coordinates[0][1] = x[0][1];
 
3047
    coordinates[0][2] = x[0][2];
 
3048
    coordinates[1][0] = x[1][0];
 
3049
    coordinates[1][1] = x[1][1];
 
3050
    coordinates[1][2] = x[1][2];
 
3051
    coordinates[2][0] = x[2][0];
 
3052
    coordinates[2][1] = x[2][1];
 
3053
    coordinates[2][2] = x[2][2];
 
3054
    coordinates[3][0] = x[3][0];
 
3055
    coordinates[3][1] = x[3][1];
 
3056
    coordinates[3][2] = x[3][2];
 
3057
  }
 
3058
 
 
3059
  /// Return the number of sub dof maps (for a mixed element)
 
3060
  virtual unsigned int num_sub_dof_maps() const
 
3061
  {
 
3062
    return 1;
 
3063
  }
 
3064
 
 
3065
  /// Create a new dof_map for sub dof map i (for a mixed element)
 
3066
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
 
3067
  {
 
3068
    return new UFC_Poisson3D_1LinearForm_dof_map_1();
 
3069
  }
 
3070
 
 
3071
};
 
3072
 
 
3073
/// This class defines the interface for the tabulation of the cell
 
3074
/// tensor corresponding to the local contribution to a form from
 
3075
/// the integral over a cell.
 
3076
 
 
3077
class UFC_Poisson3D_1LinearForm_cell_integral_0: public ufc::cell_integral
 
3078
{
 
3079
public:
 
3080
 
 
3081
  /// Constructor
 
3082
  UFC_Poisson3D_1LinearForm_cell_integral_0() : ufc::cell_integral()
 
3083
  {
 
3084
    // Do nothing
 
3085
  }
 
3086
 
 
3087
  /// Destructor
 
3088
  virtual ~UFC_Poisson3D_1LinearForm_cell_integral_0()
 
3089
  {
 
3090
    // Do nothing
 
3091
  }
 
3092
 
 
3093
  /// Tabulate the tensor for the contribution from a local cell
 
3094
  virtual void tabulate_tensor(double* A,
 
3095
                               const double * const * w,
 
3096
                               const ufc::cell& c) const
 
3097
  {
 
3098
    // Extract vertex coordinates
 
3099
    const double * const * x = c.coordinates;
 
3100
    
 
3101
    // Compute Jacobian of affine map from reference cell
 
3102
    const double J_00 = x[1][0] - x[0][0];
 
3103
    const double J_01 = x[2][0] - x[0][0];
 
3104
    const double J_02 = x[3][0] - x[0][0];
 
3105
    const double J_10 = x[1][1] - x[0][1];
 
3106
    const double J_11 = x[2][1] - x[0][1];
 
3107
    const double J_12 = x[3][1] - x[0][1];
 
3108
    const double J_20 = x[1][2] - x[0][2];
 
3109
    const double J_21 = x[2][2] - x[0][2];
 
3110
    const double J_22 = x[3][2] - x[0][2];
 
3111
      
 
3112
    // Compute sub determinants
 
3113
    const double d_00 = J_11*J_22 - J_12*J_21;
 
3114
    
 
3115
    const double d_10 = J_02*J_21 - J_01*J_22;
 
3116
    
 
3117
    const double d_20 = J_01*J_12 - J_02*J_11;
 
3118
      
 
3119
    // Compute determinant of Jacobian
 
3120
    double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
 
3121
      
 
3122
    // Compute inverse of Jacobian
 
3123
    
 
3124
    // Set scale factor
 
3125
    const double det = std::abs(detJ);
 
3126
    
 
3127
    // Compute coefficients
 
3128
    const double c0_0_0_0 = w[0][0];
 
3129
    const double c0_0_0_1 = w[0][1];
 
3130
    const double c0_0_0_2 = w[0][2];
 
3131
    const double c0_0_0_3 = w[0][3];
 
3132
    
 
3133
    // Compute geometry tensors
 
3134
    const double G0_0 = det*c0_0_0_0;
 
3135
    const double G0_1 = det*c0_0_0_1;
 
3136
    const double G0_2 = det*c0_0_0_2;
 
3137
    const double G0_3 = det*c0_0_0_3;
 
3138
    
 
3139
    // Compute element tensor
 
3140
    A[0] = 0.0166666666666666*G0_0 + 0.00833333333333331*G0_1 + 0.00833333333333331*G0_2 + 0.00833333333333331*G0_3;
 
3141
    A[1] = 0.00833333333333331*G0_0 + 0.0166666666666666*G0_1 + 0.00833333333333331*G0_2 + 0.00833333333333331*G0_3;
 
3142
    A[2] = 0.00833333333333331*G0_0 + 0.00833333333333331*G0_1 + 0.0166666666666666*G0_2 + 0.00833333333333331*G0_3;
 
3143
    A[3] = 0.00833333333333331*G0_0 + 0.00833333333333331*G0_1 + 0.00833333333333331*G0_2 + 0.0166666666666666*G0_3;
 
3144
  }
 
3145
 
 
3146
};
 
3147
 
 
3148
/// This class defines the interface for the assembly of the global
 
3149
/// tensor corresponding to a form with r + n arguments, that is, a
 
3150
/// mapping
 
3151
///
 
3152
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
 
3153
///
 
3154
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
 
3155
/// global tensor A is defined by
 
3156
///
 
3157
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
 
3158
///
 
3159
/// where each argument Vj represents the application to the
 
3160
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
 
3161
/// fixed functions (coefficients).
 
3162
 
 
3163
class UFC_Poisson3D_1LinearForm: public ufc::form
 
3164
{
 
3165
public:
 
3166
 
 
3167
  /// Constructor
 
3168
  UFC_Poisson3D_1LinearForm() : ufc::form()
 
3169
  {
 
3170
    // Do nothing
 
3171
  }
 
3172
 
 
3173
  /// Destructor
 
3174
  virtual ~UFC_Poisson3D_1LinearForm()
 
3175
  {
 
3176
    // Do nothing
 
3177
  }
 
3178
 
 
3179
  /// Return a string identifying the form
 
3180
  virtual const char* signature() const
 
3181
  {
 
3182
    return "w0_a0 | vi0*va0*dX(0)";
 
3183
  }
 
3184
 
 
3185
  /// Return the rank of the global tensor (r)
 
3186
  virtual unsigned int rank() const
 
3187
  {
 
3188
    return 1;
 
3189
  }
 
3190
 
 
3191
  /// Return the number of coefficients (n)
 
3192
  virtual unsigned int num_coefficients() const
 
3193
  {
 
3194
    return 1;
 
3195
  }
 
3196
 
 
3197
  /// Return the number of cell integrals
 
3198
  virtual unsigned int num_cell_integrals() const
 
3199
  {
 
3200
    return 1;
 
3201
  }
 
3202
  
 
3203
  /// Return the number of exterior facet integrals
 
3204
  virtual unsigned int num_exterior_facet_integrals() const
 
3205
  {
 
3206
    return 0;
 
3207
  }
 
3208
  
 
3209
  /// Return the number of interior facet integrals
 
3210
  virtual unsigned int num_interior_facet_integrals() const
 
3211
  {
 
3212
    return 0;
 
3213
  }
 
3214
    
 
3215
  /// Create a new finite element for argument function i
 
3216
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
 
3217
  {
 
3218
    switch ( i )
 
3219
    {
 
3220
    case 0:
 
3221
      return new UFC_Poisson3D_1LinearForm_finite_element_0();
 
3222
      break;
 
3223
    case 1:
 
3224
      return new UFC_Poisson3D_1LinearForm_finite_element_1();
 
3225
      break;
 
3226
    }
 
3227
    return 0;
 
3228
  }
 
3229
  
 
3230
  /// Create a new dof map for argument function i
 
3231
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
 
3232
  {
 
3233
    switch ( i )
 
3234
    {
 
3235
    case 0:
 
3236
      return new UFC_Poisson3D_1LinearForm_dof_map_0();
 
3237
      break;
 
3238
    case 1:
 
3239
      return new UFC_Poisson3D_1LinearForm_dof_map_1();
 
3240
      break;
 
3241
    }
 
3242
    return 0;
 
3243
  }
 
3244
 
 
3245
  /// Create a new cell integral on sub domain i
 
3246
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
 
3247
  {
 
3248
    return new UFC_Poisson3D_1LinearForm_cell_integral_0();
 
3249
  }
 
3250
 
 
3251
  /// Create a new exterior facet integral on sub domain i
 
3252
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
 
3253
  {
 
3254
    return 0;
 
3255
  }
 
3256
 
 
3257
  /// Create a new interior facet integral on sub domain i
 
3258
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
 
3259
  {
 
3260
    return 0;
 
3261
  }
 
3262
 
 
3263
};
 
3264
 
 
3265
// DOLFIN wrappers
 
3266
 
 
3267
#include <dolfin/Form.h>
 
3268
 
 
3269
class Poisson3D_1BilinearForm : public dolfin::Form
 
3270
{
 
3271
public:
 
3272
 
 
3273
  Poisson3D_1BilinearForm() : dolfin::Form()
 
3274
  {
 
3275
    // Do nothing
 
3276
  }
 
3277
 
 
3278
  /// Return UFC form
 
3279
  virtual const ufc::form& form() const
 
3280
  {
 
3281
    return __form;
 
3282
  }
 
3283
  
 
3284
  /// Return array of coefficients
 
3285
  virtual const dolfin::Array<dolfin::Function*>& coefficients() const
 
3286
  {
 
3287
    return __coefficients;
 
3288
  }
 
3289
 
 
3290
private:
 
3291
 
 
3292
  // UFC form
 
3293
  UFC_Poisson3D_1BilinearForm __form;
 
3294
 
 
3295
  /// Array of coefficients
 
3296
  dolfin::Array<dolfin::Function*> __coefficients;
 
3297
 
 
3298
};
 
3299
 
 
3300
class Poisson3D_1LinearForm : public dolfin::Form
 
3301
{
 
3302
public:
 
3303
 
 
3304
  Poisson3D_1LinearForm(dolfin::Function& w0) : dolfin::Form()
 
3305
  {
 
3306
    __coefficients.push_back(&w0);
 
3307
  }
 
3308
 
 
3309
  /// Return UFC form
 
3310
  virtual const ufc::form& form() const
 
3311
  {
 
3312
    return __form;
 
3313
  }
 
3314
  
 
3315
  /// Return array of coefficients
 
3316
  virtual const dolfin::Array<dolfin::Function*>& coefficients() const
 
3317
  {
 
3318
    return __coefficients;
 
3319
  }
 
3320
 
 
3321
private:
 
3322
 
 
3323
  // UFC form
 
3324
  UFC_Poisson3D_1LinearForm __form;
 
3325
 
 
3326
  /// Array of coefficients
 
3327
  dolfin::Array<dolfin::Function*> __coefficients;
 
3328
 
 
3329
};
 
3330
 
 
3331
#endif