~ubuntu-branches/ubuntu/vivid/ffc/vivid
» Revision
18
: test/regression/references/r_quadrature/PointMeasure.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_quadrature/PointMeasure.h
- browse files at revision 18
- compare with another revision
- download diff
- download tarball
- view history from revision 18
- Committer: Package Import Robot
- Author(s): Johannes Ring
- Date: 2014-01-10 13:56:45 UTC
- mfrom: (1.1.14)
- Revision ID: package-import@ubuntu.com-20140110135645-4ozcd71y1oggj44z
Tags: 1.3.0-1
* New upstream release.
* debian/watch: Update URL for move to Bitbucket.
* debian/docs: README -> README.rst and remove TODO.
* debian/control:
- Add python-numpy to Build-Depends.
- Replace python-all with python-all-dev in Build-Depends.
- Add ${shlibs:Depends} to Depends.
- Change to Architecture: any.
- Bump Standards-Version to 3.9.5 (no changes needed).
* debian/rules: Call dh_numpy in override_dh_python2.
* debian/watch: Update URL for move to Bitbucket.
* debian/docs: README -> README.rst and remove TODO.
* debian/control:
- Add python-numpy to Build-Depends.
- Replace python-all with python-all-dev in Build-Depends.
- Add ${shlibs:Depends} to Depends.
- Change to Architecture: any.
- Bump Standards-Version to 3.9.5 (no changes needed).
* debian/rules: Call dh_numpy in override_dh_python2.
- files added:
- .bzrignore
- ffc/dolfin
- ffc/ext
- test/regression/scripts
- test/unit/elements
- files removed:
- README
- doc/sphinx/source/_static
- doc/sphinx/source/_templates
- sandbox
- sandbox/features
- sandbox/jit
- sandbox/multidomain
- sandbox/swig
- test/regression/references
- test/regression/references/output
- test/regression/references/r_auto
- test/regression/references/r_quadrature
- test/regression/references/r_quadrature_O
- files modified:
- ChangeLog
added
removed
test/regression/references/r_quadrature/PointMeasure.h
1
// This code conforms with the UFC specification version 2.2.0
2
// and was automatically generated by FFC version 1.2.0.
3
//
4
// This code was generated with the following parameters:
5
//
6
// cache_dir: ''
7
// convert_exceptions_to_warnings: True
8
// cpp_optimize: False
9
// cpp_optimize_flags: '-O2'
10
// epsilon: 1e-14
11
// error_control: False
12
// form_postfix: True
13
// format: 'ufc'
14
// log_level: 20
15
// log_prefix: ''
16
// optimize: False
17
// output_dir: '.'
18
// precision: '8'
19
// quadrature_degree: 'auto'
20
// quadrature_rule: 'auto'
21
// representation: 'quadrature'
22
// split: False
23
// swig_binary: 'swig'
24
// swig_path: ''
25
26
#ifndef __POINTMEASURE_H
27
#define __POINTMEASURE_H
28
29
#include <cmath>
30
#include <stdexcept>
31
#include <fstream>
32
#include <ufc.h>
33
34
/// This class defines the interface for a finite element.
35
36
class pointmeasure_finite_element_0: public ufc::finite_element
37
{
38
public:
39
40
/// Constructor
41
pointmeasure_finite_element_0() : ufc::finite_element()
42
{
43
// Do nothing
44
}
45
46
/// Destructor
47
virtual ~pointmeasure_finite_element_0()
48
{
49
// Do nothing
50
}
51
52
/// Return a string identifying the finite element
53
virtual const char* signature() const
54
{
55
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
56
}
57
58
/// Return the cell shape
59
virtual ufc::shape cell_shape() const
60
{
61
return ufc::triangle;
62
}
63
64
/// Return the topological dimension of the cell shape
65
virtual std::size_t topological_dimension() const
66
{
67
return 2;
68
}
69
70
/// Return the geometric dimension of the cell shape
71
virtual std::size_t geometric_dimension() const
72
{
73
return 2;
74
}
75
76
/// Return the dimension of the finite element function space
77
virtual std::size_t space_dimension() const
78
{
79
return 6;
80
}
81
82
/// Return the rank of the value space
83
virtual std::size_t value_rank() const
84
{
85
return 0;
86
}
87
88
/// Return the dimension of the value space for axis i
89
virtual std::size_t value_dimension(std::size_t i) const
90
{
91
return 1;
92
}
93
94
/// Evaluate basis function i at given point x in cell
95
virtual void evaluate_basis(std::size_t i,
96
double* values,
97
const double* x,
98
const double* vertex_coordinates,
99
int cell_orientation) const
100
{
101
// Compute Jacobian
102
double J[4];
103
compute_jacobian_triangle_2d(J, vertex_coordinates);
104
105
// Compute Jacobian inverse and determinant
106
double K[4];
107
double detJ;
108
compute_jacobian_inverse_triangle_2d(K, detJ, J);
109
110
111
// Compute constants
112
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
113
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
114
115
// Get coordinates and map to the reference (FIAT) element
116
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
117
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
118
119
// Reset values
120
*values = 0.0;
121
switch (i)
122
{
123
case 0:
124
{
125
126
// Array of basisvalues
127
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
128
129
// Declare helper variables
130
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
131
double tmp1 = (1.0 - Y)/2.0;
132
double tmp2 = tmp1*tmp1;
133
134
// Compute basisvalues
135
basisvalues[0] = 1.0;
136
basisvalues[1] = tmp0;
137
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
138
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
139
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
140
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
141
basisvalues[0] *= std::sqrt(0.5);
142
basisvalues[2] *= std::sqrt(1.0);
143
basisvalues[5] *= std::sqrt(1.5);
144
basisvalues[1] *= std::sqrt(3.0);
145
basisvalues[4] *= std::sqrt(4.5);
146
basisvalues[3] *= std::sqrt(7.5);
147
148
// Table(s) of coefficients
149
static const double coefficients0[6] = \
150
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
151
152
// Compute value(s)
153
for (unsigned int r = 0; r < 6; r++)
154
{
155
*values += coefficients0[r]*basisvalues[r];
156
}// end loop over 'r'
157
break;
158
}
159
case 1:
160
{
161
162
// Array of basisvalues
163
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
164
165
// Declare helper variables
166
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
167
double tmp1 = (1.0 - Y)/2.0;
168
double tmp2 = tmp1*tmp1;
169
170
// Compute basisvalues
171
basisvalues[0] = 1.0;
172
basisvalues[1] = tmp0;
173
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
174
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
175
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
176
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
177
basisvalues[0] *= std::sqrt(0.5);
178
basisvalues[2] *= std::sqrt(1.0);
179
basisvalues[5] *= std::sqrt(1.5);
180
basisvalues[1] *= std::sqrt(3.0);
181
basisvalues[4] *= std::sqrt(4.5);
182
basisvalues[3] *= std::sqrt(7.5);
183
184
// Table(s) of coefficients
185
static const double coefficients0[6] = \
186
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
187
188
// Compute value(s)
189
for (unsigned int r = 0; r < 6; r++)
190
{
191
*values += coefficients0[r]*basisvalues[r];
192
}// end loop over 'r'
193
break;
194
}
195
case 2:
196
{
197
198
// Array of basisvalues
199
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
200
201
// Declare helper variables
202
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
203
double tmp1 = (1.0 - Y)/2.0;
204
double tmp2 = tmp1*tmp1;
205
206
// Compute basisvalues
207
basisvalues[0] = 1.0;
208
basisvalues[1] = tmp0;
209
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
210
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
211
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
212
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
213
basisvalues[0] *= std::sqrt(0.5);
214
basisvalues[2] *= std::sqrt(1.0);
215
basisvalues[5] *= std::sqrt(1.5);
216
basisvalues[1] *= std::sqrt(3.0);
217
basisvalues[4] *= std::sqrt(4.5);
218
basisvalues[3] *= std::sqrt(7.5);
219
220
// Table(s) of coefficients
221
static const double coefficients0[6] = \
222
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
223
224
// Compute value(s)
225
for (unsigned int r = 0; r < 6; r++)
226
{
227
*values += coefficients0[r]*basisvalues[r];
228
}// end loop over 'r'
229
break;
230
}
231
case 3:
232
{
233
234
// Array of basisvalues
235
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
236
237
// Declare helper variables
238
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
239
double tmp1 = (1.0 - Y)/2.0;
240
double tmp2 = tmp1*tmp1;
241
242
// Compute basisvalues
243
basisvalues[0] = 1.0;
244
basisvalues[1] = tmp0;
245
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
246
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
247
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
248
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
249
basisvalues[0] *= std::sqrt(0.5);
250
basisvalues[2] *= std::sqrt(1.0);
251
basisvalues[5] *= std::sqrt(1.5);
252
basisvalues[1] *= std::sqrt(3.0);
253
basisvalues[4] *= std::sqrt(4.5);
254
basisvalues[3] *= std::sqrt(7.5);
255
256
// Table(s) of coefficients
257
static const double coefficients0[6] = \
258
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
259
260
// Compute value(s)
261
for (unsigned int r = 0; r < 6; r++)
262
{
263
*values += coefficients0[r]*basisvalues[r];
264
}// end loop over 'r'
265
break;
266
}
267
case 4:
268
{
269
270
// Array of basisvalues
271
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
272
273
// Declare helper variables
274
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
275
double tmp1 = (1.0 - Y)/2.0;
276
double tmp2 = tmp1*tmp1;
277
278
// Compute basisvalues
279
basisvalues[0] = 1.0;
280
basisvalues[1] = tmp0;
281
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
282
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
283
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
284
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
285
basisvalues[0] *= std::sqrt(0.5);
286
basisvalues[2] *= std::sqrt(1.0);
287
basisvalues[5] *= std::sqrt(1.5);
288
basisvalues[1] *= std::sqrt(3.0);
289
basisvalues[4] *= std::sqrt(4.5);
290
basisvalues[3] *= std::sqrt(7.5);
291
292
// Table(s) of coefficients
293
static const double coefficients0[6] = \
294
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
295
296
// Compute value(s)
297
for (unsigned int r = 0; r < 6; r++)
298
{
299
*values += coefficients0[r]*basisvalues[r];
300
}// end loop over 'r'
301
break;
302
}
303
case 5:
304
{
305
306
// Array of basisvalues
307
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
308
309
// Declare helper variables
310
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
311
double tmp1 = (1.0 - Y)/2.0;
312
double tmp2 = tmp1*tmp1;
313
314
// Compute basisvalues
315
basisvalues[0] = 1.0;
316
basisvalues[1] = tmp0;
317
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
318
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
319
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
320
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
321
basisvalues[0] *= std::sqrt(0.5);
322
basisvalues[2] *= std::sqrt(1.0);
323
basisvalues[5] *= std::sqrt(1.5);
324
basisvalues[1] *= std::sqrt(3.0);
325
basisvalues[4] *= std::sqrt(4.5);
326
basisvalues[3] *= std::sqrt(7.5);
327
328
// Table(s) of coefficients
329
static const double coefficients0[6] = \
330
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
331
332
// Compute value(s)
333
for (unsigned int r = 0; r < 6; r++)
334
{
335
*values += coefficients0[r]*basisvalues[r];
336
}// end loop over 'r'
337
break;
338
}
339
}
340
341
}
342
343
/// Evaluate all basis functions at given point x in cell
344
virtual void evaluate_basis_all(double* values,
345
const double* x,
346
const double* vertex_coordinates,
347
int cell_orientation) const
348
{
349
// Helper variable to hold values of a single dof.
350
double dof_values = 0.0;
351
352
// Loop dofs and call evaluate_basis
353
for (unsigned int r = 0; r < 6; r++)
354
{
355
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
356
values[r] = dof_values;
357
}// end loop over 'r'
358
}
359
360
/// Evaluate order n derivatives of basis function i at given point x in cell
361
virtual void evaluate_basis_derivatives(std::size_t i,
362
std::size_t n,
363
double* values,
364
const double* x,
365
const double* vertex_coordinates,
366
int cell_orientation) const
367
{
368
// Compute Jacobian
369
double J[4];
370
compute_jacobian_triangle_2d(J, vertex_coordinates);
371
372
// Compute Jacobian inverse and determinant
373
double K[4];
374
double detJ;
375
compute_jacobian_inverse_triangle_2d(K, detJ, J);
376
377
378
// Compute constants
379
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
380
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
381
382
// Get coordinates and map to the reference (FIAT) element
383
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
384
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
385
386
// Compute number of derivatives.
387
unsigned int num_derivatives = 1;
388
for (unsigned int r = 0; r < n; r++)
389
{
390
num_derivatives *= 2;
391
}// end loop over 'r'
392
393
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
394
unsigned int **combinations = new unsigned int *[num_derivatives];
395
for (unsigned int row = 0; row < num_derivatives; row++)
396
{
397
combinations[row] = new unsigned int [n];
398
for (unsigned int col = 0; col < n; col++)
399
combinations[row][col] = 0;
400
}
401
402
// Generate combinations of derivatives
403
for (unsigned int row = 1; row < num_derivatives; row++)
404
{
405
for (unsigned int num = 0; num < row; num++)
406
{
407
for (unsigned int col = n-1; col+1 > 0; col--)
408
{
409
if (combinations[row][col] + 1 > 1)
410
combinations[row][col] = 0;
411
else
412
{
413
combinations[row][col] += 1;
414
break;
415
}
416
}
417
}
418
}
419
420
// Compute inverse of Jacobian
421
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
422
423
// Declare transformation matrix
424
// Declare pointer to two dimensional array and initialise
425
double **transform = new double *[num_derivatives];
426
427
for (unsigned int j = 0; j < num_derivatives; j++)
428
{
429
transform[j] = new double [num_derivatives];
430
for (unsigned int k = 0; k < num_derivatives; k++)
431
transform[j][k] = 1;
432
}
433
434
// Construct transformation matrix
435
for (unsigned int row = 0; row < num_derivatives; row++)
436
{
437
for (unsigned int col = 0; col < num_derivatives; col++)
438
{
439
for (unsigned int k = 0; k < n; k++)
440
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
441
}
442
}
443
444
// Reset values. Assuming that values is always an array.
445
for (unsigned int r = 0; r < num_derivatives; r++)
446
{
447
values[r] = 0.0;
448
}// end loop over 'r'
449
450
switch (i)
451
{
452
case 0:
453
{
454
455
// Array of basisvalues
456
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
457
458
// Declare helper variables
459
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
460
double tmp1 = (1.0 - Y)/2.0;
461
double tmp2 = tmp1*tmp1;
462
463
// Compute basisvalues
464
basisvalues[0] = 1.0;
465
basisvalues[1] = tmp0;
466
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
467
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
468
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
469
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
470
basisvalues[0] *= std::sqrt(0.5);
471
basisvalues[2] *= std::sqrt(1.0);
472
basisvalues[5] *= std::sqrt(1.5);
473
basisvalues[1] *= std::sqrt(3.0);
474
basisvalues[4] *= std::sqrt(4.5);
475
basisvalues[3] *= std::sqrt(7.5);
476
477
// Table(s) of coefficients
478
static const double coefficients0[6] = \
479
{0.0, -0.17320508, -0.1, 0.12171612, 0.094280904, 0.054433105};
480
481
// Tables of derivatives of the polynomial base (transpose).
482
static const double dmats0[6][6] = \
483
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
484
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
485
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
486
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
487
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
488
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
489
490
static const double dmats1[6][6] = \
491
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
492
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
493
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
494
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
495
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
496
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
497
498
// Compute reference derivatives.
499
// Declare pointer to array of derivatives on FIAT element.
500
double *derivatives = new double[num_derivatives];
501
for (unsigned int r = 0; r < num_derivatives; r++)
502
{
503
derivatives[r] = 0.0;
504
}// end loop over 'r'
505
506
// Declare derivative matrix (of polynomial basis).
507
double dmats[6][6] = \
508
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
509
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
510
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
511
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
512
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
513
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
514
515
// Declare (auxiliary) derivative matrix (of polynomial basis).
516
double dmats_old[6][6] = \
517
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
518
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
519
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
520
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
521
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
522
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
523
524
// Loop possible derivatives.
525
for (unsigned int r = 0; r < num_derivatives; r++)
526
{
527
// Resetting dmats values to compute next derivative.
528
for (unsigned int t = 0; t < 6; t++)
529
{
530
for (unsigned int u = 0; u < 6; u++)
531
{
532
dmats[t][u] = 0.0;
533
if (t == u)
534
{
535
dmats[t][u] = 1.0;
536
}
537
538
}// end loop over 'u'
539
}// end loop over 't'
540
541
// Looping derivative order to generate dmats.
542
for (unsigned int s = 0; s < n; s++)
543
{
544
// Updating dmats_old with new values and resetting dmats.
545
for (unsigned int t = 0; t < 6; t++)
546
{
547
for (unsigned int u = 0; u < 6; u++)
548
{
549
dmats_old[t][u] = dmats[t][u];
550
dmats[t][u] = 0.0;
551
}// end loop over 'u'
552
}// end loop over 't'
553
554
// Update dmats using an inner product.
555
if (combinations[r][s] == 0)
556
{
557
for (unsigned int t = 0; t < 6; t++)
558
{
559
for (unsigned int u = 0; u < 6; u++)
560
{
561
for (unsigned int tu = 0; tu < 6; tu++)
562
{
563
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
564
}// end loop over 'tu'
565
}// end loop over 'u'
566
}// end loop over 't'
567
}
568
569
if (combinations[r][s] == 1)
570
{
571
for (unsigned int t = 0; t < 6; t++)
572
{
573
for (unsigned int u = 0; u < 6; u++)
574
{
575
for (unsigned int tu = 0; tu < 6; tu++)
576
{
577
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
578
}// end loop over 'tu'
579
}// end loop over 'u'
580
}// end loop over 't'
581
}
582
583
}// end loop over 's'
584
for (unsigned int s = 0; s < 6; s++)
585
{
586
for (unsigned int t = 0; t < 6; t++)
587
{
588
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
589
}// end loop over 't'
590
}// end loop over 's'
591
}// end loop over 'r'
592
593
// Transform derivatives back to physical element
594
for (unsigned int r = 0; r < num_derivatives; r++)
595
{
596
for (unsigned int s = 0; s < num_derivatives; s++)
597
{
598
values[r] += transform[r][s]*derivatives[s];
599
}// end loop over 's'
600
}// end loop over 'r'
601
602
// Delete pointer to array of derivatives on FIAT element
603
delete [] derivatives;
604
605
// Delete pointer to array of combinations of derivatives and transform
606
for (unsigned int r = 0; r < num_derivatives; r++)
607
{
608
delete [] combinations[r];
609
}// end loop over 'r'
610
delete [] combinations;
611
for (unsigned int r = 0; r < num_derivatives; r++)
612
{
613
delete [] transform[r];
614
}// end loop over 'r'
615
delete [] transform;
616
break;
617
}
618
case 1:
619
{
620
621
// Array of basisvalues
622
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
623
624
// Declare helper variables
625
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
626
double tmp1 = (1.0 - Y)/2.0;
627
double tmp2 = tmp1*tmp1;
628
629
// Compute basisvalues
630
basisvalues[0] = 1.0;
631
basisvalues[1] = tmp0;
632
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
633
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
634
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
635
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
636
basisvalues[0] *= std::sqrt(0.5);
637
basisvalues[2] *= std::sqrt(1.0);
638
basisvalues[5] *= std::sqrt(1.5);
639
basisvalues[1] *= std::sqrt(3.0);
640
basisvalues[4] *= std::sqrt(4.5);
641
basisvalues[3] *= std::sqrt(7.5);
642
643
// Table(s) of coefficients
644
static const double coefficients0[6] = \
645
{0.0, 0.17320508, -0.1, 0.12171612, -0.094280904, 0.054433105};
646
647
// Tables of derivatives of the polynomial base (transpose).
648
static const double dmats0[6][6] = \
649
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
650
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
651
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
652
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
653
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
654
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
655
656
static const double dmats1[6][6] = \
657
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
658
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
659
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
660
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
661
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
662
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
663
664
// Compute reference derivatives.
665
// Declare pointer to array of derivatives on FIAT element.
666
double *derivatives = new double[num_derivatives];
667
for (unsigned int r = 0; r < num_derivatives; r++)
668
{
669
derivatives[r] = 0.0;
670
}// end loop over 'r'
671
672
// Declare derivative matrix (of polynomial basis).
673
double dmats[6][6] = \
674
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
675
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
676
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
677
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
678
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
679
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
680
681
// Declare (auxiliary) derivative matrix (of polynomial basis).
682
double dmats_old[6][6] = \
683
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
684
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
685
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
686
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
687
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
688
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
689
690
// Loop possible derivatives.
691
for (unsigned int r = 0; r < num_derivatives; r++)
692
{
693
// Resetting dmats values to compute next derivative.
694
for (unsigned int t = 0; t < 6; t++)
695
{
696
for (unsigned int u = 0; u < 6; u++)
697
{
698
dmats[t][u] = 0.0;
699
if (t == u)
700
{
701
dmats[t][u] = 1.0;
702
}
703
704
}// end loop over 'u'
705
}// end loop over 't'
706
707
// Looping derivative order to generate dmats.
708
for (unsigned int s = 0; s < n; s++)
709
{
710
// Updating dmats_old with new values and resetting dmats.
711
for (unsigned int t = 0; t < 6; t++)
712
{
713
for (unsigned int u = 0; u < 6; u++)
714
{
715
dmats_old[t][u] = dmats[t][u];
716
dmats[t][u] = 0.0;
717
}// end loop over 'u'
718
}// end loop over 't'
719
720
// Update dmats using an inner product.
721
if (combinations[r][s] == 0)
722
{
723
for (unsigned int t = 0; t < 6; t++)
724
{
725
for (unsigned int u = 0; u < 6; u++)
726
{
727
for (unsigned int tu = 0; tu < 6; tu++)
728
{
729
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
730
}// end loop over 'tu'
731
}// end loop over 'u'
732
}// end loop over 't'
733
}
734
735
if (combinations[r][s] == 1)
736
{
737
for (unsigned int t = 0; t < 6; t++)
738
{
739
for (unsigned int u = 0; u < 6; u++)
740
{
741
for (unsigned int tu = 0; tu < 6; tu++)
742
{
743
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
744
}// end loop over 'tu'
745
}// end loop over 'u'
746
}// end loop over 't'
747
}
748
749
}// end loop over 's'
750
for (unsigned int s = 0; s < 6; s++)
751
{
752
for (unsigned int t = 0; t < 6; t++)
753
{
754
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
755
}// end loop over 't'
756
}// end loop over 's'
757
}// end loop over 'r'
758
759
// Transform derivatives back to physical element
760
for (unsigned int r = 0; r < num_derivatives; r++)
761
{
762
for (unsigned int s = 0; s < num_derivatives; s++)
763
{
764
values[r] += transform[r][s]*derivatives[s];
765
}// end loop over 's'
766
}// end loop over 'r'
767
768
// Delete pointer to array of derivatives on FIAT element
769
delete [] derivatives;
770
771
// Delete pointer to array of combinations of derivatives and transform
772
for (unsigned int r = 0; r < num_derivatives; r++)
773
{
774
delete [] combinations[r];
775
}// end loop over 'r'
776
delete [] combinations;
777
for (unsigned int r = 0; r < num_derivatives; r++)
778
{
779
delete [] transform[r];
780
}// end loop over 'r'
781
delete [] transform;
782
break;
783
}
784
case 2:
785
{
786
787
// Array of basisvalues
788
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
789
790
// Declare helper variables
791
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
792
double tmp1 = (1.0 - Y)/2.0;
793
double tmp2 = tmp1*tmp1;
794
795
// Compute basisvalues
796
basisvalues[0] = 1.0;
797
basisvalues[1] = tmp0;
798
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
799
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
800
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
801
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
802
basisvalues[0] *= std::sqrt(0.5);
803
basisvalues[2] *= std::sqrt(1.0);
804
basisvalues[5] *= std::sqrt(1.5);
805
basisvalues[1] *= std::sqrt(3.0);
806
basisvalues[4] *= std::sqrt(4.5);
807
basisvalues[3] *= std::sqrt(7.5);
808
809
// Table(s) of coefficients
810
static const double coefficients0[6] = \
811
{0.0, 0.0, 0.2, 0.0, 0.0, 0.16329932};
812
813
// Tables of derivatives of the polynomial base (transpose).
814
static const double dmats0[6][6] = \
815
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
816
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
817
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
818
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
819
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
820
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
821
822
static const double dmats1[6][6] = \
823
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
824
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
825
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
826
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
827
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
828
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
829
830
// Compute reference derivatives.
831
// Declare pointer to array of derivatives on FIAT element.
832
double *derivatives = new double[num_derivatives];
833
for (unsigned int r = 0; r < num_derivatives; r++)
834
{
835
derivatives[r] = 0.0;
836
}// end loop over 'r'
837
838
// Declare derivative matrix (of polynomial basis).
839
double dmats[6][6] = \
840
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
841
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
842
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
843
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
844
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
845
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
846
847
// Declare (auxiliary) derivative matrix (of polynomial basis).
848
double dmats_old[6][6] = \
849
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
850
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
851
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
852
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
853
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
854
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
855
856
// Loop possible derivatives.
857
for (unsigned int r = 0; r < num_derivatives; r++)
858
{
859
// Resetting dmats values to compute next derivative.
860
for (unsigned int t = 0; t < 6; t++)
861
{
862
for (unsigned int u = 0; u < 6; u++)
863
{
864
dmats[t][u] = 0.0;
865
if (t == u)
866
{
867
dmats[t][u] = 1.0;
868
}
869
870
}// end loop over 'u'
871
}// end loop over 't'
872
873
// Looping derivative order to generate dmats.
874
for (unsigned int s = 0; s < n; s++)
875
{
876
// Updating dmats_old with new values and resetting dmats.
877
for (unsigned int t = 0; t < 6; t++)
878
{
879
for (unsigned int u = 0; u < 6; u++)
880
{
881
dmats_old[t][u] = dmats[t][u];
882
dmats[t][u] = 0.0;
883
}// end loop over 'u'
884
}// end loop over 't'
885
886
// Update dmats using an inner product.
887
if (combinations[r][s] == 0)
888
{
889
for (unsigned int t = 0; t < 6; t++)
890
{
891
for (unsigned int u = 0; u < 6; u++)
892
{
893
for (unsigned int tu = 0; tu < 6; tu++)
894
{
895
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
896
}// end loop over 'tu'
897
}// end loop over 'u'
898
}// end loop over 't'
899
}
900
901
if (combinations[r][s] == 1)
902
{
903
for (unsigned int t = 0; t < 6; t++)
904
{
905
for (unsigned int u = 0; u < 6; u++)
906
{
907
for (unsigned int tu = 0; tu < 6; tu++)
908
{
909
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
910
}// end loop over 'tu'
911
}// end loop over 'u'
912
}// end loop over 't'
913
}
914
915
}// end loop over 's'
916
for (unsigned int s = 0; s < 6; s++)
917
{
918
for (unsigned int t = 0; t < 6; t++)
919
{
920
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
921
}// end loop over 't'
922
}// end loop over 's'
923
}// end loop over 'r'
924
925
// Transform derivatives back to physical element
926
for (unsigned int r = 0; r < num_derivatives; r++)
927
{
928
for (unsigned int s = 0; s < num_derivatives; s++)
929
{
930
values[r] += transform[r][s]*derivatives[s];
931
}// end loop over 's'
932
}// end loop over 'r'
933
934
// Delete pointer to array of derivatives on FIAT element
935
delete [] derivatives;
936
937
// Delete pointer to array of combinations of derivatives and transform
938
for (unsigned int r = 0; r < num_derivatives; r++)
939
{
940
delete [] combinations[r];
941
}// end loop over 'r'
942
delete [] combinations;
943
for (unsigned int r = 0; r < num_derivatives; r++)
944
{
945
delete [] transform[r];
946
}// end loop over 'r'
947
delete [] transform;
948
break;
949
}
950
case 3:
951
{
952
953
// Array of basisvalues
954
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
955
956
// Declare helper variables
957
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
958
double tmp1 = (1.0 - Y)/2.0;
959
double tmp2 = tmp1*tmp1;
960
961
// Compute basisvalues
962
basisvalues[0] = 1.0;
963
basisvalues[1] = tmp0;
964
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
965
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
966
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
967
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
968
basisvalues[0] *= std::sqrt(0.5);
969
basisvalues[2] *= std::sqrt(1.0);
970
basisvalues[5] *= std::sqrt(1.5);
971
basisvalues[1] *= std::sqrt(3.0);
972
basisvalues[4] *= std::sqrt(4.5);
973
basisvalues[3] *= std::sqrt(7.5);
974
975
// Table(s) of coefficients
976
static const double coefficients0[6] = \
977
{0.47140452, 0.23094011, 0.13333333, 0.0, 0.18856181, -0.16329932};
978
979
// Tables of derivatives of the polynomial base (transpose).
980
static const double dmats0[6][6] = \
981
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
982
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
983
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
984
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
985
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
986
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
987
988
static const double dmats1[6][6] = \
989
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
990
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
991
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
992
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
993
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
994
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
995
996
// Compute reference derivatives.
997
// Declare pointer to array of derivatives on FIAT element.
998
double *derivatives = new double[num_derivatives];
999
for (unsigned int r = 0; r < num_derivatives; r++)
1000
{
1001
derivatives[r] = 0.0;
1002
}// end loop over 'r'
1003
1004
// Declare derivative matrix (of polynomial basis).
1005
double dmats[6][6] = \
1006
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1007
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1008
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1009
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1010
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1011
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1012
1013
// Declare (auxiliary) derivative matrix (of polynomial basis).
1014
double dmats_old[6][6] = \
1015
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1016
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1017
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1018
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1019
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1020
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1021
1022
// Loop possible derivatives.
1023
for (unsigned int r = 0; r < num_derivatives; r++)
1024
{
1025
// Resetting dmats values to compute next derivative.
1026
for (unsigned int t = 0; t < 6; t++)
1027
{
1028
for (unsigned int u = 0; u < 6; u++)
1029
{
1030
dmats[t][u] = 0.0;
1031
if (t == u)
1032
{
1033
dmats[t][u] = 1.0;
1034
}
1035
1036
}// end loop over 'u'
1037
}// end loop over 't'
1038
1039
// Looping derivative order to generate dmats.
1040
for (unsigned int s = 0; s < n; s++)
1041
{
1042
// Updating dmats_old with new values and resetting dmats.
1043
for (unsigned int t = 0; t < 6; t++)
1044
{
1045
for (unsigned int u = 0; u < 6; u++)
1046
{
1047
dmats_old[t][u] = dmats[t][u];
1048
dmats[t][u] = 0.0;
1049
}// end loop over 'u'
1050
}// end loop over 't'
1051
1052
// Update dmats using an inner product.
1053
if (combinations[r][s] == 0)
1054
{
1055
for (unsigned int t = 0; t < 6; t++)
1056
{
1057
for (unsigned int u = 0; u < 6; u++)
1058
{
1059
for (unsigned int tu = 0; tu < 6; tu++)
1060
{
1061
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1062
}// end loop over 'tu'
1063
}// end loop over 'u'
1064
}// end loop over 't'
1065
}
1066
1067
if (combinations[r][s] == 1)
1068
{
1069
for (unsigned int t = 0; t < 6; t++)
1070
{
1071
for (unsigned int u = 0; u < 6; u++)
1072
{
1073
for (unsigned int tu = 0; tu < 6; tu++)
1074
{
1075
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1076
}// end loop over 'tu'
1077
}// end loop over 'u'
1078
}// end loop over 't'
1079
}
1080
1081
}// end loop over 's'
1082
for (unsigned int s = 0; s < 6; s++)
1083
{
1084
for (unsigned int t = 0; t < 6; t++)
1085
{
1086
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1087
}// end loop over 't'
1088
}// end loop over 's'
1089
}// end loop over 'r'
1090
1091
// Transform derivatives back to physical element
1092
for (unsigned int r = 0; r < num_derivatives; r++)
1093
{
1094
for (unsigned int s = 0; s < num_derivatives; s++)
1095
{
1096
values[r] += transform[r][s]*derivatives[s];
1097
}// end loop over 's'
1098
}// end loop over 'r'
1099
1100
// Delete pointer to array of derivatives on FIAT element
1101
delete [] derivatives;
1102
1103
// Delete pointer to array of combinations of derivatives and transform
1104
for (unsigned int r = 0; r < num_derivatives; r++)
1105
{
1106
delete [] combinations[r];
1107
}// end loop over 'r'
1108
delete [] combinations;
1109
for (unsigned int r = 0; r < num_derivatives; r++)
1110
{
1111
delete [] transform[r];
1112
}// end loop over 'r'
1113
delete [] transform;
1114
break;
1115
}
1116
case 4:
1117
{
1118
1119
// Array of basisvalues
1120
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1121
1122
// Declare helper variables
1123
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1124
double tmp1 = (1.0 - Y)/2.0;
1125
double tmp2 = tmp1*tmp1;
1126
1127
// Compute basisvalues
1128
basisvalues[0] = 1.0;
1129
basisvalues[1] = tmp0;
1130
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1131
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1132
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1133
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
1134
basisvalues[0] *= std::sqrt(0.5);
1135
basisvalues[2] *= std::sqrt(1.0);
1136
basisvalues[5] *= std::sqrt(1.5);
1137
basisvalues[1] *= std::sqrt(3.0);
1138
basisvalues[4] *= std::sqrt(4.5);
1139
basisvalues[3] *= std::sqrt(7.5);
1140
1141
// Table(s) of coefficients
1142
static const double coefficients0[6] = \
1143
{0.47140452, -0.23094011, 0.13333333, 0.0, -0.18856181, -0.16329932};
1144
1145
// Tables of derivatives of the polynomial base (transpose).
1146
static const double dmats0[6][6] = \
1147
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1148
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
1149
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1150
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
1151
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
1152
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1153
1154
static const double dmats1[6][6] = \
1155
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1156
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
1157
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
1158
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
1159
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
1160
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
1161
1162
// Compute reference derivatives.
1163
// Declare pointer to array of derivatives on FIAT element.
1164
double *derivatives = new double[num_derivatives];
1165
for (unsigned int r = 0; r < num_derivatives; r++)
1166
{
1167
derivatives[r] = 0.0;
1168
}// end loop over 'r'
1169
1170
// Declare derivative matrix (of polynomial basis).
1171
double dmats[6][6] = \
1172
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1173
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1174
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1175
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1176
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1177
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1178
1179
// Declare (auxiliary) derivative matrix (of polynomial basis).
1180
double dmats_old[6][6] = \
1181
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1182
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1183
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1184
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1185
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1186
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1187
1188
// Loop possible derivatives.
1189
for (unsigned int r = 0; r < num_derivatives; r++)
1190
{
1191
// Resetting dmats values to compute next derivative.
1192
for (unsigned int t = 0; t < 6; t++)
1193
{
1194
for (unsigned int u = 0; u < 6; u++)
1195
{
1196
dmats[t][u] = 0.0;
1197
if (t == u)
1198
{
1199
dmats[t][u] = 1.0;
1200
}
1201
1202
}// end loop over 'u'
1203
}// end loop over 't'
1204
1205
// Looping derivative order to generate dmats.
1206
for (unsigned int s = 0; s < n; s++)
1207
{
1208
// Updating dmats_old with new values and resetting dmats.
1209
for (unsigned int t = 0; t < 6; t++)
1210
{
1211
for (unsigned int u = 0; u < 6; u++)
1212
{
1213
dmats_old[t][u] = dmats[t][u];
1214
dmats[t][u] = 0.0;
1215
}// end loop over 'u'
1216
}// end loop over 't'
1217
1218
// Update dmats using an inner product.
1219
if (combinations[r][s] == 0)
1220
{
1221
for (unsigned int t = 0; t < 6; t++)
1222
{
1223
for (unsigned int u = 0; u < 6; u++)
1224
{
1225
for (unsigned int tu = 0; tu < 6; tu++)
1226
{
1227
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1228
}// end loop over 'tu'
1229
}// end loop over 'u'
1230
}// end loop over 't'
1231
}
1232
1233
if (combinations[r][s] == 1)
1234
{
1235
for (unsigned int t = 0; t < 6; t++)
1236
{
1237
for (unsigned int u = 0; u < 6; u++)
1238
{
1239
for (unsigned int tu = 0; tu < 6; tu++)
1240
{
1241
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1242
}// end loop over 'tu'
1243
}// end loop over 'u'
1244
}// end loop over 't'
1245
}
1246
1247
}// end loop over 's'
1248
for (unsigned int s = 0; s < 6; s++)
1249
{
1250
for (unsigned int t = 0; t < 6; t++)
1251
{
1252
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1253
}// end loop over 't'
1254
}// end loop over 's'
1255
}// end loop over 'r'
1256
1257
// Transform derivatives back to physical element
1258
for (unsigned int r = 0; r < num_derivatives; r++)
1259
{
1260
for (unsigned int s = 0; s < num_derivatives; s++)
1261
{
1262
values[r] += transform[r][s]*derivatives[s];
1263
}// end loop over 's'
1264
}// end loop over 'r'
1265
1266
// Delete pointer to array of derivatives on FIAT element
1267
delete [] derivatives;
1268
1269
// Delete pointer to array of combinations of derivatives and transform
1270
for (unsigned int r = 0; r < num_derivatives; r++)
1271
{
1272
delete [] combinations[r];
1273
}// end loop over 'r'
1274
delete [] combinations;
1275
for (unsigned int r = 0; r < num_derivatives; r++)
1276
{
1277
delete [] transform[r];
1278
}// end loop over 'r'
1279
delete [] transform;
1280
break;
1281
}
1282
case 5:
1283
{
1284
1285
// Array of basisvalues
1286
double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1287
1288
// Declare helper variables
1289
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1290
double tmp1 = (1.0 - Y)/2.0;
1291
double tmp2 = tmp1*tmp1;
1292
1293
// Compute basisvalues
1294
basisvalues[0] = 1.0;
1295
basisvalues[1] = tmp0;
1296
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1297
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1298
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1299
basisvalues[5] = basisvalues[2]*(0.11111111 + Y*1.6666667) - basisvalues[0]*0.55555556;
1300
basisvalues[0] *= std::sqrt(0.5);
1301
basisvalues[2] *= std::sqrt(1.0);
1302
basisvalues[5] *= std::sqrt(1.5);
1303
basisvalues[1] *= std::sqrt(3.0);
1304
basisvalues[4] *= std::sqrt(4.5);
1305
basisvalues[3] *= std::sqrt(7.5);
1306
1307
// Table(s) of coefficients
1308
static const double coefficients0[6] = \
1309
{0.47140452, 0.0, -0.26666667, -0.24343225, 0.0, 0.054433105};
1310
1311
// Tables of derivatives of the polynomial base (transpose).
1312
static const double dmats0[6][6] = \
1313
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1314
{4.8989795, 0.0, 0.0, 0.0, 0.0, 0.0},
1315
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1316
{0.0, 9.486833, 0.0, 0.0, 0.0, 0.0},
1317
{4.0, 0.0, 7.0710678, 0.0, 0.0, 0.0},
1318
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1319
1320
static const double dmats1[6][6] = \
1321
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1322
{2.4494897, 0.0, 0.0, 0.0, 0.0, 0.0},
1323
{4.2426407, 0.0, 0.0, 0.0, 0.0, 0.0},
1324
{2.5819889, 4.7434165, -0.91287093, 0.0, 0.0, 0.0},
1325
{2.0, 6.1237244, 3.5355339, 0.0, 0.0, 0.0},
1326
{-2.3094011, 0.0, 8.1649658, 0.0, 0.0, 0.0}};
1327
1328
// Compute reference derivatives.
1329
// Declare pointer to array of derivatives on FIAT element.
1330
double *derivatives = new double[num_derivatives];
1331
for (unsigned int r = 0; r < num_derivatives; r++)
1332
{
1333
derivatives[r] = 0.0;
1334
}// end loop over 'r'
1335
1336
// Declare derivative matrix (of polynomial basis).
1337
double dmats[6][6] = \
1338
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1339
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1340
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1341
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1342
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1343
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1344
1345
// Declare (auxiliary) derivative matrix (of polynomial basis).
1346
double dmats_old[6][6] = \
1347
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1348
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1349
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1350
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1351
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1352
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1353
1354
// Loop possible derivatives.
1355
for (unsigned int r = 0; r < num_derivatives; r++)
1356
{
1357
// Resetting dmats values to compute next derivative.
1358
for (unsigned int t = 0; t < 6; t++)
1359
{
1360
for (unsigned int u = 0; u < 6; u++)
1361
{
1362
dmats[t][u] = 0.0;
1363
if (t == u)
1364
{
1365
dmats[t][u] = 1.0;
1366
}
1367
1368
}// end loop over 'u'
1369
}// end loop over 't'
1370
1371
// Looping derivative order to generate dmats.
1372
for (unsigned int s = 0; s < n; s++)
1373
{
1374
// Updating dmats_old with new values and resetting dmats.
1375
for (unsigned int t = 0; t < 6; t++)
1376
{
1377
for (unsigned int u = 0; u < 6; u++)
1378
{
1379
dmats_old[t][u] = dmats[t][u];
1380
dmats[t][u] = 0.0;
1381
}// end loop over 'u'
1382
}// end loop over 't'
1383
1384
// Update dmats using an inner product.
1385
if (combinations[r][s] == 0)
1386
{
1387
for (unsigned int t = 0; t < 6; t++)
1388
{
1389
for (unsigned int u = 0; u < 6; u++)
1390
{
1391
for (unsigned int tu = 0; tu < 6; tu++)
1392
{
1393
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1394
}// end loop over 'tu'
1395
}// end loop over 'u'
1396
}// end loop over 't'
1397
}
1398
1399
if (combinations[r][s] == 1)
1400
{
1401
for (unsigned int t = 0; t < 6; t++)
1402
{
1403
for (unsigned int u = 0; u < 6; u++)
1404
{
1405
for (unsigned int tu = 0; tu < 6; tu++)
1406
{
1407
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1408
}// end loop over 'tu'
1409
}// end loop over 'u'
1410
}// end loop over 't'
1411
}
1412
1413
}// end loop over 's'
1414
for (unsigned int s = 0; s < 6; s++)
1415
{
1416
for (unsigned int t = 0; t < 6; t++)
1417
{
1418
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1419
}// end loop over 't'
1420
}// end loop over 's'
1421
}// end loop over 'r'
1422
1423
// Transform derivatives back to physical element
1424
for (unsigned int r = 0; r < num_derivatives; r++)
1425
{
1426
for (unsigned int s = 0; s < num_derivatives; s++)
1427
{
1428
values[r] += transform[r][s]*derivatives[s];
1429
}// end loop over 's'
1430
}// end loop over 'r'
1431
1432
// Delete pointer to array of derivatives on FIAT element
1433
delete [] derivatives;
1434
1435
// Delete pointer to array of combinations of derivatives and transform
1436
for (unsigned int r = 0; r < num_derivatives; r++)
1437
{
1438
delete [] combinations[r];
1439
}// end loop over 'r'
1440
delete [] combinations;
1441
for (unsigned int r = 0; r < num_derivatives; r++)
1442
{
1443
delete [] transform[r];
1444
}// end loop over 'r'
1445
delete [] transform;
1446
break;
1447
}
1448
}
1449
1450
}
1451
1452
/// Evaluate order n derivatives of all basis functions at given point x in cell
1453
virtual void evaluate_basis_derivatives_all(std::size_t n,
1454
double* values,
1455
const double* x,
1456
const double* vertex_coordinates,
1457
int cell_orientation) const
1458
{
1459
// Compute number of derivatives.
1460
unsigned int num_derivatives = 1;
1461
for (unsigned int r = 0; r < n; r++)
1462
{
1463
num_derivatives *= 2;
1464
}// end loop over 'r'
1465
1466
// Helper variable to hold values of a single dof.
1467
double *dof_values = new double[num_derivatives];
1468
for (unsigned int r = 0; r < num_derivatives; r++)
1469
{
1470
dof_values[r] = 0.0;
1471
}// end loop over 'r'
1472
1473
// Loop dofs and call evaluate_basis_derivatives.
1474
for (unsigned int r = 0; r < 6; r++)
1475
{
1476
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1477
for (unsigned int s = 0; s < num_derivatives; s++)
1478
{
1479
values[r*num_derivatives + s] = dof_values[s];
1480
}// end loop over 's'
1481
}// end loop over 'r'
1482
1483
// Delete pointer.
1484
delete [] dof_values;
1485
}
1486
1487
/// Evaluate linear functional for dof i on the function f
1488
virtual double evaluate_dof(std::size_t i,
1489
const ufc::function& f,
1490
const double* vertex_coordinates,
1491
int cell_orientation,
1492
const ufc::cell& c) const
1493
{
1494
// Declare variables for result of evaluation
1495
double vals[1];
1496
1497
// Declare variable for physical coordinates
1498
double y[2];
1499
switch (i)
1500
{
1501
case 0:
1502
{
1503
y[0] = vertex_coordinates[0];
1504
y[1] = vertex_coordinates[1];
1505
f.evaluate(vals, y, c);
1506
return vals[0];
1507
break;
1508
}
1509
case 1:
1510
{
1511
y[0] = vertex_coordinates[2];
1512
y[1] = vertex_coordinates[3];
1513
f.evaluate(vals, y, c);
1514
return vals[0];
1515
break;
1516
}
1517
case 2:
1518
{
1519
y[0] = vertex_coordinates[4];
1520
y[1] = vertex_coordinates[5];
1521
f.evaluate(vals, y, c);
1522
return vals[0];
1523
break;
1524
}
1525
case 3:
1526
{
1527
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
1528
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1529
f.evaluate(vals, y, c);
1530
return vals[0];
1531
break;
1532
}
1533
case 4:
1534
{
1535
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
1536
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1537
f.evaluate(vals, y, c);
1538
return vals[0];
1539
break;
1540
}
1541
case 5:
1542
{
1543
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
1544
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1545
f.evaluate(vals, y, c);
1546
return vals[0];
1547
break;
1548
}
1549
}
1550
1551
return 0.0;
1552
}
1553
1554
/// Evaluate linear functionals for all dofs on the function f
1555
virtual void evaluate_dofs(double* values,
1556
const ufc::function& f,
1557
const double* vertex_coordinates,
1558
int cell_orientation,
1559
const ufc::cell& c) const
1560
{
1561
// Declare variables for result of evaluation
1562
double vals[1];
1563
1564
// Declare variable for physical coordinates
1565
double y[2];
1566
y[0] = vertex_coordinates[0];
1567
y[1] = vertex_coordinates[1];
1568
f.evaluate(vals, y, c);
1569
values[0] = vals[0];
1570
y[0] = vertex_coordinates[2];
1571
y[1] = vertex_coordinates[3];
1572
f.evaluate(vals, y, c);
1573
values[1] = vals[0];
1574
y[0] = vertex_coordinates[4];
1575
y[1] = vertex_coordinates[5];
1576
f.evaluate(vals, y, c);
1577
values[2] = vals[0];
1578
y[0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
1579
y[1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
1580
f.evaluate(vals, y, c);
1581
values[3] = vals[0];
1582
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
1583
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
1584
f.evaluate(vals, y, c);
1585
values[4] = vals[0];
1586
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
1587
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
1588
f.evaluate(vals, y, c);
1589
values[5] = vals[0];
1590
}
1591
1592
/// Interpolate vertex values from dof values
1593
virtual void interpolate_vertex_values(double* vertex_values,
1594
const double* dof_values,
1595
const double* vertex_coordinates,
1596
int cell_orientation,
1597
const ufc::cell& c) const
1598
{
1599
// Evaluate function and change variables
1600
vertex_values[0] = dof_values[0];
1601
vertex_values[1] = dof_values[1];
1602
vertex_values[2] = dof_values[2];
1603
}
1604
1605
/// Map coordinate xhat from reference cell to coordinate x in cell
1606
virtual void map_from_reference_cell(double* x,
1607
const double* xhat,
1608
const ufc::cell& c) const
1609
{
1610
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1611
}
1612
1613
/// Map from coordinate x in cell to coordinate xhat in reference cell
1614
virtual void map_to_reference_cell(double* xhat,
1615
const double* x,
1616
const ufc::cell& c) const
1617
{
1618
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1619
}
1620
1621
/// Return the number of sub elements (for a mixed element)
1622
virtual std::size_t num_sub_elements() const
1623
{
1624
return 0;
1625
}
1626
1627
/// Create a new finite element for sub element i (for a mixed element)
1628
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1629
{
1630
return 0;
1631
}
1632
1633
/// Create a new class instance
1634
virtual ufc::finite_element* create() const
1635
{
1636
return new pointmeasure_finite_element_0();
1637
}
1638
1639
};
1640
1641
/// This class defines the interface for a finite element.
1642
1643
class pointmeasure_finite_element_1: public ufc::finite_element
1644
{
1645
public:
1646
1647
/// Constructor
1648
pointmeasure_finite_element_1() : ufc::finite_element()
1649
{
1650
// Do nothing
1651
}
1652
1653
/// Destructor
1654
virtual ~pointmeasure_finite_element_1()
1655
{
1656
// Do nothing
1657
}
1658
1659
/// Return a string identifying the finite element
1660
virtual const char* signature() const
1661
{
1662
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
1663
}
1664
1665
/// Return the cell shape
1666
virtual ufc::shape cell_shape() const
1667
{
1668
return ufc::triangle;
1669
}
1670
1671
/// Return the topological dimension of the cell shape
1672
virtual std::size_t topological_dimension() const
1673
{
1674
return 2;
1675
}
1676
1677
/// Return the geometric dimension of the cell shape
1678
virtual std::size_t geometric_dimension() const
1679
{
1680
return 2;
1681
}
1682
1683
/// Return the dimension of the finite element function space
1684
virtual std::size_t space_dimension() const
1685
{
1686
return 3;
1687
}
1688
1689
/// Return the rank of the value space
1690
virtual std::size_t value_rank() const
1691
{
1692
return 0;
1693
}
1694
1695
/// Return the dimension of the value space for axis i
1696
virtual std::size_t value_dimension(std::size_t i) const
1697
{
1698
return 1;
1699
}
1700
1701
/// Evaluate basis function i at given point x in cell
1702
virtual void evaluate_basis(std::size_t i,
1703
double* values,
1704
const double* x,
1705
const double* vertex_coordinates,
1706
int cell_orientation) const
1707
{
1708
// Compute Jacobian
1709
double J[4];
1710
compute_jacobian_triangle_2d(J, vertex_coordinates);
1711
1712
// Compute Jacobian inverse and determinant
1713
double K[4];
1714
double detJ;
1715
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1716
1717
1718
// Compute constants
1719
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1720
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1721
1722
// Get coordinates and map to the reference (FIAT) element
1723
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1724
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1725
1726
// Reset values
1727
*values = 0.0;
1728
switch (i)
1729
{
1730
case 0:
1731
{
1732
1733
// Array of basisvalues
1734
double basisvalues[3] = {0.0, 0.0, 0.0};
1735
1736
// Declare helper variables
1737
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1738
1739
// Compute basisvalues
1740
basisvalues[0] = 1.0;
1741
basisvalues[1] = tmp0;
1742
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1743
basisvalues[0] *= std::sqrt(0.5);
1744
basisvalues[2] *= std::sqrt(1.0);
1745
basisvalues[1] *= std::sqrt(3.0);
1746
1747
// Table(s) of coefficients
1748
static const double coefficients0[3] = \
1749
{0.47140452, -0.28867513, -0.16666667};
1750
1751
// Compute value(s)
1752
for (unsigned int r = 0; r < 3; r++)
1753
{
1754
*values += coefficients0[r]*basisvalues[r];
1755
}// end loop over 'r'
1756
break;
1757
}
1758
case 1:
1759
{
1760
1761
// Array of basisvalues
1762
double basisvalues[3] = {0.0, 0.0, 0.0};
1763
1764
// Declare helper variables
1765
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1766
1767
// Compute basisvalues
1768
basisvalues[0] = 1.0;
1769
basisvalues[1] = tmp0;
1770
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1771
basisvalues[0] *= std::sqrt(0.5);
1772
basisvalues[2] *= std::sqrt(1.0);
1773
basisvalues[1] *= std::sqrt(3.0);
1774
1775
// Table(s) of coefficients
1776
static const double coefficients0[3] = \
1777
{0.47140452, 0.28867513, -0.16666667};
1778
1779
// Compute value(s)
1780
for (unsigned int r = 0; r < 3; r++)
1781
{
1782
*values += coefficients0[r]*basisvalues[r];
1783
}// end loop over 'r'
1784
break;
1785
}
1786
case 2:
1787
{
1788
1789
// Array of basisvalues
1790
double basisvalues[3] = {0.0, 0.0, 0.0};
1791
1792
// Declare helper variables
1793
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1794
1795
// Compute basisvalues
1796
basisvalues[0] = 1.0;
1797
basisvalues[1] = tmp0;
1798
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1799
basisvalues[0] *= std::sqrt(0.5);
1800
basisvalues[2] *= std::sqrt(1.0);
1801
basisvalues[1] *= std::sqrt(3.0);
1802
1803
// Table(s) of coefficients
1804
static const double coefficients0[3] = \
1805
{0.47140452, 0.0, 0.33333333};
1806
1807
// Compute value(s)
1808
for (unsigned int r = 0; r < 3; r++)
1809
{
1810
*values += coefficients0[r]*basisvalues[r];
1811
}// end loop over 'r'
1812
break;
1813
}
1814
}
1815
1816
}
1817
1818
/// Evaluate all basis functions at given point x in cell
1819
virtual void evaluate_basis_all(double* values,
1820
const double* x,
1821
const double* vertex_coordinates,
1822
int cell_orientation) const
1823
{
1824
// Helper variable to hold values of a single dof.
1825
double dof_values = 0.0;
1826
1827
// Loop dofs and call evaluate_basis
1828
for (unsigned int r = 0; r < 3; r++)
1829
{
1830
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
1831
values[r] = dof_values;
1832
}// end loop over 'r'
1833
}
1834
1835
/// Evaluate order n derivatives of basis function i at given point x in cell
1836
virtual void evaluate_basis_derivatives(std::size_t i,
1837
std::size_t n,
1838
double* values,
1839
const double* x,
1840
const double* vertex_coordinates,
1841
int cell_orientation) const
1842
{
1843
// Compute Jacobian
1844
double J[4];
1845
compute_jacobian_triangle_2d(J, vertex_coordinates);
1846
1847
// Compute Jacobian inverse and determinant
1848
double K[4];
1849
double detJ;
1850
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1851
1852
1853
// Compute constants
1854
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1855
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1856
1857
// Get coordinates and map to the reference (FIAT) element
1858
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1859
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1860
1861
// Compute number of derivatives.
1862
unsigned int num_derivatives = 1;
1863
for (unsigned int r = 0; r < n; r++)
1864
{
1865
num_derivatives *= 2;
1866
}// end loop over 'r'
1867
1868
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1869
unsigned int **combinations = new unsigned int *[num_derivatives];
1870
for (unsigned int row = 0; row < num_derivatives; row++)
1871
{
1872
combinations[row] = new unsigned int [n];
1873
for (unsigned int col = 0; col < n; col++)
1874
combinations[row][col] = 0;
1875
}
1876
1877
// Generate combinations of derivatives
1878
for (unsigned int row = 1; row < num_derivatives; row++)
1879
{
1880
for (unsigned int num = 0; num < row; num++)
1881
{
1882
for (unsigned int col = n-1; col+1 > 0; col--)
1883
{
1884
if (combinations[row][col] + 1 > 1)
1885
combinations[row][col] = 0;
1886
else
1887
{
1888
combinations[row][col] += 1;
1889
break;
1890
}
1891
}
1892
}
1893
}
1894
1895
// Compute inverse of Jacobian
1896
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
1897
1898
// Declare transformation matrix
1899
// Declare pointer to two dimensional array and initialise
1900
double **transform = new double *[num_derivatives];
1901
1902
for (unsigned int j = 0; j < num_derivatives; j++)
1903
{
1904
transform[j] = new double [num_derivatives];
1905
for (unsigned int k = 0; k < num_derivatives; k++)
1906
transform[j][k] = 1;
1907
}
1908
1909
// Construct transformation matrix
1910
for (unsigned int row = 0; row < num_derivatives; row++)
1911
{
1912
for (unsigned int col = 0; col < num_derivatives; col++)
1913
{
1914
for (unsigned int k = 0; k < n; k++)
1915
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1916
}
1917
}
1918
1919
// Reset values. Assuming that values is always an array.
1920
for (unsigned int r = 0; r < num_derivatives; r++)
1921
{
1922
values[r] = 0.0;
1923
}// end loop over 'r'
1924
1925
switch (i)
1926
{
1927
case 0:
1928
{
1929
1930
// Array of basisvalues
1931
double basisvalues[3] = {0.0, 0.0, 0.0};
1932
1933
// Declare helper variables
1934
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1935
1936
// Compute basisvalues
1937
basisvalues[0] = 1.0;
1938
basisvalues[1] = tmp0;
1939
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1940
basisvalues[0] *= std::sqrt(0.5);
1941
basisvalues[2] *= std::sqrt(1.0);
1942
basisvalues[1] *= std::sqrt(3.0);
1943
1944
// Table(s) of coefficients
1945
static const double coefficients0[3] = \
1946
{0.47140452, -0.28867513, -0.16666667};
1947
1948
// Tables of derivatives of the polynomial base (transpose).
1949
static const double dmats0[3][3] = \
1950
{{0.0, 0.0, 0.0},
1951
{4.8989795, 0.0, 0.0},
1952
{0.0, 0.0, 0.0}};
1953
1954
static const double dmats1[3][3] = \
1955
{{0.0, 0.0, 0.0},
1956
{2.4494897, 0.0, 0.0},
1957
{4.2426407, 0.0, 0.0}};
1958
1959
// Compute reference derivatives.
1960
// Declare pointer to array of derivatives on FIAT element.
1961
double *derivatives = new double[num_derivatives];
1962
for (unsigned int r = 0; r < num_derivatives; r++)
1963
{
1964
derivatives[r] = 0.0;
1965
}// end loop over 'r'
1966
1967
// Declare derivative matrix (of polynomial basis).
1968
double dmats[3][3] = \
1969
{{1.0, 0.0, 0.0},
1970
{0.0, 1.0, 0.0},
1971
{0.0, 0.0, 1.0}};
1972
1973
// Declare (auxiliary) derivative matrix (of polynomial basis).
1974
double dmats_old[3][3] = \
1975
{{1.0, 0.0, 0.0},
1976
{0.0, 1.0, 0.0},
1977
{0.0, 0.0, 1.0}};
1978
1979
// Loop possible derivatives.
1980
for (unsigned int r = 0; r < num_derivatives; r++)
1981
{
1982
// Resetting dmats values to compute next derivative.
1983
for (unsigned int t = 0; t < 3; t++)
1984
{
1985
for (unsigned int u = 0; u < 3; u++)
1986
{
1987
dmats[t][u] = 0.0;
1988
if (t == u)
1989
{
1990
dmats[t][u] = 1.0;
1991
}
1992
1993
}// end loop over 'u'
1994
}// end loop over 't'
1995
1996
// Looping derivative order to generate dmats.
1997
for (unsigned int s = 0; s < n; s++)
1998
{
1999
// Updating dmats_old with new values and resetting dmats.
2000
for (unsigned int t = 0; t < 3; t++)
2001
{
2002
for (unsigned int u = 0; u < 3; u++)
2003
{
2004
dmats_old[t][u] = dmats[t][u];
2005
dmats[t][u] = 0.0;
2006
}// end loop over 'u'
2007
}// end loop over 't'
2008
2009
// Update dmats using an inner product.
2010
if (combinations[r][s] == 0)
2011
{
2012
for (unsigned int t = 0; t < 3; t++)
2013
{
2014
for (unsigned int u = 0; u < 3; u++)
2015
{
2016
for (unsigned int tu = 0; tu < 3; tu++)
2017
{
2018
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2019
}// end loop over 'tu'
2020
}// end loop over 'u'
2021
}// end loop over 't'
2022
}
2023
2024
if (combinations[r][s] == 1)
2025
{
2026
for (unsigned int t = 0; t < 3; t++)
2027
{
2028
for (unsigned int u = 0; u < 3; u++)
2029
{
2030
for (unsigned int tu = 0; tu < 3; tu++)
2031
{
2032
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2033
}// end loop over 'tu'
2034
}// end loop over 'u'
2035
}// end loop over 't'
2036
}
2037
2038
}// end loop over 's'
2039
for (unsigned int s = 0; s < 3; s++)
2040
{
2041
for (unsigned int t = 0; t < 3; t++)
2042
{
2043
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2044
}// end loop over 't'
2045
}// end loop over 's'
2046
}// end loop over 'r'
2047
2048
// Transform derivatives back to physical element
2049
for (unsigned int r = 0; r < num_derivatives; r++)
2050
{
2051
for (unsigned int s = 0; s < num_derivatives; s++)
2052
{
2053
values[r] += transform[r][s]*derivatives[s];
2054
}// end loop over 's'
2055
}// end loop over 'r'
2056
2057
// Delete pointer to array of derivatives on FIAT element
2058
delete [] derivatives;
2059
2060
// Delete pointer to array of combinations of derivatives and transform
2061
for (unsigned int r = 0; r < num_derivatives; r++)
2062
{
2063
delete [] combinations[r];
2064
}// end loop over 'r'
2065
delete [] combinations;
2066
for (unsigned int r = 0; r < num_derivatives; r++)
2067
{
2068
delete [] transform[r];
2069
}// end loop over 'r'
2070
delete [] transform;
2071
break;
2072
}
2073
case 1:
2074
{
2075
2076
// Array of basisvalues
2077
double basisvalues[3] = {0.0, 0.0, 0.0};
2078
2079
// Declare helper variables
2080
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2081
2082
// Compute basisvalues
2083
basisvalues[0] = 1.0;
2084
basisvalues[1] = tmp0;
2085
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2086
basisvalues[0] *= std::sqrt(0.5);
2087
basisvalues[2] *= std::sqrt(1.0);
2088
basisvalues[1] *= std::sqrt(3.0);
2089
2090
// Table(s) of coefficients
2091
static const double coefficients0[3] = \
2092
{0.47140452, 0.28867513, -0.16666667};
2093
2094
// Tables of derivatives of the polynomial base (transpose).
2095
static const double dmats0[3][3] = \
2096
{{0.0, 0.0, 0.0},
2097
{4.8989795, 0.0, 0.0},
2098
{0.0, 0.0, 0.0}};
2099
2100
static const double dmats1[3][3] = \
2101
{{0.0, 0.0, 0.0},
2102
{2.4494897, 0.0, 0.0},
2103
{4.2426407, 0.0, 0.0}};
2104
2105
// Compute reference derivatives.
2106
// Declare pointer to array of derivatives on FIAT element.
2107
double *derivatives = new double[num_derivatives];
2108
for (unsigned int r = 0; r < num_derivatives; r++)
2109
{
2110
derivatives[r] = 0.0;
2111
}// end loop over 'r'
2112
2113
// Declare derivative matrix (of polynomial basis).
2114
double dmats[3][3] = \
2115
{{1.0, 0.0, 0.0},
2116
{0.0, 1.0, 0.0},
2117
{0.0, 0.0, 1.0}};
2118
2119
// Declare (auxiliary) derivative matrix (of polynomial basis).
2120
double dmats_old[3][3] = \
2121
{{1.0, 0.0, 0.0},
2122
{0.0, 1.0, 0.0},
2123
{0.0, 0.0, 1.0}};
2124
2125
// Loop possible derivatives.
2126
for (unsigned int r = 0; r < num_derivatives; r++)
2127
{
2128
// Resetting dmats values to compute next derivative.
2129
for (unsigned int t = 0; t < 3; t++)
2130
{
2131
for (unsigned int u = 0; u < 3; u++)
2132
{
2133
dmats[t][u] = 0.0;
2134
if (t == u)
2135
{
2136
dmats[t][u] = 1.0;
2137
}
2138
2139
}// end loop over 'u'
2140
}// end loop over 't'
2141
2142
// Looping derivative order to generate dmats.
2143
for (unsigned int s = 0; s < n; s++)
2144
{
2145
// Updating dmats_old with new values and resetting dmats.
2146
for (unsigned int t = 0; t < 3; t++)
2147
{
2148
for (unsigned int u = 0; u < 3; u++)
2149
{
2150
dmats_old[t][u] = dmats[t][u];
2151
dmats[t][u] = 0.0;
2152
}// end loop over 'u'
2153
}// end loop over 't'
2154
2155
// Update dmats using an inner product.
2156
if (combinations[r][s] == 0)
2157
{
2158
for (unsigned int t = 0; t < 3; t++)
2159
{
2160
for (unsigned int u = 0; u < 3; u++)
2161
{
2162
for (unsigned int tu = 0; tu < 3; tu++)
2163
{
2164
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2165
}// end loop over 'tu'
2166
}// end loop over 'u'
2167
}// end loop over 't'
2168
}
2169
2170
if (combinations[r][s] == 1)
2171
{
2172
for (unsigned int t = 0; t < 3; t++)
2173
{
2174
for (unsigned int u = 0; u < 3; u++)
2175
{
2176
for (unsigned int tu = 0; tu < 3; tu++)
2177
{
2178
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2179
}// end loop over 'tu'
2180
}// end loop over 'u'
2181
}// end loop over 't'
2182
}
2183
2184
}// end loop over 's'
2185
for (unsigned int s = 0; s < 3; s++)
2186
{
2187
for (unsigned int t = 0; t < 3; t++)
2188
{
2189
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2190
}// end loop over 't'
2191
}// end loop over 's'
2192
}// end loop over 'r'
2193
2194
// Transform derivatives back to physical element
2195
for (unsigned int r = 0; r < num_derivatives; r++)
2196
{
2197
for (unsigned int s = 0; s < num_derivatives; s++)
2198
{
2199
values[r] += transform[r][s]*derivatives[s];
2200
}// end loop over 's'
2201
}// end loop over 'r'
2202
2203
// Delete pointer to array of derivatives on FIAT element
2204
delete [] derivatives;
2205
2206
// Delete pointer to array of combinations of derivatives and transform
2207
for (unsigned int r = 0; r < num_derivatives; r++)
2208
{
2209
delete [] combinations[r];
2210
}// end loop over 'r'
2211
delete [] combinations;
2212
for (unsigned int r = 0; r < num_derivatives; r++)
2213
{
2214
delete [] transform[r];
2215
}// end loop over 'r'
2216
delete [] transform;
2217
break;
2218
}
2219
case 2:
2220
{
2221
2222
// Array of basisvalues
2223
double basisvalues[3] = {0.0, 0.0, 0.0};
2224
2225
// Declare helper variables
2226
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2227
2228
// Compute basisvalues
2229
basisvalues[0] = 1.0;
2230
basisvalues[1] = tmp0;
2231
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2232
basisvalues[0] *= std::sqrt(0.5);
2233
basisvalues[2] *= std::sqrt(1.0);
2234
basisvalues[1] *= std::sqrt(3.0);
2235
2236
// Table(s) of coefficients
2237
static const double coefficients0[3] = \
2238
{0.47140452, 0.0, 0.33333333};
2239
2240
// Tables of derivatives of the polynomial base (transpose).
2241
static const double dmats0[3][3] = \
2242
{{0.0, 0.0, 0.0},
2243
{4.8989795, 0.0, 0.0},
2244
{0.0, 0.0, 0.0}};
2245
2246
static const double dmats1[3][3] = \
2247
{{0.0, 0.0, 0.0},
2248
{2.4494897, 0.0, 0.0},
2249
{4.2426407, 0.0, 0.0}};
2250
2251
// Compute reference derivatives.
2252
// Declare pointer to array of derivatives on FIAT element.
2253
double *derivatives = new double[num_derivatives];
2254
for (unsigned int r = 0; r < num_derivatives; r++)
2255
{
2256
derivatives[r] = 0.0;
2257
}// end loop over 'r'
2258
2259
// Declare derivative matrix (of polynomial basis).
2260
double dmats[3][3] = \
2261
{{1.0, 0.0, 0.0},
2262
{0.0, 1.0, 0.0},
2263
{0.0, 0.0, 1.0}};
2264
2265
// Declare (auxiliary) derivative matrix (of polynomial basis).
2266
double dmats_old[3][3] = \
2267
{{1.0, 0.0, 0.0},
2268
{0.0, 1.0, 0.0},
2269
{0.0, 0.0, 1.0}};
2270
2271
// Loop possible derivatives.
2272
for (unsigned int r = 0; r < num_derivatives; r++)
2273
{
2274
// Resetting dmats values to compute next derivative.
2275
for (unsigned int t = 0; t < 3; t++)
2276
{
2277
for (unsigned int u = 0; u < 3; u++)
2278
{
2279
dmats[t][u] = 0.0;
2280
if (t == u)
2281
{
2282
dmats[t][u] = 1.0;
2283
}
2284
2285
}// end loop over 'u'
2286
}// end loop over 't'
2287
2288
// Looping derivative order to generate dmats.
2289
for (unsigned int s = 0; s < n; s++)
2290
{
2291
// Updating dmats_old with new values and resetting dmats.
2292
for (unsigned int t = 0; t < 3; t++)
2293
{
2294
for (unsigned int u = 0; u < 3; u++)
2295
{
2296
dmats_old[t][u] = dmats[t][u];
2297
dmats[t][u] = 0.0;
2298
}// end loop over 'u'
2299
}// end loop over 't'
2300
2301
// Update dmats using an inner product.
2302
if (combinations[r][s] == 0)
2303
{
2304
for (unsigned int t = 0; t < 3; t++)
2305
{
2306
for (unsigned int u = 0; u < 3; u++)
2307
{
2308
for (unsigned int tu = 0; tu < 3; tu++)
2309
{
2310
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2311
}// end loop over 'tu'
2312
}// end loop over 'u'
2313
}// end loop over 't'
2314
}
2315
2316
if (combinations[r][s] == 1)
2317
{
2318
for (unsigned int t = 0; t < 3; t++)
2319
{
2320
for (unsigned int u = 0; u < 3; u++)
2321
{
2322
for (unsigned int tu = 0; tu < 3; tu++)
2323
{
2324
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2325
}// end loop over 'tu'
2326
}// end loop over 'u'
2327
}// end loop over 't'
2328
}
2329
2330
}// end loop over 's'
2331
for (unsigned int s = 0; s < 3; s++)
2332
{
2333
for (unsigned int t = 0; t < 3; t++)
2334
{
2335
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2336
}// end loop over 't'
2337
}// end loop over 's'
2338
}// end loop over 'r'
2339
2340
// Transform derivatives back to physical element
2341
for (unsigned int r = 0; r < num_derivatives; r++)
2342
{
2343
for (unsigned int s = 0; s < num_derivatives; s++)
2344
{
2345
values[r] += transform[r][s]*derivatives[s];
2346
}// end loop over 's'
2347
}// end loop over 'r'
2348
2349
// Delete pointer to array of derivatives on FIAT element
2350
delete [] derivatives;
2351
2352
// Delete pointer to array of combinations of derivatives and transform
2353
for (unsigned int r = 0; r < num_derivatives; r++)
2354
{
2355
delete [] combinations[r];
2356
}// end loop over 'r'
2357
delete [] combinations;
2358
for (unsigned int r = 0; r < num_derivatives; r++)
2359
{
2360
delete [] transform[r];
2361
}// end loop over 'r'
2362
delete [] transform;
2363
break;
2364
}
2365
}
2366
2367
}
2368
2369
/// Evaluate order n derivatives of all basis functions at given point x in cell
2370
virtual void evaluate_basis_derivatives_all(std::size_t n,
2371
double* values,
2372
const double* x,
2373
const double* vertex_coordinates,
2374
int cell_orientation) const
2375
{
2376
// Compute number of derivatives.
2377
unsigned int num_derivatives = 1;
2378
for (unsigned int r = 0; r < n; r++)
2379
{
2380
num_derivatives *= 2;
2381
}// end loop over 'r'
2382
2383
// Helper variable to hold values of a single dof.
2384
double *dof_values = new double[num_derivatives];
2385
for (unsigned int r = 0; r < num_derivatives; r++)
2386
{
2387
dof_values[r] = 0.0;
2388
}// end loop over 'r'
2389
2390
// Loop dofs and call evaluate_basis_derivatives.
2391
for (unsigned int r = 0; r < 3; r++)
2392
{
2393
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2394
for (unsigned int s = 0; s < num_derivatives; s++)
2395
{
2396
values[r*num_derivatives + s] = dof_values[s];
2397
}// end loop over 's'
2398
}// end loop over 'r'
2399
2400
// Delete pointer.
2401
delete [] dof_values;
2402
}
2403
2404
/// Evaluate linear functional for dof i on the function f
2405
virtual double evaluate_dof(std::size_t i,
2406
const ufc::function& f,
2407
const double* vertex_coordinates,
2408
int cell_orientation,
2409
const ufc::cell& c) const
2410
{
2411
// Declare variables for result of evaluation
2412
double vals[1];
2413
2414
// Declare variable for physical coordinates
2415
double y[2];
2416
switch (i)
2417
{
2418
case 0:
2419
{
2420
y[0] = vertex_coordinates[0];
2421
y[1] = vertex_coordinates[1];
2422
f.evaluate(vals, y, c);
2423
return vals[0];
2424
break;
2425
}
2426
case 1:
2427
{
2428
y[0] = vertex_coordinates[2];
2429
y[1] = vertex_coordinates[3];
2430
f.evaluate(vals, y, c);
2431
return vals[0];
2432
break;
2433
}
2434
case 2:
2435
{
2436
y[0] = vertex_coordinates[4];
2437
y[1] = vertex_coordinates[5];
2438
f.evaluate(vals, y, c);
2439
return vals[0];
2440
break;
2441
}
2442
}
2443
2444
return 0.0;
2445
}
2446
2447
/// Evaluate linear functionals for all dofs on the function f
2448
virtual void evaluate_dofs(double* values,
2449
const ufc::function& f,
2450
const double* vertex_coordinates,
2451
int cell_orientation,
2452
const ufc::cell& c) const
2453
{
2454
// Declare variables for result of evaluation
2455
double vals[1];
2456
2457
// Declare variable for physical coordinates
2458
double y[2];
2459
y[0] = vertex_coordinates[0];
2460
y[1] = vertex_coordinates[1];
2461
f.evaluate(vals, y, c);
2462
values[0] = vals[0];
2463
y[0] = vertex_coordinates[2];
2464
y[1] = vertex_coordinates[3];
2465
f.evaluate(vals, y, c);
2466
values[1] = vals[0];
2467
y[0] = vertex_coordinates[4];
2468
y[1] = vertex_coordinates[5];
2469
f.evaluate(vals, y, c);
2470
values[2] = vals[0];
2471
}
2472
2473
/// Interpolate vertex values from dof values
2474
virtual void interpolate_vertex_values(double* vertex_values,
2475
const double* dof_values,
2476
const double* vertex_coordinates,
2477
int cell_orientation,
2478
const ufc::cell& c) const
2479
{
2480
// Evaluate function and change variables
2481
vertex_values[0] = dof_values[0];
2482
vertex_values[1] = dof_values[1];
2483
vertex_values[2] = dof_values[2];
2484
}
2485
2486
/// Map coordinate xhat from reference cell to coordinate x in cell
2487
virtual void map_from_reference_cell(double* x,
2488
const double* xhat,
2489
const ufc::cell& c) const
2490
{
2491
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
2492
}
2493
2494
/// Map from coordinate x in cell to coordinate xhat in reference cell
2495
virtual void map_to_reference_cell(double* xhat,
2496
const double* x,
2497
const ufc::cell& c) const
2498
{
2499
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
2500
}
2501
2502
/// Return the number of sub elements (for a mixed element)
2503
virtual std::size_t num_sub_elements() const
2504
{
2505
return 0;
2506
}
2507
2508
/// Create a new finite element for sub element i (for a mixed element)
2509
virtual ufc::finite_element* create_sub_element(std::size_t i) const
2510
{
2511
return 0;
2512
}
2513
2514
/// Create a new class instance
2515
virtual ufc::finite_element* create() const
2516
{
2517
return new pointmeasure_finite_element_1();
2518
}
2519
2520
};
2521
2522
/// This class defines the interface for a local-to-global mapping of
2523
/// degrees of freedom (dofs).
2524
2525
class pointmeasure_dofmap_0: public ufc::dofmap
2526
{
2527
public:
2528
2529
/// Constructor
2530
pointmeasure_dofmap_0() : ufc::dofmap()
2531
{
2532
// Do nothing
2533
}
2534
2535
/// Destructor
2536
virtual ~pointmeasure_dofmap_0()
2537
{
2538
// Do nothing
2539
}
2540
2541
/// Return a string identifying the dofmap
2542
virtual const char* signature() const
2543
{
2544
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
2545
}
2546
2547
/// Return true iff mesh entities of topological dimension d are needed
2548
virtual bool needs_mesh_entities(std::size_t d) const
2549
{
2550
switch (d)
2551
{
2552
case 0:
2553
{
2554
return true;
2555
break;
2556
}
2557
case 1:
2558
{
2559
return true;
2560
break;
2561
}
2562
case 2:
2563
{
2564
return false;
2565
break;
2566
}
2567
}
2568
2569
return false;
2570
}
2571
2572
/// Return the topological dimension of the associated cell shape
2573
virtual std::size_t topological_dimension() const
2574
{
2575
return 2;
2576
}
2577
2578
/// Return the geometric dimension of the associated cell shape
2579
virtual std::size_t geometric_dimension() const
2580
{
2581
return 2;
2582
}
2583
2584
/// Return the dimension of the global finite element function space
2585
virtual std::size_t global_dimension(const std::vector<std::size_t>&
2586
num_global_entities) const
2587
{
2588
return num_global_entities[0] + num_global_entities[1];
2589
}
2590
2591
/// Return the dimension of the local finite element function space for a cell
2592
virtual std::size_t local_dimension(const ufc::cell& c) const
2593
{
2594
return 6;
2595
}
2596
2597
/// Return the maximum dimension of the local finite element function space
2598
virtual std::size_t max_local_dimension() const
2599
{
2600
return 6;
2601
}
2602
2603
/// Return the number of dofs on each cell facet
2604
virtual std::size_t num_facet_dofs() const
2605
{
2606
return 3;
2607
}
2608
2609
/// Return the number of dofs associated with each cell entity of dimension d
2610
virtual std::size_t num_entity_dofs(std::size_t d) const
2611
{
2612
switch (d)
2613
{
2614
case 0:
2615
{
2616
return 1;
2617
break;
2618
}
2619
case 1:
2620
{
2621
return 1;
2622
break;
2623
}
2624
case 2:
2625
{
2626
return 0;
2627
break;
2628
}
2629
}
2630
2631
return 0;
2632
}
2633
2634
/// Tabulate the local-to-global mapping of dofs on a cell
2635
virtual void tabulate_dofs(std::size_t* dofs,
2636
const std::vector<std::size_t>& num_global_entities,
2637
const ufc::cell& c) const
2638
{
2639
unsigned int offset = 0;
2640
dofs[0] = offset + c.entity_indices[0][0];
2641
dofs[1] = offset + c.entity_indices[0][1];
2642
dofs[2] = offset + c.entity_indices[0][2];
2643
offset += num_global_entities[0];
2644
dofs[3] = offset + c.entity_indices[1][0];
2645
dofs[4] = offset + c.entity_indices[1][1];
2646
dofs[5] = offset + c.entity_indices[1][2];
2647
offset += num_global_entities[1];
2648
}
2649
2650
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2651
virtual void tabulate_facet_dofs(std::size_t* dofs,
2652
std::size_t facet) const
2653
{
2654
switch (facet)
2655
{
2656
case 0:
2657
{
2658
dofs[0] = 1;
2659
dofs[1] = 2;
2660
dofs[2] = 3;
2661
break;
2662
}
2663
case 1:
2664
{
2665
dofs[0] = 0;
2666
dofs[1] = 2;
2667
dofs[2] = 4;
2668
break;
2669
}
2670
case 2:
2671
{
2672
dofs[0] = 0;
2673
dofs[1] = 1;
2674
dofs[2] = 5;
2675
break;
2676
}
2677
}
2678
2679
}
2680
2681
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2682
virtual void tabulate_entity_dofs(std::size_t* dofs,
2683
std::size_t d, std::size_t i) const
2684
{
2685
if (d > 2)
2686
{
2687
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2688
}
2689
2690
switch (d)
2691
{
2692
case 0:
2693
{
2694
if (i > 2)
2695
{
2696
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2697
}
2698
2699
switch (i)
2700
{
2701
case 0:
2702
{
2703
dofs[0] = 0;
2704
break;
2705
}
2706
case 1:
2707
{
2708
dofs[0] = 1;
2709
break;
2710
}
2711
case 2:
2712
{
2713
dofs[0] = 2;
2714
break;
2715
}
2716
}
2717
2718
break;
2719
}
2720
case 1:
2721
{
2722
if (i > 2)
2723
{
2724
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2725
}
2726
2727
switch (i)
2728
{
2729
case 0:
2730
{
2731
dofs[0] = 3;
2732
break;
2733
}
2734
case 1:
2735
{
2736
dofs[0] = 4;
2737
break;
2738
}
2739
case 2:
2740
{
2741
dofs[0] = 5;
2742
break;
2743
}
2744
}
2745
2746
break;
2747
}
2748
case 2:
2749
{
2750
2751
break;
2752
}
2753
}
2754
2755
}
2756
2757
/// Tabulate the coordinates of all dofs on a cell
2758
virtual void tabulate_coordinates(double** dof_coordinates,
2759
const double* vertex_coordinates) const
2760
{
2761
dof_coordinates[0][0] = vertex_coordinates[0];
2762
dof_coordinates[0][1] = vertex_coordinates[1];
2763
dof_coordinates[1][0] = vertex_coordinates[2];
2764
dof_coordinates[1][1] = vertex_coordinates[3];
2765
dof_coordinates[2][0] = vertex_coordinates[4];
2766
dof_coordinates[2][1] = vertex_coordinates[5];
2767
dof_coordinates[3][0] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[4];
2768
dof_coordinates[3][1] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[5];
2769
dof_coordinates[4][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[4];
2770
dof_coordinates[4][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[5];
2771
dof_coordinates[5][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[2];
2772
dof_coordinates[5][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[3];
2773
}
2774
2775
/// Return the number of sub dofmaps (for a mixed element)
2776
virtual std::size_t num_sub_dofmaps() const
2777
{
2778
return 0;
2779
}
2780
2781
/// Create a new dofmap for sub dofmap i (for a mixed element)
2782
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
2783
{
2784
return 0;
2785
}
2786
2787
/// Create a new class instance
2788
virtual ufc::dofmap* create() const
2789
{
2790
return new pointmeasure_dofmap_0();
2791
}
2792
2793
};
2794
2795
/// This class defines the interface for a local-to-global mapping of
2796
/// degrees of freedom (dofs).
2797
2798
class pointmeasure_dofmap_1: public ufc::dofmap
2799
{
2800
public:
2801
2802
/// Constructor
2803
pointmeasure_dofmap_1() : ufc::dofmap()
2804
{
2805
// Do nothing
2806
}
2807
2808
/// Destructor
2809
virtual ~pointmeasure_dofmap_1()
2810
{
2811
// Do nothing
2812
}
2813
2814
/// Return a string identifying the dofmap
2815
virtual const char* signature() const
2816
{
2817
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
2818
}
2819
2820
/// Return true iff mesh entities of topological dimension d are needed
2821
virtual bool needs_mesh_entities(std::size_t d) const
2822
{
2823
switch (d)
2824
{
2825
case 0:
2826
{
2827
return true;
2828
break;
2829
}
2830
case 1:
2831
{
2832
return false;
2833
break;
2834
}
2835
case 2:
2836
{
2837
return false;
2838
break;
2839
}
2840
}
2841
2842
return false;
2843
}
2844
2845
/// Return the topological dimension of the associated cell shape
2846
virtual std::size_t topological_dimension() const
2847
{
2848
return 2;
2849
}
2850
2851
/// Return the geometric dimension of the associated cell shape
2852
virtual std::size_t geometric_dimension() const
2853
{
2854
return 2;
2855
}
2856
2857
/// Return the dimension of the global finite element function space
2858
virtual std::size_t global_dimension(const std::vector<std::size_t>&
2859
num_global_entities) const
2860
{
2861
return num_global_entities[0];
2862
}
2863
2864
/// Return the dimension of the local finite element function space for a cell
2865
virtual std::size_t local_dimension(const ufc::cell& c) const
2866
{
2867
return 3;
2868
}
2869
2870
/// Return the maximum dimension of the local finite element function space
2871
virtual std::size_t max_local_dimension() const
2872
{
2873
return 3;
2874
}
2875
2876
/// Return the number of dofs on each cell facet
2877
virtual std::size_t num_facet_dofs() const
2878
{
2879
return 2;
2880
}
2881
2882
/// Return the number of dofs associated with each cell entity of dimension d
2883
virtual std::size_t num_entity_dofs(std::size_t d) const
2884
{
2885
switch (d)
2886
{
2887
case 0:
2888
{
2889
return 1;
2890
break;
2891
}
2892
case 1:
2893
{
2894
return 0;
2895
break;
2896
}
2897
case 2:
2898
{
2899
return 0;
2900
break;
2901
}
2902
}
2903
2904
return 0;
2905
}
2906
2907
/// Tabulate the local-to-global mapping of dofs on a cell
2908
virtual void tabulate_dofs(std::size_t* dofs,
2909
const std::vector<std::size_t>& num_global_entities,
2910
const ufc::cell& c) const
2911
{
2912
dofs[0] = c.entity_indices[0][0];
2913
dofs[1] = c.entity_indices[0][1];
2914
dofs[2] = c.entity_indices[0][2];
2915
}
2916
2917
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2918
virtual void tabulate_facet_dofs(std::size_t* dofs,
2919
std::size_t facet) const
2920
{
2921
switch (facet)
2922
{
2923
case 0:
2924
{
2925
dofs[0] = 1;
2926
dofs[1] = 2;
2927
break;
2928
}
2929
case 1:
2930
{
2931
dofs[0] = 0;
2932
dofs[1] = 2;
2933
break;
2934
}
2935
case 2:
2936
{
2937
dofs[0] = 0;
2938
dofs[1] = 1;
2939
break;
2940
}
2941
}
2942
2943
}
2944
2945
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2946
virtual void tabulate_entity_dofs(std::size_t* dofs,
2947
std::size_t d, std::size_t i) const
2948
{
2949
if (d > 2)
2950
{
2951
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2952
}
2953
2954
switch (d)
2955
{
2956
case 0:
2957
{
2958
if (i > 2)
2959
{
2960
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2961
}
2962
2963
switch (i)
2964
{
2965
case 0:
2966
{
2967
dofs[0] = 0;
2968
break;
2969
}
2970
case 1:
2971
{
2972
dofs[0] = 1;
2973
break;
2974
}
2975
case 2:
2976
{
2977
dofs[0] = 2;
2978
break;
2979
}
2980
}
2981
2982
break;
2983
}
2984
case 1:
2985
{
2986
2987
break;
2988
}
2989
case 2:
2990
{
2991
2992
break;
2993
}
2994
}
2995
2996
}
2997
2998
/// Tabulate the coordinates of all dofs on a cell
2999
virtual void tabulate_coordinates(double** dof_coordinates,
3000
const double* vertex_coordinates) const
3001
{
3002
dof_coordinates[0][0] = vertex_coordinates[0];
3003
dof_coordinates[0][1] = vertex_coordinates[1];
3004
dof_coordinates[1][0] = vertex_coordinates[2];
3005
dof_coordinates[1][1] = vertex_coordinates[3];
3006
dof_coordinates[2][0] = vertex_coordinates[4];
3007
dof_coordinates[2][1] = vertex_coordinates[5];
3008
}
3009
3010
/// Return the number of sub dofmaps (for a mixed element)
3011
virtual std::size_t num_sub_dofmaps() const
3012
{
3013
return 0;
3014
}
3015
3016
/// Create a new dofmap for sub dofmap i (for a mixed element)
3017
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3018
{
3019
return 0;
3020
}
3021
3022
/// Create a new class instance
3023
virtual ufc::dofmap* create() const
3024
{
3025
return new pointmeasure_dofmap_1();
3026
}
3027
3028
};
3029
3030
/// This class defines the interface for the tabulation of the cell
3031
/// tensor corresponding to the local contribution to a form from
3032
/// the integral over a cell.
3033
3034
class pointmeasure_cell_integral_0_otherwise: public ufc::cell_integral
3035
{
3036
public:
3037
3038
/// Constructor
3039
pointmeasure_cell_integral_0_otherwise() : ufc::cell_integral()
3040
{
3041
// Do nothing
3042
}
3043
3044
/// Destructor
3045
virtual ~pointmeasure_cell_integral_0_otherwise()
3046
{
3047
// Do nothing
3048
}
3049
3050
/// Tabulate the tensor for the contribution from a local cell
3051
virtual void tabulate_tensor(double* A,
3052
const double * const * w,
3053
const double* vertex_coordinates,
3054
int cell_orientation) const
3055
{
3056
// Compute Jacobian
3057
double J[4];
3058
compute_jacobian_triangle_2d(J, vertex_coordinates);
3059
3060
// Compute Jacobian inverse and determinant
3061
double K[4];
3062
double detJ;
3063
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3064
3065
// Set scale factor
3066
const double det = std::abs(detJ);
3067
3068
// Cell volume
3069
3070
// Compute circumradius of triangle in 2D
3071
3072
3073
// Facet area
3074
3075
// Array of quadrature weights.
3076
static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
3077
// Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
3078
3079
// Value of basis functions at quadrature points.
3080
static const double FE0[3][3] = \
3081
{{0.66666667, 0.16666667, 0.16666667},
3082
{0.16666667, 0.16666667, 0.66666667},
3083
{0.16666667, 0.66666667, 0.16666667}};
3084
3085
// Reset values in the element tensor.
3086
for (unsigned int r = 0; r < 9; r++)
3087
{
3088
A[r] = 0.0;
3089
}// end loop over 'r'
3090
3091
// Compute element tensor using UFL quadrature representation
3092
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3093
3094
// Loop quadrature points for integral.
3095
// Number of operations to compute element tensor for following IP loop = 108
3096
for (unsigned int ip = 0; ip < 3; ip++)
3097
{
3098
3099
// Number of operations for primary indices: 36
3100
for (unsigned int j = 0; j < 3; j++)
3101
{
3102
for (unsigned int k = 0; k < 3; k++)
3103
{
3104
// Number of operations to compute entry: 4
3105
A[j*3 + k] += FE0[ip][j]*FE0[ip][k]*W3[ip]*det;
3106
}// end loop over 'k'
3107
}// end loop over 'j'
3108
}// end loop over 'ip'
3109
}
3110
3111
};
3112
3113
/// This class defines the interface for the tabulation of
3114
/// an expression evaluated at exactly one point.
3115
3116
class pointmeasure_point_integral_0_1: public ufc::point_integral
3117
{
3118
public:
3119
3120
/// Constructor
3121
pointmeasure_point_integral_0_1() : ufc::point_integral()
3122
{
3123
// Do nothing
3124
}
3125
3126
/// Destructor
3127
virtual ~pointmeasure_point_integral_0_1()
3128
{
3129
// Do nothing
3130
}
3131
3132
/// Tabulate the tensor for the contribution from the local vertex
3133
virtual void tabulate_tensor(double* A,
3134
const double * const * w,
3135
const double* vertex_coordinates,
3136
std::size_t vertex) const
3137
{
3138
// Compute Jacobian
3139
double J[4];
3140
compute_jacobian_triangle_2d(J, vertex_coordinates);
3141
3142
// Compute Jacobian inverse and determinant
3143
double K[4];
3144
double detJ;
3145
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3146
3147
3148
3149
// Get vertices on edge
3150
3151
// Compute scale factor (length of edge scaled by length of reference interval)
3152
3153
3154
// Array of quadrature weights.
3155
static const double W1 = 1.0;
3156
// Quadrature points on the UFC reference element: ()
3157
3158
// Value of basis functions at quadrature points.
3159
static const double FE0_v0[1][3] = \
3160
{{1.0, 0.0, 0.0}};
3161
3162
static const double FE0_v1[1][3] = \
3163
{{0.0, 1.0, 0.0}};
3164
3165
static const double FE0_v2[1][3] = \
3166
{{0.0, 0.0, 1.0}};
3167
3168
// Reset values in the element tensor.
3169
for (unsigned int r = 0; r < 9; r++)
3170
{
3171
A[r] = 0.0;
3172
}// end loop over 'r'
3173
3174
// Compute element tensor using UFL quadrature representation
3175
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3176
switch (vertex)
3177
{
3178
case 0:
3179
{
3180
// Total number of operations to compute element tensor (from this point): 69
3181
3182
// Loop quadrature points for integral.
3183
// Number of operations to compute element tensor for following IP loop = 69
3184
// Only 1 integration point, omitting IP loop.
3185
3186
// Coefficient declarations.
3187
double F0 = 0.0;
3188
3189
// Total number of operations to compute function values = 6
3190
for (unsigned int r = 0; r < 3; r++)
3191
{
3192
F0 += FE0_v0[0][r]*w[0][r];
3193
}// end loop over 'r'
3194
3195
// Number of operations for primary indices: 63
3196
for (unsigned int j = 0; j < 3; j++)
3197
{
3198
for (unsigned int k = 0; k < 3; k++)
3199
{
3200
// Number of operations to compute entry: 7
3201
A[j*3 + k] += (FE0_v0[0][j]*FE0_v0[0][k] + FE0_v0[0][j]*FE0_v0[0][k]*F0*F0)*W1;
3202
}// end loop over 'k'
3203
}// end loop over 'j'
3204
break;
3205
}
3206
case 1:
3207
{
3208
// Total number of operations to compute element tensor (from this point): 69
3209
3210
// Loop quadrature points for integral.
3211
// Number of operations to compute element tensor for following IP loop = 69
3212
// Only 1 integration point, omitting IP loop.
3213
3214
// Coefficient declarations.
3215
double F0 = 0.0;
3216
3217
// Total number of operations to compute function values = 6
3218
for (unsigned int r = 0; r < 3; r++)
3219
{
3220
F0 += FE0_v1[0][r]*w[0][r];
3221
}// end loop over 'r'
3222
3223
// Number of operations for primary indices: 63
3224
for (unsigned int j = 0; j < 3; j++)
3225
{
3226
for (unsigned int k = 0; k < 3; k++)
3227
{
3228
// Number of operations to compute entry: 7
3229
A[j*3 + k] += (FE0_v1[0][j]*FE0_v1[0][k]*F0*F0 + FE0_v1[0][j]*FE0_v1[0][k])*W1;
3230
}// end loop over 'k'
3231
}// end loop over 'j'
3232
break;
3233
}
3234
case 2:
3235
{
3236
// Total number of operations to compute element tensor (from this point): 69
3237
3238
// Loop quadrature points for integral.
3239
// Number of operations to compute element tensor for following IP loop = 69
3240
// Only 1 integration point, omitting IP loop.
3241
3242
// Coefficient declarations.
3243
double F0 = 0.0;
3244
3245
// Total number of operations to compute function values = 6
3246
for (unsigned int r = 0; r < 3; r++)
3247
{
3248
F0 += FE0_v2[0][r]*w[0][r];
3249
}// end loop over 'r'
3250
3251
// Number of operations for primary indices: 63
3252
for (unsigned int j = 0; j < 3; j++)
3253
{
3254
for (unsigned int k = 0; k < 3; k++)
3255
{
3256
// Number of operations to compute entry: 7
3257
A[j*3 + k] += (FE0_v2[0][j]*FE0_v2[0][k]*F0*F0 + FE0_v2[0][j]*FE0_v2[0][k])*W1;
3258
}// end loop over 'k'
3259
}// end loop over 'j'
3260
break;
3261
}
3262
}
3263
3264
}
3265
3266
};
3267
3268
/// This class defines the interface for the tabulation of
3269
/// an expression evaluated at exactly one point.
3270
3271
class pointmeasure_point_integral_0_otherwise: public ufc::point_integral
3272
{
3273
public:
3274
3275
/// Constructor
3276
pointmeasure_point_integral_0_otherwise() : ufc::point_integral()
3277
{
3278
// Do nothing
3279
}
3280
3281
/// Destructor
3282
virtual ~pointmeasure_point_integral_0_otherwise()
3283
{
3284
// Do nothing
3285
}
3286
3287
/// Tabulate the tensor for the contribution from the local vertex
3288
virtual void tabulate_tensor(double* A,
3289
const double * const * w,
3290
const double* vertex_coordinates,
3291
std::size_t vertex) const
3292
{
3293
// Compute Jacobian
3294
double J[4];
3295
compute_jacobian_triangle_2d(J, vertex_coordinates);
3296
3297
// Compute Jacobian inverse and determinant
3298
double K[4];
3299
double detJ;
3300
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3301
3302
3303
3304
// Get vertices on edge
3305
3306
// Compute scale factor (length of edge scaled by length of reference interval)
3307
3308
3309
// Array of quadrature weights.
3310
static const double W1 = 1.0;
3311
// Quadrature points on the UFC reference element: ()
3312
3313
// Value of basis functions at quadrature points.
3314
static const double FE0_v0[1][3] = \
3315
{{1.0, 0.0, 0.0}};
3316
3317
static const double FE0_v1[1][3] = \
3318
{{0.0, 1.0, 0.0}};
3319
3320
static const double FE0_v2[1][3] = \
3321
{{0.0, 0.0, 1.0}};
3322
3323
// Reset values in the element tensor.
3324
for (unsigned int r = 0; r < 9; r++)
3325
{
3326
A[r] = 0.0;
3327
}// end loop over 'r'
3328
3329
// Compute element tensor using UFL quadrature representation
3330
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3331
switch (vertex)
3332
{
3333
case 0:
3334
{
3335
// Total number of operations to compute element tensor (from this point): 27
3336
3337
// Loop quadrature points for integral.
3338
// Number of operations to compute element tensor for following IP loop = 27
3339
// Only 1 integration point, omitting IP loop.
3340
3341
// Number of operations for primary indices: 27
3342
for (unsigned int j = 0; j < 3; j++)
3343
{
3344
for (unsigned int k = 0; k < 3; k++)
3345
{
3346
// Number of operations to compute entry: 3
3347
A[j*3 + k] += FE0_v0[0][j]*FE0_v0[0][k]*W1;
3348
}// end loop over 'k'
3349
}// end loop over 'j'
3350
break;
3351
}
3352
case 1:
3353
{
3354
// Total number of operations to compute element tensor (from this point): 27
3355
3356
// Loop quadrature points for integral.
3357
// Number of operations to compute element tensor for following IP loop = 27
3358
// Only 1 integration point, omitting IP loop.
3359
3360
// Number of operations for primary indices: 27
3361
for (unsigned int j = 0; j < 3; j++)
3362
{
3363
for (unsigned int k = 0; k < 3; k++)
3364
{
3365
// Number of operations to compute entry: 3
3366
A[j*3 + k] += FE0_v1[0][j]*FE0_v1[0][k]*W1;
3367
}// end loop over 'k'
3368
}// end loop over 'j'
3369
break;
3370
}
3371
case 2:
3372
{
3373
// Total number of operations to compute element tensor (from this point): 27
3374
3375
// Loop quadrature points for integral.
3376
// Number of operations to compute element tensor for following IP loop = 27
3377
// Only 1 integration point, omitting IP loop.
3378
3379
// Number of operations for primary indices: 27
3380
for (unsigned int j = 0; j < 3; j++)
3381
{
3382
for (unsigned int k = 0; k < 3; k++)
3383
{
3384
// Number of operations to compute entry: 3
3385
A[j*3 + k] += FE0_v2[0][j]*FE0_v2[0][k]*W1;
3386
}// end loop over 'k'
3387
}// end loop over 'j'
3388
break;
3389
}
3390
}
3391
3392
}
3393
3394
};
3395
3396
/// This class defines the interface for the tabulation of
3397
/// an expression evaluated at exactly one point.
3398
3399
class pointmeasure_point_integral_1_otherwise: public ufc::point_integral
3400
{
3401
public:
3402
3403
/// Constructor
3404
pointmeasure_point_integral_1_otherwise() : ufc::point_integral()
3405
{
3406
// Do nothing
3407
}
3408
3409
/// Destructor
3410
virtual ~pointmeasure_point_integral_1_otherwise()
3411
{
3412
// Do nothing
3413
}
3414
3415
/// Tabulate the tensor for the contribution from the local vertex
3416
virtual void tabulate_tensor(double* A,
3417
const double * const * w,
3418
const double* vertex_coordinates,
3419
std::size_t vertex) const
3420
{
3421
// Compute Jacobian
3422
double J[4];
3423
compute_jacobian_triangle_2d(J, vertex_coordinates);
3424
3425
// Compute Jacobian inverse and determinant
3426
double K[4];
3427
double detJ;
3428
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3429
3430
3431
3432
// Get vertices on edge
3433
3434
// Compute scale factor (length of edge scaled by length of reference interval)
3435
3436
3437
// Array of quadrature weights.
3438
static const double W1 = 1.0;
3439
// Quadrature points on the UFC reference element: ()
3440
3441
// Value of basis functions at quadrature points.
3442
static const double FE0_v0[1][3] = \
3443
{{1.0, 0.0, 0.0}};
3444
3445
static const double FE0_v1[1][3] = \
3446
{{0.0, 1.0, 0.0}};
3447
3448
static const double FE0_v2[1][3] = \
3449
{{0.0, 0.0, 1.0}};
3450
3451
static const double FE1_v0[1][6] = \
3452
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
3453
3454
static const double FE1_v1[1][6] = \
3455
{{0.0, 1.0, 0.0, 0.0, 0.0, 0.0}};
3456
3457
static const double FE1_v2[1][6] = \
3458
{{0.0, 0.0, 1.0, 0.0, 0.0, 0.0}};
3459
3460
// Reset values in the element tensor.
3461
for (unsigned int r = 0; r < 3; r++)
3462
{
3463
A[r] = 0.0;
3464
}// end loop over 'r'
3465
3466
// Compute element tensor using UFL quadrature representation
3467
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3468
switch (vertex)
3469
{
3470
case 0:
3471
{
3472
// Total number of operations to compute element tensor (from this point): 21
3473
3474
// Loop quadrature points for integral.
3475
// Number of operations to compute element tensor for following IP loop = 21
3476
// Only 1 integration point, omitting IP loop.
3477
3478
// Coefficient declarations.
3479
double F0 = 0.0;
3480
3481
// Total number of operations to compute function values = 12
3482
for (unsigned int r = 0; r < 6; r++)
3483
{
3484
F0 += FE1_v0[0][r]*w[0][r];
3485
}// end loop over 'r'
3486
3487
// Number of operations for primary indices: 9
3488
for (unsigned int j = 0; j < 3; j++)
3489
{
3490
// Number of operations to compute entry: 3
3491
A[j] += FE0_v0[0][j]*F0*W1;
3492
}// end loop over 'j'
3493
break;
3494
}
3495
case 1:
3496
{
3497
// Total number of operations to compute element tensor (from this point): 21
3498
3499
// Loop quadrature points for integral.
3500
// Number of operations to compute element tensor for following IP loop = 21
3501
// Only 1 integration point, omitting IP loop.
3502
3503
// Coefficient declarations.
3504
double F0 = 0.0;
3505
3506
// Total number of operations to compute function values = 12
3507
for (unsigned int r = 0; r < 6; r++)
3508
{
3509
F0 += FE1_v1[0][r]*w[0][r];
3510
}// end loop over 'r'
3511
3512
// Number of operations for primary indices: 9
3513
for (unsigned int j = 0; j < 3; j++)
3514
{
3515
// Number of operations to compute entry: 3
3516
A[j] += FE0_v1[0][j]*F0*W1;
3517
}// end loop over 'j'
3518
break;
3519
}
3520
case 2:
3521
{
3522
// Total number of operations to compute element tensor (from this point): 21
3523
3524
// Loop quadrature points for integral.
3525
// Number of operations to compute element tensor for following IP loop = 21
3526
// Only 1 integration point, omitting IP loop.
3527
3528
// Coefficient declarations.
3529
double F0 = 0.0;
3530
3531
// Total number of operations to compute function values = 12
3532
for (unsigned int r = 0; r < 6; r++)
3533
{
3534
F0 += FE1_v2[0][r]*w[0][r];
3535
}// end loop over 'r'
3536
3537
// Number of operations for primary indices: 9
3538
for (unsigned int j = 0; j < 3; j++)
3539
{
3540
// Number of operations to compute entry: 3
3541
A[j] += FE0_v2[0][j]*F0*W1;
3542
}// end loop over 'j'
3543
break;
3544
}
3545
}
3546
3547
}
3548
3549
};
3550
3551
/// This class defines the interface for the assembly of the global
3552
/// tensor corresponding to a form with r + n arguments, that is, a
3553
/// mapping
3554
///
3555
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3556
///
3557
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3558
/// global tensor A is defined by
3559
///
3560
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3561
///
3562
/// where each argument Vj represents the application to the
3563
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3564
/// fixed functions (coefficients).
3565
3566
class pointmeasure_form_0: public ufc::form
3567
{
3568
public:
3569
3570
/// Constructor
3571
pointmeasure_form_0() : ufc::form()
3572
{
3573
// Do nothing
3574
}
3575
3576
/// Destructor
3577
virtual ~pointmeasure_form_0()
3578
{
3579
// Do nothing
3580
}
3581
3582
/// Return a string identifying the form
3583
virtual const char* signature() const
3584
{
3585
return "856f5aa1a8f7253511ba78602090dc3cf582e71ee4c859de2e80d2fd37c92108a283e09a2ae620c72f13f280a38368ed9c5e78e015492db3cbebd6d8f66016fe";
3586
}
3587
3588
/// Return the rank of the global tensor (r)
3589
virtual std::size_t rank() const
3590
{
3591
return 2;
3592
}
3593
3594
/// Return the number of coefficients (n)
3595
virtual std::size_t num_coefficients() const
3596
{
3597
return 1;
3598
}
3599
3600
/// Return the number of cell domains
3601
virtual std::size_t num_cell_domains() const
3602
{
3603
return 0;
3604
}
3605
3606
/// Return the number of exterior facet domains
3607
virtual std::size_t num_exterior_facet_domains() const
3608
{
3609
return 0;
3610
}
3611
3612
/// Return the number of interior facet domains
3613
virtual std::size_t num_interior_facet_domains() const
3614
{
3615
return 0;
3616
}
3617
3618
/// Return the number of point domains
3619
virtual std::size_t num_point_domains() const
3620
{
3621
return 2;
3622
}
3623
3624
/// Return whether the form has any cell integrals
3625
virtual bool has_cell_integrals() const
3626
{
3627
return true;
3628
}
3629
3630
/// Return whether the form has any exterior facet integrals
3631
virtual bool has_exterior_facet_integrals() const
3632
{
3633
return false;
3634
}
3635
3636
/// Return whether the form has any interior facet integrals
3637
virtual bool has_interior_facet_integrals() const
3638
{
3639
return false;
3640
}
3641
3642
/// Return whether the form has any point integrals
3643
virtual bool has_point_integrals() const
3644
{
3645
return true;
3646
}
3647
3648
/// Create a new finite element for argument function i
3649
virtual ufc::finite_element* create_finite_element(std::size_t i) const
3650
{
3651
switch (i)
3652
{
3653
case 0:
3654
{
3655
return new pointmeasure_finite_element_1();
3656
break;
3657
}
3658
case 1:
3659
{
3660
return new pointmeasure_finite_element_1();
3661
break;
3662
}
3663
case 2:
3664
{
3665
return new pointmeasure_finite_element_1();
3666
break;
3667
}
3668
}
3669
3670
return 0;
3671
}
3672
3673
/// Create a new dofmap for argument function i
3674
virtual ufc::dofmap* create_dofmap(std::size_t i) const
3675
{
3676
switch (i)
3677
{
3678
case 0:
3679
{
3680
return new pointmeasure_dofmap_1();
3681
break;
3682
}
3683
case 1:
3684
{
3685
return new pointmeasure_dofmap_1();
3686
break;
3687
}
3688
case 2:
3689
{
3690
return new pointmeasure_dofmap_1();
3691
break;
3692
}
3693
}
3694
3695
return 0;
3696
}
3697
3698
/// Create a new cell integral on sub domain i
3699
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
3700
{
3701
return 0;
3702
}
3703
3704
/// Create a new exterior facet integral on sub domain i
3705
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
3706
{
3707
return 0;
3708
}
3709
3710
/// Create a new interior facet integral on sub domain i
3711
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
3712
{
3713
return 0;
3714
}
3715
3716
/// Create a new point integral on sub domain i
3717
virtual ufc::point_integral* create_point_integral(std::size_t i) const
3718
{
3719
switch (i)
3720
{
3721
case 1:
3722
{
3723
return new pointmeasure_point_integral_0_1();
3724
break;
3725
}
3726
}
3727
3728
return 0;
3729
}
3730
3731
/// Create a new cell integral on everywhere else
3732
virtual ufc::cell_integral* create_default_cell_integral() const
3733
{
3734
return new pointmeasure_cell_integral_0_otherwise();
3735
}
3736
3737
/// Create a new exterior facet integral on everywhere else
3738
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
3739
{
3740
return 0;
3741
}
3742
3743
/// Create a new interior facet integral on everywhere else
3744
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
3745
{
3746
return 0;
3747
}
3748
3749
/// Create a new point integral on everywhere else
3750
virtual ufc::point_integral* create_default_point_integral() const
3751
{
3752
return new pointmeasure_point_integral_0_otherwise();
3753
}
3754
3755
};
3756
3757
/// This class defines the interface for the assembly of the global
3758
/// tensor corresponding to a form with r + n arguments, that is, a
3759
/// mapping
3760
///
3761
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3762
///
3763
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3764
/// global tensor A is defined by
3765
///
3766
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3767
///
3768
/// where each argument Vj represents the application to the
3769
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3770
/// fixed functions (coefficients).
3771
3772
class pointmeasure_form_1: public ufc::form
3773
{
3774
public:
3775
3776
/// Constructor
3777
pointmeasure_form_1() : ufc::form()
3778
{
3779
// Do nothing
3780
}
3781
3782
/// Destructor
3783
virtual ~pointmeasure_form_1()
3784
{
3785
// Do nothing
3786
}
3787
3788
/// Return a string identifying the form
3789
virtual const char* signature() const
3790
{
3791
return "cbc400c8cb006d3574d6a58223a2c80d90b230d781feaabb65968d87ac0550f09264f2d1f32621512f1d42a7bc048c03d38e4482a1abf0fba8fad9c1cab8cc38";
3792
}
3793
3794
/// Return the rank of the global tensor (r)
3795
virtual std::size_t rank() const
3796
{
3797
return 1;
3798
}
3799
3800
/// Return the number of coefficients (n)
3801
virtual std::size_t num_coefficients() const
3802
{
3803
return 1;
3804
}
3805
3806
/// Return the number of cell domains
3807
virtual std::size_t num_cell_domains() const
3808
{
3809
return 0;
3810
}
3811
3812
/// Return the number of exterior facet domains
3813
virtual std::size_t num_exterior_facet_domains() const
3814
{
3815
return 0;
3816
}
3817
3818
/// Return the number of interior facet domains
3819
virtual std::size_t num_interior_facet_domains() const
3820
{
3821
return 0;
3822
}
3823
3824
/// Return the number of point domains
3825
virtual std::size_t num_point_domains() const
3826
{
3827
return 0;
3828
}
3829
3830
/// Return whether the form has any cell integrals
3831
virtual bool has_cell_integrals() const
3832
{
3833
return false;
3834
}
3835
3836
/// Return whether the form has any exterior facet integrals
3837
virtual bool has_exterior_facet_integrals() const
3838
{
3839
return false;
3840
}
3841
3842
/// Return whether the form has any interior facet integrals
3843
virtual bool has_interior_facet_integrals() const
3844
{
3845
return false;
3846
}
3847
3848
/// Return whether the form has any point integrals
3849
virtual bool has_point_integrals() const
3850
{
3851
return true;
3852
}
3853
3854
/// Create a new finite element for argument function i
3855
virtual ufc::finite_element* create_finite_element(std::size_t i) const
3856
{
3857
switch (i)
3858
{
3859
case 0:
3860
{
3861
return new pointmeasure_finite_element_1();
3862
break;
3863
}
3864
case 1:
3865
{
3866
return new pointmeasure_finite_element_0();
3867
break;
3868
}
3869
}
3870
3871
return 0;
3872
}
3873
3874
/// Create a new dofmap for argument function i
3875
virtual ufc::dofmap* create_dofmap(std::size_t i) const
3876
{
3877
switch (i)
3878
{
3879
case 0:
3880
{
3881
return new pointmeasure_dofmap_1();
3882
break;
3883
}
3884
case 1:
3885
{
3886
return new pointmeasure_dofmap_0();
3887
break;
3888
}
3889
}
3890
3891
return 0;
3892
}
3893
3894
/// Create a new cell integral on sub domain i
3895
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
3896
{
3897
return 0;
3898
}
3899
3900
/// Create a new exterior facet integral on sub domain i
3901
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
3902
{
3903
return 0;
3904
}
3905
3906
/// Create a new interior facet integral on sub domain i
3907
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
3908
{
3909
return 0;
3910
}
3911
3912
/// Create a new point integral on sub domain i
3913
virtual ufc::point_integral* create_point_integral(std::size_t i) const
3914
{
3915
return 0;
3916
}
3917
3918
/// Create a new cell integral on everywhere else
3919
virtual ufc::cell_integral* create_default_cell_integral() const
3920
{
3921
return 0;
3922
}
3923
3924
/// Create a new exterior facet integral on everywhere else
3925
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
3926
{
3927
return 0;
3928
}
3929
3930
/// Create a new interior facet integral on everywhere else
3931
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
3932
{
3933
return 0;
3934
}
3935
3936
/// Create a new point integral on everywhere else
3937
virtual ufc::point_integral* create_default_point_integral() const
3938
{
3939
return new pointmeasure_point_integral_1_otherwise();
3940
}
3941
3942
};
3943
3944
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1