~ubuntu-branches/ubuntu/trusty/ffc/trusty
» Revision
4
: test/regression/reference/quadrature/Biharmonic.h
« back to all changes in this revision
Viewing changes to test/regression/reference/quadrature/Biharmonic.h
- browse files at revision 4
- compare with another revision
- download diff
- download tarball
- view history from revision 4
- Committer: Bazaar Package Importer
- Author(s): Johannes Ring
- Date: 2010-02-03 20:22:35 UTC
- mfrom: (1.1.2 upstream)
- Revision ID: james.westby@ubuntu.com-20100203202235-fe8d0kajuvgy2sqn
Tags: 0.9.0-1
* New upstream release.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.
- files added:
- demo/VectorPoisson.ufl
- ffc/quadrature
- ffc/tensor
- test/evaluate_basis
- test/evaluate_basis_derivatives
- test/regression/references
- test/unit/evaluate_basis
- test/unit/misc
- test/unit/symbolics
- files removed:
- ffc/common
- ffc/compiler
- ffc/compiler/quadrature
- ffc/compiler/tensor
- ffc/fem
- ffc/jit
- test/benchmarks
- test/regression/reference
- test/regression/reference/quadrature
- test/regression/reference/tensor
- test/simple_verify_tensors
- test/verify_tensors
- test/verify_tensors/references
- test/verify_tensors/test_cache
- files modified:
- AUTHORS
added
removed
test/regression/reference/quadrature/Biharmonic.h
1
// This code conforms with the UFC specification version 1.0
2
// and was automatically generated by FFC version 0.7.0.
3
4
#ifndef __BIHARMONIC_H
5
#define __BIHARMONIC_H
6
7
#include <cmath>
8
#include <stdexcept>
9
#include <ufc.h>
10
11
/// This class defines the interface for a finite element.
12
13
class biharmonic_0_finite_element_0: public ufc::finite_element
14
{
15
public:
16
17
/// Constructor
18
biharmonic_0_finite_element_0() : ufc::finite_element()
19
{
20
// Do nothing
21
}
22
23
/// Destructor
24
virtual ~biharmonic_0_finite_element_0()
25
{
26
// Do nothing
27
}
28
29
/// Return a string identifying the finite element
30
virtual const char* signature() const
31
{
32
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
33
}
34
35
/// Return the cell shape
36
virtual ufc::shape cell_shape() const
37
{
38
return ufc::triangle;
39
}
40
41
/// Return the dimension of the finite element function space
42
virtual unsigned int space_dimension() const
43
{
44
return 6;
45
}
46
47
/// Return the rank of the value space
48
virtual unsigned int value_rank() const
49
{
50
return 0;
51
}
52
53
/// Return the dimension of the value space for axis i
54
virtual unsigned int value_dimension(unsigned int i) const
55
{
56
return 1;
57
}
58
59
/// Evaluate basis function i at given point in cell
60
virtual void evaluate_basis(unsigned int i,
61
double* values,
62
const double* coordinates,
63
const ufc::cell& c) const
64
{
65
// Extract vertex coordinates
66
const double * const * element_coordinates = c.coordinates;
67
68
// Compute Jacobian of affine map from reference cell
69
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
70
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
71
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
72
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
73
74
// Compute determinant of Jacobian
75
const double detJ = J_00*J_11 - J_01*J_10;
76
77
// Compute inverse of Jacobian
78
79
// Get coordinates and map to the reference (UFC) element
80
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
81
element_coordinates[0][0]*element_coordinates[2][1] +\
82
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
83
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
84
element_coordinates[1][0]*element_coordinates[0][1] -\
85
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
86
87
// Map coordinates to the reference square
88
if (std::abs(y - 1.0) < 1e-08)
89
x = -1.0;
90
else
91
x = 2.0 *x/(1.0 - y) - 1.0;
92
y = 2.0*y - 1.0;
93
94
// Reset values
95
*values = 0;
96
97
// Map degree of freedom to element degree of freedom
98
const unsigned int dof = i;
99
100
// Generate scalings
101
const double scalings_y_0 = 1;
102
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
103
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
104
105
// Compute psitilde_a
106
const double psitilde_a_0 = 1;
107
const double psitilde_a_1 = x;
108
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
109
110
// Compute psitilde_bs
111
const double psitilde_bs_0_0 = 1;
112
const double psitilde_bs_0_1 = 1.5*y + 0.5;
113
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
114
const double psitilde_bs_1_0 = 1;
115
const double psitilde_bs_1_1 = 2.5*y + 1.5;
116
const double psitilde_bs_2_0 = 1;
117
118
// Compute basisvalues
119
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
120
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
121
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
122
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
123
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
124
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
125
126
// Table(s) of coefficients
127
static const double coefficients0[6][6] = \
128
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
129
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
130
{0, 0, 0.2, 0, 0, 0.163299316},
131
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
132
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
133
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
134
135
// Extract relevant coefficients
136
const double coeff0_0 = coefficients0[dof][0];
137
const double coeff0_1 = coefficients0[dof][1];
138
const double coeff0_2 = coefficients0[dof][2];
139
const double coeff0_3 = coefficients0[dof][3];
140
const double coeff0_4 = coefficients0[dof][4];
141
const double coeff0_5 = coefficients0[dof][5];
142
143
// Compute value(s)
144
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
145
}
146
147
/// Evaluate all basis functions at given point in cell
148
virtual void evaluate_basis_all(double* values,
149
const double* coordinates,
150
const ufc::cell& c) const
151
{
152
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
153
}
154
155
/// Evaluate order n derivatives of basis function i at given point in cell
156
virtual void evaluate_basis_derivatives(unsigned int i,
157
unsigned int n,
158
double* values,
159
const double* coordinates,
160
const ufc::cell& c) const
161
{
162
// Extract vertex coordinates
163
const double * const * element_coordinates = c.coordinates;
164
165
// Compute Jacobian of affine map from reference cell
166
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
167
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
168
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
169
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
170
171
// Compute determinant of Jacobian
172
const double detJ = J_00*J_11 - J_01*J_10;
173
174
// Compute inverse of Jacobian
175
176
// Get coordinates and map to the reference (UFC) element
177
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
178
element_coordinates[0][0]*element_coordinates[2][1] +\
179
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
180
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
181
element_coordinates[1][0]*element_coordinates[0][1] -\
182
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
183
184
// Map coordinates to the reference square
185
if (std::abs(y - 1.0) < 1e-08)
186
x = -1.0;
187
else
188
x = 2.0 *x/(1.0 - y) - 1.0;
189
y = 2.0*y - 1.0;
190
191
// Compute number of derivatives
192
unsigned int num_derivatives = 1;
193
194
for (unsigned int j = 0; j < n; j++)
195
num_derivatives *= 2;
196
197
198
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
199
unsigned int **combinations = new unsigned int *[num_derivatives];
200
201
for (unsigned int j = 0; j < num_derivatives; j++)
202
{
203
combinations[j] = new unsigned int [n];
204
for (unsigned int k = 0; k < n; k++)
205
combinations[j][k] = 0;
206
}
207
208
// Generate combinations of derivatives
209
for (unsigned int row = 1; row < num_derivatives; row++)
210
{
211
for (unsigned int num = 0; num < row; num++)
212
{
213
for (unsigned int col = n-1; col+1 > 0; col--)
214
{
215
if (combinations[row][col] + 1 > 1)
216
combinations[row][col] = 0;
217
else
218
{
219
combinations[row][col] += 1;
220
break;
221
}
222
}
223
}
224
}
225
226
// Compute inverse of Jacobian
227
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
228
229
// Declare transformation matrix
230
// Declare pointer to two dimensional array and initialise
231
double **transform = new double *[num_derivatives];
232
233
for (unsigned int j = 0; j < num_derivatives; j++)
234
{
235
transform[j] = new double [num_derivatives];
236
for (unsigned int k = 0; k < num_derivatives; k++)
237
transform[j][k] = 1;
238
}
239
240
// Construct transformation matrix
241
for (unsigned int row = 0; row < num_derivatives; row++)
242
{
243
for (unsigned int col = 0; col < num_derivatives; col++)
244
{
245
for (unsigned int k = 0; k < n; k++)
246
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
247
}
248
}
249
250
// Reset values
251
for (unsigned int j = 0; j < 1*num_derivatives; j++)
252
values[j] = 0;
253
254
// Map degree of freedom to element degree of freedom
255
const unsigned int dof = i;
256
257
// Generate scalings
258
const double scalings_y_0 = 1;
259
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
260
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
261
262
// Compute psitilde_a
263
const double psitilde_a_0 = 1;
264
const double psitilde_a_1 = x;
265
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
266
267
// Compute psitilde_bs
268
const double psitilde_bs_0_0 = 1;
269
const double psitilde_bs_0_1 = 1.5*y + 0.5;
270
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
271
const double psitilde_bs_1_0 = 1;
272
const double psitilde_bs_1_1 = 2.5*y + 1.5;
273
const double psitilde_bs_2_0 = 1;
274
275
// Compute basisvalues
276
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
277
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
278
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
279
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
280
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
281
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
282
283
// Table(s) of coefficients
284
static const double coefficients0[6][6] = \
285
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
286
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
287
{0, 0, 0.2, 0, 0, 0.163299316},
288
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
289
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
290
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
291
292
// Interesting (new) part
293
// Tables of derivatives of the polynomial base (transpose)
294
static const double dmats0[6][6] = \
295
{{0, 0, 0, 0, 0, 0},
296
{4.89897949, 0, 0, 0, 0, 0},
297
{0, 0, 0, 0, 0, 0},
298
{0, 9.48683298, 0, 0, 0, 0},
299
{4, 0, 7.07106781, 0, 0, 0},
300
{0, 0, 0, 0, 0, 0}};
301
302
static const double dmats1[6][6] = \
303
{{0, 0, 0, 0, 0, 0},
304
{2.44948974, 0, 0, 0, 0, 0},
305
{4.24264069, 0, 0, 0, 0, 0},
306
{2.5819889, 4.74341649, -0.912870929, 0, 0, 0},
307
{2, 6.12372436, 3.53553391, 0, 0, 0},
308
{-2.30940108, 0, 8.16496581, 0, 0, 0}};
309
310
// Compute reference derivatives
311
// Declare pointer to array of derivatives on FIAT element
312
double *derivatives = new double [num_derivatives];
313
314
// Declare coefficients
315
double coeff0_0 = 0;
316
double coeff0_1 = 0;
317
double coeff0_2 = 0;
318
double coeff0_3 = 0;
319
double coeff0_4 = 0;
320
double coeff0_5 = 0;
321
322
// Declare new coefficients
323
double new_coeff0_0 = 0;
324
double new_coeff0_1 = 0;
325
double new_coeff0_2 = 0;
326
double new_coeff0_3 = 0;
327
double new_coeff0_4 = 0;
328
double new_coeff0_5 = 0;
329
330
// Loop possible derivatives
331
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
332
{
333
// Get values from coefficients array
334
new_coeff0_0 = coefficients0[dof][0];
335
new_coeff0_1 = coefficients0[dof][1];
336
new_coeff0_2 = coefficients0[dof][2];
337
new_coeff0_3 = coefficients0[dof][3];
338
new_coeff0_4 = coefficients0[dof][4];
339
new_coeff0_5 = coefficients0[dof][5];
340
341
// Loop derivative order
342
for (unsigned int j = 0; j < n; j++)
343
{
344
// Update old coefficients
345
coeff0_0 = new_coeff0_0;
346
coeff0_1 = new_coeff0_1;
347
coeff0_2 = new_coeff0_2;
348
coeff0_3 = new_coeff0_3;
349
coeff0_4 = new_coeff0_4;
350
coeff0_5 = new_coeff0_5;
351
352
if(combinations[deriv_num][j] == 0)
353
{
354
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
355
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
356
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
357
new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
358
new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
359
new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
360
}
361
if(combinations[deriv_num][j] == 1)
362
{
363
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
364
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
365
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
366
new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
367
new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
368
new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
369
}
370
371
}
372
// Compute derivatives on reference element as dot product of coefficients and basisvalues
373
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
374
}
375
376
// Transform derivatives back to physical element
377
for (unsigned int row = 0; row < num_derivatives; row++)
378
{
379
for (unsigned int col = 0; col < num_derivatives; col++)
380
{
381
values[row] += transform[row][col]*derivatives[col];
382
}
383
}
384
// Delete pointer to array of derivatives on FIAT element
385
delete [] derivatives;
386
387
// Delete pointer to array of combinations of derivatives and transform
388
for (unsigned int row = 0; row < num_derivatives; row++)
389
{
390
delete [] combinations[row];
391
delete [] transform[row];
392
}
393
394
delete [] combinations;
395
delete [] transform;
396
}
397
398
/// Evaluate order n derivatives of all basis functions at given point in cell
399
virtual void evaluate_basis_derivatives_all(unsigned int n,
400
double* values,
401
const double* coordinates,
402
const ufc::cell& c) const
403
{
404
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
405
}
406
407
/// Evaluate linear functional for dof i on the function f
408
virtual double evaluate_dof(unsigned int i,
409
const ufc::function& f,
410
const ufc::cell& c) const
411
{
412
// The reference points, direction and weights:
413
static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
414
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
415
static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
416
417
const double * const * x = c.coordinates;
418
double result = 0.0;
419
// Iterate over the points:
420
// Evaluate basis functions for affine mapping
421
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
422
const double w1 = X[i][0][0];
423
const double w2 = X[i][0][1];
424
425
// Compute affine mapping y = F(X)
426
double y[2];
427
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
428
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
429
430
// Evaluate function at physical points
431
double values[1];
432
f.evaluate(values, y, c);
433
434
// Map function values using appropriate mapping
435
// Affine map: Do nothing
436
437
// Note that we do not map the weights (yet).
438
439
// Take directional components
440
for(int k = 0; k < 1; k++)
441
result += values[k]*D[i][0][k];
442
// Multiply by weights
443
result *= W[i][0];
444
445
return result;
446
}
447
448
/// Evaluate linear functionals for all dofs on the function f
449
virtual void evaluate_dofs(double* values,
450
const ufc::function& f,
451
const ufc::cell& c) const
452
{
453
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
454
}
455
456
/// Interpolate vertex values from dof values
457
virtual void interpolate_vertex_values(double* vertex_values,
458
const double* dof_values,
459
const ufc::cell& c) const
460
{
461
// Evaluate at vertices and use affine mapping
462
vertex_values[0] = dof_values[0];
463
vertex_values[1] = dof_values[1];
464
vertex_values[2] = dof_values[2];
465
}
466
467
/// Return the number of sub elements (for a mixed element)
468
virtual unsigned int num_sub_elements() const
469
{
470
return 1;
471
}
472
473
/// Create a new finite element for sub element i (for a mixed element)
474
virtual ufc::finite_element* create_sub_element(unsigned int i) const
475
{
476
return new biharmonic_0_finite_element_0();
477
}
478
479
};
480
481
/// This class defines the interface for a finite element.
482
483
class biharmonic_0_finite_element_1: public ufc::finite_element
484
{
485
public:
486
487
/// Constructor
488
biharmonic_0_finite_element_1() : ufc::finite_element()
489
{
490
// Do nothing
491
}
492
493
/// Destructor
494
virtual ~biharmonic_0_finite_element_1()
495
{
496
// Do nothing
497
}
498
499
/// Return a string identifying the finite element
500
virtual const char* signature() const
501
{
502
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
503
}
504
505
/// Return the cell shape
506
virtual ufc::shape cell_shape() const
507
{
508
return ufc::triangle;
509
}
510
511
/// Return the dimension of the finite element function space
512
virtual unsigned int space_dimension() const
513
{
514
return 6;
515
}
516
517
/// Return the rank of the value space
518
virtual unsigned int value_rank() const
519
{
520
return 0;
521
}
522
523
/// Return the dimension of the value space for axis i
524
virtual unsigned int value_dimension(unsigned int i) const
525
{
526
return 1;
527
}
528
529
/// Evaluate basis function i at given point in cell
530
virtual void evaluate_basis(unsigned int i,
531
double* values,
532
const double* coordinates,
533
const ufc::cell& c) const
534
{
535
// Extract vertex coordinates
536
const double * const * element_coordinates = c.coordinates;
537
538
// Compute Jacobian of affine map from reference cell
539
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
540
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
541
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
542
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
543
544
// Compute determinant of Jacobian
545
const double detJ = J_00*J_11 - J_01*J_10;
546
547
// Compute inverse of Jacobian
548
549
// Get coordinates and map to the reference (UFC) element
550
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
551
element_coordinates[0][0]*element_coordinates[2][1] +\
552
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
553
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
554
element_coordinates[1][0]*element_coordinates[0][1] -\
555
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
556
557
// Map coordinates to the reference square
558
if (std::abs(y - 1.0) < 1e-08)
559
x = -1.0;
560
else
561
x = 2.0 *x/(1.0 - y) - 1.0;
562
y = 2.0*y - 1.0;
563
564
// Reset values
565
*values = 0;
566
567
// Map degree of freedom to element degree of freedom
568
const unsigned int dof = i;
569
570
// Generate scalings
571
const double scalings_y_0 = 1;
572
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
573
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
574
575
// Compute psitilde_a
576
const double psitilde_a_0 = 1;
577
const double psitilde_a_1 = x;
578
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
579
580
// Compute psitilde_bs
581
const double psitilde_bs_0_0 = 1;
582
const double psitilde_bs_0_1 = 1.5*y + 0.5;
583
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
584
const double psitilde_bs_1_0 = 1;
585
const double psitilde_bs_1_1 = 2.5*y + 1.5;
586
const double psitilde_bs_2_0 = 1;
587
588
// Compute basisvalues
589
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
590
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
591
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
592
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
593
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
594
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
595
596
// Table(s) of coefficients
597
static const double coefficients0[6][6] = \
598
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
599
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
600
{0, 0, 0.2, 0, 0, 0.163299316},
601
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
602
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
603
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
604
605
// Extract relevant coefficients
606
const double coeff0_0 = coefficients0[dof][0];
607
const double coeff0_1 = coefficients0[dof][1];
608
const double coeff0_2 = coefficients0[dof][2];
609
const double coeff0_3 = coefficients0[dof][3];
610
const double coeff0_4 = coefficients0[dof][4];
611
const double coeff0_5 = coefficients0[dof][5];
612
613
// Compute value(s)
614
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
615
}
616
617
/// Evaluate all basis functions at given point in cell
618
virtual void evaluate_basis_all(double* values,
619
const double* coordinates,
620
const ufc::cell& c) const
621
{
622
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
623
}
624
625
/// Evaluate order n derivatives of basis function i at given point in cell
626
virtual void evaluate_basis_derivatives(unsigned int i,
627
unsigned int n,
628
double* values,
629
const double* coordinates,
630
const ufc::cell& c) const
631
{
632
// Extract vertex coordinates
633
const double * const * element_coordinates = c.coordinates;
634
635
// Compute Jacobian of affine map from reference cell
636
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
637
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
638
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
639
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
640
641
// Compute determinant of Jacobian
642
const double detJ = J_00*J_11 - J_01*J_10;
643
644
// Compute inverse of Jacobian
645
646
// Get coordinates and map to the reference (UFC) element
647
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
648
element_coordinates[0][0]*element_coordinates[2][1] +\
649
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
650
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
651
element_coordinates[1][0]*element_coordinates[0][1] -\
652
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
653
654
// Map coordinates to the reference square
655
if (std::abs(y - 1.0) < 1e-08)
656
x = -1.0;
657
else
658
x = 2.0 *x/(1.0 - y) - 1.0;
659
y = 2.0*y - 1.0;
660
661
// Compute number of derivatives
662
unsigned int num_derivatives = 1;
663
664
for (unsigned int j = 0; j < n; j++)
665
num_derivatives *= 2;
666
667
668
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
669
unsigned int **combinations = new unsigned int *[num_derivatives];
670
671
for (unsigned int j = 0; j < num_derivatives; j++)
672
{
673
combinations[j] = new unsigned int [n];
674
for (unsigned int k = 0; k < n; k++)
675
combinations[j][k] = 0;
676
}
677
678
// Generate combinations of derivatives
679
for (unsigned int row = 1; row < num_derivatives; row++)
680
{
681
for (unsigned int num = 0; num < row; num++)
682
{
683
for (unsigned int col = n-1; col+1 > 0; col--)
684
{
685
if (combinations[row][col] + 1 > 1)
686
combinations[row][col] = 0;
687
else
688
{
689
combinations[row][col] += 1;
690
break;
691
}
692
}
693
}
694
}
695
696
// Compute inverse of Jacobian
697
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
698
699
// Declare transformation matrix
700
// Declare pointer to two dimensional array and initialise
701
double **transform = new double *[num_derivatives];
702
703
for (unsigned int j = 0; j < num_derivatives; j++)
704
{
705
transform[j] = new double [num_derivatives];
706
for (unsigned int k = 0; k < num_derivatives; k++)
707
transform[j][k] = 1;
708
}
709
710
// Construct transformation matrix
711
for (unsigned int row = 0; row < num_derivatives; row++)
712
{
713
for (unsigned int col = 0; col < num_derivatives; col++)
714
{
715
for (unsigned int k = 0; k < n; k++)
716
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
717
}
718
}
719
720
// Reset values
721
for (unsigned int j = 0; j < 1*num_derivatives; j++)
722
values[j] = 0;
723
724
// Map degree of freedom to element degree of freedom
725
const unsigned int dof = i;
726
727
// Generate scalings
728
const double scalings_y_0 = 1;
729
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
730
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
731
732
// Compute psitilde_a
733
const double psitilde_a_0 = 1;
734
const double psitilde_a_1 = x;
735
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
736
737
// Compute psitilde_bs
738
const double psitilde_bs_0_0 = 1;
739
const double psitilde_bs_0_1 = 1.5*y + 0.5;
740
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
741
const double psitilde_bs_1_0 = 1;
742
const double psitilde_bs_1_1 = 2.5*y + 1.5;
743
const double psitilde_bs_2_0 = 1;
744
745
// Compute basisvalues
746
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
747
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
748
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
749
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
750
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
751
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
752
753
// Table(s) of coefficients
754
static const double coefficients0[6][6] = \
755
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
756
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
757
{0, 0, 0.2, 0, 0, 0.163299316},
758
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
759
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
760
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
761
762
// Interesting (new) part
763
// Tables of derivatives of the polynomial base (transpose)
764
static const double dmats0[6][6] = \
765
{{0, 0, 0, 0, 0, 0},
766
{4.89897949, 0, 0, 0, 0, 0},
767
{0, 0, 0, 0, 0, 0},
768
{0, 9.48683298, 0, 0, 0, 0},
769
{4, 0, 7.07106781, 0, 0, 0},
770
{0, 0, 0, 0, 0, 0}};
771
772
static const double dmats1[6][6] = \
773
{{0, 0, 0, 0, 0, 0},
774
{2.44948974, 0, 0, 0, 0, 0},
775
{4.24264069, 0, 0, 0, 0, 0},
776
{2.5819889, 4.74341649, -0.912870929, 0, 0, 0},
777
{2, 6.12372436, 3.53553391, 0, 0, 0},
778
{-2.30940108, 0, 8.16496581, 0, 0, 0}};
779
780
// Compute reference derivatives
781
// Declare pointer to array of derivatives on FIAT element
782
double *derivatives = new double [num_derivatives];
783
784
// Declare coefficients
785
double coeff0_0 = 0;
786
double coeff0_1 = 0;
787
double coeff0_2 = 0;
788
double coeff0_3 = 0;
789
double coeff0_4 = 0;
790
double coeff0_5 = 0;
791
792
// Declare new coefficients
793
double new_coeff0_0 = 0;
794
double new_coeff0_1 = 0;
795
double new_coeff0_2 = 0;
796
double new_coeff0_3 = 0;
797
double new_coeff0_4 = 0;
798
double new_coeff0_5 = 0;
799
800
// Loop possible derivatives
801
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
802
{
803
// Get values from coefficients array
804
new_coeff0_0 = coefficients0[dof][0];
805
new_coeff0_1 = coefficients0[dof][1];
806
new_coeff0_2 = coefficients0[dof][2];
807
new_coeff0_3 = coefficients0[dof][3];
808
new_coeff0_4 = coefficients0[dof][4];
809
new_coeff0_5 = coefficients0[dof][5];
810
811
// Loop derivative order
812
for (unsigned int j = 0; j < n; j++)
813
{
814
// Update old coefficients
815
coeff0_0 = new_coeff0_0;
816
coeff0_1 = new_coeff0_1;
817
coeff0_2 = new_coeff0_2;
818
coeff0_3 = new_coeff0_3;
819
coeff0_4 = new_coeff0_4;
820
coeff0_5 = new_coeff0_5;
821
822
if(combinations[deriv_num][j] == 0)
823
{
824
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
825
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
826
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
827
new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
828
new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
829
new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
830
}
831
if(combinations[deriv_num][j] == 1)
832
{
833
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
834
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
835
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
836
new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
837
new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
838
new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
839
}
840
841
}
842
// Compute derivatives on reference element as dot product of coefficients and basisvalues
843
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
844
}
845
846
// Transform derivatives back to physical element
847
for (unsigned int row = 0; row < num_derivatives; row++)
848
{
849
for (unsigned int col = 0; col < num_derivatives; col++)
850
{
851
values[row] += transform[row][col]*derivatives[col];
852
}
853
}
854
// Delete pointer to array of derivatives on FIAT element
855
delete [] derivatives;
856
857
// Delete pointer to array of combinations of derivatives and transform
858
for (unsigned int row = 0; row < num_derivatives; row++)
859
{
860
delete [] combinations[row];
861
delete [] transform[row];
862
}
863
864
delete [] combinations;
865
delete [] transform;
866
}
867
868
/// Evaluate order n derivatives of all basis functions at given point in cell
869
virtual void evaluate_basis_derivatives_all(unsigned int n,
870
double* values,
871
const double* coordinates,
872
const ufc::cell& c) const
873
{
874
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
875
}
876
877
/// Evaluate linear functional for dof i on the function f
878
virtual double evaluate_dof(unsigned int i,
879
const ufc::function& f,
880
const ufc::cell& c) const
881
{
882
// The reference points, direction and weights:
883
static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
884
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
885
static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
886
887
const double * const * x = c.coordinates;
888
double result = 0.0;
889
// Iterate over the points:
890
// Evaluate basis functions for affine mapping
891
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
892
const double w1 = X[i][0][0];
893
const double w2 = X[i][0][1];
894
895
// Compute affine mapping y = F(X)
896
double y[2];
897
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
898
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
899
900
// Evaluate function at physical points
901
double values[1];
902
f.evaluate(values, y, c);
903
904
// Map function values using appropriate mapping
905
// Affine map: Do nothing
906
907
// Note that we do not map the weights (yet).
908
909
// Take directional components
910
for(int k = 0; k < 1; k++)
911
result += values[k]*D[i][0][k];
912
// Multiply by weights
913
result *= W[i][0];
914
915
return result;
916
}
917
918
/// Evaluate linear functionals for all dofs on the function f
919
virtual void evaluate_dofs(double* values,
920
const ufc::function& f,
921
const ufc::cell& c) const
922
{
923
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
924
}
925
926
/// Interpolate vertex values from dof values
927
virtual void interpolate_vertex_values(double* vertex_values,
928
const double* dof_values,
929
const ufc::cell& c) const
930
{
931
// Evaluate at vertices and use affine mapping
932
vertex_values[0] = dof_values[0];
933
vertex_values[1] = dof_values[1];
934
vertex_values[2] = dof_values[2];
935
}
936
937
/// Return the number of sub elements (for a mixed element)
938
virtual unsigned int num_sub_elements() const
939
{
940
return 1;
941
}
942
943
/// Create a new finite element for sub element i (for a mixed element)
944
virtual ufc::finite_element* create_sub_element(unsigned int i) const
945
{
946
return new biharmonic_0_finite_element_1();
947
}
948
949
};
950
951
/// This class defines the interface for a finite element.
952
953
class biharmonic_0_finite_element_2: public ufc::finite_element
954
{
955
public:
956
957
/// Constructor
958
biharmonic_0_finite_element_2() : ufc::finite_element()
959
{
960
// Do nothing
961
}
962
963
/// Destructor
964
virtual ~biharmonic_0_finite_element_2()
965
{
966
// Do nothing
967
}
968
969
/// Return a string identifying the finite element
970
virtual const char* signature() const
971
{
972
return "FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
973
}
974
975
/// Return the cell shape
976
virtual ufc::shape cell_shape() const
977
{
978
return ufc::triangle;
979
}
980
981
/// Return the dimension of the finite element function space
982
virtual unsigned int space_dimension() const
983
{
984
return 1;
985
}
986
987
/// Return the rank of the value space
988
virtual unsigned int value_rank() const
989
{
990
return 0;
991
}
992
993
/// Return the dimension of the value space for axis i
994
virtual unsigned int value_dimension(unsigned int i) const
995
{
996
return 1;
997
}
998
999
/// Evaluate basis function i at given point in cell
1000
virtual void evaluate_basis(unsigned int i,
1001
double* values,
1002
const double* coordinates,
1003
const ufc::cell& c) const
1004
{
1005
// Extract vertex coordinates
1006
const double * const * element_coordinates = c.coordinates;
1007
1008
// Compute Jacobian of affine map from reference cell
1009
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1010
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1011
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1012
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1013
1014
// Compute determinant of Jacobian
1015
const double detJ = J_00*J_11 - J_01*J_10;
1016
1017
// Compute inverse of Jacobian
1018
1019
// Get coordinates and map to the reference (UFC) element
1020
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1021
element_coordinates[0][0]*element_coordinates[2][1] +\
1022
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1023
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1024
element_coordinates[1][0]*element_coordinates[0][1] -\
1025
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1026
1027
// Map coordinates to the reference square
1028
if (std::abs(y - 1.0) < 1e-08)
1029
x = -1.0;
1030
else
1031
x = 2.0 *x/(1.0 - y) - 1.0;
1032
y = 2.0*y - 1.0;
1033
1034
// Reset values
1035
*values = 0;
1036
1037
// Map degree of freedom to element degree of freedom
1038
const unsigned int dof = i;
1039
1040
// Generate scalings
1041
const double scalings_y_0 = 1;
1042
1043
// Compute psitilde_a
1044
const double psitilde_a_0 = 1;
1045
1046
// Compute psitilde_bs
1047
const double psitilde_bs_0_0 = 1;
1048
1049
// Compute basisvalues
1050
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1051
1052
// Table(s) of coefficients
1053
static const double coefficients0[1][1] = \
1054
{{1.41421356}};
1055
1056
// Extract relevant coefficients
1057
const double coeff0_0 = coefficients0[dof][0];
1058
1059
// Compute value(s)
1060
*values = coeff0_0*basisvalue0;
1061
}
1062
1063
/// Evaluate all basis functions at given point in cell
1064
virtual void evaluate_basis_all(double* values,
1065
const double* coordinates,
1066
const ufc::cell& c) const
1067
{
1068
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1069
}
1070
1071
/// Evaluate order n derivatives of basis function i at given point in cell
1072
virtual void evaluate_basis_derivatives(unsigned int i,
1073
unsigned int n,
1074
double* values,
1075
const double* coordinates,
1076
const ufc::cell& c) const
1077
{
1078
// Extract vertex coordinates
1079
const double * const * element_coordinates = c.coordinates;
1080
1081
// Compute Jacobian of affine map from reference cell
1082
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1083
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1084
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1085
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1086
1087
// Compute determinant of Jacobian
1088
const double detJ = J_00*J_11 - J_01*J_10;
1089
1090
// Compute inverse of Jacobian
1091
1092
// Get coordinates and map to the reference (UFC) element
1093
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1094
element_coordinates[0][0]*element_coordinates[2][1] +\
1095
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1096
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1097
element_coordinates[1][0]*element_coordinates[0][1] -\
1098
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1099
1100
// Map coordinates to the reference square
1101
if (std::abs(y - 1.0) < 1e-08)
1102
x = -1.0;
1103
else
1104
x = 2.0 *x/(1.0 - y) - 1.0;
1105
y = 2.0*y - 1.0;
1106
1107
// Compute number of derivatives
1108
unsigned int num_derivatives = 1;
1109
1110
for (unsigned int j = 0; j < n; j++)
1111
num_derivatives *= 2;
1112
1113
1114
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1115
unsigned int **combinations = new unsigned int *[num_derivatives];
1116
1117
for (unsigned int j = 0; j < num_derivatives; j++)
1118
{
1119
combinations[j] = new unsigned int [n];
1120
for (unsigned int k = 0; k < n; k++)
1121
combinations[j][k] = 0;
1122
}
1123
1124
// Generate combinations of derivatives
1125
for (unsigned int row = 1; row < num_derivatives; row++)
1126
{
1127
for (unsigned int num = 0; num < row; num++)
1128
{
1129
for (unsigned int col = n-1; col+1 > 0; col--)
1130
{
1131
if (combinations[row][col] + 1 > 1)
1132
combinations[row][col] = 0;
1133
else
1134
{
1135
combinations[row][col] += 1;
1136
break;
1137
}
1138
}
1139
}
1140
}
1141
1142
// Compute inverse of Jacobian
1143
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1144
1145
// Declare transformation matrix
1146
// Declare pointer to two dimensional array and initialise
1147
double **transform = new double *[num_derivatives];
1148
1149
for (unsigned int j = 0; j < num_derivatives; j++)
1150
{
1151
transform[j] = new double [num_derivatives];
1152
for (unsigned int k = 0; k < num_derivatives; k++)
1153
transform[j][k] = 1;
1154
}
1155
1156
// Construct transformation matrix
1157
for (unsigned int row = 0; row < num_derivatives; row++)
1158
{
1159
for (unsigned int col = 0; col < num_derivatives; col++)
1160
{
1161
for (unsigned int k = 0; k < n; k++)
1162
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1163
}
1164
}
1165
1166
// Reset values
1167
for (unsigned int j = 0; j < 1*num_derivatives; j++)
1168
values[j] = 0;
1169
1170
// Map degree of freedom to element degree of freedom
1171
const unsigned int dof = i;
1172
1173
// Generate scalings
1174
const double scalings_y_0 = 1;
1175
1176
// Compute psitilde_a
1177
const double psitilde_a_0 = 1;
1178
1179
// Compute psitilde_bs
1180
const double psitilde_bs_0_0 = 1;
1181
1182
// Compute basisvalues
1183
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1184
1185
// Table(s) of coefficients
1186
static const double coefficients0[1][1] = \
1187
{{1.41421356}};
1188
1189
// Interesting (new) part
1190
// Tables of derivatives of the polynomial base (transpose)
1191
static const double dmats0[1][1] = \
1192
{{0}};
1193
1194
static const double dmats1[1][1] = \
1195
{{0}};
1196
1197
// Compute reference derivatives
1198
// Declare pointer to array of derivatives on FIAT element
1199
double *derivatives = new double [num_derivatives];
1200
1201
// Declare coefficients
1202
double coeff0_0 = 0;
1203
1204
// Declare new coefficients
1205
double new_coeff0_0 = 0;
1206
1207
// Loop possible derivatives
1208
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1209
{
1210
// Get values from coefficients array
1211
new_coeff0_0 = coefficients0[dof][0];
1212
1213
// Loop derivative order
1214
for (unsigned int j = 0; j < n; j++)
1215
{
1216
// Update old coefficients
1217
coeff0_0 = new_coeff0_0;
1218
1219
if(combinations[deriv_num][j] == 0)
1220
{
1221
new_coeff0_0 = coeff0_0*dmats0[0][0];
1222
}
1223
if(combinations[deriv_num][j] == 1)
1224
{
1225
new_coeff0_0 = coeff0_0*dmats1[0][0];
1226
}
1227
1228
}
1229
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1230
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
1231
}
1232
1233
// Transform derivatives back to physical element
1234
for (unsigned int row = 0; row < num_derivatives; row++)
1235
{
1236
for (unsigned int col = 0; col < num_derivatives; col++)
1237
{
1238
values[row] += transform[row][col]*derivatives[col];
1239
}
1240
}
1241
// Delete pointer to array of derivatives on FIAT element
1242
delete [] derivatives;
1243
1244
// Delete pointer to array of combinations of derivatives and transform
1245
for (unsigned int row = 0; row < num_derivatives; row++)
1246
{
1247
delete [] combinations[row];
1248
delete [] transform[row];
1249
}
1250
1251
delete [] combinations;
1252
delete [] transform;
1253
}
1254
1255
/// Evaluate order n derivatives of all basis functions at given point in cell
1256
virtual void evaluate_basis_derivatives_all(unsigned int n,
1257
double* values,
1258
const double* coordinates,
1259
const ufc::cell& c) const
1260
{
1261
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1262
}
1263
1264
/// Evaluate linear functional for dof i on the function f
1265
virtual double evaluate_dof(unsigned int i,
1266
const ufc::function& f,
1267
const ufc::cell& c) const
1268
{
1269
// The reference points, direction and weights:
1270
static const double X[1][1][2] = {{{0.333333333, 0.333333333}}};
1271
static const double W[1][1] = {{1}};
1272
static const double D[1][1][1] = {{{1}}};
1273
1274
const double * const * x = c.coordinates;
1275
double result = 0.0;
1276
// Iterate over the points:
1277
// Evaluate basis functions for affine mapping
1278
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1279
const double w1 = X[i][0][0];
1280
const double w2 = X[i][0][1];
1281
1282
// Compute affine mapping y = F(X)
1283
double y[2];
1284
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1285
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1286
1287
// Evaluate function at physical points
1288
double values[1];
1289
f.evaluate(values, y, c);
1290
1291
// Map function values using appropriate mapping
1292
// Affine map: Do nothing
1293
1294
// Note that we do not map the weights (yet).
1295
1296
// Take directional components
1297
for(int k = 0; k < 1; k++)
1298
result += values[k]*D[i][0][k];
1299
// Multiply by weights
1300
result *= W[i][0];
1301
1302
return result;
1303
}
1304
1305
/// Evaluate linear functionals for all dofs on the function f
1306
virtual void evaluate_dofs(double* values,
1307
const ufc::function& f,
1308
const ufc::cell& c) const
1309
{
1310
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1311
}
1312
1313
/// Interpolate vertex values from dof values
1314
virtual void interpolate_vertex_values(double* vertex_values,
1315
const double* dof_values,
1316
const ufc::cell& c) const
1317
{
1318
// Evaluate at vertices and use affine mapping
1319
vertex_values[0] = dof_values[0];
1320
vertex_values[1] = dof_values[0];
1321
vertex_values[2] = dof_values[0];
1322
}
1323
1324
/// Return the number of sub elements (for a mixed element)
1325
virtual unsigned int num_sub_elements() const
1326
{
1327
return 1;
1328
}
1329
1330
/// Create a new finite element for sub element i (for a mixed element)
1331
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1332
{
1333
return new biharmonic_0_finite_element_2();
1334
}
1335
1336
};
1337
1338
/// This class defines the interface for a finite element.
1339
1340
class biharmonic_0_finite_element_3: public ufc::finite_element
1341
{
1342
public:
1343
1344
/// Constructor
1345
biharmonic_0_finite_element_3() : ufc::finite_element()
1346
{
1347
// Do nothing
1348
}
1349
1350
/// Destructor
1351
virtual ~biharmonic_0_finite_element_3()
1352
{
1353
// Do nothing
1354
}
1355
1356
/// Return a string identifying the finite element
1357
virtual const char* signature() const
1358
{
1359
return "FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
1360
}
1361
1362
/// Return the cell shape
1363
virtual ufc::shape cell_shape() const
1364
{
1365
return ufc::triangle;
1366
}
1367
1368
/// Return the dimension of the finite element function space
1369
virtual unsigned int space_dimension() const
1370
{
1371
return 1;
1372
}
1373
1374
/// Return the rank of the value space
1375
virtual unsigned int value_rank() const
1376
{
1377
return 0;
1378
}
1379
1380
/// Return the dimension of the value space for axis i
1381
virtual unsigned int value_dimension(unsigned int i) const
1382
{
1383
return 1;
1384
}
1385
1386
/// Evaluate basis function i at given point in cell
1387
virtual void evaluate_basis(unsigned int i,
1388
double* values,
1389
const double* coordinates,
1390
const ufc::cell& c) const
1391
{
1392
// Extract vertex coordinates
1393
const double * const * element_coordinates = c.coordinates;
1394
1395
// Compute Jacobian of affine map from reference cell
1396
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1397
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1398
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1399
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1400
1401
// Compute determinant of Jacobian
1402
const double detJ = J_00*J_11 - J_01*J_10;
1403
1404
// Compute inverse of Jacobian
1405
1406
// Get coordinates and map to the reference (UFC) element
1407
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1408
element_coordinates[0][0]*element_coordinates[2][1] +\
1409
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1410
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1411
element_coordinates[1][0]*element_coordinates[0][1] -\
1412
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1413
1414
// Map coordinates to the reference square
1415
if (std::abs(y - 1.0) < 1e-08)
1416
x = -1.0;
1417
else
1418
x = 2.0 *x/(1.0 - y) - 1.0;
1419
y = 2.0*y - 1.0;
1420
1421
// Reset values
1422
*values = 0;
1423
1424
// Map degree of freedom to element degree of freedom
1425
const unsigned int dof = i;
1426
1427
// Generate scalings
1428
const double scalings_y_0 = 1;
1429
1430
// Compute psitilde_a
1431
const double psitilde_a_0 = 1;
1432
1433
// Compute psitilde_bs
1434
const double psitilde_bs_0_0 = 1;
1435
1436
// Compute basisvalues
1437
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1438
1439
// Table(s) of coefficients
1440
static const double coefficients0[1][1] = \
1441
{{1.41421356}};
1442
1443
// Extract relevant coefficients
1444
const double coeff0_0 = coefficients0[dof][0];
1445
1446
// Compute value(s)
1447
*values = coeff0_0*basisvalue0;
1448
}
1449
1450
/// Evaluate all basis functions at given point in cell
1451
virtual void evaluate_basis_all(double* values,
1452
const double* coordinates,
1453
const ufc::cell& c) const
1454
{
1455
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1456
}
1457
1458
/// Evaluate order n derivatives of basis function i at given point in cell
1459
virtual void evaluate_basis_derivatives(unsigned int i,
1460
unsigned int n,
1461
double* values,
1462
const double* coordinates,
1463
const ufc::cell& c) const
1464
{
1465
// Extract vertex coordinates
1466
const double * const * element_coordinates = c.coordinates;
1467
1468
// Compute Jacobian of affine map from reference cell
1469
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1470
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1471
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1472
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1473
1474
// Compute determinant of Jacobian
1475
const double detJ = J_00*J_11 - J_01*J_10;
1476
1477
// Compute inverse of Jacobian
1478
1479
// Get coordinates and map to the reference (UFC) element
1480
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1481
element_coordinates[0][0]*element_coordinates[2][1] +\
1482
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1483
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1484
element_coordinates[1][0]*element_coordinates[0][1] -\
1485
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1486
1487
// Map coordinates to the reference square
1488
if (std::abs(y - 1.0) < 1e-08)
1489
x = -1.0;
1490
else
1491
x = 2.0 *x/(1.0 - y) - 1.0;
1492
y = 2.0*y - 1.0;
1493
1494
// Compute number of derivatives
1495
unsigned int num_derivatives = 1;
1496
1497
for (unsigned int j = 0; j < n; j++)
1498
num_derivatives *= 2;
1499
1500
1501
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1502
unsigned int **combinations = new unsigned int *[num_derivatives];
1503
1504
for (unsigned int j = 0; j < num_derivatives; j++)
1505
{
1506
combinations[j] = new unsigned int [n];
1507
for (unsigned int k = 0; k < n; k++)
1508
combinations[j][k] = 0;
1509
}
1510
1511
// Generate combinations of derivatives
1512
for (unsigned int row = 1; row < num_derivatives; row++)
1513
{
1514
for (unsigned int num = 0; num < row; num++)
1515
{
1516
for (unsigned int col = n-1; col+1 > 0; col--)
1517
{
1518
if (combinations[row][col] + 1 > 1)
1519
combinations[row][col] = 0;
1520
else
1521
{
1522
combinations[row][col] += 1;
1523
break;
1524
}
1525
}
1526
}
1527
}
1528
1529
// Compute inverse of Jacobian
1530
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1531
1532
// Declare transformation matrix
1533
// Declare pointer to two dimensional array and initialise
1534
double **transform = new double *[num_derivatives];
1535
1536
for (unsigned int j = 0; j < num_derivatives; j++)
1537
{
1538
transform[j] = new double [num_derivatives];
1539
for (unsigned int k = 0; k < num_derivatives; k++)
1540
transform[j][k] = 1;
1541
}
1542
1543
// Construct transformation matrix
1544
for (unsigned int row = 0; row < num_derivatives; row++)
1545
{
1546
for (unsigned int col = 0; col < num_derivatives; col++)
1547
{
1548
for (unsigned int k = 0; k < n; k++)
1549
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1550
}
1551
}
1552
1553
// Reset values
1554
for (unsigned int j = 0; j < 1*num_derivatives; j++)
1555
values[j] = 0;
1556
1557
// Map degree of freedom to element degree of freedom
1558
const unsigned int dof = i;
1559
1560
// Generate scalings
1561
const double scalings_y_0 = 1;
1562
1563
// Compute psitilde_a
1564
const double psitilde_a_0 = 1;
1565
1566
// Compute psitilde_bs
1567
const double psitilde_bs_0_0 = 1;
1568
1569
// Compute basisvalues
1570
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1571
1572
// Table(s) of coefficients
1573
static const double coefficients0[1][1] = \
1574
{{1.41421356}};
1575
1576
// Interesting (new) part
1577
// Tables of derivatives of the polynomial base (transpose)
1578
static const double dmats0[1][1] = \
1579
{{0}};
1580
1581
static const double dmats1[1][1] = \
1582
{{0}};
1583
1584
// Compute reference derivatives
1585
// Declare pointer to array of derivatives on FIAT element
1586
double *derivatives = new double [num_derivatives];
1587
1588
// Declare coefficients
1589
double coeff0_0 = 0;
1590
1591
// Declare new coefficients
1592
double new_coeff0_0 = 0;
1593
1594
// Loop possible derivatives
1595
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1596
{
1597
// Get values from coefficients array
1598
new_coeff0_0 = coefficients0[dof][0];
1599
1600
// Loop derivative order
1601
for (unsigned int j = 0; j < n; j++)
1602
{
1603
// Update old coefficients
1604
coeff0_0 = new_coeff0_0;
1605
1606
if(combinations[deriv_num][j] == 0)
1607
{
1608
new_coeff0_0 = coeff0_0*dmats0[0][0];
1609
}
1610
if(combinations[deriv_num][j] == 1)
1611
{
1612
new_coeff0_0 = coeff0_0*dmats1[0][0];
1613
}
1614
1615
}
1616
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1617
derivatives[deriv_num] = new_coeff0_0*basisvalue0;
1618
}
1619
1620
// Transform derivatives back to physical element
1621
for (unsigned int row = 0; row < num_derivatives; row++)
1622
{
1623
for (unsigned int col = 0; col < num_derivatives; col++)
1624
{
1625
values[row] += transform[row][col]*derivatives[col];
1626
}
1627
}
1628
// Delete pointer to array of derivatives on FIAT element
1629
delete [] derivatives;
1630
1631
// Delete pointer to array of combinations of derivatives and transform
1632
for (unsigned int row = 0; row < num_derivatives; row++)
1633
{
1634
delete [] combinations[row];
1635
delete [] transform[row];
1636
}
1637
1638
delete [] combinations;
1639
delete [] transform;
1640
}
1641
1642
/// Evaluate order n derivatives of all basis functions at given point in cell
1643
virtual void evaluate_basis_derivatives_all(unsigned int n,
1644
double* values,
1645
const double* coordinates,
1646
const ufc::cell& c) const
1647
{
1648
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1649
}
1650
1651
/// Evaluate linear functional for dof i on the function f
1652
virtual double evaluate_dof(unsigned int i,
1653
const ufc::function& f,
1654
const ufc::cell& c) const
1655
{
1656
// The reference points, direction and weights:
1657
static const double X[1][1][2] = {{{0.333333333, 0.333333333}}};
1658
static const double W[1][1] = {{1}};
1659
static const double D[1][1][1] = {{{1}}};
1660
1661
const double * const * x = c.coordinates;
1662
double result = 0.0;
1663
// Iterate over the points:
1664
// Evaluate basis functions for affine mapping
1665
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1666
const double w1 = X[i][0][0];
1667
const double w2 = X[i][0][1];
1668
1669
// Compute affine mapping y = F(X)
1670
double y[2];
1671
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1672
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1673
1674
// Evaluate function at physical points
1675
double values[1];
1676
f.evaluate(values, y, c);
1677
1678
// Map function values using appropriate mapping
1679
// Affine map: Do nothing
1680
1681
// Note that we do not map the weights (yet).
1682
1683
// Take directional components
1684
for(int k = 0; k < 1; k++)
1685
result += values[k]*D[i][0][k];
1686
// Multiply by weights
1687
result *= W[i][0];
1688
1689
return result;
1690
}
1691
1692
/// Evaluate linear functionals for all dofs on the function f
1693
virtual void evaluate_dofs(double* values,
1694
const ufc::function& f,
1695
const ufc::cell& c) const
1696
{
1697
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1698
}
1699
1700
/// Interpolate vertex values from dof values
1701
virtual void interpolate_vertex_values(double* vertex_values,
1702
const double* dof_values,
1703
const ufc::cell& c) const
1704
{
1705
// Evaluate at vertices and use affine mapping
1706
vertex_values[0] = dof_values[0];
1707
vertex_values[1] = dof_values[0];
1708
vertex_values[2] = dof_values[0];
1709
}
1710
1711
/// Return the number of sub elements (for a mixed element)
1712
virtual unsigned int num_sub_elements() const
1713
{
1714
return 1;
1715
}
1716
1717
/// Create a new finite element for sub element i (for a mixed element)
1718
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1719
{
1720
return new biharmonic_0_finite_element_3();
1721
}
1722
1723
};
1724
1725
/// This class defines the interface for a local-to-global mapping of
1726
/// degrees of freedom (dofs).
1727
1728
class biharmonic_0_dof_map_0: public ufc::dof_map
1729
{
1730
private:
1731
1732
unsigned int __global_dimension;
1733
1734
public:
1735
1736
/// Constructor
1737
biharmonic_0_dof_map_0() : ufc::dof_map()
1738
{
1739
__global_dimension = 0;
1740
}
1741
1742
/// Destructor
1743
virtual ~biharmonic_0_dof_map_0()
1744
{
1745
// Do nothing
1746
}
1747
1748
/// Return a string identifying the dof map
1749
virtual const char* signature() const
1750
{
1751
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
1752
}
1753
1754
/// Return true iff mesh entities of topological dimension d are needed
1755
virtual bool needs_mesh_entities(unsigned int d) const
1756
{
1757
switch ( d )
1758
{
1759
case 0:
1760
return true;
1761
break;
1762
case 1:
1763
return true;
1764
break;
1765
case 2:
1766
return false;
1767
break;
1768
}
1769
return false;
1770
}
1771
1772
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1773
virtual bool init_mesh(const ufc::mesh& m)
1774
{
1775
__global_dimension = m.num_entities[0] + m.num_entities[1];
1776
return false;
1777
}
1778
1779
/// Initialize dof map for given cell
1780
virtual void init_cell(const ufc::mesh& m,
1781
const ufc::cell& c)
1782
{
1783
// Do nothing
1784
}
1785
1786
/// Finish initialization of dof map for cells
1787
virtual void init_cell_finalize()
1788
{
1789
// Do nothing
1790
}
1791
1792
/// Return the dimension of the global finite element function space
1793
virtual unsigned int global_dimension() const
1794
{
1795
return __global_dimension;
1796
}
1797
1798
/// Return the dimension of the local finite element function space for a cell
1799
virtual unsigned int local_dimension(const ufc::cell& c) const
1800
{
1801
return 6;
1802
}
1803
1804
/// Return the maximum dimension of the local finite element function space
1805
virtual unsigned int max_local_dimension() const
1806
{
1807
return 6;
1808
}
1809
1810
// Return the geometric dimension of the coordinates this dof map provides
1811
virtual unsigned int geometric_dimension() const
1812
{
1813
return 2;
1814
}
1815
1816
/// Return the number of dofs on each cell facet
1817
virtual unsigned int num_facet_dofs() const
1818
{
1819
return 3;
1820
}
1821
1822
/// Return the number of dofs associated with each cell entity of dimension d
1823
virtual unsigned int num_entity_dofs(unsigned int d) const
1824
{
1825
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1826
}
1827
1828
/// Tabulate the local-to-global mapping of dofs on a cell
1829
virtual void tabulate_dofs(unsigned int* dofs,
1830
const ufc::mesh& m,
1831
const ufc::cell& c) const
1832
{
1833
dofs[0] = c.entity_indices[0][0];
1834
dofs[1] = c.entity_indices[0][1];
1835
dofs[2] = c.entity_indices[0][2];
1836
unsigned int offset = m.num_entities[0];
1837
dofs[3] = offset + c.entity_indices[1][0];
1838
dofs[4] = offset + c.entity_indices[1][1];
1839
dofs[5] = offset + c.entity_indices[1][2];
1840
}
1841
1842
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1843
virtual void tabulate_facet_dofs(unsigned int* dofs,
1844
unsigned int facet) const
1845
{
1846
switch ( facet )
1847
{
1848
case 0:
1849
dofs[0] = 1;
1850
dofs[1] = 2;
1851
dofs[2] = 3;
1852
break;
1853
case 1:
1854
dofs[0] = 0;
1855
dofs[1] = 2;
1856
dofs[2] = 4;
1857
break;
1858
case 2:
1859
dofs[0] = 0;
1860
dofs[1] = 1;
1861
dofs[2] = 5;
1862
break;
1863
}
1864
}
1865
1866
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1867
virtual void tabulate_entity_dofs(unsigned int* dofs,
1868
unsigned int d, unsigned int i) const
1869
{
1870
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1871
}
1872
1873
/// Tabulate the coordinates of all dofs on a cell
1874
virtual void tabulate_coordinates(double** coordinates,
1875
const ufc::cell& c) const
1876
{
1877
const double * const * x = c.coordinates;
1878
coordinates[0][0] = x[0][0];
1879
coordinates[0][1] = x[0][1];
1880
coordinates[1][0] = x[1][0];
1881
coordinates[1][1] = x[1][1];
1882
coordinates[2][0] = x[2][0];
1883
coordinates[2][1] = x[2][1];
1884
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
1885
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
1886
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
1887
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
1888
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
1889
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
1890
}
1891
1892
/// Return the number of sub dof maps (for a mixed element)
1893
virtual unsigned int num_sub_dof_maps() const
1894
{
1895
return 1;
1896
}
1897
1898
/// Create a new dof_map for sub dof map i (for a mixed element)
1899
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1900
{
1901
return new biharmonic_0_dof_map_0();
1902
}
1903
1904
};
1905
1906
/// This class defines the interface for a local-to-global mapping of
1907
/// degrees of freedom (dofs).
1908
1909
class biharmonic_0_dof_map_1: public ufc::dof_map
1910
{
1911
private:
1912
1913
unsigned int __global_dimension;
1914
1915
public:
1916
1917
/// Constructor
1918
biharmonic_0_dof_map_1() : ufc::dof_map()
1919
{
1920
__global_dimension = 0;
1921
}
1922
1923
/// Destructor
1924
virtual ~biharmonic_0_dof_map_1()
1925
{
1926
// Do nothing
1927
}
1928
1929
/// Return a string identifying the dof map
1930
virtual const char* signature() const
1931
{
1932
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
1933
}
1934
1935
/// Return true iff mesh entities of topological dimension d are needed
1936
virtual bool needs_mesh_entities(unsigned int d) const
1937
{
1938
switch ( d )
1939
{
1940
case 0:
1941
return true;
1942
break;
1943
case 1:
1944
return true;
1945
break;
1946
case 2:
1947
return false;
1948
break;
1949
}
1950
return false;
1951
}
1952
1953
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1954
virtual bool init_mesh(const ufc::mesh& m)
1955
{
1956
__global_dimension = m.num_entities[0] + m.num_entities[1];
1957
return false;
1958
}
1959
1960
/// Initialize dof map for given cell
1961
virtual void init_cell(const ufc::mesh& m,
1962
const ufc::cell& c)
1963
{
1964
// Do nothing
1965
}
1966
1967
/// Finish initialization of dof map for cells
1968
virtual void init_cell_finalize()
1969
{
1970
// Do nothing
1971
}
1972
1973
/// Return the dimension of the global finite element function space
1974
virtual unsigned int global_dimension() const
1975
{
1976
return __global_dimension;
1977
}
1978
1979
/// Return the dimension of the local finite element function space for a cell
1980
virtual unsigned int local_dimension(const ufc::cell& c) const
1981
{
1982
return 6;
1983
}
1984
1985
/// Return the maximum dimension of the local finite element function space
1986
virtual unsigned int max_local_dimension() const
1987
{
1988
return 6;
1989
}
1990
1991
// Return the geometric dimension of the coordinates this dof map provides
1992
virtual unsigned int geometric_dimension() const
1993
{
1994
return 2;
1995
}
1996
1997
/// Return the number of dofs on each cell facet
1998
virtual unsigned int num_facet_dofs() const
1999
{
2000
return 3;
2001
}
2002
2003
/// Return the number of dofs associated with each cell entity of dimension d
2004
virtual unsigned int num_entity_dofs(unsigned int d) const
2005
{
2006
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2007
}
2008
2009
/// Tabulate the local-to-global mapping of dofs on a cell
2010
virtual void tabulate_dofs(unsigned int* dofs,
2011
const ufc::mesh& m,
2012
const ufc::cell& c) const
2013
{
2014
dofs[0] = c.entity_indices[0][0];
2015
dofs[1] = c.entity_indices[0][1];
2016
dofs[2] = c.entity_indices[0][2];
2017
unsigned int offset = m.num_entities[0];
2018
dofs[3] = offset + c.entity_indices[1][0];
2019
dofs[4] = offset + c.entity_indices[1][1];
2020
dofs[5] = offset + c.entity_indices[1][2];
2021
}
2022
2023
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2024
virtual void tabulate_facet_dofs(unsigned int* dofs,
2025
unsigned int facet) const
2026
{
2027
switch ( facet )
2028
{
2029
case 0:
2030
dofs[0] = 1;
2031
dofs[1] = 2;
2032
dofs[2] = 3;
2033
break;
2034
case 1:
2035
dofs[0] = 0;
2036
dofs[1] = 2;
2037
dofs[2] = 4;
2038
break;
2039
case 2:
2040
dofs[0] = 0;
2041
dofs[1] = 1;
2042
dofs[2] = 5;
2043
break;
2044
}
2045
}
2046
2047
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2048
virtual void tabulate_entity_dofs(unsigned int* dofs,
2049
unsigned int d, unsigned int i) const
2050
{
2051
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2052
}
2053
2054
/// Tabulate the coordinates of all dofs on a cell
2055
virtual void tabulate_coordinates(double** coordinates,
2056
const ufc::cell& c) const
2057
{
2058
const double * const * x = c.coordinates;
2059
coordinates[0][0] = x[0][0];
2060
coordinates[0][1] = x[0][1];
2061
coordinates[1][0] = x[1][0];
2062
coordinates[1][1] = x[1][1];
2063
coordinates[2][0] = x[2][0];
2064
coordinates[2][1] = x[2][1];
2065
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
2066
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
2067
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
2068
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
2069
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
2070
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
2071
}
2072
2073
/// Return the number of sub dof maps (for a mixed element)
2074
virtual unsigned int num_sub_dof_maps() const
2075
{
2076
return 1;
2077
}
2078
2079
/// Create a new dof_map for sub dof map i (for a mixed element)
2080
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2081
{
2082
return new biharmonic_0_dof_map_1();
2083
}
2084
2085
};
2086
2087
/// This class defines the interface for a local-to-global mapping of
2088
/// degrees of freedom (dofs).
2089
2090
class biharmonic_0_dof_map_2: public ufc::dof_map
2091
{
2092
private:
2093
2094
unsigned int __global_dimension;
2095
2096
public:
2097
2098
/// Constructor
2099
biharmonic_0_dof_map_2() : ufc::dof_map()
2100
{
2101
__global_dimension = 0;
2102
}
2103
2104
/// Destructor
2105
virtual ~biharmonic_0_dof_map_2()
2106
{
2107
// Do nothing
2108
}
2109
2110
/// Return a string identifying the dof map
2111
virtual const char* signature() const
2112
{
2113
return "FFC dof map for FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
2114
}
2115
2116
/// Return true iff mesh entities of topological dimension d are needed
2117
virtual bool needs_mesh_entities(unsigned int d) const
2118
{
2119
switch ( d )
2120
{
2121
case 0:
2122
return false;
2123
break;
2124
case 1:
2125
return false;
2126
break;
2127
case 2:
2128
return true;
2129
break;
2130
}
2131
return false;
2132
}
2133
2134
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2135
virtual bool init_mesh(const ufc::mesh& m)
2136
{
2137
__global_dimension = m.num_entities[2];
2138
return false;
2139
}
2140
2141
/// Initialize dof map for given cell
2142
virtual void init_cell(const ufc::mesh& m,
2143
const ufc::cell& c)
2144
{
2145
// Do nothing
2146
}
2147
2148
/// Finish initialization of dof map for cells
2149
virtual void init_cell_finalize()
2150
{
2151
// Do nothing
2152
}
2153
2154
/// Return the dimension of the global finite element function space
2155
virtual unsigned int global_dimension() const
2156
{
2157
return __global_dimension;
2158
}
2159
2160
/// Return the dimension of the local finite element function space for a cell
2161
virtual unsigned int local_dimension(const ufc::cell& c) const
2162
{
2163
return 1;
2164
}
2165
2166
/// Return the maximum dimension of the local finite element function space
2167
virtual unsigned int max_local_dimension() const
2168
{
2169
return 1;
2170
}
2171
2172
// Return the geometric dimension of the coordinates this dof map provides
2173
virtual unsigned int geometric_dimension() const
2174
{
2175
return 2;
2176
}
2177
2178
/// Return the number of dofs on each cell facet
2179
virtual unsigned int num_facet_dofs() const
2180
{
2181
return 0;
2182
}
2183
2184
/// Return the number of dofs associated with each cell entity of dimension d
2185
virtual unsigned int num_entity_dofs(unsigned int d) const
2186
{
2187
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2188
}
2189
2190
/// Tabulate the local-to-global mapping of dofs on a cell
2191
virtual void tabulate_dofs(unsigned int* dofs,
2192
const ufc::mesh& m,
2193
const ufc::cell& c) const
2194
{
2195
dofs[0] = c.entity_indices[2][0];
2196
}
2197
2198
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2199
virtual void tabulate_facet_dofs(unsigned int* dofs,
2200
unsigned int facet) const
2201
{
2202
switch ( facet )
2203
{
2204
case 0:
2205
2206
break;
2207
case 1:
2208
2209
break;
2210
case 2:
2211
2212
break;
2213
}
2214
}
2215
2216
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2217
virtual void tabulate_entity_dofs(unsigned int* dofs,
2218
unsigned int d, unsigned int i) const
2219
{
2220
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2221
}
2222
2223
/// Tabulate the coordinates of all dofs on a cell
2224
virtual void tabulate_coordinates(double** coordinates,
2225
const ufc::cell& c) const
2226
{
2227
const double * const * x = c.coordinates;
2228
coordinates[0][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
2229
coordinates[0][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
2230
}
2231
2232
/// Return the number of sub dof maps (for a mixed element)
2233
virtual unsigned int num_sub_dof_maps() const
2234
{
2235
return 1;
2236
}
2237
2238
/// Create a new dof_map for sub dof map i (for a mixed element)
2239
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2240
{
2241
return new biharmonic_0_dof_map_2();
2242
}
2243
2244
};
2245
2246
/// This class defines the interface for a local-to-global mapping of
2247
/// degrees of freedom (dofs).
2248
2249
class biharmonic_0_dof_map_3: public ufc::dof_map
2250
{
2251
private:
2252
2253
unsigned int __global_dimension;
2254
2255
public:
2256
2257
/// Constructor
2258
biharmonic_0_dof_map_3() : ufc::dof_map()
2259
{
2260
__global_dimension = 0;
2261
}
2262
2263
/// Destructor
2264
virtual ~biharmonic_0_dof_map_3()
2265
{
2266
// Do nothing
2267
}
2268
2269
/// Return a string identifying the dof map
2270
virtual const char* signature() const
2271
{
2272
return "FFC dof map for FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
2273
}
2274
2275
/// Return true iff mesh entities of topological dimension d are needed
2276
virtual bool needs_mesh_entities(unsigned int d) const
2277
{
2278
switch ( d )
2279
{
2280
case 0:
2281
return false;
2282
break;
2283
case 1:
2284
return false;
2285
break;
2286
case 2:
2287
return true;
2288
break;
2289
}
2290
return false;
2291
}
2292
2293
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2294
virtual bool init_mesh(const ufc::mesh& m)
2295
{
2296
__global_dimension = m.num_entities[2];
2297
return false;
2298
}
2299
2300
/// Initialize dof map for given cell
2301
virtual void init_cell(const ufc::mesh& m,
2302
const ufc::cell& c)
2303
{
2304
// Do nothing
2305
}
2306
2307
/// Finish initialization of dof map for cells
2308
virtual void init_cell_finalize()
2309
{
2310
// Do nothing
2311
}
2312
2313
/// Return the dimension of the global finite element function space
2314
virtual unsigned int global_dimension() const
2315
{
2316
return __global_dimension;
2317
}
2318
2319
/// Return the dimension of the local finite element function space for a cell
2320
virtual unsigned int local_dimension(const ufc::cell& c) const
2321
{
2322
return 1;
2323
}
2324
2325
/// Return the maximum dimension of the local finite element function space
2326
virtual unsigned int max_local_dimension() const
2327
{
2328
return 1;
2329
}
2330
2331
// Return the geometric dimension of the coordinates this dof map provides
2332
virtual unsigned int geometric_dimension() const
2333
{
2334
return 2;
2335
}
2336
2337
/// Return the number of dofs on each cell facet
2338
virtual unsigned int num_facet_dofs() const
2339
{
2340
return 0;
2341
}
2342
2343
/// Return the number of dofs associated with each cell entity of dimension d
2344
virtual unsigned int num_entity_dofs(unsigned int d) const
2345
{
2346
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2347
}
2348
2349
/// Tabulate the local-to-global mapping of dofs on a cell
2350
virtual void tabulate_dofs(unsigned int* dofs,
2351
const ufc::mesh& m,
2352
const ufc::cell& c) const
2353
{
2354
dofs[0] = c.entity_indices[2][0];
2355
}
2356
2357
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2358
virtual void tabulate_facet_dofs(unsigned int* dofs,
2359
unsigned int facet) const
2360
{
2361
switch ( facet )
2362
{
2363
case 0:
2364
2365
break;
2366
case 1:
2367
2368
break;
2369
case 2:
2370
2371
break;
2372
}
2373
}
2374
2375
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2376
virtual void tabulate_entity_dofs(unsigned int* dofs,
2377
unsigned int d, unsigned int i) const
2378
{
2379
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2380
}
2381
2382
/// Tabulate the coordinates of all dofs on a cell
2383
virtual void tabulate_coordinates(double** coordinates,
2384
const ufc::cell& c) const
2385
{
2386
const double * const * x = c.coordinates;
2387
coordinates[0][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];
2388
coordinates[0][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];
2389
}
2390
2391
/// Return the number of sub dof maps (for a mixed element)
2392
virtual unsigned int num_sub_dof_maps() const
2393
{
2394
return 1;
2395
}
2396
2397
/// Create a new dof_map for sub dof map i (for a mixed element)
2398
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2399
{
2400
return new biharmonic_0_dof_map_3();
2401
}
2402
2403
};
2404
2405
/// This class defines the interface for the tabulation of the cell
2406
/// tensor corresponding to the local contribution to a form from
2407
/// the integral over a cell.
2408
2409
class biharmonic_0_cell_integral_0_quadrature: public ufc::cell_integral
2410
{
2411
public:
2412
2413
/// Constructor
2414
biharmonic_0_cell_integral_0_quadrature() : ufc::cell_integral()
2415
{
2416
// Do nothing
2417
}
2418
2419
/// Destructor
2420
virtual ~biharmonic_0_cell_integral_0_quadrature()
2421
{
2422
// Do nothing
2423
}
2424
2425
/// Tabulate the tensor for the contribution from a local cell
2426
virtual void tabulate_tensor(double* A,
2427
const double * const * w,
2428
const ufc::cell& c) const
2429
{
2430
// Extract vertex coordinates
2431
const double * const * x = c.coordinates;
2432
2433
// Compute Jacobian of affine map from reference cell
2434
const double J_00 = x[1][0] - x[0][0];
2435
const double J_01 = x[2][0] - x[0][0];
2436
const double J_10 = x[1][1] - x[0][1];
2437
const double J_11 = x[2][1] - x[0][1];
2438
2439
// Compute determinant of Jacobian
2440
double detJ = J_00*J_11 - J_01*J_10;
2441
2442
// Compute inverse of Jacobian
2443
const double Jinv_00 = J_11 / detJ;
2444
const double Jinv_01 = -J_01 / detJ;
2445
const double Jinv_10 = -J_10 / detJ;
2446
const double Jinv_11 = J_00 / detJ;
2447
2448
// Set scale factor
2449
const double det = std::abs(detJ);
2450
2451
2452
// Array of quadrature weights
2453
static const double W1 = 0.5;
2454
// Quadrature points on the UFC reference element: (0.333333333, 0.333333333)
2455
2456
// Value of basis functions at quadrature points.
2457
static const double FE0_D02[1][6] = \
2458
{{4, 0, 4, 0, -8, 0}};
2459
2460
static const double FE0_D11[1][6] = \
2461
{{4, 0, 0, 4, -4, -4}};
2462
2463
static const double FE0_D20[1][6] = \
2464
{{4, 4, 0, 0, 0, -8}};
2465
2466
2467
// Compute element tensor using UFL quadrature representation
2468
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
2469
// Total number of operations to compute element tensor: 1800
2470
2471
// Loop quadrature points for integral
2472
// Number of operations to compute element tensor for following IP loop = 1800
2473
// Only 1 integration point, omitting IP loop.
2474
2475
// Number of operations for primary indices: 1800
2476
for (unsigned int j = 0; j < 6; j++)
2477
{
2478
for (unsigned int k = 0; k < 6; k++)
2479
{
2480
// Number of operations to compute entry: 50
2481
A[j*6 + k] += ((Jinv_01*Jinv_01*FE0_D20[0][j] + Jinv_01*Jinv_11*FE0_D11[0][j] + Jinv_11*Jinv_01*FE0_D11[0][j] + Jinv_11*Jinv_11*FE0_D02[0][j]) + (Jinv_00*Jinv_00*FE0_D20[0][j] + Jinv_00*Jinv_10*FE0_D11[0][j] + Jinv_10*Jinv_00*FE0_D11[0][j] + Jinv_10*Jinv_10*FE0_D02[0][j]))*((Jinv_01*Jinv_01*FE0_D20[0][k] + Jinv_01*Jinv_11*FE0_D11[0][k] + Jinv_11*Jinv_01*FE0_D11[0][k] + Jinv_11*Jinv_11*FE0_D02[0][k]) + (Jinv_00*Jinv_00*FE0_D20[0][k] + Jinv_00*Jinv_10*FE0_D11[0][k] + Jinv_10*Jinv_00*FE0_D11[0][k] + Jinv_10*Jinv_10*FE0_D02[0][k]))*W1*det;
2482
}// end loop over 'k'
2483
}// end loop over 'j'
2484
}
2485
2486
};
2487
2488
/// This class defines the interface for the tabulation of the cell
2489
/// tensor corresponding to the local contribution to a form from
2490
/// the integral over a cell.
2491
2492
class biharmonic_0_cell_integral_0: public ufc::cell_integral
2493
{
2494
private:
2495
2496
biharmonic_0_cell_integral_0_quadrature integral_0_quadrature;
2497
2498
public:
2499
2500
/// Constructor
2501
biharmonic_0_cell_integral_0() : ufc::cell_integral()
2502
{
2503
// Do nothing
2504
}
2505
2506
/// Destructor
2507
virtual ~biharmonic_0_cell_integral_0()
2508
{
2509
// Do nothing
2510
}
2511
2512
/// Tabulate the tensor for the contribution from a local cell
2513
virtual void tabulate_tensor(double* A,
2514
const double * const * w,
2515
const ufc::cell& c) const
2516
{
2517
// Reset values of the element tensor block
2518
for (unsigned int j = 0; j < 36; j++)
2519
A[j] = 0;
2520
2521
// Add all contributions to element tensor
2522
integral_0_quadrature.tabulate_tensor(A, w, c);
2523
}
2524
2525
};
2526
2527
/// This class defines the interface for the tabulation of the
2528
/// interior facet tensor corresponding to the local contribution to
2529
/// a form from the integral over an interior facet.
2530
2531
class biharmonic_0_interior_facet_integral_0_quadrature: public ufc::interior_facet_integral
2532
{
2533
public:
2534
2535
/// Constructor
2536
biharmonic_0_interior_facet_integral_0_quadrature() : ufc::interior_facet_integral()
2537
{
2538
// Do nothing
2539
}
2540
2541
/// Destructor
2542
virtual ~biharmonic_0_interior_facet_integral_0_quadrature()
2543
{
2544
// Do nothing
2545
}
2546
2547
/// Tabulate the tensor for the contribution from a local interior facet
2548
virtual void tabulate_tensor(double* A,
2549
const double * const * w,
2550
const ufc::cell& c0,
2551
const ufc::cell& c1,
2552
unsigned int facet0,
2553
unsigned int facet1) const
2554
{
2555
// Extract vertex coordinates
2556
const double * const * x0 = c0.coordinates;
2557
2558
// Compute Jacobian of affine map from reference cell
2559
const double J0_00 = x0[1][0] - x0[0][0];
2560
const double J0_01 = x0[2][0] - x0[0][0];
2561
const double J0_10 = x0[1][1] - x0[0][1];
2562
const double J0_11 = x0[2][1] - x0[0][1];
2563
2564
// Compute determinant of Jacobian
2565
double detJ0 = J0_00*J0_11 - J0_01*J0_10;
2566
2567
// Compute inverse of Jacobian
2568
const double Jinv0_00 = J0_11 / detJ0;
2569
const double Jinv0_01 = -J0_01 / detJ0;
2570
const double Jinv0_10 = -J0_10 / detJ0;
2571
const double Jinv0_11 = J0_00 / detJ0;
2572
2573
// Extract vertex coordinates
2574
const double * const * x1 = c1.coordinates;
2575
2576
// Compute Jacobian of affine map from reference cell
2577
const double J1_00 = x1[1][0] - x1[0][0];
2578
const double J1_01 = x1[2][0] - x1[0][0];
2579
const double J1_10 = x1[1][1] - x1[0][1];
2580
const double J1_11 = x1[2][1] - x1[0][1];
2581
2582
// Compute determinant of Jacobian
2583
double detJ1 = J1_00*J1_11 - J1_01*J1_10;
2584
2585
// Compute inverse of Jacobian
2586
const double Jinv1_00 = J1_11 / detJ1;
2587
const double Jinv1_01 = -J1_01 / detJ1;
2588
const double Jinv1_10 = -J1_10 / detJ1;
2589
const double Jinv1_11 = J1_00 / detJ1;
2590
2591
// Vertices on edges
2592
static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
2593
2594
// Get vertices
2595
const unsigned int v0 = edge_vertices[facet0][0];
2596
const unsigned int v1 = edge_vertices[facet0][1];
2597
2598
// Compute scale factor (length of edge scaled by length of reference interval)
2599
const double dx0 = x0[v1][0] - x0[v0][0];
2600
const double dx1 = x0[v1][1] - x0[v0][1];
2601
const double det = std::sqrt(dx0*dx0 + dx1*dx1);
2602
2603
const bool direction = dx1*(x0[facet0][0] - x0[v0][0]) - dx0*(x0[facet0][1] - x0[v0][1]) < 0;
2604
2605
// Compute facet normals from the facet scale factor constants
2606
const double n00 = direction ? dx1 / det : -dx1 / det;
2607
const double n01 = direction ? -dx0 / det : dx0 / det;
2608
2609
// Compute facet normals from the facet scale factor constants
2610
const double n10 = !direction ? dx1 / det : -dx1 / det;
2611
const double n11 = !direction ? -dx0 / det : dx0 / det;
2612
2613
2614
// Array of quadrature weights
2615
static const double W2[2] = {0.5, 0.5};
2616
// Quadrature points on the UFC reference element: (0.211324865), (0.788675135)
2617
2618
// Value of basis functions at quadrature points.
2619
static const double FE0_f0_D01[2][6] = \
2620
{{1, 0, -0.154700538, 3.15470054, -0.845299462, -3.15470054},
2621
{1, 0, 2.15470054, 0.845299462, -3.15470054, -0.845299462}};
2622
2623
static const double FE0_f0_D02[2][6] = \
2624
{{4, 0, 4, 0, -8, 0},
2625
{4, 0, 4, 0, -8, 0}};
2626
2627
static const double FE0_f0_D10[2][6] = \
2628
{{1, 2.15470054, 0, 0.845299462, -0.845299462, -3.15470054},
2629
{1, -0.154700538, 0, 3.15470054, -3.15470054, -0.845299462}};
2630
2631
static const double FE0_f0_D11[2][6] = \
2632
{{4, 0, 0, 4, -4, -4},
2633
{4, 0, 0, 4, -4, -4}};
2634
2635
static const double FE0_f0_D20[2][6] = \
2636
{{4, 4, 0, 0, 0, -8},
2637
{4, 4, 0, 0, 0, -8}};
2638
2639
static const double FE0_f1_D01[2][6] = \
2640
{{-2.15470054, 0, -0.154700538, 0, 2.30940108, 0},
2641
{0.154700538, 0, 2.15470054, 0, -2.30940108, 0}};
2642
2643
static const double FE0_f1_D10[2][6] = \
2644
{{-2.15470054, -1, 0, 0.845299462, -0.845299462, 3.15470054},
2645
{0.154700538, -1, 0, 3.15470054, -3.15470054, 0.845299462}};
2646
2647
static const double FE0_f2_D01[2][6] = \
2648
{{-2.15470054, 0, -1, 0.845299462, 3.15470054, -0.845299462},
2649
{0.154700538, 0, -1, 3.15470054, 0.845299462, -3.15470054}};
2650
2651
static const double FE0_f2_D10[2][6] = \
2652
{{-2.15470054, -0.154700538, 0, 0, 0, 2.30940108},
2653
{0.154700538, 2.15470054, 0, 0, 0, -2.30940108}};
2654
2655
2656
// Compute element tensor using UFL quadrature representation
2657
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
2658
switch ( facet0 )
2659
{
2660
case 0:
2661
switch ( facet1 )
2662
{
2663
case 0:
2664
{
2665
// Total number of operations to compute element tensor (from this point): 27360
2666
2667
// Loop quadrature points for integral
2668
// Number of operations to compute element tensor for following IP loop = 27360
2669
for (unsigned int ip = 0; ip < 2; ip++)
2670
{
2671
2672
// Number of operations for primary indices: 13680
2673
for (unsigned int j = 0; j < 6; j++)
2674
{
2675
for (unsigned int k = 0; k < 6; k++)
2676
{
2677
// Number of operations to compute entry: 95
2678
A[(j + 6)*12 + k] += (((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*w[1][0]/(w[0][0]) + (((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*-1))*W2[ip]*det;
2679
// Number of operations to compute entry: 95
2680
A[(j + 6)*12 + (k + 6)] += ((((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*-1) + ((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*w[1][0]/(w[0][0]))*W2[ip]*det;
2681
// Number of operations to compute entry: 95
2682
A[j*12 + (k + 6)] += (((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*w[1][0]/(w[0][0]) + (((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*-1))*W2[ip]*det;
2683
// Number of operations to compute entry: 95
2684
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*-1) + ((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2685
}// end loop over 'k'
2686
}// end loop over 'j'
2687
}// end loop over 'ip'
2688
}
2689
break;
2690
case 1:
2691
{
2692
// Total number of operations to compute element tensor (from this point): 27360
2693
2694
// Loop quadrature points for integral
2695
// Number of operations to compute element tensor for following IP loop = 27360
2696
for (unsigned int ip = 0; ip < 2; ip++)
2697
{
2698
2699
// Number of operations for primary indices: 13680
2700
for (unsigned int j = 0; j < 6; j++)
2701
{
2702
for (unsigned int k = 0; k < 6; k++)
2703
{
2704
// Number of operations to compute entry: 95
2705
A[(j + 6)*12 + k] += ((((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*-1) + ((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2706
// Number of operations to compute entry: 95
2707
A[(j + 6)*12 + (k + 6)] += (((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*-1))*W2[ip]*det;
2708
// Number of operations to compute entry: 95
2709
A[j*12 + (k + 6)] += (((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*-1))*W2[ip]*det;
2710
// Number of operations to compute entry: 95
2711
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*-1) + ((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2712
}// end loop over 'k'
2713
}// end loop over 'j'
2714
}// end loop over 'ip'
2715
}
2716
break;
2717
case 2:
2718
{
2719
// Total number of operations to compute element tensor (from this point): 27360
2720
2721
// Loop quadrature points for integral
2722
// Number of operations to compute element tensor for following IP loop = 27360
2723
for (unsigned int ip = 0; ip < 2; ip++)
2724
{
2725
2726
// Number of operations for primary indices: 13680
2727
for (unsigned int j = 0; j < 6; j++)
2728
{
2729
for (unsigned int k = 0; k < 6; k++)
2730
{
2731
// Number of operations to compute entry: 95
2732
A[(j + 6)*12 + k] += (((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*w[1][0]/(w[0][0]) + (((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*-1))*W2[ip]*det;
2733
// Number of operations to compute entry: 95
2734
A[(j + 6)*12 + (k + 6)] += (((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*-1 + ((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*-1))*W2[ip]*det;
2735
// Number of operations to compute entry: 95
2736
A[j*12 + (k + 6)] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*-1 + ((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*-1) + ((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*w[1][0]/(w[0][0]))*W2[ip]*det;
2737
// Number of operations to compute entry: 95
2738
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*-1) + ((Jinv0_01*FE0_f0_D10[ip][j] + Jinv0_11*FE0_f0_D01[ip][j])*n01 + (Jinv0_00*FE0_f0_D10[ip][j] + Jinv0_10*FE0_f0_D01[ip][j])*n00)*((Jinv0_01*FE0_f0_D10[ip][k] + Jinv0_11*FE0_f0_D01[ip][k])*n01 + (Jinv0_00*FE0_f0_D10[ip][k] + Jinv0_10*FE0_f0_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2739
}// end loop over 'k'
2740
}// end loop over 'j'
2741
}// end loop over 'ip'
2742
}
2743
break;
2744
}
2745
break;
2746
case 1:
2747
switch ( facet1 )
2748
{
2749
case 0:
2750
{
2751
// Total number of operations to compute element tensor (from this point): 27360
2752
2753
// Loop quadrature points for integral
2754
// Number of operations to compute element tensor for following IP loop = 27360
2755
for (unsigned int ip = 0; ip < 2; ip++)
2756
{
2757
2758
// Number of operations for primary indices: 13680
2759
for (unsigned int j = 0; j < 6; j++)
2760
{
2761
for (unsigned int k = 0; k < 6; k++)
2762
{
2763
// Number of operations to compute entry: 95
2764
A[(j + 6)*12 + k] += (((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*w[1][0]/(w[0][0]) + (((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*-1))*W2[ip]*det;
2765
// Number of operations to compute entry: 95
2766
A[(j + 6)*12 + (k + 6)] += ((((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*-1) + ((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*w[1][0]/(w[0][0]))*W2[ip]*det;
2767
// Number of operations to compute entry: 95
2768
A[j*12 + (k + 6)] += (((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*w[1][0]/(w[0][0]) + (((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*-1))*W2[ip]*det;
2769
// Number of operations to compute entry: 95
2770
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*-1) + ((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2771
}// end loop over 'k'
2772
}// end loop over 'j'
2773
}// end loop over 'ip'
2774
}
2775
break;
2776
case 1:
2777
{
2778
// Total number of operations to compute element tensor (from this point): 27360
2779
2780
// Loop quadrature points for integral
2781
// Number of operations to compute element tensor for following IP loop = 27360
2782
for (unsigned int ip = 0; ip < 2; ip++)
2783
{
2784
2785
// Number of operations for primary indices: 13680
2786
for (unsigned int j = 0; j < 6; j++)
2787
{
2788
for (unsigned int k = 0; k < 6; k++)
2789
{
2790
// Number of operations to compute entry: 95
2791
A[(j + 6)*12 + k] += ((((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*-1) + ((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2792
// Number of operations to compute entry: 95
2793
A[(j + 6)*12 + (k + 6)] += (((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*-1))*W2[ip]*det;
2794
// Number of operations to compute entry: 95
2795
A[j*12 + (k + 6)] += (((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*-1))*W2[ip]*det;
2796
// Number of operations to compute entry: 95
2797
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*-1) + ((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2798
}// end loop over 'k'
2799
}// end loop over 'j'
2800
}// end loop over 'ip'
2801
}
2802
break;
2803
case 2:
2804
{
2805
// Total number of operations to compute element tensor (from this point): 27360
2806
2807
// Loop quadrature points for integral
2808
// Number of operations to compute element tensor for following IP loop = 27360
2809
for (unsigned int ip = 0; ip < 2; ip++)
2810
{
2811
2812
// Number of operations for primary indices: 13680
2813
for (unsigned int j = 0; j < 6; j++)
2814
{
2815
for (unsigned int k = 0; k < 6; k++)
2816
{
2817
// Number of operations to compute entry: 95
2818
A[(j + 6)*12 + k] += (((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*w[1][0]/(w[0][0]) + (((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*-1))*W2[ip]*det;
2819
// Number of operations to compute entry: 95
2820
A[(j + 6)*12 + (k + 6)] += (((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*-1 + ((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*-1))*W2[ip]*det;
2821
// Number of operations to compute entry: 95
2822
A[j*12 + (k + 6)] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*-1 + ((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*-1) + ((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*w[1][0]/(w[0][0]))*W2[ip]*det;
2823
// Number of operations to compute entry: 95
2824
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*-1) + ((Jinv0_00*FE0_f1_D10[ip][j] + Jinv0_10*FE0_f1_D01[ip][j])*n00 + (Jinv0_01*FE0_f1_D10[ip][j] + Jinv0_11*FE0_f1_D01[ip][j])*n01)*((Jinv0_01*FE0_f1_D10[ip][k] + Jinv0_11*FE0_f1_D01[ip][k])*n01 + (Jinv0_00*FE0_f1_D10[ip][k] + Jinv0_10*FE0_f1_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2825
}// end loop over 'k'
2826
}// end loop over 'j'
2827
}// end loop over 'ip'
2828
}
2829
break;
2830
}
2831
break;
2832
case 2:
2833
switch ( facet1 )
2834
{
2835
case 0:
2836
{
2837
// Total number of operations to compute element tensor (from this point): 27360
2838
2839
// Loop quadrature points for integral
2840
// Number of operations to compute element tensor for following IP loop = 27360
2841
for (unsigned int ip = 0; ip < 2; ip++)
2842
{
2843
2844
// Number of operations for primary indices: 13680
2845
for (unsigned int j = 0; j < 6; j++)
2846
{
2847
for (unsigned int k = 0; k < 6; k++)
2848
{
2849
// Number of operations to compute entry: 95
2850
A[(j + 6)*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*-1) + ((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2851
// Number of operations to compute entry: 95
2852
A[(j + 6)*12 + (k + 6)] += ((((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*-1) + ((Jinv1_01*FE0_f0_D10[ip][j] + Jinv1_11*FE0_f0_D01[ip][j])*n11 + (Jinv1_00*FE0_f0_D10[ip][j] + Jinv1_10*FE0_f0_D01[ip][j])*n10)*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*w[1][0]/(w[0][0]))*W2[ip]*det;
2853
// Number of operations to compute entry: 95
2854
A[j*12 + (k + 6)] += ((((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*-1) + ((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*((Jinv1_00*FE0_f0_D10[ip][k] + Jinv1_10*FE0_f0_D01[ip][k])*n10 + (Jinv1_01*FE0_f0_D10[ip][k] + Jinv1_11*FE0_f0_D01[ip][k])*n11)*w[1][0]/(w[0][0]))*W2[ip]*det;
2855
// Number of operations to compute entry: 95
2856
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*-1) + ((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2857
}// end loop over 'k'
2858
}// end loop over 'j'
2859
}// end loop over 'ip'
2860
}
2861
break;
2862
case 1:
2863
{
2864
// Total number of operations to compute element tensor (from this point): 27360
2865
2866
// Loop quadrature points for integral
2867
// Number of operations to compute element tensor for following IP loop = 27360
2868
for (unsigned int ip = 0; ip < 2; ip++)
2869
{
2870
2871
// Number of operations for primary indices: 13680
2872
for (unsigned int j = 0; j < 6; j++)
2873
{
2874
for (unsigned int k = 0; k < 6; k++)
2875
{
2876
// Number of operations to compute entry: 95
2877
A[(j + 6)*12 + k] += (((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*w[1][0]/(w[0][0]) + (((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*-1))*W2[ip]*det;
2878
// Number of operations to compute entry: 95
2879
A[(j + 6)*12 + (k + 6)] += (((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f1_D10[ip][j] + Jinv1_11*FE0_f1_D01[ip][j])*n11 + (Jinv1_00*FE0_f1_D10[ip][j] + Jinv1_10*FE0_f1_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*-1))*W2[ip]*det;
2880
// Number of operations to compute entry: 95
2881
A[j*12 + (k + 6)] += ((((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*-1) + ((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*((Jinv1_01*FE0_f1_D10[ip][k] + Jinv1_11*FE0_f1_D01[ip][k])*n11 + (Jinv1_00*FE0_f1_D10[ip][k] + Jinv1_10*FE0_f1_D01[ip][k])*n10)*w[1][0]/(w[0][0]))*W2[ip]*det;
2882
// Number of operations to compute entry: 95
2883
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*-1) + ((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2884
}// end loop over 'k'
2885
}// end loop over 'j'
2886
}// end loop over 'ip'
2887
}
2888
break;
2889
case 2:
2890
{
2891
// Total number of operations to compute element tensor (from this point): 27360
2892
2893
// Loop quadrature points for integral
2894
// Number of operations to compute element tensor for following IP loop = 27360
2895
for (unsigned int ip = 0; ip < 2; ip++)
2896
{
2897
2898
// Number of operations for primary indices: 13680
2899
for (unsigned int j = 0; j < 6; j++)
2900
{
2901
for (unsigned int k = 0; k < 6; k++)
2902
{
2903
// Number of operations to compute entry: 95
2904
A[(j + 6)*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*-1 + ((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*-1) + ((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2905
// Number of operations to compute entry: 95
2906
A[(j + 6)*12 + (k + 6)] += (((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv1_00*Jinv1_00*FE0_f0_D20[ip][j] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][j] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][j]) + (Jinv1_01*Jinv1_01*FE0_f0_D20[ip][j] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][j] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*-1 + ((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv1_01*FE0_f2_D10[ip][j] + Jinv1_11*FE0_f2_D01[ip][j])*n11 + (Jinv1_00*FE0_f2_D10[ip][j] + Jinv1_10*FE0_f2_D01[ip][j])*n10)*-1))*W2[ip]*det;
2907
// Number of operations to compute entry: 95
2908
A[j*12 + (k + 6)] += (((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*w[1][0]/(w[0][0]) + (((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv1_01*FE0_f2_D10[ip][k] + Jinv1_11*FE0_f2_D01[ip][k])*n11 + (Jinv1_00*FE0_f2_D10[ip][k] + Jinv1_10*FE0_f2_D01[ip][k])*n10)*-1 + ((Jinv1_01*Jinv1_01*FE0_f0_D20[ip][k] + Jinv1_01*Jinv1_11*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_01*FE0_f0_D11[ip][k] + Jinv1_11*Jinv1_11*FE0_f0_D02[ip][k]) + (Jinv1_00*Jinv1_00*FE0_f0_D20[ip][k] + Jinv1_00*Jinv1_10*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_00*FE0_f0_D11[ip][k] + Jinv1_10*Jinv1_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*-1))*W2[ip]*det;
2909
// Number of operations to compute entry: 95
2910
A[j*12 + k] += ((((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][k] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][k] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][k]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][k] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][k] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][k]))*0.5*((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*-1 + ((Jinv0_01*Jinv0_01*FE0_f0_D20[ip][j] + Jinv0_01*Jinv0_11*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_01*FE0_f0_D11[ip][j] + Jinv0_11*Jinv0_11*FE0_f0_D02[ip][j]) + (Jinv0_00*Jinv0_00*FE0_f0_D20[ip][j] + Jinv0_00*Jinv0_10*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_00*FE0_f0_D11[ip][j] + Jinv0_10*Jinv0_10*FE0_f0_D02[ip][j]))*0.5*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*-1) + ((Jinv0_01*FE0_f2_D10[ip][j] + Jinv0_11*FE0_f2_D01[ip][j])*n01 + (Jinv0_00*FE0_f2_D10[ip][j] + Jinv0_10*FE0_f2_D01[ip][j])*n00)*((Jinv0_01*FE0_f2_D10[ip][k] + Jinv0_11*FE0_f2_D01[ip][k])*n01 + (Jinv0_00*FE0_f2_D10[ip][k] + Jinv0_10*FE0_f2_D01[ip][k])*n00)*w[1][0]/(w[0][0]))*W2[ip]*det;
2911
}// end loop over 'k'
2912
}// end loop over 'j'
2913
}// end loop over 'ip'
2914
}
2915
break;
2916
}
2917
break;
2918
}
2919
}
2920
2921
};
2922
2923
/// This class defines the interface for the tabulation of the
2924
/// interior facet tensor corresponding to the local contribution to
2925
/// a form from the integral over an interior facet.
2926
2927
class biharmonic_0_interior_facet_integral_0: public ufc::interior_facet_integral
2928
{
2929
private:
2930
2931
biharmonic_0_interior_facet_integral_0_quadrature integral_0_quadrature;
2932
2933
public:
2934
2935
/// Constructor
2936
biharmonic_0_interior_facet_integral_0() : ufc::interior_facet_integral()
2937
{
2938
// Do nothing
2939
}
2940
2941
/// Destructor
2942
virtual ~biharmonic_0_interior_facet_integral_0()
2943
{
2944
// Do nothing
2945
}
2946
2947
/// Tabulate the tensor for the contribution from a local interior facet
2948
virtual void tabulate_tensor(double* A,
2949
const double * const * w,
2950
const ufc::cell& c0,
2951
const ufc::cell& c1,
2952
unsigned int facet0,
2953
unsigned int facet1) const
2954
{
2955
// Reset values of the element tensor block
2956
for (unsigned int j = 0; j < 144; j++)
2957
A[j] = 0;
2958
2959
// Add all contributions to element tensor
2960
integral_0_quadrature.tabulate_tensor(A, w, c0, c1, facet0, facet1);
2961
}
2962
2963
};
2964
2965
/// This class defines the interface for the assembly of the global
2966
/// tensor corresponding to a form with r + n arguments, that is, a
2967
/// mapping
2968
///
2969
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
2970
///
2971
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
2972
/// global tensor A is defined by
2973
///
2974
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
2975
///
2976
/// where each argument Vj represents the application to the
2977
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
2978
/// fixed functions (coefficients).
2979
2980
class biharmonic_form_0: public ufc::form
2981
{
2982
public:
2983
2984
/// Constructor
2985
biharmonic_form_0() : ufc::form()
2986
{
2987
// Do nothing
2988
}
2989
2990
/// Destructor
2991
virtual ~biharmonic_form_0()
2992
{
2993
// Do nothing
2994
}
2995
2996
/// Return a string identifying the form
2997
virtual const char* signature() const
2998
{
2999
return "Form([Integral(Product(IndexSum(Indexed(ComponentTensor(SpatialDerivative(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), IndexSum(Indexed(ComponentTensor(SpatialDerivative(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((Index(3),), {Index(3): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2}))), Measure('cell', 0, None)), Integral(Sum(Product(Division(PositiveRestricted(Constant(Cell('triangle', 1, Space(2)), 1)), PositiveRestricted(Constant(Cell('triangle', 1, Space(2)), 0))), Product(Sum(IndexSum(Product(Indexed(NegativeRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(4),), {Index(4): 2})), MultiIndex((Index(4),), {Index(4): 2}))), MultiIndex((Index(5),), {Index(5): 2})), Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(5),), {Index(5): 2}))), MultiIndex((Index(5),), {Index(5): 2})), IndexSum(Product(Indexed(PositiveRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(6),), {Index(6): 2})), MultiIndex((Index(6),), {Index(6): 2}))), MultiIndex((Index(7),), {Index(7): 2})), Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(7),), {Index(7): 2}))), MultiIndex((Index(7),), {Index(7): 2}))), Sum(IndexSum(Product(Indexed(NegativeRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(8),), {Index(8): 2})), MultiIndex((Index(8),), {Index(8): 2}))), MultiIndex((Index(9),), {Index(9): 2})), Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(9),), {Index(9): 2}))), MultiIndex((Index(9),), {Index(9): 2})), IndexSum(Product(Indexed(PositiveRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((Index(10),), {Index(10): 2}))), MultiIndex((Index(11),), {Index(11): 2})), Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(11),), {Index(11): 2}))), MultiIndex((Index(11),), {Index(11): 2}))))), Sum(Product(IntValue(-1, (), (), {}), Product(Product(FloatValue(0.5, (), (), {}), Sum(NegativeRestricted(IndexSum(Indexed(ComponentTensor(SpatialDerivative(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(12),), {Index(12): 2})), MultiIndex((Index(13),), {Index(13): 2})), MultiIndex((Index(13),), {Index(13): 2})), MultiIndex((Index(12),), {Index(12): 2})), MultiIndex((Index(12),), {Index(12): 2}))), PositiveRestricted(IndexSum(Indexed(ComponentTensor(SpatialDerivative(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(14),), {Index(14): 2})), MultiIndex((Index(15),), {Index(15): 2})), MultiIndex((Index(15),), {Index(15): 2})), MultiIndex((Index(14),), {Index(14): 2})), MultiIndex((Index(14),), {Index(14): 2}))))), Sum(IndexSum(Product(Indexed(NegativeRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(16),), {Index(16): 2})), MultiIndex((Index(16),), {Index(16): 2}))), MultiIndex((Index(17),), {Index(17): 2})), Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(17),), {Index(17): 2}))), MultiIndex((Index(17),), {Index(17): 2})), IndexSum(Product(Indexed(PositiveRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(18),), {Index(18): 2})), MultiIndex((Index(18),), {Index(18): 2}))), MultiIndex((Index(19),), {Index(19): 2})), Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(19),), {Index(19): 2}))), MultiIndex((Index(19),), {Index(19): 2}))))), Product(IntValue(-1, (), (), {}), Product(Product(FloatValue(0.5, (), (), {}), Sum(NegativeRestricted(IndexSum(Indexed(ComponentTensor(SpatialDerivative(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(20),), {Index(20): 2})), MultiIndex((Index(21),), {Index(21): 2})), MultiIndex((Index(21),), {Index(21): 2})), MultiIndex((Index(20),), {Index(20): 2})), MultiIndex((Index(20),), {Index(20): 2}))), PositiveRestricted(IndexSum(Indexed(ComponentTensor(SpatialDerivative(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 1), MultiIndex((Index(22),), {Index(22): 2})), MultiIndex((Index(23),), {Index(23): 2})), MultiIndex((Index(23),), {Index(23): 2})), MultiIndex((Index(22),), {Index(22): 2})), MultiIndex((Index(22),), {Index(22): 2}))))), Sum(IndexSum(Product(Indexed(NegativeRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(24),), {Index(24): 2})), MultiIndex((Index(24),), {Index(24): 2}))), MultiIndex((Index(25),), {Index(25): 2})), Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(25),), {Index(25): 2}))), MultiIndex((Index(25),), {Index(25): 2})), IndexSum(Product(Indexed(PositiveRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), MultiIndex((Index(26),), {Index(26): 2})), MultiIndex((Index(26),), {Index(26): 2}))), MultiIndex((Index(27),), {Index(27): 2})), Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(27),), {Index(27): 2}))), MultiIndex((Index(27),), {Index(27): 2}))))))), Measure('interior_facet', 0, None))])";
3000
}
3001
3002
/// Return the rank of the global tensor (r)
3003
virtual unsigned int rank() const
3004
{
3005
return 2;
3006
}
3007
3008
/// Return the number of coefficients (n)
3009
virtual unsigned int num_coefficients() const
3010
{
3011
return 2;
3012
}
3013
3014
/// Return the number of cell integrals
3015
virtual unsigned int num_cell_integrals() const
3016
{
3017
return 1;
3018
}
3019
3020
/// Return the number of exterior facet integrals
3021
virtual unsigned int num_exterior_facet_integrals() const
3022
{
3023
return 0;
3024
}
3025
3026
/// Return the number of interior facet integrals
3027
virtual unsigned int num_interior_facet_integrals() const
3028
{
3029
return 1;
3030
}
3031
3032
/// Create a new finite element for argument function i
3033
virtual ufc::finite_element* create_finite_element(unsigned int i) const
3034
{
3035
switch ( i )
3036
{
3037
case 0:
3038
return new biharmonic_0_finite_element_0();
3039
break;
3040
case 1:
3041
return new biharmonic_0_finite_element_1();
3042
break;
3043
case 2:
3044
return new biharmonic_0_finite_element_2();
3045
break;
3046
case 3:
3047
return new biharmonic_0_finite_element_3();
3048
break;
3049
}
3050
return 0;
3051
}
3052
3053
/// Create a new dof map for argument function i
3054
virtual ufc::dof_map* create_dof_map(unsigned int i) const
3055
{
3056
switch ( i )
3057
{
3058
case 0:
3059
return new biharmonic_0_dof_map_0();
3060
break;
3061
case 1:
3062
return new biharmonic_0_dof_map_1();
3063
break;
3064
case 2:
3065
return new biharmonic_0_dof_map_2();
3066
break;
3067
case 3:
3068
return new biharmonic_0_dof_map_3();
3069
break;
3070
}
3071
return 0;
3072
}
3073
3074
/// Create a new cell integral on sub domain i
3075
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
3076
{
3077
return new biharmonic_0_cell_integral_0();
3078
}
3079
3080
/// Create a new exterior facet integral on sub domain i
3081
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
3082
{
3083
return 0;
3084
}
3085
3086
/// Create a new interior facet integral on sub domain i
3087
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
3088
{
3089
return new biharmonic_0_interior_facet_integral_0();
3090
}
3091
3092
};
3093
3094
/// This class defines the interface for a finite element.
3095
3096
class biharmonic_1_finite_element_0: public ufc::finite_element
3097
{
3098
public:
3099
3100
/// Constructor
3101
biharmonic_1_finite_element_0() : ufc::finite_element()
3102
{
3103
// Do nothing
3104
}
3105
3106
/// Destructor
3107
virtual ~biharmonic_1_finite_element_0()
3108
{
3109
// Do nothing
3110
}
3111
3112
/// Return a string identifying the finite element
3113
virtual const char* signature() const
3114
{
3115
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
3116
}
3117
3118
/// Return the cell shape
3119
virtual ufc::shape cell_shape() const
3120
{
3121
return ufc::triangle;
3122
}
3123
3124
/// Return the dimension of the finite element function space
3125
virtual unsigned int space_dimension() const
3126
{
3127
return 6;
3128
}
3129
3130
/// Return the rank of the value space
3131
virtual unsigned int value_rank() const
3132
{
3133
return 0;
3134
}
3135
3136
/// Return the dimension of the value space for axis i
3137
virtual unsigned int value_dimension(unsigned int i) const
3138
{
3139
return 1;
3140
}
3141
3142
/// Evaluate basis function i at given point in cell
3143
virtual void evaluate_basis(unsigned int i,
3144
double* values,
3145
const double* coordinates,
3146
const ufc::cell& c) const
3147
{
3148
// Extract vertex coordinates
3149
const double * const * element_coordinates = c.coordinates;
3150
3151
// Compute Jacobian of affine map from reference cell
3152
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
3153
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
3154
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
3155
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
3156
3157
// Compute determinant of Jacobian
3158
const double detJ = J_00*J_11 - J_01*J_10;
3159
3160
// Compute inverse of Jacobian
3161
3162
// Get coordinates and map to the reference (UFC) element
3163
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
3164
element_coordinates[0][0]*element_coordinates[2][1] +\
3165
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
3166
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
3167
element_coordinates[1][0]*element_coordinates[0][1] -\
3168
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
3169
3170
// Map coordinates to the reference square
3171
if (std::abs(y - 1.0) < 1e-08)
3172
x = -1.0;
3173
else
3174
x = 2.0 *x/(1.0 - y) - 1.0;
3175
y = 2.0*y - 1.0;
3176
3177
// Reset values
3178
*values = 0;
3179
3180
// Map degree of freedom to element degree of freedom
3181
const unsigned int dof = i;
3182
3183
// Generate scalings
3184
const double scalings_y_0 = 1;
3185
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
3186
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
3187
3188
// Compute psitilde_a
3189
const double psitilde_a_0 = 1;
3190
const double psitilde_a_1 = x;
3191
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
3192
3193
// Compute psitilde_bs
3194
const double psitilde_bs_0_0 = 1;
3195
const double psitilde_bs_0_1 = 1.5*y + 0.5;
3196
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
3197
const double psitilde_bs_1_0 = 1;
3198
const double psitilde_bs_1_1 = 2.5*y + 1.5;
3199
const double psitilde_bs_2_0 = 1;
3200
3201
// Compute basisvalues
3202
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
3203
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
3204
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
3205
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
3206
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
3207
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
3208
3209
// Table(s) of coefficients
3210
static const double coefficients0[6][6] = \
3211
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
3212
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
3213
{0, 0, 0.2, 0, 0, 0.163299316},
3214
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
3215
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
3216
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
3217
3218
// Extract relevant coefficients
3219
const double coeff0_0 = coefficients0[dof][0];
3220
const double coeff0_1 = coefficients0[dof][1];
3221
const double coeff0_2 = coefficients0[dof][2];
3222
const double coeff0_3 = coefficients0[dof][3];
3223
const double coeff0_4 = coefficients0[dof][4];
3224
const double coeff0_5 = coefficients0[dof][5];
3225
3226
// Compute value(s)
3227
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
3228
}
3229
3230
/// Evaluate all basis functions at given point in cell
3231
virtual void evaluate_basis_all(double* values,
3232
const double* coordinates,
3233
const ufc::cell& c) const
3234
{
3235
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
3236
}
3237
3238
/// Evaluate order n derivatives of basis function i at given point in cell
3239
virtual void evaluate_basis_derivatives(unsigned int i,
3240
unsigned int n,
3241
double* values,
3242
const double* coordinates,
3243
const ufc::cell& c) const
3244
{
3245
// Extract vertex coordinates
3246
const double * const * element_coordinates = c.coordinates;
3247
3248
// Compute Jacobian of affine map from reference cell
3249
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
3250
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
3251
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
3252
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
3253
3254
// Compute determinant of Jacobian
3255
const double detJ = J_00*J_11 - J_01*J_10;
3256
3257
// Compute inverse of Jacobian
3258
3259
// Get coordinates and map to the reference (UFC) element
3260
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
3261
element_coordinates[0][0]*element_coordinates[2][1] +\
3262
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
3263
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
3264
element_coordinates[1][0]*element_coordinates[0][1] -\
3265
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
3266
3267
// Map coordinates to the reference square
3268
if (std::abs(y - 1.0) < 1e-08)
3269
x = -1.0;
3270
else
3271
x = 2.0 *x/(1.0 - y) - 1.0;
3272
y = 2.0*y - 1.0;
3273
3274
// Compute number of derivatives
3275
unsigned int num_derivatives = 1;
3276
3277
for (unsigned int j = 0; j < n; j++)
3278
num_derivatives *= 2;
3279
3280
3281
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
3282
unsigned int **combinations = new unsigned int *[num_derivatives];
3283
3284
for (unsigned int j = 0; j < num_derivatives; j++)
3285
{
3286
combinations[j] = new unsigned int [n];
3287
for (unsigned int k = 0; k < n; k++)
3288
combinations[j][k] = 0;
3289
}
3290
3291
// Generate combinations of derivatives
3292
for (unsigned int row = 1; row < num_derivatives; row++)
3293
{
3294
for (unsigned int num = 0; num < row; num++)
3295
{
3296
for (unsigned int col = n-1; col+1 > 0; col--)
3297
{
3298
if (combinations[row][col] + 1 > 1)
3299
combinations[row][col] = 0;
3300
else
3301
{
3302
combinations[row][col] += 1;
3303
break;
3304
}
3305
}
3306
}
3307
}
3308
3309
// Compute inverse of Jacobian
3310
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
3311
3312
// Declare transformation matrix
3313
// Declare pointer to two dimensional array and initialise
3314
double **transform = new double *[num_derivatives];
3315
3316
for (unsigned int j = 0; j < num_derivatives; j++)
3317
{
3318
transform[j] = new double [num_derivatives];
3319
for (unsigned int k = 0; k < num_derivatives; k++)
3320
transform[j][k] = 1;
3321
}
3322
3323
// Construct transformation matrix
3324
for (unsigned int row = 0; row < num_derivatives; row++)
3325
{
3326
for (unsigned int col = 0; col < num_derivatives; col++)
3327
{
3328
for (unsigned int k = 0; k < n; k++)
3329
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
3330
}
3331
}
3332
3333
// Reset values
3334
for (unsigned int j = 0; j < 1*num_derivatives; j++)
3335
values[j] = 0;
3336
3337
// Map degree of freedom to element degree of freedom
3338
const unsigned int dof = i;
3339
3340
// Generate scalings
3341
const double scalings_y_0 = 1;
3342
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
3343
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
3344
3345
// Compute psitilde_a
3346
const double psitilde_a_0 = 1;
3347
const double psitilde_a_1 = x;
3348
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
3349
3350
// Compute psitilde_bs
3351
const double psitilde_bs_0_0 = 1;
3352
const double psitilde_bs_0_1 = 1.5*y + 0.5;
3353
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
3354
const double psitilde_bs_1_0 = 1;
3355
const double psitilde_bs_1_1 = 2.5*y + 1.5;
3356
const double psitilde_bs_2_0 = 1;
3357
3358
// Compute basisvalues
3359
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
3360
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
3361
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
3362
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
3363
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
3364
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
3365
3366
// Table(s) of coefficients
3367
static const double coefficients0[6][6] = \
3368
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
3369
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
3370
{0, 0, 0.2, 0, 0, 0.163299316},
3371
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
3372
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
3373
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
3374
3375
// Interesting (new) part
3376
// Tables of derivatives of the polynomial base (transpose)
3377
static const double dmats0[6][6] = \
3378
{{0, 0, 0, 0, 0, 0},
3379
{4.89897949, 0, 0, 0, 0, 0},
3380
{0, 0, 0, 0, 0, 0},
3381
{0, 9.48683298, 0, 0, 0, 0},
3382
{4, 0, 7.07106781, 0, 0, 0},
3383
{0, 0, 0, 0, 0, 0}};
3384
3385
static const double dmats1[6][6] = \
3386
{{0, 0, 0, 0, 0, 0},
3387
{2.44948974, 0, 0, 0, 0, 0},
3388
{4.24264069, 0, 0, 0, 0, 0},
3389
{2.5819889, 4.74341649, -0.912870929, 0, 0, 0},
3390
{2, 6.12372436, 3.53553391, 0, 0, 0},
3391
{-2.30940108, 0, 8.16496581, 0, 0, 0}};
3392
3393
// Compute reference derivatives
3394
// Declare pointer to array of derivatives on FIAT element
3395
double *derivatives = new double [num_derivatives];
3396
3397
// Declare coefficients
3398
double coeff0_0 = 0;
3399
double coeff0_1 = 0;
3400
double coeff0_2 = 0;
3401
double coeff0_3 = 0;
3402
double coeff0_4 = 0;
3403
double coeff0_5 = 0;
3404
3405
// Declare new coefficients
3406
double new_coeff0_0 = 0;
3407
double new_coeff0_1 = 0;
3408
double new_coeff0_2 = 0;
3409
double new_coeff0_3 = 0;
3410
double new_coeff0_4 = 0;
3411
double new_coeff0_5 = 0;
3412
3413
// Loop possible derivatives
3414
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
3415
{
3416
// Get values from coefficients array
3417
new_coeff0_0 = coefficients0[dof][0];
3418
new_coeff0_1 = coefficients0[dof][1];
3419
new_coeff0_2 = coefficients0[dof][2];
3420
new_coeff0_3 = coefficients0[dof][3];
3421
new_coeff0_4 = coefficients0[dof][4];
3422
new_coeff0_5 = coefficients0[dof][5];
3423
3424
// Loop derivative order
3425
for (unsigned int j = 0; j < n; j++)
3426
{
3427
// Update old coefficients
3428
coeff0_0 = new_coeff0_0;
3429
coeff0_1 = new_coeff0_1;
3430
coeff0_2 = new_coeff0_2;
3431
coeff0_3 = new_coeff0_3;
3432
coeff0_4 = new_coeff0_4;
3433
coeff0_5 = new_coeff0_5;
3434
3435
if(combinations[deriv_num][j] == 0)
3436
{
3437
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
3438
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
3439
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
3440
new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
3441
new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
3442
new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
3443
}
3444
if(combinations[deriv_num][j] == 1)
3445
{
3446
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
3447
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
3448
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
3449
new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
3450
new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
3451
new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
3452
}
3453
3454
}
3455
// Compute derivatives on reference element as dot product of coefficients and basisvalues
3456
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
3457
}
3458
3459
// Transform derivatives back to physical element
3460
for (unsigned int row = 0; row < num_derivatives; row++)
3461
{
3462
for (unsigned int col = 0; col < num_derivatives; col++)
3463
{
3464
values[row] += transform[row][col]*derivatives[col];
3465
}
3466
}
3467
// Delete pointer to array of derivatives on FIAT element
3468
delete [] derivatives;
3469
3470
// Delete pointer to array of combinations of derivatives and transform
3471
for (unsigned int row = 0; row < num_derivatives; row++)
3472
{
3473
delete [] combinations[row];
3474
delete [] transform[row];
3475
}
3476
3477
delete [] combinations;
3478
delete [] transform;
3479
}
3480
3481
/// Evaluate order n derivatives of all basis functions at given point in cell
3482
virtual void evaluate_basis_derivatives_all(unsigned int n,
3483
double* values,
3484
const double* coordinates,
3485
const ufc::cell& c) const
3486
{
3487
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
3488
}
3489
3490
/// Evaluate linear functional for dof i on the function f
3491
virtual double evaluate_dof(unsigned int i,
3492
const ufc::function& f,
3493
const ufc::cell& c) const
3494
{
3495
// The reference points, direction and weights:
3496
static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
3497
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
3498
static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
3499
3500
const double * const * x = c.coordinates;
3501
double result = 0.0;
3502
// Iterate over the points:
3503
// Evaluate basis functions for affine mapping
3504
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
3505
const double w1 = X[i][0][0];
3506
const double w2 = X[i][0][1];
3507
3508
// Compute affine mapping y = F(X)
3509
double y[2];
3510
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
3511
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
3512
3513
// Evaluate function at physical points
3514
double values[1];
3515
f.evaluate(values, y, c);
3516
3517
// Map function values using appropriate mapping
3518
// Affine map: Do nothing
3519
3520
// Note that we do not map the weights (yet).
3521
3522
// Take directional components
3523
for(int k = 0; k < 1; k++)
3524
result += values[k]*D[i][0][k];
3525
// Multiply by weights
3526
result *= W[i][0];
3527
3528
return result;
3529
}
3530
3531
/// Evaluate linear functionals for all dofs on the function f
3532
virtual void evaluate_dofs(double* values,
3533
const ufc::function& f,
3534
const ufc::cell& c) const
3535
{
3536
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
3537
}
3538
3539
/// Interpolate vertex values from dof values
3540
virtual void interpolate_vertex_values(double* vertex_values,
3541
const double* dof_values,
3542
const ufc::cell& c) const
3543
{
3544
// Evaluate at vertices and use affine mapping
3545
vertex_values[0] = dof_values[0];
3546
vertex_values[1] = dof_values[1];
3547
vertex_values[2] = dof_values[2];
3548
}
3549
3550
/// Return the number of sub elements (for a mixed element)
3551
virtual unsigned int num_sub_elements() const
3552
{
3553
return 1;
3554
}
3555
3556
/// Create a new finite element for sub element i (for a mixed element)
3557
virtual ufc::finite_element* create_sub_element(unsigned int i) const
3558
{
3559
return new biharmonic_1_finite_element_0();
3560
}
3561
3562
};
3563
3564
/// This class defines the interface for a finite element.
3565
3566
class biharmonic_1_finite_element_1: public ufc::finite_element
3567
{
3568
public:
3569
3570
/// Constructor
3571
biharmonic_1_finite_element_1() : ufc::finite_element()
3572
{
3573
// Do nothing
3574
}
3575
3576
/// Destructor
3577
virtual ~biharmonic_1_finite_element_1()
3578
{
3579
// Do nothing
3580
}
3581
3582
/// Return a string identifying the finite element
3583
virtual const char* signature() const
3584
{
3585
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
3586
}
3587
3588
/// Return the cell shape
3589
virtual ufc::shape cell_shape() const
3590
{
3591
return ufc::triangle;
3592
}
3593
3594
/// Return the dimension of the finite element function space
3595
virtual unsigned int space_dimension() const
3596
{
3597
return 6;
3598
}
3599
3600
/// Return the rank of the value space
3601
virtual unsigned int value_rank() const
3602
{
3603
return 0;
3604
}
3605
3606
/// Return the dimension of the value space for axis i
3607
virtual unsigned int value_dimension(unsigned int i) const
3608
{
3609
return 1;
3610
}
3611
3612
/// Evaluate basis function i at given point in cell
3613
virtual void evaluate_basis(unsigned int i,
3614
double* values,
3615
const double* coordinates,
3616
const ufc::cell& c) const
3617
{
3618
// Extract vertex coordinates
3619
const double * const * element_coordinates = c.coordinates;
3620
3621
// Compute Jacobian of affine map from reference cell
3622
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
3623
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
3624
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
3625
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
3626
3627
// Compute determinant of Jacobian
3628
const double detJ = J_00*J_11 - J_01*J_10;
3629
3630
// Compute inverse of Jacobian
3631
3632
// Get coordinates and map to the reference (UFC) element
3633
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
3634
element_coordinates[0][0]*element_coordinates[2][1] +\
3635
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
3636
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
3637
element_coordinates[1][0]*element_coordinates[0][1] -\
3638
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
3639
3640
// Map coordinates to the reference square
3641
if (std::abs(y - 1.0) < 1e-08)
3642
x = -1.0;
3643
else
3644
x = 2.0 *x/(1.0 - y) - 1.0;
3645
y = 2.0*y - 1.0;
3646
3647
// Reset values
3648
*values = 0;
3649
3650
// Map degree of freedom to element degree of freedom
3651
const unsigned int dof = i;
3652
3653
// Generate scalings
3654
const double scalings_y_0 = 1;
3655
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
3656
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
3657
3658
// Compute psitilde_a
3659
const double psitilde_a_0 = 1;
3660
const double psitilde_a_1 = x;
3661
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
3662
3663
// Compute psitilde_bs
3664
const double psitilde_bs_0_0 = 1;
3665
const double psitilde_bs_0_1 = 1.5*y + 0.5;
3666
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
3667
const double psitilde_bs_1_0 = 1;
3668
const double psitilde_bs_1_1 = 2.5*y + 1.5;
3669
const double psitilde_bs_2_0 = 1;
3670
3671
// Compute basisvalues
3672
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
3673
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
3674
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
3675
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
3676
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
3677
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
3678
3679
// Table(s) of coefficients
3680
static const double coefficients0[6][6] = \
3681
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
3682
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
3683
{0, 0, 0.2, 0, 0, 0.163299316},
3684
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
3685
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
3686
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
3687
3688
// Extract relevant coefficients
3689
const double coeff0_0 = coefficients0[dof][0];
3690
const double coeff0_1 = coefficients0[dof][1];
3691
const double coeff0_2 = coefficients0[dof][2];
3692
const double coeff0_3 = coefficients0[dof][3];
3693
const double coeff0_4 = coefficients0[dof][4];
3694
const double coeff0_5 = coefficients0[dof][5];
3695
3696
// Compute value(s)
3697
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;
3698
}
3699
3700
/// Evaluate all basis functions at given point in cell
3701
virtual void evaluate_basis_all(double* values,
3702
const double* coordinates,
3703
const ufc::cell& c) const
3704
{
3705
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
3706
}
3707
3708
/// Evaluate order n derivatives of basis function i at given point in cell
3709
virtual void evaluate_basis_derivatives(unsigned int i,
3710
unsigned int n,
3711
double* values,
3712
const double* coordinates,
3713
const ufc::cell& c) const
3714
{
3715
// Extract vertex coordinates
3716
const double * const * element_coordinates = c.coordinates;
3717
3718
// Compute Jacobian of affine map from reference cell
3719
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
3720
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
3721
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
3722
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
3723
3724
// Compute determinant of Jacobian
3725
const double detJ = J_00*J_11 - J_01*J_10;
3726
3727
// Compute inverse of Jacobian
3728
3729
// Get coordinates and map to the reference (UFC) element
3730
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
3731
element_coordinates[0][0]*element_coordinates[2][1] +\
3732
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
3733
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
3734
element_coordinates[1][0]*element_coordinates[0][1] -\
3735
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
3736
3737
// Map coordinates to the reference square
3738
if (std::abs(y - 1.0) < 1e-08)
3739
x = -1.0;
3740
else
3741
x = 2.0 *x/(1.0 - y) - 1.0;
3742
y = 2.0*y - 1.0;
3743
3744
// Compute number of derivatives
3745
unsigned int num_derivatives = 1;
3746
3747
for (unsigned int j = 0; j < n; j++)
3748
num_derivatives *= 2;
3749
3750
3751
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
3752
unsigned int **combinations = new unsigned int *[num_derivatives];
3753
3754
for (unsigned int j = 0; j < num_derivatives; j++)
3755
{
3756
combinations[j] = new unsigned int [n];
3757
for (unsigned int k = 0; k < n; k++)
3758
combinations[j][k] = 0;
3759
}
3760
3761
// Generate combinations of derivatives
3762
for (unsigned int row = 1; row < num_derivatives; row++)
3763
{
3764
for (unsigned int num = 0; num < row; num++)
3765
{
3766
for (unsigned int col = n-1; col+1 > 0; col--)
3767
{
3768
if (combinations[row][col] + 1 > 1)
3769
combinations[row][col] = 0;
3770
else
3771
{
3772
combinations[row][col] += 1;
3773
break;
3774
}
3775
}
3776
}
3777
}
3778
3779
// Compute inverse of Jacobian
3780
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
3781
3782
// Declare transformation matrix
3783
// Declare pointer to two dimensional array and initialise
3784
double **transform = new double *[num_derivatives];
3785
3786
for (unsigned int j = 0; j < num_derivatives; j++)
3787
{
3788
transform[j] = new double [num_derivatives];
3789
for (unsigned int k = 0; k < num_derivatives; k++)
3790
transform[j][k] = 1;
3791
}
3792
3793
// Construct transformation matrix
3794
for (unsigned int row = 0; row < num_derivatives; row++)
3795
{
3796
for (unsigned int col = 0; col < num_derivatives; col++)
3797
{
3798
for (unsigned int k = 0; k < n; k++)
3799
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
3800
}
3801
}
3802
3803
// Reset values
3804
for (unsigned int j = 0; j < 1*num_derivatives; j++)
3805
values[j] = 0;
3806
3807
// Map degree of freedom to element degree of freedom
3808
const unsigned int dof = i;
3809
3810
// Generate scalings
3811
const double scalings_y_0 = 1;
3812
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
3813
const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
3814
3815
// Compute psitilde_a
3816
const double psitilde_a_0 = 1;
3817
const double psitilde_a_1 = x;
3818
const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
3819
3820
// Compute psitilde_bs
3821
const double psitilde_bs_0_0 = 1;
3822
const double psitilde_bs_0_1 = 1.5*y + 0.5;
3823
const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
3824
const double psitilde_bs_1_0 = 1;
3825
const double psitilde_bs_1_1 = 2.5*y + 1.5;
3826
const double psitilde_bs_2_0 = 1;
3827
3828
// Compute basisvalues
3829
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
3830
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
3831
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
3832
const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
3833
const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
3834
const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
3835
3836
// Table(s) of coefficients
3837
static const double coefficients0[6][6] = \
3838
{{0, -0.173205081, -0.1, 0.121716124, 0.0942809042, 0.0544331054},
3839
{0, 0.173205081, -0.1, 0.121716124, -0.0942809042, 0.0544331054},
3840
{0, 0, 0.2, 0, 0, 0.163299316},
3841
{0.471404521, 0.230940108, 0.133333333, 0, 0.188561808, -0.163299316},
3842
{0.471404521, -0.230940108, 0.133333333, 0, -0.188561808, -0.163299316},
3843
{0.471404521, 0, -0.266666667, -0.243432248, 0, 0.0544331054}};
3844
3845
// Interesting (new) part
3846
// Tables of derivatives of the polynomial base (transpose)
3847
static const double dmats0[6][6] = \
3848
{{0, 0, 0, 0, 0, 0},
3849
{4.89897949, 0, 0, 0, 0, 0},
3850
{0, 0, 0, 0, 0, 0},
3851
{0, 9.48683298, 0, 0, 0, 0},
3852
{4, 0, 7.07106781, 0, 0, 0},
3853
{0, 0, 0, 0, 0, 0}};
3854
3855
static const double dmats1[6][6] = \
3856
{{0, 0, 0, 0, 0, 0},
3857
{2.44948974, 0, 0, 0, 0, 0},
3858
{4.24264069, 0, 0, 0, 0, 0},
3859
{2.5819889, 4.74341649, -0.912870929, 0, 0, 0},
3860
{2, 6.12372436, 3.53553391, 0, 0, 0},
3861
{-2.30940108, 0, 8.16496581, 0, 0, 0}};
3862
3863
// Compute reference derivatives
3864
// Declare pointer to array of derivatives on FIAT element
3865
double *derivatives = new double [num_derivatives];
3866
3867
// Declare coefficients
3868
double coeff0_0 = 0;
3869
double coeff0_1 = 0;
3870
double coeff0_2 = 0;
3871
double coeff0_3 = 0;
3872
double coeff0_4 = 0;
3873
double coeff0_5 = 0;
3874
3875
// Declare new coefficients
3876
double new_coeff0_0 = 0;
3877
double new_coeff0_1 = 0;
3878
double new_coeff0_2 = 0;
3879
double new_coeff0_3 = 0;
3880
double new_coeff0_4 = 0;
3881
double new_coeff0_5 = 0;
3882
3883
// Loop possible derivatives
3884
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
3885
{
3886
// Get values from coefficients array
3887
new_coeff0_0 = coefficients0[dof][0];
3888
new_coeff0_1 = coefficients0[dof][1];
3889
new_coeff0_2 = coefficients0[dof][2];
3890
new_coeff0_3 = coefficients0[dof][3];
3891
new_coeff0_4 = coefficients0[dof][4];
3892
new_coeff0_5 = coefficients0[dof][5];
3893
3894
// Loop derivative order
3895
for (unsigned int j = 0; j < n; j++)
3896
{
3897
// Update old coefficients
3898
coeff0_0 = new_coeff0_0;
3899
coeff0_1 = new_coeff0_1;
3900
coeff0_2 = new_coeff0_2;
3901
coeff0_3 = new_coeff0_3;
3902
coeff0_4 = new_coeff0_4;
3903
coeff0_5 = new_coeff0_5;
3904
3905
if(combinations[deriv_num][j] == 0)
3906
{
3907
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];
3908
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];
3909
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];
3910
new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];
3911
new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];
3912
new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];
3913
}
3914
if(combinations[deriv_num][j] == 1)
3915
{
3916
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];
3917
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];
3918
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];
3919
new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];
3920
new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];
3921
new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];
3922
}
3923
3924
}
3925
// Compute derivatives on reference element as dot product of coefficients and basisvalues
3926
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;
3927
}
3928
3929
// Transform derivatives back to physical element
3930
for (unsigned int row = 0; row < num_derivatives; row++)
3931
{
3932
for (unsigned int col = 0; col < num_derivatives; col++)
3933
{
3934
values[row] += transform[row][col]*derivatives[col];
3935
}
3936
}
3937
// Delete pointer to array of derivatives on FIAT element
3938
delete [] derivatives;
3939
3940
// Delete pointer to array of combinations of derivatives and transform
3941
for (unsigned int row = 0; row < num_derivatives; row++)
3942
{
3943
delete [] combinations[row];
3944
delete [] transform[row];
3945
}
3946
3947
delete [] combinations;
3948
delete [] transform;
3949
}
3950
3951
/// Evaluate order n derivatives of all basis functions at given point in cell
3952
virtual void evaluate_basis_derivatives_all(unsigned int n,
3953
double* values,
3954
const double* coordinates,
3955
const ufc::cell& c) const
3956
{
3957
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
3958
}
3959
3960
/// Evaluate linear functional for dof i on the function f
3961
virtual double evaluate_dof(unsigned int i,
3962
const ufc::function& f,
3963
const ufc::cell& c) const
3964
{
3965
// The reference points, direction and weights:
3966
static const double X[6][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.5, 0.5}}, {{0, 0.5}}, {{0.5, 0}}};
3967
static const double W[6][1] = {{1}, {1}, {1}, {1}, {1}, {1}};
3968
static const double D[6][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
3969
3970
const double * const * x = c.coordinates;
3971
double result = 0.0;
3972
// Iterate over the points:
3973
// Evaluate basis functions for affine mapping
3974
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
3975
const double w1 = X[i][0][0];
3976
const double w2 = X[i][0][1];
3977
3978
// Compute affine mapping y = F(X)
3979
double y[2];
3980
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
3981
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
3982
3983
// Evaluate function at physical points
3984
double values[1];
3985
f.evaluate(values, y, c);
3986
3987
// Map function values using appropriate mapping
3988
// Affine map: Do nothing
3989
3990
// Note that we do not map the weights (yet).
3991
3992
// Take directional components
3993
for(int k = 0; k < 1; k++)
3994
result += values[k]*D[i][0][k];
3995
// Multiply by weights
3996
result *= W[i][0];
3997
3998
return result;
3999
}
4000
4001
/// Evaluate linear functionals for all dofs on the function f
4002
virtual void evaluate_dofs(double* values,
4003
const ufc::function& f,
4004
const ufc::cell& c) const
4005
{
4006
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4007
}
4008
4009
/// Interpolate vertex values from dof values
4010
virtual void interpolate_vertex_values(double* vertex_values,
4011
const double* dof_values,
4012
const ufc::cell& c) const
4013
{
4014
// Evaluate at vertices and use affine mapping
4015
vertex_values[0] = dof_values[0];
4016
vertex_values[1] = dof_values[1];
4017
vertex_values[2] = dof_values[2];
4018
}
4019
4020
/// Return the number of sub elements (for a mixed element)
4021
virtual unsigned int num_sub_elements() const
4022
{
4023
return 1;
4024
}
4025
4026
/// Create a new finite element for sub element i (for a mixed element)
4027
virtual ufc::finite_element* create_sub_element(unsigned int i) const
4028
{
4029
return new biharmonic_1_finite_element_1();
4030
}
4031
4032
};
4033
4034
/// This class defines the interface for a local-to-global mapping of
4035
/// degrees of freedom (dofs).
4036
4037
class biharmonic_1_dof_map_0: public ufc::dof_map
4038
{
4039
private:
4040
4041
unsigned int __global_dimension;
4042
4043
public:
4044
4045
/// Constructor
4046
biharmonic_1_dof_map_0() : ufc::dof_map()
4047
{
4048
__global_dimension = 0;
4049
}
4050
4051
/// Destructor
4052
virtual ~biharmonic_1_dof_map_0()
4053
{
4054
// Do nothing
4055
}
4056
4057
/// Return a string identifying the dof map
4058
virtual const char* signature() const
4059
{
4060
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
4061
}
4062
4063
/// Return true iff mesh entities of topological dimension d are needed
4064
virtual bool needs_mesh_entities(unsigned int d) const
4065
{
4066
switch ( d )
4067
{
4068
case 0:
4069
return true;
4070
break;
4071
case 1:
4072
return true;
4073
break;
4074
case 2:
4075
return false;
4076
break;
4077
}
4078
return false;
4079
}
4080
4081
/// Initialize dof map for mesh (return true iff init_cell() is needed)
4082
virtual bool init_mesh(const ufc::mesh& m)
4083
{
4084
__global_dimension = m.num_entities[0] + m.num_entities[1];
4085
return false;
4086
}
4087
4088
/// Initialize dof map for given cell
4089
virtual void init_cell(const ufc::mesh& m,
4090
const ufc::cell& c)
4091
{
4092
// Do nothing
4093
}
4094
4095
/// Finish initialization of dof map for cells
4096
virtual void init_cell_finalize()
4097
{
4098
// Do nothing
4099
}
4100
4101
/// Return the dimension of the global finite element function space
4102
virtual unsigned int global_dimension() const
4103
{
4104
return __global_dimension;
4105
}
4106
4107
/// Return the dimension of the local finite element function space for a cell
4108
virtual unsigned int local_dimension(const ufc::cell& c) const
4109
{
4110
return 6;
4111
}
4112
4113
/// Return the maximum dimension of the local finite element function space
4114
virtual unsigned int max_local_dimension() const
4115
{
4116
return 6;
4117
}
4118
4119
// Return the geometric dimension of the coordinates this dof map provides
4120
virtual unsigned int geometric_dimension() const
4121
{
4122
return 2;
4123
}
4124
4125
/// Return the number of dofs on each cell facet
4126
virtual unsigned int num_facet_dofs() const
4127
{
4128
return 3;
4129
}
4130
4131
/// Return the number of dofs associated with each cell entity of dimension d
4132
virtual unsigned int num_entity_dofs(unsigned int d) const
4133
{
4134
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4135
}
4136
4137
/// Tabulate the local-to-global mapping of dofs on a cell
4138
virtual void tabulate_dofs(unsigned int* dofs,
4139
const ufc::mesh& m,
4140
const ufc::cell& c) const
4141
{
4142
dofs[0] = c.entity_indices[0][0];
4143
dofs[1] = c.entity_indices[0][1];
4144
dofs[2] = c.entity_indices[0][2];
4145
unsigned int offset = m.num_entities[0];
4146
dofs[3] = offset + c.entity_indices[1][0];
4147
dofs[4] = offset + c.entity_indices[1][1];
4148
dofs[5] = offset + c.entity_indices[1][2];
4149
}
4150
4151
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4152
virtual void tabulate_facet_dofs(unsigned int* dofs,
4153
unsigned int facet) const
4154
{
4155
switch ( facet )
4156
{
4157
case 0:
4158
dofs[0] = 1;
4159
dofs[1] = 2;
4160
dofs[2] = 3;
4161
break;
4162
case 1:
4163
dofs[0] = 0;
4164
dofs[1] = 2;
4165
dofs[2] = 4;
4166
break;
4167
case 2:
4168
dofs[0] = 0;
4169
dofs[1] = 1;
4170
dofs[2] = 5;
4171
break;
4172
}
4173
}
4174
4175
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4176
virtual void tabulate_entity_dofs(unsigned int* dofs,
4177
unsigned int d, unsigned int i) const
4178
{
4179
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4180
}
4181
4182
/// Tabulate the coordinates of all dofs on a cell
4183
virtual void tabulate_coordinates(double** coordinates,
4184
const ufc::cell& c) const
4185
{
4186
const double * const * x = c.coordinates;
4187
coordinates[0][0] = x[0][0];
4188
coordinates[0][1] = x[0][1];
4189
coordinates[1][0] = x[1][0];
4190
coordinates[1][1] = x[1][1];
4191
coordinates[2][0] = x[2][0];
4192
coordinates[2][1] = x[2][1];
4193
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
4194
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
4195
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
4196
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
4197
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
4198
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
4199
}
4200
4201
/// Return the number of sub dof maps (for a mixed element)
4202
virtual unsigned int num_sub_dof_maps() const
4203
{
4204
return 1;
4205
}
4206
4207
/// Create a new dof_map for sub dof map i (for a mixed element)
4208
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
4209
{
4210
return new biharmonic_1_dof_map_0();
4211
}
4212
4213
};
4214
4215
/// This class defines the interface for a local-to-global mapping of
4216
/// degrees of freedom (dofs).
4217
4218
class biharmonic_1_dof_map_1: public ufc::dof_map
4219
{
4220
private:
4221
4222
unsigned int __global_dimension;
4223
4224
public:
4225
4226
/// Constructor
4227
biharmonic_1_dof_map_1() : ufc::dof_map()
4228
{
4229
__global_dimension = 0;
4230
}
4231
4232
/// Destructor
4233
virtual ~biharmonic_1_dof_map_1()
4234
{
4235
// Do nothing
4236
}
4237
4238
/// Return a string identifying the dof map
4239
virtual const char* signature() const
4240
{
4241
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2)";
4242
}
4243
4244
/// Return true iff mesh entities of topological dimension d are needed
4245
virtual bool needs_mesh_entities(unsigned int d) const
4246
{
4247
switch ( d )
4248
{
4249
case 0:
4250
return true;
4251
break;
4252
case 1:
4253
return true;
4254
break;
4255
case 2:
4256
return false;
4257
break;
4258
}
4259
return false;
4260
}
4261
4262
/// Initialize dof map for mesh (return true iff init_cell() is needed)
4263
virtual bool init_mesh(const ufc::mesh& m)
4264
{
4265
__global_dimension = m.num_entities[0] + m.num_entities[1];
4266
return false;
4267
}
4268
4269
/// Initialize dof map for given cell
4270
virtual void init_cell(const ufc::mesh& m,
4271
const ufc::cell& c)
4272
{
4273
// Do nothing
4274
}
4275
4276
/// Finish initialization of dof map for cells
4277
virtual void init_cell_finalize()
4278
{
4279
// Do nothing
4280
}
4281
4282
/// Return the dimension of the global finite element function space
4283
virtual unsigned int global_dimension() const
4284
{
4285
return __global_dimension;
4286
}
4287
4288
/// Return the dimension of the local finite element function space for a cell
4289
virtual unsigned int local_dimension(const ufc::cell& c) const
4290
{
4291
return 6;
4292
}
4293
4294
/// Return the maximum dimension of the local finite element function space
4295
virtual unsigned int max_local_dimension() const
4296
{
4297
return 6;
4298
}
4299
4300
// Return the geometric dimension of the coordinates this dof map provides
4301
virtual unsigned int geometric_dimension() const
4302
{
4303
return 2;
4304
}
4305
4306
/// Return the number of dofs on each cell facet
4307
virtual unsigned int num_facet_dofs() const
4308
{
4309
return 3;
4310
}
4311
4312
/// Return the number of dofs associated with each cell entity of dimension d
4313
virtual unsigned int num_entity_dofs(unsigned int d) const
4314
{
4315
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4316
}
4317
4318
/// Tabulate the local-to-global mapping of dofs on a cell
4319
virtual void tabulate_dofs(unsigned int* dofs,
4320
const ufc::mesh& m,
4321
const ufc::cell& c) const
4322
{
4323
dofs[0] = c.entity_indices[0][0];
4324
dofs[1] = c.entity_indices[0][1];
4325
dofs[2] = c.entity_indices[0][2];
4326
unsigned int offset = m.num_entities[0];
4327
dofs[3] = offset + c.entity_indices[1][0];
4328
dofs[4] = offset + c.entity_indices[1][1];
4329
dofs[5] = offset + c.entity_indices[1][2];
4330
}
4331
4332
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4333
virtual void tabulate_facet_dofs(unsigned int* dofs,
4334
unsigned int facet) const
4335
{
4336
switch ( facet )
4337
{
4338
case 0:
4339
dofs[0] = 1;
4340
dofs[1] = 2;
4341
dofs[2] = 3;
4342
break;
4343
case 1:
4344
dofs[0] = 0;
4345
dofs[1] = 2;
4346
dofs[2] = 4;
4347
break;
4348
case 2:
4349
dofs[0] = 0;
4350
dofs[1] = 1;
4351
dofs[2] = 5;
4352
break;
4353
}
4354
}
4355
4356
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4357
virtual void tabulate_entity_dofs(unsigned int* dofs,
4358
unsigned int d, unsigned int i) const
4359
{
4360
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
4361
}
4362
4363
/// Tabulate the coordinates of all dofs on a cell
4364
virtual void tabulate_coordinates(double** coordinates,
4365
const ufc::cell& c) const
4366
{
4367
const double * const * x = c.coordinates;
4368
coordinates[0][0] = x[0][0];
4369
coordinates[0][1] = x[0][1];
4370
coordinates[1][0] = x[1][0];
4371
coordinates[1][1] = x[1][1];
4372
coordinates[2][0] = x[2][0];
4373
coordinates[2][1] = x[2][1];
4374
coordinates[3][0] = 0.5*x[1][0] + 0.5*x[2][0];
4375
coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];
4376
coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];
4377
coordinates[4][1] = 0.5*x[0][1] + 0.5*x[2][1];
4378
coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];
4379
coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];
4380
}
4381
4382
/// Return the number of sub dof maps (for a mixed element)
4383
virtual unsigned int num_sub_dof_maps() const
4384
{
4385
return 1;
4386
}
4387
4388
/// Create a new dof_map for sub dof map i (for a mixed element)
4389
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
4390
{
4391
return new biharmonic_1_dof_map_1();
4392
}
4393
4394
};
4395
4396
/// This class defines the interface for the tabulation of the cell
4397
/// tensor corresponding to the local contribution to a form from
4398
/// the integral over a cell.
4399
4400
class biharmonic_1_cell_integral_0_quadrature: public ufc::cell_integral
4401
{
4402
public:
4403
4404
/// Constructor
4405
biharmonic_1_cell_integral_0_quadrature() : ufc::cell_integral()
4406
{
4407
// Do nothing
4408
}
4409
4410
/// Destructor
4411
virtual ~biharmonic_1_cell_integral_0_quadrature()
4412
{
4413
// Do nothing
4414
}
4415
4416
/// Tabulate the tensor for the contribution from a local cell
4417
virtual void tabulate_tensor(double* A,
4418
const double * const * w,
4419
const ufc::cell& c) const
4420
{
4421
// Extract vertex coordinates
4422
const double * const * x = c.coordinates;
4423
4424
// Compute Jacobian of affine map from reference cell
4425
const double J_00 = x[1][0] - x[0][0];
4426
const double J_01 = x[2][0] - x[0][0];
4427
const double J_10 = x[1][1] - x[0][1];
4428
const double J_11 = x[2][1] - x[0][1];
4429
4430
// Compute determinant of Jacobian
4431
double detJ = J_00*J_11 - J_01*J_10;
4432
4433
// Compute inverse of Jacobian
4434
4435
// Set scale factor
4436
const double det = std::abs(detJ);
4437
4438
4439
// Array of quadrature weights
4440
static const double W9[9] = {0.0558144205, 0.0636780851, 0.0193963833, 0.0893030728, 0.101884936, 0.0310342133, 0.0558144205, 0.0636780851, 0.0193963833};
4441
// Quadrature points on the UFC reference element: (0.102717655, 0.0885879595), (0.0665540678, 0.409466864), (0.0239311323, 0.787659462), (0.45570602, 0.0885879595), (0.295266568, 0.409466864), (0.106170269, 0.787659462), (0.808694386, 0.0885879595), (0.523979068, 0.409466864), (0.188409406, 0.787659462)
4442
4443
// Value of basis functions at quadrature points.
4444
static const double FE0[9][6] = \
4445
{{0.499278833, -0.0816158216, -0.0728923064, 0.0363981898, 0.286562342, 0.332268763},
4446
{0.0251290591, -0.0576951799, -0.0741406383, 0.109006742, 0.858208264, 0.139491754},
4447
{-0.117413197, -0.0227857341, 0.453155394, 0.0753983311, 0.593609805, 0.0180354017},
4448
{-0.0403700665, -0.0403700665, -0.0728923064, 0.161480266, 0.161480266, 0.830671908},
4449
{-0.120901876, -0.120901876, -0.0741406383, 0.483607503, 0.483607503, 0.348729384},
4450
{-0.083626017, -0.083626017, 0.453155394, 0.334504068, 0.334504068, 0.0450885042},
4451
{-0.0816158216, 0.499278833, -0.0728923064, 0.286562342, 0.0363981898, 0.332268763},
4452
{-0.0576951799, 0.0251290591, -0.0741406383, 0.858208264, 0.109006742, 0.139491754},
4453
{-0.0227857341, -0.117413197, 0.453155394, 0.593609805, 0.0753983311, 0.0180354017}};
4454
4455
4456
// Compute element tensor using UFL quadrature representation
4457
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
4458
// Total number of operations to compute element tensor: 324
4459
4460
// Loop quadrature points for integral
4461
// Number of operations to compute element tensor for following IP loop = 324
4462
for (unsigned int ip = 0; ip < 9; ip++)
4463
{
4464
4465
// Function declarations
4466
double F0 = 0;
4467
4468
// Total number of operations to compute function values = 12
4469
for (unsigned int r = 0; r < 6; r++)
4470
{
4471
F0 += FE0[ip][r]*w[0][r];
4472
}// end loop over 'r'
4473
4474
// Number of operations for primary indices: 24
4475
for (unsigned int j = 0; j < 6; j++)
4476
{
4477
// Number of operations to compute entry: 4
4478
A[j] += FE0[ip][j]*F0*W9[ip]*det;
4479
}// end loop over 'j'
4480
}// end loop over 'ip'
4481
}
4482
4483
};
4484
4485
/// This class defines the interface for the tabulation of the cell
4486
/// tensor corresponding to the local contribution to a form from
4487
/// the integral over a cell.
4488
4489
class biharmonic_1_cell_integral_0: public ufc::cell_integral
4490
{
4491
private:
4492
4493
biharmonic_1_cell_integral_0_quadrature integral_0_quadrature;
4494
4495
public:
4496
4497
/// Constructor
4498
biharmonic_1_cell_integral_0() : ufc::cell_integral()
4499
{
4500
// Do nothing
4501
}
4502
4503
/// Destructor
4504
virtual ~biharmonic_1_cell_integral_0()
4505
{
4506
// Do nothing
4507
}
4508
4509
/// Tabulate the tensor for the contribution from a local cell
4510
virtual void tabulate_tensor(double* A,
4511
const double * const * w,
4512
const ufc::cell& c) const
4513
{
4514
// Reset values of the element tensor block
4515
for (unsigned int j = 0; j < 6; j++)
4516
A[j] = 0;
4517
4518
// Add all contributions to element tensor
4519
integral_0_quadrature.tabulate_tensor(A, w, c);
4520
}
4521
4522
};
4523
4524
/// This class defines the interface for the assembly of the global
4525
/// tensor corresponding to a form with r + n arguments, that is, a
4526
/// mapping
4527
///
4528
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
4529
///
4530
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
4531
/// global tensor A is defined by
4532
///
4533
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
4534
///
4535
/// where each argument Vj represents the application to the
4536
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
4537
/// fixed functions (coefficients).
4538
4539
class biharmonic_form_1: public ufc::form
4540
{
4541
public:
4542
4543
/// Constructor
4544
biharmonic_form_1() : ufc::form()
4545
{
4546
// Do nothing
4547
}
4548
4549
/// Destructor
4550
virtual ~biharmonic_form_1()
4551
{
4552
// Do nothing
4553
}
4554
4555
/// Return a string identifying the form
4556
virtual const char* signature() const
4557
{
4558
return "Form([Integral(Product(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0), Function(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 2), 0)), Measure('cell', 0, None))])";
4559
}
4560
4561
/// Return the rank of the global tensor (r)
4562
virtual unsigned int rank() const
4563
{
4564
return 1;
4565
}
4566
4567
/// Return the number of coefficients (n)
4568
virtual unsigned int num_coefficients() const
4569
{
4570
return 1;
4571
}
4572
4573
/// Return the number of cell integrals
4574
virtual unsigned int num_cell_integrals() const
4575
{
4576
return 1;
4577
}
4578
4579
/// Return the number of exterior facet integrals
4580
virtual unsigned int num_exterior_facet_integrals() const
4581
{
4582
return 0;
4583
}
4584
4585
/// Return the number of interior facet integrals
4586
virtual unsigned int num_interior_facet_integrals() const
4587
{
4588
return 0;
4589
}
4590
4591
/// Create a new finite element for argument function i
4592
virtual ufc::finite_element* create_finite_element(unsigned int i) const
4593
{
4594
switch ( i )
4595
{
4596
case 0:
4597
return new biharmonic_1_finite_element_0();
4598
break;
4599
case 1:
4600
return new biharmonic_1_finite_element_1();
4601
break;
4602
}
4603
return 0;
4604
}
4605
4606
/// Create a new dof map for argument function i
4607
virtual ufc::dof_map* create_dof_map(unsigned int i) const
4608
{
4609
switch ( i )
4610
{
4611
case 0:
4612
return new biharmonic_1_dof_map_0();
4613
break;
4614
case 1:
4615
return new biharmonic_1_dof_map_1();
4616
break;
4617
}
4618
return 0;
4619
}
4620
4621
/// Create a new cell integral on sub domain i
4622
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
4623
{
4624
return new biharmonic_1_cell_integral_0();
4625
}
4626
4627
/// Create a new exterior facet integral on sub domain i
4628
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
4629
{
4630
return 0;
4631
}
4632
4633
/// Create a new interior facet integral on sub domain i
4634
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
4635
{
4636
return 0;
4637
}
4638
4639
};
4640
4641
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1