~ubuntu-branches/ubuntu/trusty/ffc/trusty
» Revision
4
: test/regression/reference/quadrature/Equation.h
« back to all changes in this revision
Viewing changes to test/regression/reference/quadrature/Equation.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/Equation.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 __EQUATION_H
5
#define __EQUATION_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 equation_0_finite_element_0: public ufc::finite_element
14
{
15
public:
16
17
/// Constructor
18
equation_0_finite_element_0() : ufc::finite_element()
19
{
20
// Do nothing
21
}
22
23
/// Destructor
24
virtual ~equation_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)), 1)";
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 3;
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
104
// Compute psitilde_a
105
const double psitilde_a_0 = 1;
106
const double psitilde_a_1 = x;
107
108
// Compute psitilde_bs
109
const double psitilde_bs_0_0 = 1;
110
const double psitilde_bs_0_1 = 1.5*y + 0.5;
111
const double psitilde_bs_1_0 = 1;
112
113
// Compute basisvalues
114
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
115
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
116
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
117
118
// Table(s) of coefficients
119
static const double coefficients0[3][3] = \
120
{{0.471404521, -0.288675135, -0.166666667},
121
{0.471404521, 0.288675135, -0.166666667},
122
{0.471404521, 0, 0.333333333}};
123
124
// Extract relevant coefficients
125
const double coeff0_0 = coefficients0[dof][0];
126
const double coeff0_1 = coefficients0[dof][1];
127
const double coeff0_2 = coefficients0[dof][2];
128
129
// Compute value(s)
130
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
131
}
132
133
/// Evaluate all basis functions at given point in cell
134
virtual void evaluate_basis_all(double* values,
135
const double* coordinates,
136
const ufc::cell& c) const
137
{
138
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
139
}
140
141
/// Evaluate order n derivatives of basis function i at given point in cell
142
virtual void evaluate_basis_derivatives(unsigned int i,
143
unsigned int n,
144
double* values,
145
const double* coordinates,
146
const ufc::cell& c) const
147
{
148
// Extract vertex coordinates
149
const double * const * element_coordinates = c.coordinates;
150
151
// Compute Jacobian of affine map from reference cell
152
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
153
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
154
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
155
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
156
157
// Compute determinant of Jacobian
158
const double detJ = J_00*J_11 - J_01*J_10;
159
160
// Compute inverse of Jacobian
161
162
// Get coordinates and map to the reference (UFC) element
163
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
164
element_coordinates[0][0]*element_coordinates[2][1] +\
165
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
166
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
167
element_coordinates[1][0]*element_coordinates[0][1] -\
168
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
169
170
// Map coordinates to the reference square
171
if (std::abs(y - 1.0) < 1e-08)
172
x = -1.0;
173
else
174
x = 2.0 *x/(1.0 - y) - 1.0;
175
y = 2.0*y - 1.0;
176
177
// Compute number of derivatives
178
unsigned int num_derivatives = 1;
179
180
for (unsigned int j = 0; j < n; j++)
181
num_derivatives *= 2;
182
183
184
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
185
unsigned int **combinations = new unsigned int *[num_derivatives];
186
187
for (unsigned int j = 0; j < num_derivatives; j++)
188
{
189
combinations[j] = new unsigned int [n];
190
for (unsigned int k = 0; k < n; k++)
191
combinations[j][k] = 0;
192
}
193
194
// Generate combinations of derivatives
195
for (unsigned int row = 1; row < num_derivatives; row++)
196
{
197
for (unsigned int num = 0; num < row; num++)
198
{
199
for (unsigned int col = n-1; col+1 > 0; col--)
200
{
201
if (combinations[row][col] + 1 > 1)
202
combinations[row][col] = 0;
203
else
204
{
205
combinations[row][col] += 1;
206
break;
207
}
208
}
209
}
210
}
211
212
// Compute inverse of Jacobian
213
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
214
215
// Declare transformation matrix
216
// Declare pointer to two dimensional array and initialise
217
double **transform = new double *[num_derivatives];
218
219
for (unsigned int j = 0; j < num_derivatives; j++)
220
{
221
transform[j] = new double [num_derivatives];
222
for (unsigned int k = 0; k < num_derivatives; k++)
223
transform[j][k] = 1;
224
}
225
226
// Construct transformation matrix
227
for (unsigned int row = 0; row < num_derivatives; row++)
228
{
229
for (unsigned int col = 0; col < num_derivatives; col++)
230
{
231
for (unsigned int k = 0; k < n; k++)
232
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
233
}
234
}
235
236
// Reset values
237
for (unsigned int j = 0; j < 1*num_derivatives; j++)
238
values[j] = 0;
239
240
// Map degree of freedom to element degree of freedom
241
const unsigned int dof = i;
242
243
// Generate scalings
244
const double scalings_y_0 = 1;
245
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
246
247
// Compute psitilde_a
248
const double psitilde_a_0 = 1;
249
const double psitilde_a_1 = x;
250
251
// Compute psitilde_bs
252
const double psitilde_bs_0_0 = 1;
253
const double psitilde_bs_0_1 = 1.5*y + 0.5;
254
const double psitilde_bs_1_0 = 1;
255
256
// Compute basisvalues
257
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
258
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
259
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
260
261
// Table(s) of coefficients
262
static const double coefficients0[3][3] = \
263
{{0.471404521, -0.288675135, -0.166666667},
264
{0.471404521, 0.288675135, -0.166666667},
265
{0.471404521, 0, 0.333333333}};
266
267
// Interesting (new) part
268
// Tables of derivatives of the polynomial base (transpose)
269
static const double dmats0[3][3] = \
270
{{0, 0, 0},
271
{4.89897949, 0, 0},
272
{0, 0, 0}};
273
274
static const double dmats1[3][3] = \
275
{{0, 0, 0},
276
{2.44948974, 0, 0},
277
{4.24264069, 0, 0}};
278
279
// Compute reference derivatives
280
// Declare pointer to array of derivatives on FIAT element
281
double *derivatives = new double [num_derivatives];
282
283
// Declare coefficients
284
double coeff0_0 = 0;
285
double coeff0_1 = 0;
286
double coeff0_2 = 0;
287
288
// Declare new coefficients
289
double new_coeff0_0 = 0;
290
double new_coeff0_1 = 0;
291
double new_coeff0_2 = 0;
292
293
// Loop possible derivatives
294
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
295
{
296
// Get values from coefficients array
297
new_coeff0_0 = coefficients0[dof][0];
298
new_coeff0_1 = coefficients0[dof][1];
299
new_coeff0_2 = coefficients0[dof][2];
300
301
// Loop derivative order
302
for (unsigned int j = 0; j < n; j++)
303
{
304
// Update old coefficients
305
coeff0_0 = new_coeff0_0;
306
coeff0_1 = new_coeff0_1;
307
coeff0_2 = new_coeff0_2;
308
309
if(combinations[deriv_num][j] == 0)
310
{
311
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
312
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
313
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
314
}
315
if(combinations[deriv_num][j] == 1)
316
{
317
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
318
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
319
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
320
}
321
322
}
323
// Compute derivatives on reference element as dot product of coefficients and basisvalues
324
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
325
}
326
327
// Transform derivatives back to physical element
328
for (unsigned int row = 0; row < num_derivatives; row++)
329
{
330
for (unsigned int col = 0; col < num_derivatives; col++)
331
{
332
values[row] += transform[row][col]*derivatives[col];
333
}
334
}
335
// Delete pointer to array of derivatives on FIAT element
336
delete [] derivatives;
337
338
// Delete pointer to array of combinations of derivatives and transform
339
for (unsigned int row = 0; row < num_derivatives; row++)
340
{
341
delete [] combinations[row];
342
delete [] transform[row];
343
}
344
345
delete [] combinations;
346
delete [] transform;
347
}
348
349
/// Evaluate order n derivatives of all basis functions at given point in cell
350
virtual void evaluate_basis_derivatives_all(unsigned int n,
351
double* values,
352
const double* coordinates,
353
const ufc::cell& c) const
354
{
355
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
356
}
357
358
/// Evaluate linear functional for dof i on the function f
359
virtual double evaluate_dof(unsigned int i,
360
const ufc::function& f,
361
const ufc::cell& c) const
362
{
363
// The reference points, direction and weights:
364
static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
365
static const double W[3][1] = {{1}, {1}, {1}};
366
static const double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
367
368
const double * const * x = c.coordinates;
369
double result = 0.0;
370
// Iterate over the points:
371
// Evaluate basis functions for affine mapping
372
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
373
const double w1 = X[i][0][0];
374
const double w2 = X[i][0][1];
375
376
// Compute affine mapping y = F(X)
377
double y[2];
378
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
379
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
380
381
// Evaluate function at physical points
382
double values[1];
383
f.evaluate(values, y, c);
384
385
// Map function values using appropriate mapping
386
// Affine map: Do nothing
387
388
// Note that we do not map the weights (yet).
389
390
// Take directional components
391
for(int k = 0; k < 1; k++)
392
result += values[k]*D[i][0][k];
393
// Multiply by weights
394
result *= W[i][0];
395
396
return result;
397
}
398
399
/// Evaluate linear functionals for all dofs on the function f
400
virtual void evaluate_dofs(double* values,
401
const ufc::function& f,
402
const ufc::cell& c) const
403
{
404
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
405
}
406
407
/// Interpolate vertex values from dof values
408
virtual void interpolate_vertex_values(double* vertex_values,
409
const double* dof_values,
410
const ufc::cell& c) const
411
{
412
// Evaluate at vertices and use affine mapping
413
vertex_values[0] = dof_values[0];
414
vertex_values[1] = dof_values[1];
415
vertex_values[2] = dof_values[2];
416
}
417
418
/// Return the number of sub elements (for a mixed element)
419
virtual unsigned int num_sub_elements() const
420
{
421
return 1;
422
}
423
424
/// Create a new finite element for sub element i (for a mixed element)
425
virtual ufc::finite_element* create_sub_element(unsigned int i) const
426
{
427
return new equation_0_finite_element_0();
428
}
429
430
};
431
432
/// This class defines the interface for a finite element.
433
434
class equation_0_finite_element_1: public ufc::finite_element
435
{
436
public:
437
438
/// Constructor
439
equation_0_finite_element_1() : ufc::finite_element()
440
{
441
// Do nothing
442
}
443
444
/// Destructor
445
virtual ~equation_0_finite_element_1()
446
{
447
// Do nothing
448
}
449
450
/// Return a string identifying the finite element
451
virtual const char* signature() const
452
{
453
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
454
}
455
456
/// Return the cell shape
457
virtual ufc::shape cell_shape() const
458
{
459
return ufc::triangle;
460
}
461
462
/// Return the dimension of the finite element function space
463
virtual unsigned int space_dimension() const
464
{
465
return 3;
466
}
467
468
/// Return the rank of the value space
469
virtual unsigned int value_rank() const
470
{
471
return 0;
472
}
473
474
/// Return the dimension of the value space for axis i
475
virtual unsigned int value_dimension(unsigned int i) const
476
{
477
return 1;
478
}
479
480
/// Evaluate basis function i at given point in cell
481
virtual void evaluate_basis(unsigned int i,
482
double* values,
483
const double* coordinates,
484
const ufc::cell& c) const
485
{
486
// Extract vertex coordinates
487
const double * const * element_coordinates = c.coordinates;
488
489
// Compute Jacobian of affine map from reference cell
490
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
491
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
492
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
493
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
494
495
// Compute determinant of Jacobian
496
const double detJ = J_00*J_11 - J_01*J_10;
497
498
// Compute inverse of Jacobian
499
500
// Get coordinates and map to the reference (UFC) element
501
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
502
element_coordinates[0][0]*element_coordinates[2][1] +\
503
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
504
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
505
element_coordinates[1][0]*element_coordinates[0][1] -\
506
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
507
508
// Map coordinates to the reference square
509
if (std::abs(y - 1.0) < 1e-08)
510
x = -1.0;
511
else
512
x = 2.0 *x/(1.0 - y) - 1.0;
513
y = 2.0*y - 1.0;
514
515
// Reset values
516
*values = 0;
517
518
// Map degree of freedom to element degree of freedom
519
const unsigned int dof = i;
520
521
// Generate scalings
522
const double scalings_y_0 = 1;
523
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
524
525
// Compute psitilde_a
526
const double psitilde_a_0 = 1;
527
const double psitilde_a_1 = x;
528
529
// Compute psitilde_bs
530
const double psitilde_bs_0_0 = 1;
531
const double psitilde_bs_0_1 = 1.5*y + 0.5;
532
const double psitilde_bs_1_0 = 1;
533
534
// Compute basisvalues
535
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
536
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
537
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
538
539
// Table(s) of coefficients
540
static const double coefficients0[3][3] = \
541
{{0.471404521, -0.288675135, -0.166666667},
542
{0.471404521, 0.288675135, -0.166666667},
543
{0.471404521, 0, 0.333333333}};
544
545
// Extract relevant coefficients
546
const double coeff0_0 = coefficients0[dof][0];
547
const double coeff0_1 = coefficients0[dof][1];
548
const double coeff0_2 = coefficients0[dof][2];
549
550
// Compute value(s)
551
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
552
}
553
554
/// Evaluate all basis functions at given point in cell
555
virtual void evaluate_basis_all(double* values,
556
const double* coordinates,
557
const ufc::cell& c) const
558
{
559
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
560
}
561
562
/// Evaluate order n derivatives of basis function i at given point in cell
563
virtual void evaluate_basis_derivatives(unsigned int i,
564
unsigned int n,
565
double* values,
566
const double* coordinates,
567
const ufc::cell& c) const
568
{
569
// Extract vertex coordinates
570
const double * const * element_coordinates = c.coordinates;
571
572
// Compute Jacobian of affine map from reference cell
573
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
574
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
575
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
576
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
577
578
// Compute determinant of Jacobian
579
const double detJ = J_00*J_11 - J_01*J_10;
580
581
// Compute inverse of Jacobian
582
583
// Get coordinates and map to the reference (UFC) element
584
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
585
element_coordinates[0][0]*element_coordinates[2][1] +\
586
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
587
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
588
element_coordinates[1][0]*element_coordinates[0][1] -\
589
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
590
591
// Map coordinates to the reference square
592
if (std::abs(y - 1.0) < 1e-08)
593
x = -1.0;
594
else
595
x = 2.0 *x/(1.0 - y) - 1.0;
596
y = 2.0*y - 1.0;
597
598
// Compute number of derivatives
599
unsigned int num_derivatives = 1;
600
601
for (unsigned int j = 0; j < n; j++)
602
num_derivatives *= 2;
603
604
605
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
606
unsigned int **combinations = new unsigned int *[num_derivatives];
607
608
for (unsigned int j = 0; j < num_derivatives; j++)
609
{
610
combinations[j] = new unsigned int [n];
611
for (unsigned int k = 0; k < n; k++)
612
combinations[j][k] = 0;
613
}
614
615
// Generate combinations of derivatives
616
for (unsigned int row = 1; row < num_derivatives; row++)
617
{
618
for (unsigned int num = 0; num < row; num++)
619
{
620
for (unsigned int col = n-1; col+1 > 0; col--)
621
{
622
if (combinations[row][col] + 1 > 1)
623
combinations[row][col] = 0;
624
else
625
{
626
combinations[row][col] += 1;
627
break;
628
}
629
}
630
}
631
}
632
633
// Compute inverse of Jacobian
634
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
635
636
// Declare transformation matrix
637
// Declare pointer to two dimensional array and initialise
638
double **transform = new double *[num_derivatives];
639
640
for (unsigned int j = 0; j < num_derivatives; j++)
641
{
642
transform[j] = new double [num_derivatives];
643
for (unsigned int k = 0; k < num_derivatives; k++)
644
transform[j][k] = 1;
645
}
646
647
// Construct transformation matrix
648
for (unsigned int row = 0; row < num_derivatives; row++)
649
{
650
for (unsigned int col = 0; col < num_derivatives; col++)
651
{
652
for (unsigned int k = 0; k < n; k++)
653
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
654
}
655
}
656
657
// Reset values
658
for (unsigned int j = 0; j < 1*num_derivatives; j++)
659
values[j] = 0;
660
661
// Map degree of freedom to element degree of freedom
662
const unsigned int dof = i;
663
664
// Generate scalings
665
const double scalings_y_0 = 1;
666
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
667
668
// Compute psitilde_a
669
const double psitilde_a_0 = 1;
670
const double psitilde_a_1 = x;
671
672
// Compute psitilde_bs
673
const double psitilde_bs_0_0 = 1;
674
const double psitilde_bs_0_1 = 1.5*y + 0.5;
675
const double psitilde_bs_1_0 = 1;
676
677
// Compute basisvalues
678
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
679
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
680
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
681
682
// Table(s) of coefficients
683
static const double coefficients0[3][3] = \
684
{{0.471404521, -0.288675135, -0.166666667},
685
{0.471404521, 0.288675135, -0.166666667},
686
{0.471404521, 0, 0.333333333}};
687
688
// Interesting (new) part
689
// Tables of derivatives of the polynomial base (transpose)
690
static const double dmats0[3][3] = \
691
{{0, 0, 0},
692
{4.89897949, 0, 0},
693
{0, 0, 0}};
694
695
static const double dmats1[3][3] = \
696
{{0, 0, 0},
697
{2.44948974, 0, 0},
698
{4.24264069, 0, 0}};
699
700
// Compute reference derivatives
701
// Declare pointer to array of derivatives on FIAT element
702
double *derivatives = new double [num_derivatives];
703
704
// Declare coefficients
705
double coeff0_0 = 0;
706
double coeff0_1 = 0;
707
double coeff0_2 = 0;
708
709
// Declare new coefficients
710
double new_coeff0_0 = 0;
711
double new_coeff0_1 = 0;
712
double new_coeff0_2 = 0;
713
714
// Loop possible derivatives
715
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
716
{
717
// Get values from coefficients array
718
new_coeff0_0 = coefficients0[dof][0];
719
new_coeff0_1 = coefficients0[dof][1];
720
new_coeff0_2 = coefficients0[dof][2];
721
722
// Loop derivative order
723
for (unsigned int j = 0; j < n; j++)
724
{
725
// Update old coefficients
726
coeff0_0 = new_coeff0_0;
727
coeff0_1 = new_coeff0_1;
728
coeff0_2 = new_coeff0_2;
729
730
if(combinations[deriv_num][j] == 0)
731
{
732
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
733
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
734
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
735
}
736
if(combinations[deriv_num][j] == 1)
737
{
738
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
739
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
740
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
741
}
742
743
}
744
// Compute derivatives on reference element as dot product of coefficients and basisvalues
745
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
746
}
747
748
// Transform derivatives back to physical element
749
for (unsigned int row = 0; row < num_derivatives; row++)
750
{
751
for (unsigned int col = 0; col < num_derivatives; col++)
752
{
753
values[row] += transform[row][col]*derivatives[col];
754
}
755
}
756
// Delete pointer to array of derivatives on FIAT element
757
delete [] derivatives;
758
759
// Delete pointer to array of combinations of derivatives and transform
760
for (unsigned int row = 0; row < num_derivatives; row++)
761
{
762
delete [] combinations[row];
763
delete [] transform[row];
764
}
765
766
delete [] combinations;
767
delete [] transform;
768
}
769
770
/// Evaluate order n derivatives of all basis functions at given point in cell
771
virtual void evaluate_basis_derivatives_all(unsigned int n,
772
double* values,
773
const double* coordinates,
774
const ufc::cell& c) const
775
{
776
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
777
}
778
779
/// Evaluate linear functional for dof i on the function f
780
virtual double evaluate_dof(unsigned int i,
781
const ufc::function& f,
782
const ufc::cell& c) const
783
{
784
// The reference points, direction and weights:
785
static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
786
static const double W[3][1] = {{1}, {1}, {1}};
787
static const double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
788
789
const double * const * x = c.coordinates;
790
double result = 0.0;
791
// Iterate over the points:
792
// Evaluate basis functions for affine mapping
793
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
794
const double w1 = X[i][0][0];
795
const double w2 = X[i][0][1];
796
797
// Compute affine mapping y = F(X)
798
double y[2];
799
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
800
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
801
802
// Evaluate function at physical points
803
double values[1];
804
f.evaluate(values, y, c);
805
806
// Map function values using appropriate mapping
807
// Affine map: Do nothing
808
809
// Note that we do not map the weights (yet).
810
811
// Take directional components
812
for(int k = 0; k < 1; k++)
813
result += values[k]*D[i][0][k];
814
// Multiply by weights
815
result *= W[i][0];
816
817
return result;
818
}
819
820
/// Evaluate linear functionals for all dofs on the function f
821
virtual void evaluate_dofs(double* values,
822
const ufc::function& f,
823
const ufc::cell& c) const
824
{
825
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
826
}
827
828
/// Interpolate vertex values from dof values
829
virtual void interpolate_vertex_values(double* vertex_values,
830
const double* dof_values,
831
const ufc::cell& c) const
832
{
833
// Evaluate at vertices and use affine mapping
834
vertex_values[0] = dof_values[0];
835
vertex_values[1] = dof_values[1];
836
vertex_values[2] = dof_values[2];
837
}
838
839
/// Return the number of sub elements (for a mixed element)
840
virtual unsigned int num_sub_elements() const
841
{
842
return 1;
843
}
844
845
/// Create a new finite element for sub element i (for a mixed element)
846
virtual ufc::finite_element* create_sub_element(unsigned int i) const
847
{
848
return new equation_0_finite_element_1();
849
}
850
851
};
852
853
/// This class defines the interface for a local-to-global mapping of
854
/// degrees of freedom (dofs).
855
856
class equation_0_dof_map_0: public ufc::dof_map
857
{
858
private:
859
860
unsigned int __global_dimension;
861
862
public:
863
864
/// Constructor
865
equation_0_dof_map_0() : ufc::dof_map()
866
{
867
__global_dimension = 0;
868
}
869
870
/// Destructor
871
virtual ~equation_0_dof_map_0()
872
{
873
// Do nothing
874
}
875
876
/// Return a string identifying the dof map
877
virtual const char* signature() const
878
{
879
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
880
}
881
882
/// Return true iff mesh entities of topological dimension d are needed
883
virtual bool needs_mesh_entities(unsigned int d) const
884
{
885
switch ( d )
886
{
887
case 0:
888
return true;
889
break;
890
case 1:
891
return false;
892
break;
893
case 2:
894
return false;
895
break;
896
}
897
return false;
898
}
899
900
/// Initialize dof map for mesh (return true iff init_cell() is needed)
901
virtual bool init_mesh(const ufc::mesh& m)
902
{
903
__global_dimension = m.num_entities[0];
904
return false;
905
}
906
907
/// Initialize dof map for given cell
908
virtual void init_cell(const ufc::mesh& m,
909
const ufc::cell& c)
910
{
911
// Do nothing
912
}
913
914
/// Finish initialization of dof map for cells
915
virtual void init_cell_finalize()
916
{
917
// Do nothing
918
}
919
920
/// Return the dimension of the global finite element function space
921
virtual unsigned int global_dimension() const
922
{
923
return __global_dimension;
924
}
925
926
/// Return the dimension of the local finite element function space for a cell
927
virtual unsigned int local_dimension(const ufc::cell& c) const
928
{
929
return 3;
930
}
931
932
/// Return the maximum dimension of the local finite element function space
933
virtual unsigned int max_local_dimension() const
934
{
935
return 3;
936
}
937
938
// Return the geometric dimension of the coordinates this dof map provides
939
virtual unsigned int geometric_dimension() const
940
{
941
return 2;
942
}
943
944
/// Return the number of dofs on each cell facet
945
virtual unsigned int num_facet_dofs() const
946
{
947
return 2;
948
}
949
950
/// Return the number of dofs associated with each cell entity of dimension d
951
virtual unsigned int num_entity_dofs(unsigned int d) const
952
{
953
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
954
}
955
956
/// Tabulate the local-to-global mapping of dofs on a cell
957
virtual void tabulate_dofs(unsigned int* dofs,
958
const ufc::mesh& m,
959
const ufc::cell& c) const
960
{
961
dofs[0] = c.entity_indices[0][0];
962
dofs[1] = c.entity_indices[0][1];
963
dofs[2] = c.entity_indices[0][2];
964
}
965
966
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
967
virtual void tabulate_facet_dofs(unsigned int* dofs,
968
unsigned int facet) const
969
{
970
switch ( facet )
971
{
972
case 0:
973
dofs[0] = 1;
974
dofs[1] = 2;
975
break;
976
case 1:
977
dofs[0] = 0;
978
dofs[1] = 2;
979
break;
980
case 2:
981
dofs[0] = 0;
982
dofs[1] = 1;
983
break;
984
}
985
}
986
987
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
988
virtual void tabulate_entity_dofs(unsigned int* dofs,
989
unsigned int d, unsigned int i) const
990
{
991
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
992
}
993
994
/// Tabulate the coordinates of all dofs on a cell
995
virtual void tabulate_coordinates(double** coordinates,
996
const ufc::cell& c) const
997
{
998
const double * const * x = c.coordinates;
999
coordinates[0][0] = x[0][0];
1000
coordinates[0][1] = x[0][1];
1001
coordinates[1][0] = x[1][0];
1002
coordinates[1][1] = x[1][1];
1003
coordinates[2][0] = x[2][0];
1004
coordinates[2][1] = x[2][1];
1005
}
1006
1007
/// Return the number of sub dof maps (for a mixed element)
1008
virtual unsigned int num_sub_dof_maps() const
1009
{
1010
return 1;
1011
}
1012
1013
/// Create a new dof_map for sub dof map i (for a mixed element)
1014
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1015
{
1016
return new equation_0_dof_map_0();
1017
}
1018
1019
};
1020
1021
/// This class defines the interface for a local-to-global mapping of
1022
/// degrees of freedom (dofs).
1023
1024
class equation_0_dof_map_1: public ufc::dof_map
1025
{
1026
private:
1027
1028
unsigned int __global_dimension;
1029
1030
public:
1031
1032
/// Constructor
1033
equation_0_dof_map_1() : ufc::dof_map()
1034
{
1035
__global_dimension = 0;
1036
}
1037
1038
/// Destructor
1039
virtual ~equation_0_dof_map_1()
1040
{
1041
// Do nothing
1042
}
1043
1044
/// Return a string identifying the dof map
1045
virtual const char* signature() const
1046
{
1047
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
1048
}
1049
1050
/// Return true iff mesh entities of topological dimension d are needed
1051
virtual bool needs_mesh_entities(unsigned int d) const
1052
{
1053
switch ( d )
1054
{
1055
case 0:
1056
return true;
1057
break;
1058
case 1:
1059
return false;
1060
break;
1061
case 2:
1062
return false;
1063
break;
1064
}
1065
return false;
1066
}
1067
1068
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1069
virtual bool init_mesh(const ufc::mesh& m)
1070
{
1071
__global_dimension = m.num_entities[0];
1072
return false;
1073
}
1074
1075
/// Initialize dof map for given cell
1076
virtual void init_cell(const ufc::mesh& m,
1077
const ufc::cell& c)
1078
{
1079
// Do nothing
1080
}
1081
1082
/// Finish initialization of dof map for cells
1083
virtual void init_cell_finalize()
1084
{
1085
// Do nothing
1086
}
1087
1088
/// Return the dimension of the global finite element function space
1089
virtual unsigned int global_dimension() const
1090
{
1091
return __global_dimension;
1092
}
1093
1094
/// Return the dimension of the local finite element function space for a cell
1095
virtual unsigned int local_dimension(const ufc::cell& c) const
1096
{
1097
return 3;
1098
}
1099
1100
/// Return the maximum dimension of the local finite element function space
1101
virtual unsigned int max_local_dimension() const
1102
{
1103
return 3;
1104
}
1105
1106
// Return the geometric dimension of the coordinates this dof map provides
1107
virtual unsigned int geometric_dimension() const
1108
{
1109
return 2;
1110
}
1111
1112
/// Return the number of dofs on each cell facet
1113
virtual unsigned int num_facet_dofs() const
1114
{
1115
return 2;
1116
}
1117
1118
/// Return the number of dofs associated with each cell entity of dimension d
1119
virtual unsigned int num_entity_dofs(unsigned int d) const
1120
{
1121
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1122
}
1123
1124
/// Tabulate the local-to-global mapping of dofs on a cell
1125
virtual void tabulate_dofs(unsigned int* dofs,
1126
const ufc::mesh& m,
1127
const ufc::cell& c) const
1128
{
1129
dofs[0] = c.entity_indices[0][0];
1130
dofs[1] = c.entity_indices[0][1];
1131
dofs[2] = c.entity_indices[0][2];
1132
}
1133
1134
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
1135
virtual void tabulate_facet_dofs(unsigned int* dofs,
1136
unsigned int facet) const
1137
{
1138
switch ( facet )
1139
{
1140
case 0:
1141
dofs[0] = 1;
1142
dofs[1] = 2;
1143
break;
1144
case 1:
1145
dofs[0] = 0;
1146
dofs[1] = 2;
1147
break;
1148
case 2:
1149
dofs[0] = 0;
1150
dofs[1] = 1;
1151
break;
1152
}
1153
}
1154
1155
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
1156
virtual void tabulate_entity_dofs(unsigned int* dofs,
1157
unsigned int d, unsigned int i) const
1158
{
1159
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1160
}
1161
1162
/// Tabulate the coordinates of all dofs on a cell
1163
virtual void tabulate_coordinates(double** coordinates,
1164
const ufc::cell& c) const
1165
{
1166
const double * const * x = c.coordinates;
1167
coordinates[0][0] = x[0][0];
1168
coordinates[0][1] = x[0][1];
1169
coordinates[1][0] = x[1][0];
1170
coordinates[1][1] = x[1][1];
1171
coordinates[2][0] = x[2][0];
1172
coordinates[2][1] = x[2][1];
1173
}
1174
1175
/// Return the number of sub dof maps (for a mixed element)
1176
virtual unsigned int num_sub_dof_maps() const
1177
{
1178
return 1;
1179
}
1180
1181
/// Create a new dof_map for sub dof map i (for a mixed element)
1182
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
1183
{
1184
return new equation_0_dof_map_1();
1185
}
1186
1187
};
1188
1189
/// This class defines the interface for the tabulation of the cell
1190
/// tensor corresponding to the local contribution to a form from
1191
/// the integral over a cell.
1192
1193
class equation_0_cell_integral_0_quadrature: public ufc::cell_integral
1194
{
1195
public:
1196
1197
/// Constructor
1198
equation_0_cell_integral_0_quadrature() : ufc::cell_integral()
1199
{
1200
// Do nothing
1201
}
1202
1203
/// Destructor
1204
virtual ~equation_0_cell_integral_0_quadrature()
1205
{
1206
// Do nothing
1207
}
1208
1209
/// Tabulate the tensor for the contribution from a local cell
1210
virtual void tabulate_tensor(double* A,
1211
const double * const * w,
1212
const ufc::cell& c) const
1213
{
1214
// Extract vertex coordinates
1215
const double * const * x = c.coordinates;
1216
1217
// Compute Jacobian of affine map from reference cell
1218
const double J_00 = x[1][0] - x[0][0];
1219
const double J_01 = x[2][0] - x[0][0];
1220
const double J_10 = x[1][1] - x[0][1];
1221
const double J_11 = x[2][1] - x[0][1];
1222
1223
// Compute determinant of Jacobian
1224
double detJ = J_00*J_11 - J_01*J_10;
1225
1226
// Compute inverse of Jacobian
1227
const double Jinv_00 = J_11 / detJ;
1228
const double Jinv_01 = -J_01 / detJ;
1229
const double Jinv_10 = -J_10 / detJ;
1230
const double Jinv_11 = J_00 / detJ;
1231
1232
// Set scale factor
1233
const double det = std::abs(detJ);
1234
1235
1236
// Array of quadrature weights
1237
static const double W4[4] = {0.159020691, 0.0909793091, 0.159020691, 0.0909793091};
1238
// Quadrature points on the UFC reference element: (0.178558728, 0.155051026), (0.0750311102, 0.644948974), (0.666390246, 0.155051026), (0.280019915, 0.644948974)
1239
1240
// Value of basis functions at quadrature points.
1241
static const double FE0[4][3] = \
1242
{{0.666390246, 0.178558728, 0.155051026},
1243
{0.280019915, 0.0750311102, 0.644948974},
1244
{0.178558728, 0.666390246, 0.155051026},
1245
{0.0750311102, 0.280019915, 0.644948974}};
1246
1247
static const double FE0_D01[4][3] = \
1248
{{-1, 0, 1},
1249
{-1, 0, 1},
1250
{-1, 0, 1},
1251
{-1, 0, 1}};
1252
1253
static const double FE0_D10[4][3] = \
1254
{{-1, 1, 0},
1255
{-1, 1, 0},
1256
{-1, 1, 0},
1257
{-1, 1, 0}};
1258
1259
1260
// Compute element tensor using UFL quadrature representation
1261
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
1262
// Total number of operations to compute element tensor: 828
1263
1264
// Loop quadrature points for integral
1265
// Number of operations to compute element tensor for following IP loop = 828
1266
for (unsigned int ip = 0; ip < 4; ip++)
1267
{
1268
1269
// Number of operations for primary indices: 207
1270
for (unsigned int j = 0; j < 3; j++)
1271
{
1272
for (unsigned int k = 0; k < 3; k++)
1273
{
1274
// Number of operations to compute entry: 23
1275
A[j*3 + k] += (FE0[ip][j]*FE0[ip][k] + ((Jinv_00*FE0_D10[ip][k] + Jinv_10*FE0_D01[ip][k])*0.5*(Jinv_00*FE0_D10[ip][j] + Jinv_10*FE0_D01[ip][j]) + (Jinv_01*FE0_D10[ip][k] + Jinv_11*FE0_D01[ip][k])*0.5*(Jinv_01*FE0_D10[ip][j] + Jinv_11*FE0_D01[ip][j]))*0.1)*W4[ip]*det;
1276
}// end loop over 'k'
1277
}// end loop over 'j'
1278
}// end loop over 'ip'
1279
}
1280
1281
};
1282
1283
/// This class defines the interface for the tabulation of the cell
1284
/// tensor corresponding to the local contribution to a form from
1285
/// the integral over a cell.
1286
1287
class equation_0_cell_integral_0: public ufc::cell_integral
1288
{
1289
private:
1290
1291
equation_0_cell_integral_0_quadrature integral_0_quadrature;
1292
1293
public:
1294
1295
/// Constructor
1296
equation_0_cell_integral_0() : ufc::cell_integral()
1297
{
1298
// Do nothing
1299
}
1300
1301
/// Destructor
1302
virtual ~equation_0_cell_integral_0()
1303
{
1304
// Do nothing
1305
}
1306
1307
/// Tabulate the tensor for the contribution from a local cell
1308
virtual void tabulate_tensor(double* A,
1309
const double * const * w,
1310
const ufc::cell& c) const
1311
{
1312
// Reset values of the element tensor block
1313
for (unsigned int j = 0; j < 9; j++)
1314
A[j] = 0;
1315
1316
// Add all contributions to element tensor
1317
integral_0_quadrature.tabulate_tensor(A, w, c);
1318
}
1319
1320
};
1321
1322
/// This class defines the interface for the assembly of the global
1323
/// tensor corresponding to a form with r + n arguments, that is, a
1324
/// mapping
1325
///
1326
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
1327
///
1328
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
1329
/// global tensor A is defined by
1330
///
1331
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
1332
///
1333
/// where each argument Vj represents the application to the
1334
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
1335
/// fixed functions (coefficients).
1336
1337
class equation_form_0: public ufc::form
1338
{
1339
public:
1340
1341
/// Constructor
1342
equation_form_0() : ufc::form()
1343
{
1344
// Do nothing
1345
}
1346
1347
/// Destructor
1348
virtual ~equation_form_0()
1349
{
1350
// Do nothing
1351
}
1352
1353
/// Return a string identifying the form
1354
virtual const char* signature() const
1355
{
1356
return "Form([Integral(Sum(Product(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1)), Product(FloatValue(0.10000000000000001, (), (), {}), IndexSum(Product(Indexed(ComponentTensor(Product(FloatValue(0.5, (), (), {}), SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1), MultiIndex((Index(0),), {Index(0): 2}))), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})))), Measure('cell', 0, None))])";
1357
}
1358
1359
/// Return the rank of the global tensor (r)
1360
virtual unsigned int rank() const
1361
{
1362
return 2;
1363
}
1364
1365
/// Return the number of coefficients (n)
1366
virtual unsigned int num_coefficients() const
1367
{
1368
return 0;
1369
}
1370
1371
/// Return the number of cell integrals
1372
virtual unsigned int num_cell_integrals() const
1373
{
1374
return 1;
1375
}
1376
1377
/// Return the number of exterior facet integrals
1378
virtual unsigned int num_exterior_facet_integrals() const
1379
{
1380
return 0;
1381
}
1382
1383
/// Return the number of interior facet integrals
1384
virtual unsigned int num_interior_facet_integrals() const
1385
{
1386
return 0;
1387
}
1388
1389
/// Create a new finite element for argument function i
1390
virtual ufc::finite_element* create_finite_element(unsigned int i) const
1391
{
1392
switch ( i )
1393
{
1394
case 0:
1395
return new equation_0_finite_element_0();
1396
break;
1397
case 1:
1398
return new equation_0_finite_element_1();
1399
break;
1400
}
1401
return 0;
1402
}
1403
1404
/// Create a new dof map for argument function i
1405
virtual ufc::dof_map* create_dof_map(unsigned int i) const
1406
{
1407
switch ( i )
1408
{
1409
case 0:
1410
return new equation_0_dof_map_0();
1411
break;
1412
case 1:
1413
return new equation_0_dof_map_1();
1414
break;
1415
}
1416
return 0;
1417
}
1418
1419
/// Create a new cell integral on sub domain i
1420
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
1421
{
1422
return new equation_0_cell_integral_0();
1423
}
1424
1425
/// Create a new exterior facet integral on sub domain i
1426
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
1427
{
1428
return 0;
1429
}
1430
1431
/// Create a new interior facet integral on sub domain i
1432
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
1433
{
1434
return 0;
1435
}
1436
1437
};
1438
1439
/// This class defines the interface for a finite element.
1440
1441
class equation_1_finite_element_0: public ufc::finite_element
1442
{
1443
public:
1444
1445
/// Constructor
1446
equation_1_finite_element_0() : ufc::finite_element()
1447
{
1448
// Do nothing
1449
}
1450
1451
/// Destructor
1452
virtual ~equation_1_finite_element_0()
1453
{
1454
// Do nothing
1455
}
1456
1457
/// Return a string identifying the finite element
1458
virtual const char* signature() const
1459
{
1460
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
1461
}
1462
1463
/// Return the cell shape
1464
virtual ufc::shape cell_shape() const
1465
{
1466
return ufc::triangle;
1467
}
1468
1469
/// Return the dimension of the finite element function space
1470
virtual unsigned int space_dimension() const
1471
{
1472
return 3;
1473
}
1474
1475
/// Return the rank of the value space
1476
virtual unsigned int value_rank() const
1477
{
1478
return 0;
1479
}
1480
1481
/// Return the dimension of the value space for axis i
1482
virtual unsigned int value_dimension(unsigned int i) const
1483
{
1484
return 1;
1485
}
1486
1487
/// Evaluate basis function i at given point in cell
1488
virtual void evaluate_basis(unsigned int i,
1489
double* values,
1490
const double* coordinates,
1491
const ufc::cell& c) const
1492
{
1493
// Extract vertex coordinates
1494
const double * const * element_coordinates = c.coordinates;
1495
1496
// Compute Jacobian of affine map from reference cell
1497
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1498
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1499
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1500
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1501
1502
// Compute determinant of Jacobian
1503
const double detJ = J_00*J_11 - J_01*J_10;
1504
1505
// Compute inverse of Jacobian
1506
1507
// Get coordinates and map to the reference (UFC) element
1508
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1509
element_coordinates[0][0]*element_coordinates[2][1] +\
1510
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1511
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1512
element_coordinates[1][0]*element_coordinates[0][1] -\
1513
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1514
1515
// Map coordinates to the reference square
1516
if (std::abs(y - 1.0) < 1e-08)
1517
x = -1.0;
1518
else
1519
x = 2.0 *x/(1.0 - y) - 1.0;
1520
y = 2.0*y - 1.0;
1521
1522
// Reset values
1523
*values = 0;
1524
1525
// Map degree of freedom to element degree of freedom
1526
const unsigned int dof = i;
1527
1528
// Generate scalings
1529
const double scalings_y_0 = 1;
1530
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1531
1532
// Compute psitilde_a
1533
const double psitilde_a_0 = 1;
1534
const double psitilde_a_1 = x;
1535
1536
// Compute psitilde_bs
1537
const double psitilde_bs_0_0 = 1;
1538
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1539
const double psitilde_bs_1_0 = 1;
1540
1541
// Compute basisvalues
1542
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1543
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1544
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1545
1546
// Table(s) of coefficients
1547
static const double coefficients0[3][3] = \
1548
{{0.471404521, -0.288675135, -0.166666667},
1549
{0.471404521, 0.288675135, -0.166666667},
1550
{0.471404521, 0, 0.333333333}};
1551
1552
// Extract relevant coefficients
1553
const double coeff0_0 = coefficients0[dof][0];
1554
const double coeff0_1 = coefficients0[dof][1];
1555
const double coeff0_2 = coefficients0[dof][2];
1556
1557
// Compute value(s)
1558
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
1559
}
1560
1561
/// Evaluate all basis functions at given point in cell
1562
virtual void evaluate_basis_all(double* values,
1563
const double* coordinates,
1564
const ufc::cell& c) const
1565
{
1566
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1567
}
1568
1569
/// Evaluate order n derivatives of basis function i at given point in cell
1570
virtual void evaluate_basis_derivatives(unsigned int i,
1571
unsigned int n,
1572
double* values,
1573
const double* coordinates,
1574
const ufc::cell& c) const
1575
{
1576
// Extract vertex coordinates
1577
const double * const * element_coordinates = c.coordinates;
1578
1579
// Compute Jacobian of affine map from reference cell
1580
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1581
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1582
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1583
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1584
1585
// Compute determinant of Jacobian
1586
const double detJ = J_00*J_11 - J_01*J_10;
1587
1588
// Compute inverse of Jacobian
1589
1590
// Get coordinates and map to the reference (UFC) element
1591
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1592
element_coordinates[0][0]*element_coordinates[2][1] +\
1593
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1594
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1595
element_coordinates[1][0]*element_coordinates[0][1] -\
1596
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1597
1598
// Map coordinates to the reference square
1599
if (std::abs(y - 1.0) < 1e-08)
1600
x = -1.0;
1601
else
1602
x = 2.0 *x/(1.0 - y) - 1.0;
1603
y = 2.0*y - 1.0;
1604
1605
// Compute number of derivatives
1606
unsigned int num_derivatives = 1;
1607
1608
for (unsigned int j = 0; j < n; j++)
1609
num_derivatives *= 2;
1610
1611
1612
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1613
unsigned int **combinations = new unsigned int *[num_derivatives];
1614
1615
for (unsigned int j = 0; j < num_derivatives; j++)
1616
{
1617
combinations[j] = new unsigned int [n];
1618
for (unsigned int k = 0; k < n; k++)
1619
combinations[j][k] = 0;
1620
}
1621
1622
// Generate combinations of derivatives
1623
for (unsigned int row = 1; row < num_derivatives; row++)
1624
{
1625
for (unsigned int num = 0; num < row; num++)
1626
{
1627
for (unsigned int col = n-1; col+1 > 0; col--)
1628
{
1629
if (combinations[row][col] + 1 > 1)
1630
combinations[row][col] = 0;
1631
else
1632
{
1633
combinations[row][col] += 1;
1634
break;
1635
}
1636
}
1637
}
1638
}
1639
1640
// Compute inverse of Jacobian
1641
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1642
1643
// Declare transformation matrix
1644
// Declare pointer to two dimensional array and initialise
1645
double **transform = new double *[num_derivatives];
1646
1647
for (unsigned int j = 0; j < num_derivatives; j++)
1648
{
1649
transform[j] = new double [num_derivatives];
1650
for (unsigned int k = 0; k < num_derivatives; k++)
1651
transform[j][k] = 1;
1652
}
1653
1654
// Construct transformation matrix
1655
for (unsigned int row = 0; row < num_derivatives; row++)
1656
{
1657
for (unsigned int col = 0; col < num_derivatives; col++)
1658
{
1659
for (unsigned int k = 0; k < n; k++)
1660
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1661
}
1662
}
1663
1664
// Reset values
1665
for (unsigned int j = 0; j < 1*num_derivatives; j++)
1666
values[j] = 0;
1667
1668
// Map degree of freedom to element degree of freedom
1669
const unsigned int dof = i;
1670
1671
// Generate scalings
1672
const double scalings_y_0 = 1;
1673
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1674
1675
// Compute psitilde_a
1676
const double psitilde_a_0 = 1;
1677
const double psitilde_a_1 = x;
1678
1679
// Compute psitilde_bs
1680
const double psitilde_bs_0_0 = 1;
1681
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1682
const double psitilde_bs_1_0 = 1;
1683
1684
// Compute basisvalues
1685
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1686
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1687
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1688
1689
// Table(s) of coefficients
1690
static const double coefficients0[3][3] = \
1691
{{0.471404521, -0.288675135, -0.166666667},
1692
{0.471404521, 0.288675135, -0.166666667},
1693
{0.471404521, 0, 0.333333333}};
1694
1695
// Interesting (new) part
1696
// Tables of derivatives of the polynomial base (transpose)
1697
static const double dmats0[3][3] = \
1698
{{0, 0, 0},
1699
{4.89897949, 0, 0},
1700
{0, 0, 0}};
1701
1702
static const double dmats1[3][3] = \
1703
{{0, 0, 0},
1704
{2.44948974, 0, 0},
1705
{4.24264069, 0, 0}};
1706
1707
// Compute reference derivatives
1708
// Declare pointer to array of derivatives on FIAT element
1709
double *derivatives = new double [num_derivatives];
1710
1711
// Declare coefficients
1712
double coeff0_0 = 0;
1713
double coeff0_1 = 0;
1714
double coeff0_2 = 0;
1715
1716
// Declare new coefficients
1717
double new_coeff0_0 = 0;
1718
double new_coeff0_1 = 0;
1719
double new_coeff0_2 = 0;
1720
1721
// Loop possible derivatives
1722
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1723
{
1724
// Get values from coefficients array
1725
new_coeff0_0 = coefficients0[dof][0];
1726
new_coeff0_1 = coefficients0[dof][1];
1727
new_coeff0_2 = coefficients0[dof][2];
1728
1729
// Loop derivative order
1730
for (unsigned int j = 0; j < n; j++)
1731
{
1732
// Update old coefficients
1733
coeff0_0 = new_coeff0_0;
1734
coeff0_1 = new_coeff0_1;
1735
coeff0_2 = new_coeff0_2;
1736
1737
if(combinations[deriv_num][j] == 0)
1738
{
1739
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
1740
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
1741
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
1742
}
1743
if(combinations[deriv_num][j] == 1)
1744
{
1745
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
1746
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
1747
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
1748
}
1749
1750
}
1751
// Compute derivatives on reference element as dot product of coefficients and basisvalues
1752
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
1753
}
1754
1755
// Transform derivatives back to physical element
1756
for (unsigned int row = 0; row < num_derivatives; row++)
1757
{
1758
for (unsigned int col = 0; col < num_derivatives; col++)
1759
{
1760
values[row] += transform[row][col]*derivatives[col];
1761
}
1762
}
1763
// Delete pointer to array of derivatives on FIAT element
1764
delete [] derivatives;
1765
1766
// Delete pointer to array of combinations of derivatives and transform
1767
for (unsigned int row = 0; row < num_derivatives; row++)
1768
{
1769
delete [] combinations[row];
1770
delete [] transform[row];
1771
}
1772
1773
delete [] combinations;
1774
delete [] transform;
1775
}
1776
1777
/// Evaluate order n derivatives of all basis functions at given point in cell
1778
virtual void evaluate_basis_derivatives_all(unsigned int n,
1779
double* values,
1780
const double* coordinates,
1781
const ufc::cell& c) const
1782
{
1783
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1784
}
1785
1786
/// Evaluate linear functional for dof i on the function f
1787
virtual double evaluate_dof(unsigned int i,
1788
const ufc::function& f,
1789
const ufc::cell& c) const
1790
{
1791
// The reference points, direction and weights:
1792
static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
1793
static const double W[3][1] = {{1}, {1}, {1}};
1794
static const double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
1795
1796
const double * const * x = c.coordinates;
1797
double result = 0.0;
1798
// Iterate over the points:
1799
// Evaluate basis functions for affine mapping
1800
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1801
const double w1 = X[i][0][0];
1802
const double w2 = X[i][0][1];
1803
1804
// Compute affine mapping y = F(X)
1805
double y[2];
1806
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1807
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1808
1809
// Evaluate function at physical points
1810
double values[1];
1811
f.evaluate(values, y, c);
1812
1813
// Map function values using appropriate mapping
1814
// Affine map: Do nothing
1815
1816
// Note that we do not map the weights (yet).
1817
1818
// Take directional components
1819
for(int k = 0; k < 1; k++)
1820
result += values[k]*D[i][0][k];
1821
// Multiply by weights
1822
result *= W[i][0];
1823
1824
return result;
1825
}
1826
1827
/// Evaluate linear functionals for all dofs on the function f
1828
virtual void evaluate_dofs(double* values,
1829
const ufc::function& f,
1830
const ufc::cell& c) const
1831
{
1832
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1833
}
1834
1835
/// Interpolate vertex values from dof values
1836
virtual void interpolate_vertex_values(double* vertex_values,
1837
const double* dof_values,
1838
const ufc::cell& c) const
1839
{
1840
// Evaluate at vertices and use affine mapping
1841
vertex_values[0] = dof_values[0];
1842
vertex_values[1] = dof_values[1];
1843
vertex_values[2] = dof_values[2];
1844
}
1845
1846
/// Return the number of sub elements (for a mixed element)
1847
virtual unsigned int num_sub_elements() const
1848
{
1849
return 1;
1850
}
1851
1852
/// Create a new finite element for sub element i (for a mixed element)
1853
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1854
{
1855
return new equation_1_finite_element_0();
1856
}
1857
1858
};
1859
1860
/// This class defines the interface for a finite element.
1861
1862
class equation_1_finite_element_1: public ufc::finite_element
1863
{
1864
public:
1865
1866
/// Constructor
1867
equation_1_finite_element_1() : ufc::finite_element()
1868
{
1869
// Do nothing
1870
}
1871
1872
/// Destructor
1873
virtual ~equation_1_finite_element_1()
1874
{
1875
// Do nothing
1876
}
1877
1878
/// Return a string identifying the finite element
1879
virtual const char* signature() const
1880
{
1881
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
1882
}
1883
1884
/// Return the cell shape
1885
virtual ufc::shape cell_shape() const
1886
{
1887
return ufc::triangle;
1888
}
1889
1890
/// Return the dimension of the finite element function space
1891
virtual unsigned int space_dimension() const
1892
{
1893
return 3;
1894
}
1895
1896
/// Return the rank of the value space
1897
virtual unsigned int value_rank() const
1898
{
1899
return 0;
1900
}
1901
1902
/// Return the dimension of the value space for axis i
1903
virtual unsigned int value_dimension(unsigned int i) const
1904
{
1905
return 1;
1906
}
1907
1908
/// Evaluate basis function i at given point in cell
1909
virtual void evaluate_basis(unsigned int i,
1910
double* values,
1911
const double* coordinates,
1912
const ufc::cell& c) const
1913
{
1914
// Extract vertex coordinates
1915
const double * const * element_coordinates = c.coordinates;
1916
1917
// Compute Jacobian of affine map from reference cell
1918
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1919
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1920
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1921
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1922
1923
// Compute determinant of Jacobian
1924
const double detJ = J_00*J_11 - J_01*J_10;
1925
1926
// Compute inverse of Jacobian
1927
1928
// Get coordinates and map to the reference (UFC) element
1929
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1930
element_coordinates[0][0]*element_coordinates[2][1] +\
1931
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1932
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1933
element_coordinates[1][0]*element_coordinates[0][1] -\
1934
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1935
1936
// Map coordinates to the reference square
1937
if (std::abs(y - 1.0) < 1e-08)
1938
x = -1.0;
1939
else
1940
x = 2.0 *x/(1.0 - y) - 1.0;
1941
y = 2.0*y - 1.0;
1942
1943
// Reset values
1944
*values = 0;
1945
1946
// Map degree of freedom to element degree of freedom
1947
const unsigned int dof = i;
1948
1949
// Generate scalings
1950
const double scalings_y_0 = 1;
1951
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1952
1953
// Compute psitilde_a
1954
const double psitilde_a_0 = 1;
1955
const double psitilde_a_1 = x;
1956
1957
// Compute psitilde_bs
1958
const double psitilde_bs_0_0 = 1;
1959
const double psitilde_bs_0_1 = 1.5*y + 0.5;
1960
const double psitilde_bs_1_0 = 1;
1961
1962
// Compute basisvalues
1963
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1964
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1965
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1966
1967
// Table(s) of coefficients
1968
static const double coefficients0[3][3] = \
1969
{{0.471404521, -0.288675135, -0.166666667},
1970
{0.471404521, 0.288675135, -0.166666667},
1971
{0.471404521, 0, 0.333333333}};
1972
1973
// Extract relevant coefficients
1974
const double coeff0_0 = coefficients0[dof][0];
1975
const double coeff0_1 = coefficients0[dof][1];
1976
const double coeff0_2 = coefficients0[dof][2];
1977
1978
// Compute value(s)
1979
*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
1980
}
1981
1982
/// Evaluate all basis functions at given point in cell
1983
virtual void evaluate_basis_all(double* values,
1984
const double* coordinates,
1985
const ufc::cell& c) const
1986
{
1987
throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1988
}
1989
1990
/// Evaluate order n derivatives of basis function i at given point in cell
1991
virtual void evaluate_basis_derivatives(unsigned int i,
1992
unsigned int n,
1993
double* values,
1994
const double* coordinates,
1995
const ufc::cell& c) const
1996
{
1997
// Extract vertex coordinates
1998
const double * const * element_coordinates = c.coordinates;
1999
2000
// Compute Jacobian of affine map from reference cell
2001
const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
2002
const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
2003
const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
2004
const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
2005
2006
// Compute determinant of Jacobian
2007
const double detJ = J_00*J_11 - J_01*J_10;
2008
2009
// Compute inverse of Jacobian
2010
2011
// Get coordinates and map to the reference (UFC) element
2012
double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
2013
element_coordinates[0][0]*element_coordinates[2][1] +\
2014
J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
2015
double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
2016
element_coordinates[1][0]*element_coordinates[0][1] -\
2017
J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
2018
2019
// Map coordinates to the reference square
2020
if (std::abs(y - 1.0) < 1e-08)
2021
x = -1.0;
2022
else
2023
x = 2.0 *x/(1.0 - y) - 1.0;
2024
y = 2.0*y - 1.0;
2025
2026
// Compute number of derivatives
2027
unsigned int num_derivatives = 1;
2028
2029
for (unsigned int j = 0; j < n; j++)
2030
num_derivatives *= 2;
2031
2032
2033
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
2034
unsigned int **combinations = new unsigned int *[num_derivatives];
2035
2036
for (unsigned int j = 0; j < num_derivatives; j++)
2037
{
2038
combinations[j] = new unsigned int [n];
2039
for (unsigned int k = 0; k < n; k++)
2040
combinations[j][k] = 0;
2041
}
2042
2043
// Generate combinations of derivatives
2044
for (unsigned int row = 1; row < num_derivatives; row++)
2045
{
2046
for (unsigned int num = 0; num < row; num++)
2047
{
2048
for (unsigned int col = n-1; col+1 > 0; col--)
2049
{
2050
if (combinations[row][col] + 1 > 1)
2051
combinations[row][col] = 0;
2052
else
2053
{
2054
combinations[row][col] += 1;
2055
break;
2056
}
2057
}
2058
}
2059
}
2060
2061
// Compute inverse of Jacobian
2062
const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
2063
2064
// Declare transformation matrix
2065
// Declare pointer to two dimensional array and initialise
2066
double **transform = new double *[num_derivatives];
2067
2068
for (unsigned int j = 0; j < num_derivatives; j++)
2069
{
2070
transform[j] = new double [num_derivatives];
2071
for (unsigned int k = 0; k < num_derivatives; k++)
2072
transform[j][k] = 1;
2073
}
2074
2075
// Construct transformation matrix
2076
for (unsigned int row = 0; row < num_derivatives; row++)
2077
{
2078
for (unsigned int col = 0; col < num_derivatives; col++)
2079
{
2080
for (unsigned int k = 0; k < n; k++)
2081
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2082
}
2083
}
2084
2085
// Reset values
2086
for (unsigned int j = 0; j < 1*num_derivatives; j++)
2087
values[j] = 0;
2088
2089
// Map degree of freedom to element degree of freedom
2090
const unsigned int dof = i;
2091
2092
// Generate scalings
2093
const double scalings_y_0 = 1;
2094
const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
2095
2096
// Compute psitilde_a
2097
const double psitilde_a_0 = 1;
2098
const double psitilde_a_1 = x;
2099
2100
// Compute psitilde_bs
2101
const double psitilde_bs_0_0 = 1;
2102
const double psitilde_bs_0_1 = 1.5*y + 0.5;
2103
const double psitilde_bs_1_0 = 1;
2104
2105
// Compute basisvalues
2106
const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
2107
const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
2108
const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
2109
2110
// Table(s) of coefficients
2111
static const double coefficients0[3][3] = \
2112
{{0.471404521, -0.288675135, -0.166666667},
2113
{0.471404521, 0.288675135, -0.166666667},
2114
{0.471404521, 0, 0.333333333}};
2115
2116
// Interesting (new) part
2117
// Tables of derivatives of the polynomial base (transpose)
2118
static const double dmats0[3][3] = \
2119
{{0, 0, 0},
2120
{4.89897949, 0, 0},
2121
{0, 0, 0}};
2122
2123
static const double dmats1[3][3] = \
2124
{{0, 0, 0},
2125
{2.44948974, 0, 0},
2126
{4.24264069, 0, 0}};
2127
2128
// Compute reference derivatives
2129
// Declare pointer to array of derivatives on FIAT element
2130
double *derivatives = new double [num_derivatives];
2131
2132
// Declare coefficients
2133
double coeff0_0 = 0;
2134
double coeff0_1 = 0;
2135
double coeff0_2 = 0;
2136
2137
// Declare new coefficients
2138
double new_coeff0_0 = 0;
2139
double new_coeff0_1 = 0;
2140
double new_coeff0_2 = 0;
2141
2142
// Loop possible derivatives
2143
for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
2144
{
2145
// Get values from coefficients array
2146
new_coeff0_0 = coefficients0[dof][0];
2147
new_coeff0_1 = coefficients0[dof][1];
2148
new_coeff0_2 = coefficients0[dof][2];
2149
2150
// Loop derivative order
2151
for (unsigned int j = 0; j < n; j++)
2152
{
2153
// Update old coefficients
2154
coeff0_0 = new_coeff0_0;
2155
coeff0_1 = new_coeff0_1;
2156
coeff0_2 = new_coeff0_2;
2157
2158
if(combinations[deriv_num][j] == 0)
2159
{
2160
new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
2161
new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
2162
new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
2163
}
2164
if(combinations[deriv_num][j] == 1)
2165
{
2166
new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
2167
new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
2168
new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
2169
}
2170
2171
}
2172
// Compute derivatives on reference element as dot product of coefficients and basisvalues
2173
derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
2174
}
2175
2176
// Transform derivatives back to physical element
2177
for (unsigned int row = 0; row < num_derivatives; row++)
2178
{
2179
for (unsigned int col = 0; col < num_derivatives; col++)
2180
{
2181
values[row] += transform[row][col]*derivatives[col];
2182
}
2183
}
2184
// Delete pointer to array of derivatives on FIAT element
2185
delete [] derivatives;
2186
2187
// Delete pointer to array of combinations of derivatives and transform
2188
for (unsigned int row = 0; row < num_derivatives; row++)
2189
{
2190
delete [] combinations[row];
2191
delete [] transform[row];
2192
}
2193
2194
delete [] combinations;
2195
delete [] transform;
2196
}
2197
2198
/// Evaluate order n derivatives of all basis functions at given point in cell
2199
virtual void evaluate_basis_derivatives_all(unsigned int n,
2200
double* values,
2201
const double* coordinates,
2202
const ufc::cell& c) const
2203
{
2204
throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
2205
}
2206
2207
/// Evaluate linear functional for dof i on the function f
2208
virtual double evaluate_dof(unsigned int i,
2209
const ufc::function& f,
2210
const ufc::cell& c) const
2211
{
2212
// The reference points, direction and weights:
2213
static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
2214
static const double W[3][1] = {{1}, {1}, {1}};
2215
static const double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
2216
2217
const double * const * x = c.coordinates;
2218
double result = 0.0;
2219
// Iterate over the points:
2220
// Evaluate basis functions for affine mapping
2221
const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
2222
const double w1 = X[i][0][0];
2223
const double w2 = X[i][0][1];
2224
2225
// Compute affine mapping y = F(X)
2226
double y[2];
2227
y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
2228
y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
2229
2230
// Evaluate function at physical points
2231
double values[1];
2232
f.evaluate(values, y, c);
2233
2234
// Map function values using appropriate mapping
2235
// Affine map: Do nothing
2236
2237
// Note that we do not map the weights (yet).
2238
2239
// Take directional components
2240
for(int k = 0; k < 1; k++)
2241
result += values[k]*D[i][0][k];
2242
// Multiply by weights
2243
result *= W[i][0];
2244
2245
return result;
2246
}
2247
2248
/// Evaluate linear functionals for all dofs on the function f
2249
virtual void evaluate_dofs(double* values,
2250
const ufc::function& f,
2251
const ufc::cell& c) const
2252
{
2253
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2254
}
2255
2256
/// Interpolate vertex values from dof values
2257
virtual void interpolate_vertex_values(double* vertex_values,
2258
const double* dof_values,
2259
const ufc::cell& c) const
2260
{
2261
// Evaluate at vertices and use affine mapping
2262
vertex_values[0] = dof_values[0];
2263
vertex_values[1] = dof_values[1];
2264
vertex_values[2] = dof_values[2];
2265
}
2266
2267
/// Return the number of sub elements (for a mixed element)
2268
virtual unsigned int num_sub_elements() const
2269
{
2270
return 1;
2271
}
2272
2273
/// Create a new finite element for sub element i (for a mixed element)
2274
virtual ufc::finite_element* create_sub_element(unsigned int i) const
2275
{
2276
return new equation_1_finite_element_1();
2277
}
2278
2279
};
2280
2281
/// This class defines the interface for a local-to-global mapping of
2282
/// degrees of freedom (dofs).
2283
2284
class equation_1_dof_map_0: public ufc::dof_map
2285
{
2286
private:
2287
2288
unsigned int __global_dimension;
2289
2290
public:
2291
2292
/// Constructor
2293
equation_1_dof_map_0() : ufc::dof_map()
2294
{
2295
__global_dimension = 0;
2296
}
2297
2298
/// Destructor
2299
virtual ~equation_1_dof_map_0()
2300
{
2301
// Do nothing
2302
}
2303
2304
/// Return a string identifying the dof map
2305
virtual const char* signature() const
2306
{
2307
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
2308
}
2309
2310
/// Return true iff mesh entities of topological dimension d are needed
2311
virtual bool needs_mesh_entities(unsigned int d) const
2312
{
2313
switch ( d )
2314
{
2315
case 0:
2316
return true;
2317
break;
2318
case 1:
2319
return false;
2320
break;
2321
case 2:
2322
return false;
2323
break;
2324
}
2325
return false;
2326
}
2327
2328
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2329
virtual bool init_mesh(const ufc::mesh& m)
2330
{
2331
__global_dimension = m.num_entities[0];
2332
return false;
2333
}
2334
2335
/// Initialize dof map for given cell
2336
virtual void init_cell(const ufc::mesh& m,
2337
const ufc::cell& c)
2338
{
2339
// Do nothing
2340
}
2341
2342
/// Finish initialization of dof map for cells
2343
virtual void init_cell_finalize()
2344
{
2345
// Do nothing
2346
}
2347
2348
/// Return the dimension of the global finite element function space
2349
virtual unsigned int global_dimension() const
2350
{
2351
return __global_dimension;
2352
}
2353
2354
/// Return the dimension of the local finite element function space for a cell
2355
virtual unsigned int local_dimension(const ufc::cell& c) const
2356
{
2357
return 3;
2358
}
2359
2360
/// Return the maximum dimension of the local finite element function space
2361
virtual unsigned int max_local_dimension() const
2362
{
2363
return 3;
2364
}
2365
2366
// Return the geometric dimension of the coordinates this dof map provides
2367
virtual unsigned int geometric_dimension() const
2368
{
2369
return 2;
2370
}
2371
2372
/// Return the number of dofs on each cell facet
2373
virtual unsigned int num_facet_dofs() const
2374
{
2375
return 2;
2376
}
2377
2378
/// Return the number of dofs associated with each cell entity of dimension d
2379
virtual unsigned int num_entity_dofs(unsigned int d) const
2380
{
2381
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2382
}
2383
2384
/// Tabulate the local-to-global mapping of dofs on a cell
2385
virtual void tabulate_dofs(unsigned int* dofs,
2386
const ufc::mesh& m,
2387
const ufc::cell& c) const
2388
{
2389
dofs[0] = c.entity_indices[0][0];
2390
dofs[1] = c.entity_indices[0][1];
2391
dofs[2] = c.entity_indices[0][2];
2392
}
2393
2394
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2395
virtual void tabulate_facet_dofs(unsigned int* dofs,
2396
unsigned int facet) const
2397
{
2398
switch ( facet )
2399
{
2400
case 0:
2401
dofs[0] = 1;
2402
dofs[1] = 2;
2403
break;
2404
case 1:
2405
dofs[0] = 0;
2406
dofs[1] = 2;
2407
break;
2408
case 2:
2409
dofs[0] = 0;
2410
dofs[1] = 1;
2411
break;
2412
}
2413
}
2414
2415
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2416
virtual void tabulate_entity_dofs(unsigned int* dofs,
2417
unsigned int d, unsigned int i) const
2418
{
2419
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2420
}
2421
2422
/// Tabulate the coordinates of all dofs on a cell
2423
virtual void tabulate_coordinates(double** coordinates,
2424
const ufc::cell& c) const
2425
{
2426
const double * const * x = c.coordinates;
2427
coordinates[0][0] = x[0][0];
2428
coordinates[0][1] = x[0][1];
2429
coordinates[1][0] = x[1][0];
2430
coordinates[1][1] = x[1][1];
2431
coordinates[2][0] = x[2][0];
2432
coordinates[2][1] = x[2][1];
2433
}
2434
2435
/// Return the number of sub dof maps (for a mixed element)
2436
virtual unsigned int num_sub_dof_maps() const
2437
{
2438
return 1;
2439
}
2440
2441
/// Create a new dof_map for sub dof map i (for a mixed element)
2442
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2443
{
2444
return new equation_1_dof_map_0();
2445
}
2446
2447
};
2448
2449
/// This class defines the interface for a local-to-global mapping of
2450
/// degrees of freedom (dofs).
2451
2452
class equation_1_dof_map_1: public ufc::dof_map
2453
{
2454
private:
2455
2456
unsigned int __global_dimension;
2457
2458
public:
2459
2460
/// Constructor
2461
equation_1_dof_map_1() : ufc::dof_map()
2462
{
2463
__global_dimension = 0;
2464
}
2465
2466
/// Destructor
2467
virtual ~equation_1_dof_map_1()
2468
{
2469
// Do nothing
2470
}
2471
2472
/// Return a string identifying the dof map
2473
virtual const char* signature() const
2474
{
2475
return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
2476
}
2477
2478
/// Return true iff mesh entities of topological dimension d are needed
2479
virtual bool needs_mesh_entities(unsigned int d) const
2480
{
2481
switch ( d )
2482
{
2483
case 0:
2484
return true;
2485
break;
2486
case 1:
2487
return false;
2488
break;
2489
case 2:
2490
return false;
2491
break;
2492
}
2493
return false;
2494
}
2495
2496
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2497
virtual bool init_mesh(const ufc::mesh& m)
2498
{
2499
__global_dimension = m.num_entities[0];
2500
return false;
2501
}
2502
2503
/// Initialize dof map for given cell
2504
virtual void init_cell(const ufc::mesh& m,
2505
const ufc::cell& c)
2506
{
2507
// Do nothing
2508
}
2509
2510
/// Finish initialization of dof map for cells
2511
virtual void init_cell_finalize()
2512
{
2513
// Do nothing
2514
}
2515
2516
/// Return the dimension of the global finite element function space
2517
virtual unsigned int global_dimension() const
2518
{
2519
return __global_dimension;
2520
}
2521
2522
/// Return the dimension of the local finite element function space for a cell
2523
virtual unsigned int local_dimension(const ufc::cell& c) const
2524
{
2525
return 3;
2526
}
2527
2528
/// Return the maximum dimension of the local finite element function space
2529
virtual unsigned int max_local_dimension() const
2530
{
2531
return 3;
2532
}
2533
2534
// Return the geometric dimension of the coordinates this dof map provides
2535
virtual unsigned int geometric_dimension() const
2536
{
2537
return 2;
2538
}
2539
2540
/// Return the number of dofs on each cell facet
2541
virtual unsigned int num_facet_dofs() const
2542
{
2543
return 2;
2544
}
2545
2546
/// Return the number of dofs associated with each cell entity of dimension d
2547
virtual unsigned int num_entity_dofs(unsigned int d) const
2548
{
2549
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2550
}
2551
2552
/// Tabulate the local-to-global mapping of dofs on a cell
2553
virtual void tabulate_dofs(unsigned int* dofs,
2554
const ufc::mesh& m,
2555
const ufc::cell& c) const
2556
{
2557
dofs[0] = c.entity_indices[0][0];
2558
dofs[1] = c.entity_indices[0][1];
2559
dofs[2] = c.entity_indices[0][2];
2560
}
2561
2562
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2563
virtual void tabulate_facet_dofs(unsigned int* dofs,
2564
unsigned int facet) const
2565
{
2566
switch ( facet )
2567
{
2568
case 0:
2569
dofs[0] = 1;
2570
dofs[1] = 2;
2571
break;
2572
case 1:
2573
dofs[0] = 0;
2574
dofs[1] = 2;
2575
break;
2576
case 2:
2577
dofs[0] = 0;
2578
dofs[1] = 1;
2579
break;
2580
}
2581
}
2582
2583
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2584
virtual void tabulate_entity_dofs(unsigned int* dofs,
2585
unsigned int d, unsigned int i) const
2586
{
2587
throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2588
}
2589
2590
/// Tabulate the coordinates of all dofs on a cell
2591
virtual void tabulate_coordinates(double** coordinates,
2592
const ufc::cell& c) const
2593
{
2594
const double * const * x = c.coordinates;
2595
coordinates[0][0] = x[0][0];
2596
coordinates[0][1] = x[0][1];
2597
coordinates[1][0] = x[1][0];
2598
coordinates[1][1] = x[1][1];
2599
coordinates[2][0] = x[2][0];
2600
coordinates[2][1] = x[2][1];
2601
}
2602
2603
/// Return the number of sub dof maps (for a mixed element)
2604
virtual unsigned int num_sub_dof_maps() const
2605
{
2606
return 1;
2607
}
2608
2609
/// Create a new dof_map for sub dof map i (for a mixed element)
2610
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2611
{
2612
return new equation_1_dof_map_1();
2613
}
2614
2615
};
2616
2617
/// This class defines the interface for the tabulation of the cell
2618
/// tensor corresponding to the local contribution to a form from
2619
/// the integral over a cell.
2620
2621
class equation_1_cell_integral_0_quadrature: public ufc::cell_integral
2622
{
2623
public:
2624
2625
/// Constructor
2626
equation_1_cell_integral_0_quadrature() : ufc::cell_integral()
2627
{
2628
// Do nothing
2629
}
2630
2631
/// Destructor
2632
virtual ~equation_1_cell_integral_0_quadrature()
2633
{
2634
// Do nothing
2635
}
2636
2637
/// Tabulate the tensor for the contribution from a local cell
2638
virtual void tabulate_tensor(double* A,
2639
const double * const * w,
2640
const ufc::cell& c) const
2641
{
2642
// Extract vertex coordinates
2643
const double * const * x = c.coordinates;
2644
2645
// Compute Jacobian of affine map from reference cell
2646
const double J_00 = x[1][0] - x[0][0];
2647
const double J_01 = x[2][0] - x[0][0];
2648
const double J_10 = x[1][1] - x[0][1];
2649
const double J_11 = x[2][1] - x[0][1];
2650
2651
// Compute determinant of Jacobian
2652
double detJ = J_00*J_11 - J_01*J_10;
2653
2654
// Compute inverse of Jacobian
2655
const double Jinv_00 = J_11 / detJ;
2656
const double Jinv_01 = -J_01 / detJ;
2657
const double Jinv_10 = -J_10 / detJ;
2658
const double Jinv_11 = J_00 / detJ;
2659
2660
// Set scale factor
2661
const double det = std::abs(detJ);
2662
2663
2664
// Array of quadrature weights
2665
static const double W4[4] = {0.159020691, 0.0909793091, 0.159020691, 0.0909793091};
2666
// Quadrature points on the UFC reference element: (0.178558728, 0.155051026), (0.0750311102, 0.644948974), (0.666390246, 0.155051026), (0.280019915, 0.644948974)
2667
2668
// Value of basis functions at quadrature points.
2669
static const double FE0[4][3] = \
2670
{{0.666390246, 0.178558728, 0.155051026},
2671
{0.280019915, 0.0750311102, 0.644948974},
2672
{0.178558728, 0.666390246, 0.155051026},
2673
{0.0750311102, 0.280019915, 0.644948974}};
2674
2675
static const double FE0_D01[4][3] = \
2676
{{-1, 0, 1},
2677
{-1, 0, 1},
2678
{-1, 0, 1},
2679
{-1, 0, 1}};
2680
2681
static const double FE0_D10[4][3] = \
2682
{{-1, 1, 0},
2683
{-1, 1, 0},
2684
{-1, 1, 0},
2685
{-1, 1, 0}};
2686
2687
2688
// Compute element tensor using UFL quadrature representation
2689
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
2690
// Total number of operations to compute element tensor: 372
2691
2692
// Loop quadrature points for integral
2693
// Number of operations to compute element tensor for following IP loop = 372
2694
for (unsigned int ip = 0; ip < 4; ip++)
2695
{
2696
2697
// Function declarations
2698
double F0 = 0;
2699
double F1 = 0;
2700
double F2 = 0;
2701
2702
// Total number of operations to compute function values = 18
2703
for (unsigned int r = 0; r < 3; r++)
2704
{
2705
F0 += FE0[ip][r]*w[0][r];
2706
F1 += FE0_D10[ip][r]*w[0][r];
2707
F2 += FE0_D01[ip][r]*w[0][r];
2708
}// end loop over 'r'
2709
2710
// Number of operations for primary indices: 75
2711
for (unsigned int j = 0; j < 3; j++)
2712
{
2713
// Number of operations to compute entry: 25
2714
A[j] += (FE0[ip][j]*-1*F0 + ((Jinv_00*FE0_D10[ip][j] + Jinv_10*FE0_D01[ip][j])*0.5*(Jinv_00*F1 + Jinv_10*F2) + (Jinv_01*FE0_D10[ip][j] + Jinv_11*FE0_D01[ip][j])*0.5*(Jinv_01*F1 + Jinv_11*F2))*0.1)*-1*W4[ip]*det;
2715
}// end loop over 'j'
2716
}// end loop over 'ip'
2717
}
2718
2719
};
2720
2721
/// This class defines the interface for the tabulation of the cell
2722
/// tensor corresponding to the local contribution to a form from
2723
/// the integral over a cell.
2724
2725
class equation_1_cell_integral_0: public ufc::cell_integral
2726
{
2727
private:
2728
2729
equation_1_cell_integral_0_quadrature integral_0_quadrature;
2730
2731
public:
2732
2733
/// Constructor
2734
equation_1_cell_integral_0() : ufc::cell_integral()
2735
{
2736
// Do nothing
2737
}
2738
2739
/// Destructor
2740
virtual ~equation_1_cell_integral_0()
2741
{
2742
// Do nothing
2743
}
2744
2745
/// Tabulate the tensor for the contribution from a local cell
2746
virtual void tabulate_tensor(double* A,
2747
const double * const * w,
2748
const ufc::cell& c) const
2749
{
2750
// Reset values of the element tensor block
2751
for (unsigned int j = 0; j < 3; j++)
2752
A[j] = 0;
2753
2754
// Add all contributions to element tensor
2755
integral_0_quadrature.tabulate_tensor(A, w, c);
2756
}
2757
2758
};
2759
2760
/// This class defines the interface for the assembly of the global
2761
/// tensor corresponding to a form with r + n arguments, that is, a
2762
/// mapping
2763
///
2764
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
2765
///
2766
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
2767
/// global tensor A is defined by
2768
///
2769
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
2770
///
2771
/// where each argument Vj represents the application to the
2772
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
2773
/// fixed functions (coefficients).
2774
2775
class equation_form_1: public ufc::form
2776
{
2777
public:
2778
2779
/// Constructor
2780
equation_form_1() : ufc::form()
2781
{
2782
// Do nothing
2783
}
2784
2785
/// Destructor
2786
virtual ~equation_form_1()
2787
{
2788
// Do nothing
2789
}
2790
2791
/// Return a string identifying the form
2792
virtual const char* signature() const
2793
{
2794
return "Form([Integral(Product(IntValue(-1, (), (), {}), Sum(Product(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Product(IntValue(-1, (), (), {}), Function(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0))), Product(FloatValue(0.10000000000000001, (), (), {}), IndexSum(Product(Indexed(ComponentTensor(Product(FloatValue(0.5, (), (), {}), SpatialDerivative(Function(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(0),), {Index(0): 2}))), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2}))))), Measure('cell', 0, None))])";
2795
}
2796
2797
/// Return the rank of the global tensor (r)
2798
virtual unsigned int rank() const
2799
{
2800
return 1;
2801
}
2802
2803
/// Return the number of coefficients (n)
2804
virtual unsigned int num_coefficients() const
2805
{
2806
return 1;
2807
}
2808
2809
/// Return the number of cell integrals
2810
virtual unsigned int num_cell_integrals() const
2811
{
2812
return 1;
2813
}
2814
2815
/// Return the number of exterior facet integrals
2816
virtual unsigned int num_exterior_facet_integrals() const
2817
{
2818
return 0;
2819
}
2820
2821
/// Return the number of interior facet integrals
2822
virtual unsigned int num_interior_facet_integrals() const
2823
{
2824
return 0;
2825
}
2826
2827
/// Create a new finite element for argument function i
2828
virtual ufc::finite_element* create_finite_element(unsigned int i) const
2829
{
2830
switch ( i )
2831
{
2832
case 0:
2833
return new equation_1_finite_element_0();
2834
break;
2835
case 1:
2836
return new equation_1_finite_element_1();
2837
break;
2838
}
2839
return 0;
2840
}
2841
2842
/// Create a new dof map for argument function i
2843
virtual ufc::dof_map* create_dof_map(unsigned int i) const
2844
{
2845
switch ( i )
2846
{
2847
case 0:
2848
return new equation_1_dof_map_0();
2849
break;
2850
case 1:
2851
return new equation_1_dof_map_1();
2852
break;
2853
}
2854
return 0;
2855
}
2856
2857
/// Create a new cell integral on sub domain i
2858
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
2859
{
2860
return new equation_1_cell_integral_0();
2861
}
2862
2863
/// Create a new exterior facet integral on sub domain i
2864
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
2865
{
2866
return 0;
2867
}
2868
2869
/// Create a new interior facet integral on sub domain i
2870
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
2871
{
2872
return 0;
2873
}
2874
2875
};
2876
2877
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1