~ubuntu-branches/ubuntu/vivid/ffc/vivid
» Revision
18
: test/regression/references/r_quadrature/MixedPoisson.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_quadrature/MixedPoisson.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/MixedPoisson.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 __MIXEDPOISSON_H
27
#define __MIXEDPOISSON_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 mixedpoisson_finite_element_0: public ufc::finite_element
37
{
38
public:
39
40
/// Constructor
41
mixedpoisson_finite_element_0() : ufc::finite_element()
42
{
43
// Do nothing
44
}
45
46
/// Destructor
47
virtual ~mixedpoisson_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('Brezzi-Douglas-Marini', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, 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 1;
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
switch (i)
92
{
93
case 0:
94
{
95
return 2;
96
break;
97
}
98
}
99
100
return 0;
101
}
102
103
/// Evaluate basis function i at given point x in cell
104
virtual void evaluate_basis(std::size_t i,
105
double* values,
106
const double* x,
107
const double* vertex_coordinates,
108
int cell_orientation) const
109
{
110
// Compute Jacobian
111
double J[4];
112
compute_jacobian_triangle_2d(J, vertex_coordinates);
113
114
// Compute Jacobian inverse and determinant
115
double K[4];
116
double detJ;
117
compute_jacobian_inverse_triangle_2d(K, detJ, J);
118
119
120
121
// Compute constants
122
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
123
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
124
125
// Get coordinates and map to the reference (FIAT) element
126
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
127
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
128
129
// Reset values
130
values[0] = 0.0;
131
values[1] = 0.0;
132
switch (i)
133
{
134
case 0:
135
{
136
137
// Array of basisvalues
138
double basisvalues[3] = {0.0, 0.0, 0.0};
139
140
// Declare helper variables
141
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
142
143
// Compute basisvalues
144
basisvalues[0] = 1.0;
145
basisvalues[1] = tmp0;
146
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
147
basisvalues[0] *= std::sqrt(0.5);
148
basisvalues[2] *= std::sqrt(1.0);
149
basisvalues[1] *= std::sqrt(3.0);
150
151
// Table(s) of coefficients
152
static const double coefficients0[3] = \
153
{0.94280904, 0.57735027, -0.33333333};
154
155
static const double coefficients1[3] = \
156
{-0.47140452, 0.0, -0.33333333};
157
158
// Compute value(s)
159
for (unsigned int r = 0; r < 3; r++)
160
{
161
values[0] += coefficients0[r]*basisvalues[r];
162
values[1] += coefficients1[r]*basisvalues[r];
163
}// end loop over 'r'
164
165
// Using contravariant Piola transform to map values back to the physical element
166
const double tmp_ref0 = values[0];
167
const double tmp_ref1 = values[1];
168
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
169
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
170
break;
171
}
172
case 1:
173
{
174
175
// Array of basisvalues
176
double basisvalues[3] = {0.0, 0.0, 0.0};
177
178
// Declare helper variables
179
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
180
181
// Compute basisvalues
182
basisvalues[0] = 1.0;
183
basisvalues[1] = tmp0;
184
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
185
basisvalues[0] *= std::sqrt(0.5);
186
basisvalues[2] *= std::sqrt(1.0);
187
basisvalues[1] *= std::sqrt(3.0);
188
189
// Table(s) of coefficients
190
static const double coefficients0[3] = \
191
{-0.47140452, -0.28867513, 0.16666667};
192
193
static const double coefficients1[3] = \
194
{0.94280904, 0.0, 0.66666667};
195
196
// Compute value(s)
197
for (unsigned int r = 0; r < 3; r++)
198
{
199
values[0] += coefficients0[r]*basisvalues[r];
200
values[1] += coefficients1[r]*basisvalues[r];
201
}// end loop over 'r'
202
203
// Using contravariant Piola transform to map values back to the physical element
204
const double tmp_ref0 = values[0];
205
const double tmp_ref1 = values[1];
206
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
207
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
208
break;
209
}
210
case 2:
211
{
212
213
// Array of basisvalues
214
double basisvalues[3] = {0.0, 0.0, 0.0};
215
216
// Declare helper variables
217
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
218
219
// Compute basisvalues
220
basisvalues[0] = 1.0;
221
basisvalues[1] = tmp0;
222
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
223
basisvalues[0] *= std::sqrt(0.5);
224
basisvalues[2] *= std::sqrt(1.0);
225
basisvalues[1] *= std::sqrt(3.0);
226
227
// Table(s) of coefficients
228
static const double coefficients0[3] = \
229
{0.47140452, -0.57735027, -0.66666667};
230
231
static const double coefficients1[3] = \
232
{0.47140452, 0.0, 0.33333333};
233
234
// Compute value(s)
235
for (unsigned int r = 0; r < 3; r++)
236
{
237
values[0] += coefficients0[r]*basisvalues[r];
238
values[1] += coefficients1[r]*basisvalues[r];
239
}// end loop over 'r'
240
241
// Using contravariant Piola transform to map values back to the physical element
242
const double tmp_ref0 = values[0];
243
const double tmp_ref1 = values[1];
244
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
245
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
246
break;
247
}
248
case 3:
249
{
250
251
// Array of basisvalues
252
double basisvalues[3] = {0.0, 0.0, 0.0};
253
254
// Declare helper variables
255
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
256
257
// Compute basisvalues
258
basisvalues[0] = 1.0;
259
basisvalues[1] = tmp0;
260
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
261
basisvalues[0] *= std::sqrt(0.5);
262
basisvalues[2] *= std::sqrt(1.0);
263
basisvalues[1] *= std::sqrt(3.0);
264
265
// Table(s) of coefficients
266
static const double coefficients0[3] = \
267
{0.47140452, 0.28867513, 0.83333333};
268
269
static const double coefficients1[3] = \
270
{-0.94280904, 0.0, -0.66666667};
271
272
// Compute value(s)
273
for (unsigned int r = 0; r < 3; r++)
274
{
275
values[0] += coefficients0[r]*basisvalues[r];
276
values[1] += coefficients1[r]*basisvalues[r];
277
}// end loop over 'r'
278
279
// Using contravariant Piola transform to map values back to the physical element
280
const double tmp_ref0 = values[0];
281
const double tmp_ref1 = values[1];
282
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
283
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
284
break;
285
}
286
case 4:
287
{
288
289
// Array of basisvalues
290
double basisvalues[3] = {0.0, 0.0, 0.0};
291
292
// Declare helper variables
293
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
294
295
// Compute basisvalues
296
basisvalues[0] = 1.0;
297
basisvalues[1] = tmp0;
298
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
299
basisvalues[0] *= std::sqrt(0.5);
300
basisvalues[2] *= std::sqrt(1.0);
301
basisvalues[1] *= std::sqrt(3.0);
302
303
// Table(s) of coefficients
304
static const double coefficients0[3] = \
305
{-0.47140452, -0.28867513, 0.16666667};
306
307
static const double coefficients1[3] = \
308
{-0.47140452, 0.8660254, 0.16666667};
309
310
// Compute value(s)
311
for (unsigned int r = 0; r < 3; r++)
312
{
313
values[0] += coefficients0[r]*basisvalues[r];
314
values[1] += coefficients1[r]*basisvalues[r];
315
}// end loop over 'r'
316
317
// Using contravariant Piola transform to map values back to the physical element
318
const double tmp_ref0 = values[0];
319
const double tmp_ref1 = values[1];
320
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
321
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
322
break;
323
}
324
case 5:
325
{
326
327
// Array of basisvalues
328
double basisvalues[3] = {0.0, 0.0, 0.0};
329
330
// Declare helper variables
331
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
332
333
// Compute basisvalues
334
basisvalues[0] = 1.0;
335
basisvalues[1] = tmp0;
336
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
337
basisvalues[0] *= std::sqrt(0.5);
338
basisvalues[2] *= std::sqrt(1.0);
339
basisvalues[1] *= std::sqrt(3.0);
340
341
// Table(s) of coefficients
342
static const double coefficients0[3] = \
343
{0.94280904, 0.57735027, -0.33333333};
344
345
static const double coefficients1[3] = \
346
{-0.47140452, -0.8660254, 0.16666667};
347
348
// Compute value(s)
349
for (unsigned int r = 0; r < 3; r++)
350
{
351
values[0] += coefficients0[r]*basisvalues[r];
352
values[1] += coefficients1[r]*basisvalues[r];
353
}// end loop over 'r'
354
355
// Using contravariant Piola transform to map values back to the physical element
356
const double tmp_ref0 = values[0];
357
const double tmp_ref1 = values[1];
358
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
359
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
360
break;
361
}
362
}
363
364
}
365
366
/// Evaluate all basis functions at given point x in cell
367
virtual void evaluate_basis_all(double* values,
368
const double* x,
369
const double* vertex_coordinates,
370
int cell_orientation) const
371
{
372
// Helper variable to hold values of a single dof.
373
double dof_values[2] = {0.0, 0.0};
374
375
// Loop dofs and call evaluate_basis
376
for (unsigned int r = 0; r < 6; r++)
377
{
378
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
379
for (unsigned int s = 0; s < 2; s++)
380
{
381
values[r*2 + s] = dof_values[s];
382
}// end loop over 's'
383
}// end loop over 'r'
384
}
385
386
/// Evaluate order n derivatives of basis function i at given point x in cell
387
virtual void evaluate_basis_derivatives(std::size_t i,
388
std::size_t n,
389
double* values,
390
const double* x,
391
const double* vertex_coordinates,
392
int cell_orientation) const
393
{
394
// Compute Jacobian
395
double J[4];
396
compute_jacobian_triangle_2d(J, vertex_coordinates);
397
398
// Compute Jacobian inverse and determinant
399
double K[4];
400
double detJ;
401
compute_jacobian_inverse_triangle_2d(K, detJ, J);
402
403
404
405
// Compute constants
406
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
407
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
408
409
// Get coordinates and map to the reference (FIAT) element
410
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
411
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
412
413
// Compute number of derivatives.
414
unsigned int num_derivatives = 1;
415
for (unsigned int r = 0; r < n; r++)
416
{
417
num_derivatives *= 2;
418
}// end loop over 'r'
419
420
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
421
unsigned int **combinations = new unsigned int *[num_derivatives];
422
for (unsigned int row = 0; row < num_derivatives; row++)
423
{
424
combinations[row] = new unsigned int [n];
425
for (unsigned int col = 0; col < n; col++)
426
combinations[row][col] = 0;
427
}
428
429
// Generate combinations of derivatives
430
for (unsigned int row = 1; row < num_derivatives; row++)
431
{
432
for (unsigned int num = 0; num < row; num++)
433
{
434
for (unsigned int col = n-1; col+1 > 0; col--)
435
{
436
if (combinations[row][col] + 1 > 1)
437
combinations[row][col] = 0;
438
else
439
{
440
combinations[row][col] += 1;
441
break;
442
}
443
}
444
}
445
}
446
447
// Compute inverse of Jacobian
448
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
449
450
// Declare transformation matrix
451
// Declare pointer to two dimensional array and initialise
452
double **transform = new double *[num_derivatives];
453
454
for (unsigned int j = 0; j < num_derivatives; j++)
455
{
456
transform[j] = new double [num_derivatives];
457
for (unsigned int k = 0; k < num_derivatives; k++)
458
transform[j][k] = 1;
459
}
460
461
// Construct transformation matrix
462
for (unsigned int row = 0; row < num_derivatives; row++)
463
{
464
for (unsigned int col = 0; col < num_derivatives; col++)
465
{
466
for (unsigned int k = 0; k < n; k++)
467
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
468
}
469
}
470
471
// Reset values. Assuming that values is always an array.
472
for (unsigned int r = 0; r < 2*num_derivatives; r++)
473
{
474
values[r] = 0.0;
475
}// end loop over 'r'
476
477
switch (i)
478
{
479
case 0:
480
{
481
482
// Array of basisvalues
483
double basisvalues[3] = {0.0, 0.0, 0.0};
484
485
// Declare helper variables
486
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
487
488
// Compute basisvalues
489
basisvalues[0] = 1.0;
490
basisvalues[1] = tmp0;
491
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
492
basisvalues[0] *= std::sqrt(0.5);
493
basisvalues[2] *= std::sqrt(1.0);
494
basisvalues[1] *= std::sqrt(3.0);
495
496
// Table(s) of coefficients
497
static const double coefficients0[3] = \
498
{0.94280904, 0.57735027, -0.33333333};
499
500
static const double coefficients1[3] = \
501
{-0.47140452, 0.0, -0.33333333};
502
503
// Tables of derivatives of the polynomial base (transpose).
504
static const double dmats0[3][3] = \
505
{{0.0, 0.0, 0.0},
506
{4.8989795, 0.0, 0.0},
507
{0.0, 0.0, 0.0}};
508
509
static const double dmats1[3][3] = \
510
{{0.0, 0.0, 0.0},
511
{2.4494897, 0.0, 0.0},
512
{4.2426407, 0.0, 0.0}};
513
514
// Compute reference derivatives.
515
// Declare pointer to array of derivatives on FIAT element.
516
double *derivatives = new double[2*num_derivatives];
517
for (unsigned int r = 0; r < 2*num_derivatives; r++)
518
{
519
derivatives[r] = 0.0;
520
}// end loop over 'r'
521
522
// Declare pointer to array of reference derivatives on physical element.
523
double *derivatives_p = new double[2*num_derivatives];
524
for (unsigned int r = 0; r < 2*num_derivatives; r++)
525
{
526
derivatives_p[r] = 0.0;
527
}// end loop over 'r'
528
529
// Declare derivative matrix (of polynomial basis).
530
double dmats[3][3] = \
531
{{1.0, 0.0, 0.0},
532
{0.0, 1.0, 0.0},
533
{0.0, 0.0, 1.0}};
534
535
// Declare (auxiliary) derivative matrix (of polynomial basis).
536
double dmats_old[3][3] = \
537
{{1.0, 0.0, 0.0},
538
{0.0, 1.0, 0.0},
539
{0.0, 0.0, 1.0}};
540
541
// Loop possible derivatives.
542
for (unsigned int r = 0; r < num_derivatives; r++)
543
{
544
// Resetting dmats values to compute next derivative.
545
for (unsigned int t = 0; t < 3; t++)
546
{
547
for (unsigned int u = 0; u < 3; u++)
548
{
549
dmats[t][u] = 0.0;
550
if (t == u)
551
{
552
dmats[t][u] = 1.0;
553
}
554
555
}// end loop over 'u'
556
}// end loop over 't'
557
558
// Looping derivative order to generate dmats.
559
for (unsigned int s = 0; s < n; s++)
560
{
561
// Updating dmats_old with new values and resetting dmats.
562
for (unsigned int t = 0; t < 3; t++)
563
{
564
for (unsigned int u = 0; u < 3; u++)
565
{
566
dmats_old[t][u] = dmats[t][u];
567
dmats[t][u] = 0.0;
568
}// end loop over 'u'
569
}// end loop over 't'
570
571
// Update dmats using an inner product.
572
if (combinations[r][s] == 0)
573
{
574
for (unsigned int t = 0; t < 3; t++)
575
{
576
for (unsigned int u = 0; u < 3; u++)
577
{
578
for (unsigned int tu = 0; tu < 3; tu++)
579
{
580
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
581
}// end loop over 'tu'
582
}// end loop over 'u'
583
}// end loop over 't'
584
}
585
586
if (combinations[r][s] == 1)
587
{
588
for (unsigned int t = 0; t < 3; t++)
589
{
590
for (unsigned int u = 0; u < 3; u++)
591
{
592
for (unsigned int tu = 0; tu < 3; tu++)
593
{
594
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
595
}// end loop over 'tu'
596
}// end loop over 'u'
597
}// end loop over 't'
598
}
599
600
}// end loop over 's'
601
for (unsigned int s = 0; s < 3; s++)
602
{
603
for (unsigned int t = 0; t < 3; t++)
604
{
605
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
606
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
607
}// end loop over 't'
608
}// end loop over 's'
609
610
// Using contravariant Piola transform to map values back to the physical element.
611
const double tmp_ref0 = derivatives[r];
612
const double tmp_ref1 = derivatives[num_derivatives + r];
613
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
614
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
615
}// end loop over 'r'
616
617
// Transform derivatives back to physical element
618
for (unsigned int r = 0; r < num_derivatives; r++)
619
{
620
for (unsigned int s = 0; s < num_derivatives; s++)
621
{
622
values[r] += transform[r][s]*derivatives_p[s];
623
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
624
}// end loop over 's'
625
}// end loop over 'r'
626
627
// Delete pointer to array of derivatives on FIAT element
628
delete [] derivatives;
629
630
// Delete pointer to array of reference derivatives on physical element.
631
delete [] derivatives_p;
632
633
// Delete pointer to array of combinations of derivatives and transform
634
for (unsigned int r = 0; r < num_derivatives; r++)
635
{
636
delete [] combinations[r];
637
}// end loop over 'r'
638
delete [] combinations;
639
for (unsigned int r = 0; r < num_derivatives; r++)
640
{
641
delete [] transform[r];
642
}// end loop over 'r'
643
delete [] transform;
644
break;
645
}
646
case 1:
647
{
648
649
// Array of basisvalues
650
double basisvalues[3] = {0.0, 0.0, 0.0};
651
652
// Declare helper variables
653
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
654
655
// Compute basisvalues
656
basisvalues[0] = 1.0;
657
basisvalues[1] = tmp0;
658
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
659
basisvalues[0] *= std::sqrt(0.5);
660
basisvalues[2] *= std::sqrt(1.0);
661
basisvalues[1] *= std::sqrt(3.0);
662
663
// Table(s) of coefficients
664
static const double coefficients0[3] = \
665
{-0.47140452, -0.28867513, 0.16666667};
666
667
static const double coefficients1[3] = \
668
{0.94280904, 0.0, 0.66666667};
669
670
// Tables of derivatives of the polynomial base (transpose).
671
static const double dmats0[3][3] = \
672
{{0.0, 0.0, 0.0},
673
{4.8989795, 0.0, 0.0},
674
{0.0, 0.0, 0.0}};
675
676
static const double dmats1[3][3] = \
677
{{0.0, 0.0, 0.0},
678
{2.4494897, 0.0, 0.0},
679
{4.2426407, 0.0, 0.0}};
680
681
// Compute reference derivatives.
682
// Declare pointer to array of derivatives on FIAT element.
683
double *derivatives = new double[2*num_derivatives];
684
for (unsigned int r = 0; r < 2*num_derivatives; r++)
685
{
686
derivatives[r] = 0.0;
687
}// end loop over 'r'
688
689
// Declare pointer to array of reference derivatives on physical element.
690
double *derivatives_p = new double[2*num_derivatives];
691
for (unsigned int r = 0; r < 2*num_derivatives; r++)
692
{
693
derivatives_p[r] = 0.0;
694
}// end loop over 'r'
695
696
// Declare derivative matrix (of polynomial basis).
697
double dmats[3][3] = \
698
{{1.0, 0.0, 0.0},
699
{0.0, 1.0, 0.0},
700
{0.0, 0.0, 1.0}};
701
702
// Declare (auxiliary) derivative matrix (of polynomial basis).
703
double dmats_old[3][3] = \
704
{{1.0, 0.0, 0.0},
705
{0.0, 1.0, 0.0},
706
{0.0, 0.0, 1.0}};
707
708
// Loop possible derivatives.
709
for (unsigned int r = 0; r < num_derivatives; r++)
710
{
711
// Resetting dmats values to compute next derivative.
712
for (unsigned int t = 0; t < 3; t++)
713
{
714
for (unsigned int u = 0; u < 3; u++)
715
{
716
dmats[t][u] = 0.0;
717
if (t == u)
718
{
719
dmats[t][u] = 1.0;
720
}
721
722
}// end loop over 'u'
723
}// end loop over 't'
724
725
// Looping derivative order to generate dmats.
726
for (unsigned int s = 0; s < n; s++)
727
{
728
// Updating dmats_old with new values and resetting dmats.
729
for (unsigned int t = 0; t < 3; t++)
730
{
731
for (unsigned int u = 0; u < 3; u++)
732
{
733
dmats_old[t][u] = dmats[t][u];
734
dmats[t][u] = 0.0;
735
}// end loop over 'u'
736
}// end loop over 't'
737
738
// Update dmats using an inner product.
739
if (combinations[r][s] == 0)
740
{
741
for (unsigned int t = 0; t < 3; t++)
742
{
743
for (unsigned int u = 0; u < 3; u++)
744
{
745
for (unsigned int tu = 0; tu < 3; tu++)
746
{
747
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
748
}// end loop over 'tu'
749
}// end loop over 'u'
750
}// end loop over 't'
751
}
752
753
if (combinations[r][s] == 1)
754
{
755
for (unsigned int t = 0; t < 3; t++)
756
{
757
for (unsigned int u = 0; u < 3; u++)
758
{
759
for (unsigned int tu = 0; tu < 3; tu++)
760
{
761
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
762
}// end loop over 'tu'
763
}// end loop over 'u'
764
}// end loop over 't'
765
}
766
767
}// end loop over 's'
768
for (unsigned int s = 0; s < 3; s++)
769
{
770
for (unsigned int t = 0; t < 3; t++)
771
{
772
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
773
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
774
}// end loop over 't'
775
}// end loop over 's'
776
777
// Using contravariant Piola transform to map values back to the physical element.
778
const double tmp_ref0 = derivatives[r];
779
const double tmp_ref1 = derivatives[num_derivatives + r];
780
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
781
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
782
}// end loop over 'r'
783
784
// Transform derivatives back to physical element
785
for (unsigned int r = 0; r < num_derivatives; r++)
786
{
787
for (unsigned int s = 0; s < num_derivatives; s++)
788
{
789
values[r] += transform[r][s]*derivatives_p[s];
790
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
791
}// end loop over 's'
792
}// end loop over 'r'
793
794
// Delete pointer to array of derivatives on FIAT element
795
delete [] derivatives;
796
797
// Delete pointer to array of reference derivatives on physical element.
798
delete [] derivatives_p;
799
800
// Delete pointer to array of combinations of derivatives and transform
801
for (unsigned int r = 0; r < num_derivatives; r++)
802
{
803
delete [] combinations[r];
804
}// end loop over 'r'
805
delete [] combinations;
806
for (unsigned int r = 0; r < num_derivatives; r++)
807
{
808
delete [] transform[r];
809
}// end loop over 'r'
810
delete [] transform;
811
break;
812
}
813
case 2:
814
{
815
816
// Array of basisvalues
817
double basisvalues[3] = {0.0, 0.0, 0.0};
818
819
// Declare helper variables
820
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
821
822
// Compute basisvalues
823
basisvalues[0] = 1.0;
824
basisvalues[1] = tmp0;
825
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
826
basisvalues[0] *= std::sqrt(0.5);
827
basisvalues[2] *= std::sqrt(1.0);
828
basisvalues[1] *= std::sqrt(3.0);
829
830
// Table(s) of coefficients
831
static const double coefficients0[3] = \
832
{0.47140452, -0.57735027, -0.66666667};
833
834
static const double coefficients1[3] = \
835
{0.47140452, 0.0, 0.33333333};
836
837
// Tables of derivatives of the polynomial base (transpose).
838
static const double dmats0[3][3] = \
839
{{0.0, 0.0, 0.0},
840
{4.8989795, 0.0, 0.0},
841
{0.0, 0.0, 0.0}};
842
843
static const double dmats1[3][3] = \
844
{{0.0, 0.0, 0.0},
845
{2.4494897, 0.0, 0.0},
846
{4.2426407, 0.0, 0.0}};
847
848
// Compute reference derivatives.
849
// Declare pointer to array of derivatives on FIAT element.
850
double *derivatives = new double[2*num_derivatives];
851
for (unsigned int r = 0; r < 2*num_derivatives; r++)
852
{
853
derivatives[r] = 0.0;
854
}// end loop over 'r'
855
856
// Declare pointer to array of reference derivatives on physical element.
857
double *derivatives_p = new double[2*num_derivatives];
858
for (unsigned int r = 0; r < 2*num_derivatives; r++)
859
{
860
derivatives_p[r] = 0.0;
861
}// end loop over 'r'
862
863
// Declare derivative matrix (of polynomial basis).
864
double dmats[3][3] = \
865
{{1.0, 0.0, 0.0},
866
{0.0, 1.0, 0.0},
867
{0.0, 0.0, 1.0}};
868
869
// Declare (auxiliary) derivative matrix (of polynomial basis).
870
double dmats_old[3][3] = \
871
{{1.0, 0.0, 0.0},
872
{0.0, 1.0, 0.0},
873
{0.0, 0.0, 1.0}};
874
875
// Loop possible derivatives.
876
for (unsigned int r = 0; r < num_derivatives; r++)
877
{
878
// Resetting dmats values to compute next derivative.
879
for (unsigned int t = 0; t < 3; t++)
880
{
881
for (unsigned int u = 0; u < 3; u++)
882
{
883
dmats[t][u] = 0.0;
884
if (t == u)
885
{
886
dmats[t][u] = 1.0;
887
}
888
889
}// end loop over 'u'
890
}// end loop over 't'
891
892
// Looping derivative order to generate dmats.
893
for (unsigned int s = 0; s < n; s++)
894
{
895
// Updating dmats_old with new values and resetting dmats.
896
for (unsigned int t = 0; t < 3; t++)
897
{
898
for (unsigned int u = 0; u < 3; u++)
899
{
900
dmats_old[t][u] = dmats[t][u];
901
dmats[t][u] = 0.0;
902
}// end loop over 'u'
903
}// end loop over 't'
904
905
// Update dmats using an inner product.
906
if (combinations[r][s] == 0)
907
{
908
for (unsigned int t = 0; t < 3; t++)
909
{
910
for (unsigned int u = 0; u < 3; u++)
911
{
912
for (unsigned int tu = 0; tu < 3; tu++)
913
{
914
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
915
}// end loop over 'tu'
916
}// end loop over 'u'
917
}// end loop over 't'
918
}
919
920
if (combinations[r][s] == 1)
921
{
922
for (unsigned int t = 0; t < 3; t++)
923
{
924
for (unsigned int u = 0; u < 3; u++)
925
{
926
for (unsigned int tu = 0; tu < 3; tu++)
927
{
928
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
929
}// end loop over 'tu'
930
}// end loop over 'u'
931
}// end loop over 't'
932
}
933
934
}// end loop over 's'
935
for (unsigned int s = 0; s < 3; s++)
936
{
937
for (unsigned int t = 0; t < 3; t++)
938
{
939
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
940
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
941
}// end loop over 't'
942
}// end loop over 's'
943
944
// Using contravariant Piola transform to map values back to the physical element.
945
const double tmp_ref0 = derivatives[r];
946
const double tmp_ref1 = derivatives[num_derivatives + r];
947
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
948
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
949
}// end loop over 'r'
950
951
// Transform derivatives back to physical element
952
for (unsigned int r = 0; r < num_derivatives; r++)
953
{
954
for (unsigned int s = 0; s < num_derivatives; s++)
955
{
956
values[r] += transform[r][s]*derivatives_p[s];
957
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
958
}// end loop over 's'
959
}// end loop over 'r'
960
961
// Delete pointer to array of derivatives on FIAT element
962
delete [] derivatives;
963
964
// Delete pointer to array of reference derivatives on physical element.
965
delete [] derivatives_p;
966
967
// Delete pointer to array of combinations of derivatives and transform
968
for (unsigned int r = 0; r < num_derivatives; r++)
969
{
970
delete [] combinations[r];
971
}// end loop over 'r'
972
delete [] combinations;
973
for (unsigned int r = 0; r < num_derivatives; r++)
974
{
975
delete [] transform[r];
976
}// end loop over 'r'
977
delete [] transform;
978
break;
979
}
980
case 3:
981
{
982
983
// Array of basisvalues
984
double basisvalues[3] = {0.0, 0.0, 0.0};
985
986
// Declare helper variables
987
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
988
989
// Compute basisvalues
990
basisvalues[0] = 1.0;
991
basisvalues[1] = tmp0;
992
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
993
basisvalues[0] *= std::sqrt(0.5);
994
basisvalues[2] *= std::sqrt(1.0);
995
basisvalues[1] *= std::sqrt(3.0);
996
997
// Table(s) of coefficients
998
static const double coefficients0[3] = \
999
{0.47140452, 0.28867513, 0.83333333};
1000
1001
static const double coefficients1[3] = \
1002
{-0.94280904, 0.0, -0.66666667};
1003
1004
// Tables of derivatives of the polynomial base (transpose).
1005
static const double dmats0[3][3] = \
1006
{{0.0, 0.0, 0.0},
1007
{4.8989795, 0.0, 0.0},
1008
{0.0, 0.0, 0.0}};
1009
1010
static const double dmats1[3][3] = \
1011
{{0.0, 0.0, 0.0},
1012
{2.4494897, 0.0, 0.0},
1013
{4.2426407, 0.0, 0.0}};
1014
1015
// Compute reference derivatives.
1016
// Declare pointer to array of derivatives on FIAT element.
1017
double *derivatives = new double[2*num_derivatives];
1018
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1019
{
1020
derivatives[r] = 0.0;
1021
}// end loop over 'r'
1022
1023
// Declare pointer to array of reference derivatives on physical element.
1024
double *derivatives_p = new double[2*num_derivatives];
1025
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1026
{
1027
derivatives_p[r] = 0.0;
1028
}// end loop over 'r'
1029
1030
// Declare derivative matrix (of polynomial basis).
1031
double dmats[3][3] = \
1032
{{1.0, 0.0, 0.0},
1033
{0.0, 1.0, 0.0},
1034
{0.0, 0.0, 1.0}};
1035
1036
// Declare (auxiliary) derivative matrix (of polynomial basis).
1037
double dmats_old[3][3] = \
1038
{{1.0, 0.0, 0.0},
1039
{0.0, 1.0, 0.0},
1040
{0.0, 0.0, 1.0}};
1041
1042
// Loop possible derivatives.
1043
for (unsigned int r = 0; r < num_derivatives; r++)
1044
{
1045
// Resetting dmats values to compute next derivative.
1046
for (unsigned int t = 0; t < 3; t++)
1047
{
1048
for (unsigned int u = 0; u < 3; u++)
1049
{
1050
dmats[t][u] = 0.0;
1051
if (t == u)
1052
{
1053
dmats[t][u] = 1.0;
1054
}
1055
1056
}// end loop over 'u'
1057
}// end loop over 't'
1058
1059
// Looping derivative order to generate dmats.
1060
for (unsigned int s = 0; s < n; s++)
1061
{
1062
// Updating dmats_old with new values and resetting dmats.
1063
for (unsigned int t = 0; t < 3; t++)
1064
{
1065
for (unsigned int u = 0; u < 3; u++)
1066
{
1067
dmats_old[t][u] = dmats[t][u];
1068
dmats[t][u] = 0.0;
1069
}// end loop over 'u'
1070
}// end loop over 't'
1071
1072
// Update dmats using an inner product.
1073
if (combinations[r][s] == 0)
1074
{
1075
for (unsigned int t = 0; t < 3; t++)
1076
{
1077
for (unsigned int u = 0; u < 3; u++)
1078
{
1079
for (unsigned int tu = 0; tu < 3; tu++)
1080
{
1081
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1082
}// end loop over 'tu'
1083
}// end loop over 'u'
1084
}// end loop over 't'
1085
}
1086
1087
if (combinations[r][s] == 1)
1088
{
1089
for (unsigned int t = 0; t < 3; t++)
1090
{
1091
for (unsigned int u = 0; u < 3; u++)
1092
{
1093
for (unsigned int tu = 0; tu < 3; tu++)
1094
{
1095
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1096
}// end loop over 'tu'
1097
}// end loop over 'u'
1098
}// end loop over 't'
1099
}
1100
1101
}// end loop over 's'
1102
for (unsigned int s = 0; s < 3; s++)
1103
{
1104
for (unsigned int t = 0; t < 3; t++)
1105
{
1106
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1107
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
1108
}// end loop over 't'
1109
}// end loop over 's'
1110
1111
// Using contravariant Piola transform to map values back to the physical element.
1112
const double tmp_ref0 = derivatives[r];
1113
const double tmp_ref1 = derivatives[num_derivatives + r];
1114
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
1115
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
1116
}// end loop over 'r'
1117
1118
// Transform derivatives back to physical element
1119
for (unsigned int r = 0; r < num_derivatives; r++)
1120
{
1121
for (unsigned int s = 0; s < num_derivatives; s++)
1122
{
1123
values[r] += transform[r][s]*derivatives_p[s];
1124
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
1125
}// end loop over 's'
1126
}// end loop over 'r'
1127
1128
// Delete pointer to array of derivatives on FIAT element
1129
delete [] derivatives;
1130
1131
// Delete pointer to array of reference derivatives on physical element.
1132
delete [] derivatives_p;
1133
1134
// Delete pointer to array of combinations of derivatives and transform
1135
for (unsigned int r = 0; r < num_derivatives; r++)
1136
{
1137
delete [] combinations[r];
1138
}// end loop over 'r'
1139
delete [] combinations;
1140
for (unsigned int r = 0; r < num_derivatives; r++)
1141
{
1142
delete [] transform[r];
1143
}// end loop over 'r'
1144
delete [] transform;
1145
break;
1146
}
1147
case 4:
1148
{
1149
1150
// Array of basisvalues
1151
double basisvalues[3] = {0.0, 0.0, 0.0};
1152
1153
// Declare helper variables
1154
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1155
1156
// Compute basisvalues
1157
basisvalues[0] = 1.0;
1158
basisvalues[1] = tmp0;
1159
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1160
basisvalues[0] *= std::sqrt(0.5);
1161
basisvalues[2] *= std::sqrt(1.0);
1162
basisvalues[1] *= std::sqrt(3.0);
1163
1164
// Table(s) of coefficients
1165
static const double coefficients0[3] = \
1166
{-0.47140452, -0.28867513, 0.16666667};
1167
1168
static const double coefficients1[3] = \
1169
{-0.47140452, 0.8660254, 0.16666667};
1170
1171
// Tables of derivatives of the polynomial base (transpose).
1172
static const double dmats0[3][3] = \
1173
{{0.0, 0.0, 0.0},
1174
{4.8989795, 0.0, 0.0},
1175
{0.0, 0.0, 0.0}};
1176
1177
static const double dmats1[3][3] = \
1178
{{0.0, 0.0, 0.0},
1179
{2.4494897, 0.0, 0.0},
1180
{4.2426407, 0.0, 0.0}};
1181
1182
// Compute reference derivatives.
1183
// Declare pointer to array of derivatives on FIAT element.
1184
double *derivatives = new double[2*num_derivatives];
1185
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1186
{
1187
derivatives[r] = 0.0;
1188
}// end loop over 'r'
1189
1190
// Declare pointer to array of reference derivatives on physical element.
1191
double *derivatives_p = new double[2*num_derivatives];
1192
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1193
{
1194
derivatives_p[r] = 0.0;
1195
}// end loop over 'r'
1196
1197
// Declare derivative matrix (of polynomial basis).
1198
double dmats[3][3] = \
1199
{{1.0, 0.0, 0.0},
1200
{0.0, 1.0, 0.0},
1201
{0.0, 0.0, 1.0}};
1202
1203
// Declare (auxiliary) derivative matrix (of polynomial basis).
1204
double dmats_old[3][3] = \
1205
{{1.0, 0.0, 0.0},
1206
{0.0, 1.0, 0.0},
1207
{0.0, 0.0, 1.0}};
1208
1209
// Loop possible derivatives.
1210
for (unsigned int r = 0; r < num_derivatives; r++)
1211
{
1212
// Resetting dmats values to compute next derivative.
1213
for (unsigned int t = 0; t < 3; t++)
1214
{
1215
for (unsigned int u = 0; u < 3; u++)
1216
{
1217
dmats[t][u] = 0.0;
1218
if (t == u)
1219
{
1220
dmats[t][u] = 1.0;
1221
}
1222
1223
}// end loop over 'u'
1224
}// end loop over 't'
1225
1226
// Looping derivative order to generate dmats.
1227
for (unsigned int s = 0; s < n; s++)
1228
{
1229
// Updating dmats_old with new values and resetting dmats.
1230
for (unsigned int t = 0; t < 3; t++)
1231
{
1232
for (unsigned int u = 0; u < 3; u++)
1233
{
1234
dmats_old[t][u] = dmats[t][u];
1235
dmats[t][u] = 0.0;
1236
}// end loop over 'u'
1237
}// end loop over 't'
1238
1239
// Update dmats using an inner product.
1240
if (combinations[r][s] == 0)
1241
{
1242
for (unsigned int t = 0; t < 3; t++)
1243
{
1244
for (unsigned int u = 0; u < 3; u++)
1245
{
1246
for (unsigned int tu = 0; tu < 3; tu++)
1247
{
1248
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1249
}// end loop over 'tu'
1250
}// end loop over 'u'
1251
}// end loop over 't'
1252
}
1253
1254
if (combinations[r][s] == 1)
1255
{
1256
for (unsigned int t = 0; t < 3; t++)
1257
{
1258
for (unsigned int u = 0; u < 3; u++)
1259
{
1260
for (unsigned int tu = 0; tu < 3; tu++)
1261
{
1262
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1263
}// end loop over 'tu'
1264
}// end loop over 'u'
1265
}// end loop over 't'
1266
}
1267
1268
}// end loop over 's'
1269
for (unsigned int s = 0; s < 3; s++)
1270
{
1271
for (unsigned int t = 0; t < 3; t++)
1272
{
1273
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1274
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
1275
}// end loop over 't'
1276
}// end loop over 's'
1277
1278
// Using contravariant Piola transform to map values back to the physical element.
1279
const double tmp_ref0 = derivatives[r];
1280
const double tmp_ref1 = derivatives[num_derivatives + r];
1281
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
1282
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
1283
}// end loop over 'r'
1284
1285
// Transform derivatives back to physical element
1286
for (unsigned int r = 0; r < num_derivatives; r++)
1287
{
1288
for (unsigned int s = 0; s < num_derivatives; s++)
1289
{
1290
values[r] += transform[r][s]*derivatives_p[s];
1291
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
1292
}// end loop over 's'
1293
}// end loop over 'r'
1294
1295
// Delete pointer to array of derivatives on FIAT element
1296
delete [] derivatives;
1297
1298
// Delete pointer to array of reference derivatives on physical element.
1299
delete [] derivatives_p;
1300
1301
// Delete pointer to array of combinations of derivatives and transform
1302
for (unsigned int r = 0; r < num_derivatives; r++)
1303
{
1304
delete [] combinations[r];
1305
}// end loop over 'r'
1306
delete [] combinations;
1307
for (unsigned int r = 0; r < num_derivatives; r++)
1308
{
1309
delete [] transform[r];
1310
}// end loop over 'r'
1311
delete [] transform;
1312
break;
1313
}
1314
case 5:
1315
{
1316
1317
// Array of basisvalues
1318
double basisvalues[3] = {0.0, 0.0, 0.0};
1319
1320
// Declare helper variables
1321
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1322
1323
// Compute basisvalues
1324
basisvalues[0] = 1.0;
1325
basisvalues[1] = tmp0;
1326
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1327
basisvalues[0] *= std::sqrt(0.5);
1328
basisvalues[2] *= std::sqrt(1.0);
1329
basisvalues[1] *= std::sqrt(3.0);
1330
1331
// Table(s) of coefficients
1332
static const double coefficients0[3] = \
1333
{0.94280904, 0.57735027, -0.33333333};
1334
1335
static const double coefficients1[3] = \
1336
{-0.47140452, -0.8660254, 0.16666667};
1337
1338
// Tables of derivatives of the polynomial base (transpose).
1339
static const double dmats0[3][3] = \
1340
{{0.0, 0.0, 0.0},
1341
{4.8989795, 0.0, 0.0},
1342
{0.0, 0.0, 0.0}};
1343
1344
static const double dmats1[3][3] = \
1345
{{0.0, 0.0, 0.0},
1346
{2.4494897, 0.0, 0.0},
1347
{4.2426407, 0.0, 0.0}};
1348
1349
// Compute reference derivatives.
1350
// Declare pointer to array of derivatives on FIAT element.
1351
double *derivatives = new double[2*num_derivatives];
1352
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1353
{
1354
derivatives[r] = 0.0;
1355
}// end loop over 'r'
1356
1357
// Declare pointer to array of reference derivatives on physical element.
1358
double *derivatives_p = new double[2*num_derivatives];
1359
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1360
{
1361
derivatives_p[r] = 0.0;
1362
}// end loop over 'r'
1363
1364
// Declare derivative matrix (of polynomial basis).
1365
double dmats[3][3] = \
1366
{{1.0, 0.0, 0.0},
1367
{0.0, 1.0, 0.0},
1368
{0.0, 0.0, 1.0}};
1369
1370
// Declare (auxiliary) derivative matrix (of polynomial basis).
1371
double dmats_old[3][3] = \
1372
{{1.0, 0.0, 0.0},
1373
{0.0, 1.0, 0.0},
1374
{0.0, 0.0, 1.0}};
1375
1376
// Loop possible derivatives.
1377
for (unsigned int r = 0; r < num_derivatives; r++)
1378
{
1379
// Resetting dmats values to compute next derivative.
1380
for (unsigned int t = 0; t < 3; t++)
1381
{
1382
for (unsigned int u = 0; u < 3; u++)
1383
{
1384
dmats[t][u] = 0.0;
1385
if (t == u)
1386
{
1387
dmats[t][u] = 1.0;
1388
}
1389
1390
}// end loop over 'u'
1391
}// end loop over 't'
1392
1393
// Looping derivative order to generate dmats.
1394
for (unsigned int s = 0; s < n; s++)
1395
{
1396
// Updating dmats_old with new values and resetting dmats.
1397
for (unsigned int t = 0; t < 3; t++)
1398
{
1399
for (unsigned int u = 0; u < 3; u++)
1400
{
1401
dmats_old[t][u] = dmats[t][u];
1402
dmats[t][u] = 0.0;
1403
}// end loop over 'u'
1404
}// end loop over 't'
1405
1406
// Update dmats using an inner product.
1407
if (combinations[r][s] == 0)
1408
{
1409
for (unsigned int t = 0; t < 3; t++)
1410
{
1411
for (unsigned int u = 0; u < 3; u++)
1412
{
1413
for (unsigned int tu = 0; tu < 3; tu++)
1414
{
1415
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1416
}// end loop over 'tu'
1417
}// end loop over 'u'
1418
}// end loop over 't'
1419
}
1420
1421
if (combinations[r][s] == 1)
1422
{
1423
for (unsigned int t = 0; t < 3; t++)
1424
{
1425
for (unsigned int u = 0; u < 3; u++)
1426
{
1427
for (unsigned int tu = 0; tu < 3; tu++)
1428
{
1429
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1430
}// end loop over 'tu'
1431
}// end loop over 'u'
1432
}// end loop over 't'
1433
}
1434
1435
}// end loop over 's'
1436
for (unsigned int s = 0; s < 3; s++)
1437
{
1438
for (unsigned int t = 0; t < 3; t++)
1439
{
1440
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1441
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
1442
}// end loop over 't'
1443
}// end loop over 's'
1444
1445
// Using contravariant Piola transform to map values back to the physical element.
1446
const double tmp_ref0 = derivatives[r];
1447
const double tmp_ref1 = derivatives[num_derivatives + r];
1448
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
1449
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
1450
}// end loop over 'r'
1451
1452
// Transform derivatives back to physical element
1453
for (unsigned int r = 0; r < num_derivatives; r++)
1454
{
1455
for (unsigned int s = 0; s < num_derivatives; s++)
1456
{
1457
values[r] += transform[r][s]*derivatives_p[s];
1458
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
1459
}// end loop over 's'
1460
}// end loop over 'r'
1461
1462
// Delete pointer to array of derivatives on FIAT element
1463
delete [] derivatives;
1464
1465
// Delete pointer to array of reference derivatives on physical element.
1466
delete [] derivatives_p;
1467
1468
// Delete pointer to array of combinations of derivatives and transform
1469
for (unsigned int r = 0; r < num_derivatives; r++)
1470
{
1471
delete [] combinations[r];
1472
}// end loop over 'r'
1473
delete [] combinations;
1474
for (unsigned int r = 0; r < num_derivatives; r++)
1475
{
1476
delete [] transform[r];
1477
}// end loop over 'r'
1478
delete [] transform;
1479
break;
1480
}
1481
}
1482
1483
}
1484
1485
/// Evaluate order n derivatives of all basis functions at given point x in cell
1486
virtual void evaluate_basis_derivatives_all(std::size_t n,
1487
double* values,
1488
const double* x,
1489
const double* vertex_coordinates,
1490
int cell_orientation) const
1491
{
1492
// Compute number of derivatives.
1493
unsigned int num_derivatives = 1;
1494
for (unsigned int r = 0; r < n; r++)
1495
{
1496
num_derivatives *= 2;
1497
}// end loop over 'r'
1498
1499
// Helper variable to hold values of a single dof.
1500
double *dof_values = new double[2*num_derivatives];
1501
for (unsigned int r = 0; r < 2*num_derivatives; r++)
1502
{
1503
dof_values[r] = 0.0;
1504
}// end loop over 'r'
1505
1506
// Loop dofs and call evaluate_basis_derivatives.
1507
for (unsigned int r = 0; r < 6; r++)
1508
{
1509
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1510
for (unsigned int s = 0; s < 2*num_derivatives; s++)
1511
{
1512
values[r*2*num_derivatives + s] = dof_values[s];
1513
}// end loop over 's'
1514
}// end loop over 'r'
1515
1516
// Delete pointer.
1517
delete [] dof_values;
1518
}
1519
1520
/// Evaluate linear functional for dof i on the function f
1521
virtual double evaluate_dof(std::size_t i,
1522
const ufc::function& f,
1523
const double* vertex_coordinates,
1524
int cell_orientation,
1525
const ufc::cell& c) const
1526
{
1527
// Declare variables for result of evaluation
1528
double vals[2];
1529
1530
// Declare variable for physical coordinates
1531
double y[2];
1532
1533
double result;
1534
// Compute Jacobian
1535
double J[4];
1536
compute_jacobian_triangle_2d(J, vertex_coordinates);
1537
1538
1539
// Compute Jacobian inverse and determinant
1540
double K[4];
1541
double detJ;
1542
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1543
1544
1545
switch (i)
1546
{
1547
case 0:
1548
{
1549
y[0] = 0.66666667*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
1550
y[1] = 0.66666667*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
1551
f.evaluate(vals, y, c);
1552
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
1553
return result;
1554
break;
1555
}
1556
case 1:
1557
{
1558
y[0] = 0.33333333*vertex_coordinates[2] + 0.66666667*vertex_coordinates[4];
1559
y[1] = 0.33333333*vertex_coordinates[3] + 0.66666667*vertex_coordinates[5];
1560
f.evaluate(vals, y, c);
1561
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
1562
return result;
1563
break;
1564
}
1565
case 2:
1566
{
1567
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[4];
1568
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[5];
1569
f.evaluate(vals, y, c);
1570
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
1571
return result;
1572
break;
1573
}
1574
case 3:
1575
{
1576
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[4];
1577
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[5];
1578
f.evaluate(vals, y, c);
1579
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
1580
return result;
1581
break;
1582
}
1583
case 4:
1584
{
1585
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2];
1586
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3];
1587
f.evaluate(vals, y, c);
1588
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
1589
return result;
1590
break;
1591
}
1592
case 5:
1593
{
1594
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[2];
1595
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[3];
1596
f.evaluate(vals, y, c);
1597
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
1598
return result;
1599
break;
1600
}
1601
}
1602
1603
return 0.0;
1604
}
1605
1606
/// Evaluate linear functionals for all dofs on the function f
1607
virtual void evaluate_dofs(double* values,
1608
const ufc::function& f,
1609
const double* vertex_coordinates,
1610
int cell_orientation,
1611
const ufc::cell& c) const
1612
{
1613
// Declare variables for result of evaluation
1614
double vals[2];
1615
1616
// Declare variable for physical coordinates
1617
double y[2];
1618
1619
double result;
1620
// Compute Jacobian
1621
double J[4];
1622
compute_jacobian_triangle_2d(J, vertex_coordinates);
1623
1624
1625
// Compute Jacobian inverse and determinant
1626
double K[4];
1627
double detJ;
1628
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1629
1630
1631
y[0] = 0.66666667*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
1632
y[1] = 0.66666667*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
1633
f.evaluate(vals, y, c);
1634
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
1635
values[0] = result;
1636
y[0] = 0.33333333*vertex_coordinates[2] + 0.66666667*vertex_coordinates[4];
1637
y[1] = 0.33333333*vertex_coordinates[3] + 0.66666667*vertex_coordinates[5];
1638
f.evaluate(vals, y, c);
1639
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
1640
values[1] = result;
1641
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[4];
1642
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[5];
1643
f.evaluate(vals, y, c);
1644
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
1645
values[2] = result;
1646
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[4];
1647
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[5];
1648
f.evaluate(vals, y, c);
1649
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
1650
values[3] = result;
1651
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2];
1652
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3];
1653
f.evaluate(vals, y, c);
1654
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
1655
values[4] = result;
1656
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[2];
1657
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[3];
1658
f.evaluate(vals, y, c);
1659
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
1660
values[5] = result;
1661
}
1662
1663
/// Interpolate vertex values from dof values
1664
virtual void interpolate_vertex_values(double* vertex_values,
1665
const double* dof_values,
1666
const double* vertex_coordinates,
1667
int cell_orientation,
1668
const ufc::cell& c) const
1669
{
1670
// Compute Jacobian
1671
double J[4];
1672
compute_jacobian_triangle_2d(J, vertex_coordinates);
1673
1674
1675
// Compute Jacobian inverse and determinant
1676
double K[4];
1677
double detJ;
1678
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1679
1680
1681
1682
// Evaluate function and change variables
1683
vertex_values[0] = dof_values[2]*(1.0/detJ)*J[0]*2.0 + dof_values[3]*((1.0/detJ)*(J[0]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[1]*(-2.0))) + dof_values[5]*(1.0/detJ)*J[1];
1684
vertex_values[2] = dof_values[0]*(1.0/detJ)*J[0]*2.0 + dof_values[1]*((1.0/detJ)*(J[0]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[0]*(-1.0) + J[1])) + dof_values[5]*((1.0/detJ)*(J[0]*2.0 + J[1]*(-2.0)));
1685
vertex_values[4] = dof_values[0]*((1.0/detJ)*(J[1]*(-1.0))) + dof_values[1]*(1.0/detJ)*J[1]*2.0 + dof_values[2]*((1.0/detJ)*(J[0]*(-1.0) + J[1])) + dof_values[3]*((1.0/detJ)*(J[0]*2.0 + J[1]*(-2.0)));
1686
vertex_values[1] = dof_values[2]*(1.0/detJ)*J[2]*2.0 + dof_values[3]*((1.0/detJ)*(J[2]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[3]*(-2.0))) + dof_values[5]*(1.0/detJ)*J[3];
1687
vertex_values[3] = dof_values[0]*(1.0/detJ)*J[2]*2.0 + dof_values[1]*((1.0/detJ)*(J[2]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[2]*(-1.0) + J[3])) + dof_values[5]*((1.0/detJ)*(J[2]*2.0 + J[3]*(-2.0)));
1688
vertex_values[5] = dof_values[0]*((1.0/detJ)*(J[3]*(-1.0))) + dof_values[1]*(1.0/detJ)*J[3]*2.0 + dof_values[2]*((1.0/detJ)*(J[2]*(-1.0) + J[3])) + dof_values[3]*((1.0/detJ)*(J[2]*2.0 + J[3]*(-2.0)));
1689
}
1690
1691
/// Map coordinate xhat from reference cell to coordinate x in cell
1692
virtual void map_from_reference_cell(double* x,
1693
const double* xhat,
1694
const ufc::cell& c) const
1695
{
1696
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1697
}
1698
1699
/// Map from coordinate x in cell to coordinate xhat in reference cell
1700
virtual void map_to_reference_cell(double* xhat,
1701
const double* x,
1702
const ufc::cell& c) const
1703
{
1704
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1705
}
1706
1707
/// Return the number of sub elements (for a mixed element)
1708
virtual std::size_t num_sub_elements() const
1709
{
1710
return 0;
1711
}
1712
1713
/// Create a new finite element for sub element i (for a mixed element)
1714
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1715
{
1716
return 0;
1717
}
1718
1719
/// Create a new class instance
1720
virtual ufc::finite_element* create() const
1721
{
1722
return new mixedpoisson_finite_element_0();
1723
}
1724
1725
};
1726
1727
/// This class defines the interface for a finite element.
1728
1729
class mixedpoisson_finite_element_1: public ufc::finite_element
1730
{
1731
public:
1732
1733
/// Constructor
1734
mixedpoisson_finite_element_1() : ufc::finite_element()
1735
{
1736
// Do nothing
1737
}
1738
1739
/// Destructor
1740
virtual ~mixedpoisson_finite_element_1()
1741
{
1742
// Do nothing
1743
}
1744
1745
/// Return a string identifying the finite element
1746
virtual const char* signature() const
1747
{
1748
return "FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
1749
}
1750
1751
/// Return the cell shape
1752
virtual ufc::shape cell_shape() const
1753
{
1754
return ufc::triangle;
1755
}
1756
1757
/// Return the topological dimension of the cell shape
1758
virtual std::size_t topological_dimension() const
1759
{
1760
return 2;
1761
}
1762
1763
/// Return the geometric dimension of the cell shape
1764
virtual std::size_t geometric_dimension() const
1765
{
1766
return 2;
1767
}
1768
1769
/// Return the dimension of the finite element function space
1770
virtual std::size_t space_dimension() const
1771
{
1772
return 1;
1773
}
1774
1775
/// Return the rank of the value space
1776
virtual std::size_t value_rank() const
1777
{
1778
return 0;
1779
}
1780
1781
/// Return the dimension of the value space for axis i
1782
virtual std::size_t value_dimension(std::size_t i) const
1783
{
1784
return 1;
1785
}
1786
1787
/// Evaluate basis function i at given point x in cell
1788
virtual void evaluate_basis(std::size_t i,
1789
double* values,
1790
const double* x,
1791
const double* vertex_coordinates,
1792
int cell_orientation) const
1793
{
1794
// Compute Jacobian
1795
double J[4];
1796
compute_jacobian_triangle_2d(J, vertex_coordinates);
1797
1798
// Compute Jacobian inverse and determinant
1799
double K[4];
1800
double detJ;
1801
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1802
1803
1804
// Compute constants
1805
1806
// Get coordinates and map to the reference (FIAT) element
1807
1808
// Reset values
1809
*values = 0.0;
1810
1811
// Array of basisvalues
1812
double basisvalues[1] = {0.0};
1813
1814
// Declare helper variables
1815
1816
// Compute basisvalues
1817
basisvalues[0] = 1.0;
1818
1819
// Table(s) of coefficients
1820
static const double coefficients0[1] = \
1821
{1.0};
1822
1823
// Compute value(s)
1824
for (unsigned int r = 0; r < 1; r++)
1825
{
1826
*values += coefficients0[r]*basisvalues[r];
1827
}// end loop over 'r'
1828
}
1829
1830
/// Evaluate all basis functions at given point x in cell
1831
virtual void evaluate_basis_all(double* values,
1832
const double* x,
1833
const double* vertex_coordinates,
1834
int cell_orientation) const
1835
{
1836
// Element is constant, calling evaluate_basis.
1837
evaluate_basis(0, values, x, vertex_coordinates, cell_orientation);
1838
}
1839
1840
/// Evaluate order n derivatives of basis function i at given point x in cell
1841
virtual void evaluate_basis_derivatives(std::size_t i,
1842
std::size_t n,
1843
double* values,
1844
const double* x,
1845
const double* vertex_coordinates,
1846
int cell_orientation) const
1847
{
1848
// Compute Jacobian
1849
double J[4];
1850
compute_jacobian_triangle_2d(J, vertex_coordinates);
1851
1852
// Compute Jacobian inverse and determinant
1853
double K[4];
1854
double detJ;
1855
compute_jacobian_inverse_triangle_2d(K, detJ, J);
1856
1857
1858
// Compute constants
1859
1860
// Get coordinates and map to the reference (FIAT) element
1861
1862
// Compute number of derivatives.
1863
unsigned int num_derivatives = 1;
1864
for (unsigned int r = 0; r < n; r++)
1865
{
1866
num_derivatives *= 2;
1867
}// end loop over 'r'
1868
1869
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1870
unsigned int **combinations = new unsigned int *[num_derivatives];
1871
for (unsigned int row = 0; row < num_derivatives; row++)
1872
{
1873
combinations[row] = new unsigned int [n];
1874
for (unsigned int col = 0; col < n; col++)
1875
combinations[row][col] = 0;
1876
}
1877
1878
// Generate combinations of derivatives
1879
for (unsigned int row = 1; row < num_derivatives; row++)
1880
{
1881
for (unsigned int num = 0; num < row; num++)
1882
{
1883
for (unsigned int col = n-1; col+1 > 0; col--)
1884
{
1885
if (combinations[row][col] + 1 > 1)
1886
combinations[row][col] = 0;
1887
else
1888
{
1889
combinations[row][col] += 1;
1890
break;
1891
}
1892
}
1893
}
1894
}
1895
1896
// Compute inverse of Jacobian
1897
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
1898
1899
// Declare transformation matrix
1900
// Declare pointer to two dimensional array and initialise
1901
double **transform = new double *[num_derivatives];
1902
1903
for (unsigned int j = 0; j < num_derivatives; j++)
1904
{
1905
transform[j] = new double [num_derivatives];
1906
for (unsigned int k = 0; k < num_derivatives; k++)
1907
transform[j][k] = 1;
1908
}
1909
1910
// Construct transformation matrix
1911
for (unsigned int row = 0; row < num_derivatives; row++)
1912
{
1913
for (unsigned int col = 0; col < num_derivatives; col++)
1914
{
1915
for (unsigned int k = 0; k < n; k++)
1916
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1917
}
1918
}
1919
1920
// Reset values. Assuming that values is always an array.
1921
for (unsigned int r = 0; r < num_derivatives; r++)
1922
{
1923
values[r] = 0.0;
1924
}// end loop over 'r'
1925
1926
1927
// Array of basisvalues
1928
double basisvalues[1] = {0.0};
1929
1930
// Declare helper variables
1931
1932
// Compute basisvalues
1933
basisvalues[0] = 1.0;
1934
1935
// Table(s) of coefficients
1936
static const double coefficients0[1] = \
1937
{1.0};
1938
1939
// Tables of derivatives of the polynomial base (transpose).
1940
static const double dmats0[1][1] = \
1941
{{0.0}};
1942
1943
static const double dmats1[1][1] = \
1944
{{0.0}};
1945
1946
// Compute reference derivatives.
1947
// Declare pointer to array of derivatives on FIAT element.
1948
double *derivatives = new double[num_derivatives];
1949
for (unsigned int r = 0; r < num_derivatives; r++)
1950
{
1951
derivatives[r] = 0.0;
1952
}// end loop over 'r'
1953
1954
// Declare derivative matrix (of polynomial basis).
1955
double dmats[1][1] = \
1956
{{1.0}};
1957
1958
// Declare (auxiliary) derivative matrix (of polynomial basis).
1959
double dmats_old[1][1] = \
1960
{{1.0}};
1961
1962
// Loop possible derivatives.
1963
for (unsigned int r = 0; r < num_derivatives; r++)
1964
{
1965
// Resetting dmats values to compute next derivative.
1966
for (unsigned int t = 0; t < 1; t++)
1967
{
1968
for (unsigned int u = 0; u < 1; u++)
1969
{
1970
dmats[t][u] = 0.0;
1971
if (t == u)
1972
{
1973
dmats[t][u] = 1.0;
1974
}
1975
1976
}// end loop over 'u'
1977
}// end loop over 't'
1978
1979
// Looping derivative order to generate dmats.
1980
for (unsigned int s = 0; s < n; s++)
1981
{
1982
// Updating dmats_old with new values and resetting dmats.
1983
for (unsigned int t = 0; t < 1; t++)
1984
{
1985
for (unsigned int u = 0; u < 1; u++)
1986
{
1987
dmats_old[t][u] = dmats[t][u];
1988
dmats[t][u] = 0.0;
1989
}// end loop over 'u'
1990
}// end loop over 't'
1991
1992
// Update dmats using an inner product.
1993
if (combinations[r][s] == 0)
1994
{
1995
for (unsigned int t = 0; t < 1; t++)
1996
{
1997
for (unsigned int u = 0; u < 1; u++)
1998
{
1999
for (unsigned int tu = 0; tu < 1; tu++)
2000
{
2001
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2002
}// end loop over 'tu'
2003
}// end loop over 'u'
2004
}// end loop over 't'
2005
}
2006
2007
if (combinations[r][s] == 1)
2008
{
2009
for (unsigned int t = 0; t < 1; t++)
2010
{
2011
for (unsigned int u = 0; u < 1; u++)
2012
{
2013
for (unsigned int tu = 0; tu < 1; tu++)
2014
{
2015
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2016
}// end loop over 'tu'
2017
}// end loop over 'u'
2018
}// end loop over 't'
2019
}
2020
2021
}// end loop over 's'
2022
for (unsigned int s = 0; s < 1; s++)
2023
{
2024
for (unsigned int t = 0; t < 1; t++)
2025
{
2026
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2027
}// end loop over 't'
2028
}// end loop over 's'
2029
}// end loop over 'r'
2030
2031
// Transform derivatives back to physical element
2032
for (unsigned int r = 0; r < num_derivatives; r++)
2033
{
2034
for (unsigned int s = 0; s < num_derivatives; s++)
2035
{
2036
values[r] += transform[r][s]*derivatives[s];
2037
}// end loop over 's'
2038
}// end loop over 'r'
2039
2040
// Delete pointer to array of derivatives on FIAT element
2041
delete [] derivatives;
2042
2043
// Delete pointer to array of combinations of derivatives and transform
2044
for (unsigned int r = 0; r < num_derivatives; r++)
2045
{
2046
delete [] combinations[r];
2047
}// end loop over 'r'
2048
delete [] combinations;
2049
for (unsigned int r = 0; r < num_derivatives; r++)
2050
{
2051
delete [] transform[r];
2052
}// end loop over 'r'
2053
delete [] transform;
2054
}
2055
2056
/// Evaluate order n derivatives of all basis functions at given point x in cell
2057
virtual void evaluate_basis_derivatives_all(std::size_t n,
2058
double* values,
2059
const double* x,
2060
const double* vertex_coordinates,
2061
int cell_orientation) const
2062
{
2063
// Element is constant, calling evaluate_basis_derivatives.
2064
evaluate_basis_derivatives(0, n, values, x, vertex_coordinates, cell_orientation);
2065
}
2066
2067
/// Evaluate linear functional for dof i on the function f
2068
virtual double evaluate_dof(std::size_t i,
2069
const ufc::function& f,
2070
const double* vertex_coordinates,
2071
int cell_orientation,
2072
const ufc::cell& c) const
2073
{
2074
// Declare variables for result of evaluation
2075
double vals[1];
2076
2077
// Declare variable for physical coordinates
2078
double y[2];
2079
switch (i)
2080
{
2081
case 0:
2082
{
2083
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
2084
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
2085
f.evaluate(vals, y, c);
2086
return vals[0];
2087
break;
2088
}
2089
}
2090
2091
return 0.0;
2092
}
2093
2094
/// Evaluate linear functionals for all dofs on the function f
2095
virtual void evaluate_dofs(double* values,
2096
const ufc::function& f,
2097
const double* vertex_coordinates,
2098
int cell_orientation,
2099
const ufc::cell& c) const
2100
{
2101
// Declare variables for result of evaluation
2102
double vals[1];
2103
2104
// Declare variable for physical coordinates
2105
double y[2];
2106
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
2107
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
2108
f.evaluate(vals, y, c);
2109
values[0] = vals[0];
2110
}
2111
2112
/// Interpolate vertex values from dof values
2113
virtual void interpolate_vertex_values(double* vertex_values,
2114
const double* dof_values,
2115
const double* vertex_coordinates,
2116
int cell_orientation,
2117
const ufc::cell& c) const
2118
{
2119
// Evaluate function and change variables
2120
vertex_values[0] = dof_values[0];
2121
vertex_values[1] = dof_values[0];
2122
vertex_values[2] = dof_values[0];
2123
}
2124
2125
/// Map coordinate xhat from reference cell to coordinate x in cell
2126
virtual void map_from_reference_cell(double* x,
2127
const double* xhat,
2128
const ufc::cell& c) const
2129
{
2130
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
2131
}
2132
2133
/// Map from coordinate x in cell to coordinate xhat in reference cell
2134
virtual void map_to_reference_cell(double* xhat,
2135
const double* x,
2136
const ufc::cell& c) const
2137
{
2138
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
2139
}
2140
2141
/// Return the number of sub elements (for a mixed element)
2142
virtual std::size_t num_sub_elements() const
2143
{
2144
return 0;
2145
}
2146
2147
/// Create a new finite element for sub element i (for a mixed element)
2148
virtual ufc::finite_element* create_sub_element(std::size_t i) const
2149
{
2150
return 0;
2151
}
2152
2153
/// Create a new class instance
2154
virtual ufc::finite_element* create() const
2155
{
2156
return new mixedpoisson_finite_element_1();
2157
}
2158
2159
};
2160
2161
/// This class defines the interface for a finite element.
2162
2163
class mixedpoisson_finite_element_2: public ufc::finite_element
2164
{
2165
public:
2166
2167
/// Constructor
2168
mixedpoisson_finite_element_2() : ufc::finite_element()
2169
{
2170
// Do nothing
2171
}
2172
2173
/// Destructor
2174
virtual ~mixedpoisson_finite_element_2()
2175
{
2176
// Do nothing
2177
}
2178
2179
/// Return a string identifying the finite element
2180
virtual const char* signature() const
2181
{
2182
return "MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None), FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)], **{'value_shape': (3,) })";
2183
}
2184
2185
/// Return the cell shape
2186
virtual ufc::shape cell_shape() const
2187
{
2188
return ufc::triangle;
2189
}
2190
2191
/// Return the topological dimension of the cell shape
2192
virtual std::size_t topological_dimension() const
2193
{
2194
return 2;
2195
}
2196
2197
/// Return the geometric dimension of the cell shape
2198
virtual std::size_t geometric_dimension() const
2199
{
2200
return 2;
2201
}
2202
2203
/// Return the dimension of the finite element function space
2204
virtual std::size_t space_dimension() const
2205
{
2206
return 7;
2207
}
2208
2209
/// Return the rank of the value space
2210
virtual std::size_t value_rank() const
2211
{
2212
return 1;
2213
}
2214
2215
/// Return the dimension of the value space for axis i
2216
virtual std::size_t value_dimension(std::size_t i) const
2217
{
2218
switch (i)
2219
{
2220
case 0:
2221
{
2222
return 3;
2223
break;
2224
}
2225
}
2226
2227
return 0;
2228
}
2229
2230
/// Evaluate basis function i at given point x in cell
2231
virtual void evaluate_basis(std::size_t i,
2232
double* values,
2233
const double* x,
2234
const double* vertex_coordinates,
2235
int cell_orientation) const
2236
{
2237
// Compute Jacobian
2238
double J[4];
2239
compute_jacobian_triangle_2d(J, vertex_coordinates);
2240
2241
// Compute Jacobian inverse and determinant
2242
double K[4];
2243
double detJ;
2244
compute_jacobian_inverse_triangle_2d(K, detJ, J);
2245
2246
2247
2248
// Compute constants
2249
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
2250
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
2251
2252
// Get coordinates and map to the reference (FIAT) element
2253
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
2254
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
2255
2256
// Reset values
2257
values[0] = 0.0;
2258
values[1] = 0.0;
2259
values[2] = 0.0;
2260
switch (i)
2261
{
2262
case 0:
2263
{
2264
2265
// Array of basisvalues
2266
double basisvalues[3] = {0.0, 0.0, 0.0};
2267
2268
// Declare helper variables
2269
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2270
2271
// Compute basisvalues
2272
basisvalues[0] = 1.0;
2273
basisvalues[1] = tmp0;
2274
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2275
basisvalues[0] *= std::sqrt(0.5);
2276
basisvalues[2] *= std::sqrt(1.0);
2277
basisvalues[1] *= std::sqrt(3.0);
2278
2279
// Table(s) of coefficients
2280
static const double coefficients0[3] = \
2281
{0.94280904, 0.57735027, -0.33333333};
2282
2283
static const double coefficients1[3] = \
2284
{-0.47140452, 0.0, -0.33333333};
2285
2286
// Compute value(s)
2287
for (unsigned int r = 0; r < 3; r++)
2288
{
2289
values[0] += coefficients0[r]*basisvalues[r];
2290
values[1] += coefficients1[r]*basisvalues[r];
2291
}// end loop over 'r'
2292
2293
// Using contravariant Piola transform to map values back to the physical element
2294
const double tmp_ref0 = values[0];
2295
const double tmp_ref1 = values[1];
2296
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2297
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2298
break;
2299
}
2300
case 1:
2301
{
2302
2303
// Array of basisvalues
2304
double basisvalues[3] = {0.0, 0.0, 0.0};
2305
2306
// Declare helper variables
2307
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2308
2309
// Compute basisvalues
2310
basisvalues[0] = 1.0;
2311
basisvalues[1] = tmp0;
2312
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2313
basisvalues[0] *= std::sqrt(0.5);
2314
basisvalues[2] *= std::sqrt(1.0);
2315
basisvalues[1] *= std::sqrt(3.0);
2316
2317
// Table(s) of coefficients
2318
static const double coefficients0[3] = \
2319
{-0.47140452, -0.28867513, 0.16666667};
2320
2321
static const double coefficients1[3] = \
2322
{0.94280904, 0.0, 0.66666667};
2323
2324
// Compute value(s)
2325
for (unsigned int r = 0; r < 3; r++)
2326
{
2327
values[0] += coefficients0[r]*basisvalues[r];
2328
values[1] += coefficients1[r]*basisvalues[r];
2329
}// end loop over 'r'
2330
2331
// Using contravariant Piola transform to map values back to the physical element
2332
const double tmp_ref0 = values[0];
2333
const double tmp_ref1 = values[1];
2334
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2335
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2336
break;
2337
}
2338
case 2:
2339
{
2340
2341
// Array of basisvalues
2342
double basisvalues[3] = {0.0, 0.0, 0.0};
2343
2344
// Declare helper variables
2345
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2346
2347
// Compute basisvalues
2348
basisvalues[0] = 1.0;
2349
basisvalues[1] = tmp0;
2350
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2351
basisvalues[0] *= std::sqrt(0.5);
2352
basisvalues[2] *= std::sqrt(1.0);
2353
basisvalues[1] *= std::sqrt(3.0);
2354
2355
// Table(s) of coefficients
2356
static const double coefficients0[3] = \
2357
{0.47140452, -0.57735027, -0.66666667};
2358
2359
static const double coefficients1[3] = \
2360
{0.47140452, 0.0, 0.33333333};
2361
2362
// Compute value(s)
2363
for (unsigned int r = 0; r < 3; r++)
2364
{
2365
values[0] += coefficients0[r]*basisvalues[r];
2366
values[1] += coefficients1[r]*basisvalues[r];
2367
}// end loop over 'r'
2368
2369
// Using contravariant Piola transform to map values back to the physical element
2370
const double tmp_ref0 = values[0];
2371
const double tmp_ref1 = values[1];
2372
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2373
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2374
break;
2375
}
2376
case 3:
2377
{
2378
2379
// Array of basisvalues
2380
double basisvalues[3] = {0.0, 0.0, 0.0};
2381
2382
// Declare helper variables
2383
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2384
2385
// Compute basisvalues
2386
basisvalues[0] = 1.0;
2387
basisvalues[1] = tmp0;
2388
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2389
basisvalues[0] *= std::sqrt(0.5);
2390
basisvalues[2] *= std::sqrt(1.0);
2391
basisvalues[1] *= std::sqrt(3.0);
2392
2393
// Table(s) of coefficients
2394
static const double coefficients0[3] = \
2395
{0.47140452, 0.28867513, 0.83333333};
2396
2397
static const double coefficients1[3] = \
2398
{-0.94280904, 0.0, -0.66666667};
2399
2400
// Compute value(s)
2401
for (unsigned int r = 0; r < 3; r++)
2402
{
2403
values[0] += coefficients0[r]*basisvalues[r];
2404
values[1] += coefficients1[r]*basisvalues[r];
2405
}// end loop over 'r'
2406
2407
// Using contravariant Piola transform to map values back to the physical element
2408
const double tmp_ref0 = values[0];
2409
const double tmp_ref1 = values[1];
2410
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2411
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2412
break;
2413
}
2414
case 4:
2415
{
2416
2417
// Array of basisvalues
2418
double basisvalues[3] = {0.0, 0.0, 0.0};
2419
2420
// Declare helper variables
2421
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2422
2423
// Compute basisvalues
2424
basisvalues[0] = 1.0;
2425
basisvalues[1] = tmp0;
2426
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2427
basisvalues[0] *= std::sqrt(0.5);
2428
basisvalues[2] *= std::sqrt(1.0);
2429
basisvalues[1] *= std::sqrt(3.0);
2430
2431
// Table(s) of coefficients
2432
static const double coefficients0[3] = \
2433
{-0.47140452, -0.28867513, 0.16666667};
2434
2435
static const double coefficients1[3] = \
2436
{-0.47140452, 0.8660254, 0.16666667};
2437
2438
// Compute value(s)
2439
for (unsigned int r = 0; r < 3; r++)
2440
{
2441
values[0] += coefficients0[r]*basisvalues[r];
2442
values[1] += coefficients1[r]*basisvalues[r];
2443
}// end loop over 'r'
2444
2445
// Using contravariant Piola transform to map values back to the physical element
2446
const double tmp_ref0 = values[0];
2447
const double tmp_ref1 = values[1];
2448
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2449
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2450
break;
2451
}
2452
case 5:
2453
{
2454
2455
// Array of basisvalues
2456
double basisvalues[3] = {0.0, 0.0, 0.0};
2457
2458
// Declare helper variables
2459
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2460
2461
// Compute basisvalues
2462
basisvalues[0] = 1.0;
2463
basisvalues[1] = tmp0;
2464
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2465
basisvalues[0] *= std::sqrt(0.5);
2466
basisvalues[2] *= std::sqrt(1.0);
2467
basisvalues[1] *= std::sqrt(3.0);
2468
2469
// Table(s) of coefficients
2470
static const double coefficients0[3] = \
2471
{0.94280904, 0.57735027, -0.33333333};
2472
2473
static const double coefficients1[3] = \
2474
{-0.47140452, -0.8660254, 0.16666667};
2475
2476
// Compute value(s)
2477
for (unsigned int r = 0; r < 3; r++)
2478
{
2479
values[0] += coefficients0[r]*basisvalues[r];
2480
values[1] += coefficients1[r]*basisvalues[r];
2481
}// end loop over 'r'
2482
2483
// Using contravariant Piola transform to map values back to the physical element
2484
const double tmp_ref0 = values[0];
2485
const double tmp_ref1 = values[1];
2486
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2487
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2488
break;
2489
}
2490
case 6:
2491
{
2492
2493
// Array of basisvalues
2494
double basisvalues[1] = {0.0};
2495
2496
// Declare helper variables
2497
2498
// Compute basisvalues
2499
basisvalues[0] = 1.0;
2500
2501
// Table(s) of coefficients
2502
static const double coefficients0[1] = \
2503
{1.0};
2504
2505
// Compute value(s)
2506
for (unsigned int r = 0; r < 1; r++)
2507
{
2508
values[2] += coefficients0[r]*basisvalues[r];
2509
}// end loop over 'r'
2510
break;
2511
}
2512
}
2513
2514
}
2515
2516
/// Evaluate all basis functions at given point x in cell
2517
virtual void evaluate_basis_all(double* values,
2518
const double* x,
2519
const double* vertex_coordinates,
2520
int cell_orientation) const
2521
{
2522
// Helper variable to hold values of a single dof.
2523
double dof_values[3] = {0.0, 0.0, 0.0};
2524
2525
// Loop dofs and call evaluate_basis
2526
for (unsigned int r = 0; r < 7; r++)
2527
{
2528
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
2529
for (unsigned int s = 0; s < 3; s++)
2530
{
2531
values[r*3 + s] = dof_values[s];
2532
}// end loop over 's'
2533
}// end loop over 'r'
2534
}
2535
2536
/// Evaluate order n derivatives of basis function i at given point x in cell
2537
virtual void evaluate_basis_derivatives(std::size_t i,
2538
std::size_t n,
2539
double* values,
2540
const double* x,
2541
const double* vertex_coordinates,
2542
int cell_orientation) const
2543
{
2544
// Compute Jacobian
2545
double J[4];
2546
compute_jacobian_triangle_2d(J, vertex_coordinates);
2547
2548
// Compute Jacobian inverse and determinant
2549
double K[4];
2550
double detJ;
2551
compute_jacobian_inverse_triangle_2d(K, detJ, J);
2552
2553
2554
2555
// Compute constants
2556
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
2557
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
2558
2559
// Get coordinates and map to the reference (FIAT) element
2560
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
2561
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
2562
2563
// Compute number of derivatives.
2564
unsigned int num_derivatives = 1;
2565
for (unsigned int r = 0; r < n; r++)
2566
{
2567
num_derivatives *= 2;
2568
}// end loop over 'r'
2569
2570
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
2571
unsigned int **combinations = new unsigned int *[num_derivatives];
2572
for (unsigned int row = 0; row < num_derivatives; row++)
2573
{
2574
combinations[row] = new unsigned int [n];
2575
for (unsigned int col = 0; col < n; col++)
2576
combinations[row][col] = 0;
2577
}
2578
2579
// Generate combinations of derivatives
2580
for (unsigned int row = 1; row < num_derivatives; row++)
2581
{
2582
for (unsigned int num = 0; num < row; num++)
2583
{
2584
for (unsigned int col = n-1; col+1 > 0; col--)
2585
{
2586
if (combinations[row][col] + 1 > 1)
2587
combinations[row][col] = 0;
2588
else
2589
{
2590
combinations[row][col] += 1;
2591
break;
2592
}
2593
}
2594
}
2595
}
2596
2597
// Compute inverse of Jacobian
2598
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
2599
2600
// Declare transformation matrix
2601
// Declare pointer to two dimensional array and initialise
2602
double **transform = new double *[num_derivatives];
2603
2604
for (unsigned int j = 0; j < num_derivatives; j++)
2605
{
2606
transform[j] = new double [num_derivatives];
2607
for (unsigned int k = 0; k < num_derivatives; k++)
2608
transform[j][k] = 1;
2609
}
2610
2611
// Construct transformation matrix
2612
for (unsigned int row = 0; row < num_derivatives; row++)
2613
{
2614
for (unsigned int col = 0; col < num_derivatives; col++)
2615
{
2616
for (unsigned int k = 0; k < n; k++)
2617
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2618
}
2619
}
2620
2621
// Reset values. Assuming that values is always an array.
2622
for (unsigned int r = 0; r < 3*num_derivatives; r++)
2623
{
2624
values[r] = 0.0;
2625
}// end loop over 'r'
2626
2627
switch (i)
2628
{
2629
case 0:
2630
{
2631
2632
// Array of basisvalues
2633
double basisvalues[3] = {0.0, 0.0, 0.0};
2634
2635
// Declare helper variables
2636
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2637
2638
// Compute basisvalues
2639
basisvalues[0] = 1.0;
2640
basisvalues[1] = tmp0;
2641
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2642
basisvalues[0] *= std::sqrt(0.5);
2643
basisvalues[2] *= std::sqrt(1.0);
2644
basisvalues[1] *= std::sqrt(3.0);
2645
2646
// Table(s) of coefficients
2647
static const double coefficients0[3] = \
2648
{0.94280904, 0.57735027, -0.33333333};
2649
2650
static const double coefficients1[3] = \
2651
{-0.47140452, 0.0, -0.33333333};
2652
2653
// Tables of derivatives of the polynomial base (transpose).
2654
static const double dmats0[3][3] = \
2655
{{0.0, 0.0, 0.0},
2656
{4.8989795, 0.0, 0.0},
2657
{0.0, 0.0, 0.0}};
2658
2659
static const double dmats1[3][3] = \
2660
{{0.0, 0.0, 0.0},
2661
{2.4494897, 0.0, 0.0},
2662
{4.2426407, 0.0, 0.0}};
2663
2664
// Compute reference derivatives.
2665
// Declare pointer to array of derivatives on FIAT element.
2666
double *derivatives = new double[2*num_derivatives];
2667
for (unsigned int r = 0; r < 2*num_derivatives; r++)
2668
{
2669
derivatives[r] = 0.0;
2670
}// end loop over 'r'
2671
2672
// Declare pointer to array of reference derivatives on physical element.
2673
double *derivatives_p = new double[2*num_derivatives];
2674
for (unsigned int r = 0; r < 2*num_derivatives; r++)
2675
{
2676
derivatives_p[r] = 0.0;
2677
}// end loop over 'r'
2678
2679
// Declare derivative matrix (of polynomial basis).
2680
double dmats[3][3] = \
2681
{{1.0, 0.0, 0.0},
2682
{0.0, 1.0, 0.0},
2683
{0.0, 0.0, 1.0}};
2684
2685
// Declare (auxiliary) derivative matrix (of polynomial basis).
2686
double dmats_old[3][3] = \
2687
{{1.0, 0.0, 0.0},
2688
{0.0, 1.0, 0.0},
2689
{0.0, 0.0, 1.0}};
2690
2691
// Loop possible derivatives.
2692
for (unsigned int r = 0; r < num_derivatives; r++)
2693
{
2694
// Resetting dmats values to compute next derivative.
2695
for (unsigned int t = 0; t < 3; t++)
2696
{
2697
for (unsigned int u = 0; u < 3; u++)
2698
{
2699
dmats[t][u] = 0.0;
2700
if (t == u)
2701
{
2702
dmats[t][u] = 1.0;
2703
}
2704
2705
}// end loop over 'u'
2706
}// end loop over 't'
2707
2708
// Looping derivative order to generate dmats.
2709
for (unsigned int s = 0; s < n; s++)
2710
{
2711
// Updating dmats_old with new values and resetting dmats.
2712
for (unsigned int t = 0; t < 3; t++)
2713
{
2714
for (unsigned int u = 0; u < 3; u++)
2715
{
2716
dmats_old[t][u] = dmats[t][u];
2717
dmats[t][u] = 0.0;
2718
}// end loop over 'u'
2719
}// end loop over 't'
2720
2721
// Update dmats using an inner product.
2722
if (combinations[r][s] == 0)
2723
{
2724
for (unsigned int t = 0; t < 3; t++)
2725
{
2726
for (unsigned int u = 0; u < 3; u++)
2727
{
2728
for (unsigned int tu = 0; tu < 3; tu++)
2729
{
2730
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2731
}// end loop over 'tu'
2732
}// end loop over 'u'
2733
}// end loop over 't'
2734
}
2735
2736
if (combinations[r][s] == 1)
2737
{
2738
for (unsigned int t = 0; t < 3; t++)
2739
{
2740
for (unsigned int u = 0; u < 3; u++)
2741
{
2742
for (unsigned int tu = 0; tu < 3; tu++)
2743
{
2744
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2745
}// end loop over 'tu'
2746
}// end loop over 'u'
2747
}// end loop over 't'
2748
}
2749
2750
}// end loop over 's'
2751
for (unsigned int s = 0; s < 3; s++)
2752
{
2753
for (unsigned int t = 0; t < 3; t++)
2754
{
2755
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2756
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
2757
}// end loop over 't'
2758
}// end loop over 's'
2759
2760
// Using contravariant Piola transform to map values back to the physical element.
2761
const double tmp_ref0 = derivatives[r];
2762
const double tmp_ref1 = derivatives[num_derivatives + r];
2763
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2764
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2765
}// end loop over 'r'
2766
2767
// Transform derivatives back to physical element
2768
for (unsigned int r = 0; r < num_derivatives; r++)
2769
{
2770
for (unsigned int s = 0; s < num_derivatives; s++)
2771
{
2772
values[r] += transform[r][s]*derivatives_p[s];
2773
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
2774
}// end loop over 's'
2775
}// end loop over 'r'
2776
2777
// Delete pointer to array of derivatives on FIAT element
2778
delete [] derivatives;
2779
2780
// Delete pointer to array of reference derivatives on physical element.
2781
delete [] derivatives_p;
2782
2783
// Delete pointer to array of combinations of derivatives and transform
2784
for (unsigned int r = 0; r < num_derivatives; r++)
2785
{
2786
delete [] combinations[r];
2787
}// end loop over 'r'
2788
delete [] combinations;
2789
for (unsigned int r = 0; r < num_derivatives; r++)
2790
{
2791
delete [] transform[r];
2792
}// end loop over 'r'
2793
delete [] transform;
2794
break;
2795
}
2796
case 1:
2797
{
2798
2799
// Array of basisvalues
2800
double basisvalues[3] = {0.0, 0.0, 0.0};
2801
2802
// Declare helper variables
2803
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2804
2805
// Compute basisvalues
2806
basisvalues[0] = 1.0;
2807
basisvalues[1] = tmp0;
2808
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2809
basisvalues[0] *= std::sqrt(0.5);
2810
basisvalues[2] *= std::sqrt(1.0);
2811
basisvalues[1] *= std::sqrt(3.0);
2812
2813
// Table(s) of coefficients
2814
static const double coefficients0[3] = \
2815
{-0.47140452, -0.28867513, 0.16666667};
2816
2817
static const double coefficients1[3] = \
2818
{0.94280904, 0.0, 0.66666667};
2819
2820
// Tables of derivatives of the polynomial base (transpose).
2821
static const double dmats0[3][3] = \
2822
{{0.0, 0.0, 0.0},
2823
{4.8989795, 0.0, 0.0},
2824
{0.0, 0.0, 0.0}};
2825
2826
static const double dmats1[3][3] = \
2827
{{0.0, 0.0, 0.0},
2828
{2.4494897, 0.0, 0.0},
2829
{4.2426407, 0.0, 0.0}};
2830
2831
// Compute reference derivatives.
2832
// Declare pointer to array of derivatives on FIAT element.
2833
double *derivatives = new double[2*num_derivatives];
2834
for (unsigned int r = 0; r < 2*num_derivatives; r++)
2835
{
2836
derivatives[r] = 0.0;
2837
}// end loop over 'r'
2838
2839
// Declare pointer to array of reference derivatives on physical element.
2840
double *derivatives_p = new double[2*num_derivatives];
2841
for (unsigned int r = 0; r < 2*num_derivatives; r++)
2842
{
2843
derivatives_p[r] = 0.0;
2844
}// end loop over 'r'
2845
2846
// Declare derivative matrix (of polynomial basis).
2847
double dmats[3][3] = \
2848
{{1.0, 0.0, 0.0},
2849
{0.0, 1.0, 0.0},
2850
{0.0, 0.0, 1.0}};
2851
2852
// Declare (auxiliary) derivative matrix (of polynomial basis).
2853
double dmats_old[3][3] = \
2854
{{1.0, 0.0, 0.0},
2855
{0.0, 1.0, 0.0},
2856
{0.0, 0.0, 1.0}};
2857
2858
// Loop possible derivatives.
2859
for (unsigned int r = 0; r < num_derivatives; r++)
2860
{
2861
// Resetting dmats values to compute next derivative.
2862
for (unsigned int t = 0; t < 3; t++)
2863
{
2864
for (unsigned int u = 0; u < 3; u++)
2865
{
2866
dmats[t][u] = 0.0;
2867
if (t == u)
2868
{
2869
dmats[t][u] = 1.0;
2870
}
2871
2872
}// end loop over 'u'
2873
}// end loop over 't'
2874
2875
// Looping derivative order to generate dmats.
2876
for (unsigned int s = 0; s < n; s++)
2877
{
2878
// Updating dmats_old with new values and resetting dmats.
2879
for (unsigned int t = 0; t < 3; t++)
2880
{
2881
for (unsigned int u = 0; u < 3; u++)
2882
{
2883
dmats_old[t][u] = dmats[t][u];
2884
dmats[t][u] = 0.0;
2885
}// end loop over 'u'
2886
}// end loop over 't'
2887
2888
// Update dmats using an inner product.
2889
if (combinations[r][s] == 0)
2890
{
2891
for (unsigned int t = 0; t < 3; t++)
2892
{
2893
for (unsigned int u = 0; u < 3; u++)
2894
{
2895
for (unsigned int tu = 0; tu < 3; tu++)
2896
{
2897
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2898
}// end loop over 'tu'
2899
}// end loop over 'u'
2900
}// end loop over 't'
2901
}
2902
2903
if (combinations[r][s] == 1)
2904
{
2905
for (unsigned int t = 0; t < 3; t++)
2906
{
2907
for (unsigned int u = 0; u < 3; u++)
2908
{
2909
for (unsigned int tu = 0; tu < 3; tu++)
2910
{
2911
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2912
}// end loop over 'tu'
2913
}// end loop over 'u'
2914
}// end loop over 't'
2915
}
2916
2917
}// end loop over 's'
2918
for (unsigned int s = 0; s < 3; s++)
2919
{
2920
for (unsigned int t = 0; t < 3; t++)
2921
{
2922
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2923
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
2924
}// end loop over 't'
2925
}// end loop over 's'
2926
2927
// Using contravariant Piola transform to map values back to the physical element.
2928
const double tmp_ref0 = derivatives[r];
2929
const double tmp_ref1 = derivatives[num_derivatives + r];
2930
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2931
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2932
}// end loop over 'r'
2933
2934
// Transform derivatives back to physical element
2935
for (unsigned int r = 0; r < num_derivatives; r++)
2936
{
2937
for (unsigned int s = 0; s < num_derivatives; s++)
2938
{
2939
values[r] += transform[r][s]*derivatives_p[s];
2940
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
2941
}// end loop over 's'
2942
}// end loop over 'r'
2943
2944
// Delete pointer to array of derivatives on FIAT element
2945
delete [] derivatives;
2946
2947
// Delete pointer to array of reference derivatives on physical element.
2948
delete [] derivatives_p;
2949
2950
// Delete pointer to array of combinations of derivatives and transform
2951
for (unsigned int r = 0; r < num_derivatives; r++)
2952
{
2953
delete [] combinations[r];
2954
}// end loop over 'r'
2955
delete [] combinations;
2956
for (unsigned int r = 0; r < num_derivatives; r++)
2957
{
2958
delete [] transform[r];
2959
}// end loop over 'r'
2960
delete [] transform;
2961
break;
2962
}
2963
case 2:
2964
{
2965
2966
// Array of basisvalues
2967
double basisvalues[3] = {0.0, 0.0, 0.0};
2968
2969
// Declare helper variables
2970
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2971
2972
// Compute basisvalues
2973
basisvalues[0] = 1.0;
2974
basisvalues[1] = tmp0;
2975
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2976
basisvalues[0] *= std::sqrt(0.5);
2977
basisvalues[2] *= std::sqrt(1.0);
2978
basisvalues[1] *= std::sqrt(3.0);
2979
2980
// Table(s) of coefficients
2981
static const double coefficients0[3] = \
2982
{0.47140452, -0.57735027, -0.66666667};
2983
2984
static const double coefficients1[3] = \
2985
{0.47140452, 0.0, 0.33333333};
2986
2987
// Tables of derivatives of the polynomial base (transpose).
2988
static const double dmats0[3][3] = \
2989
{{0.0, 0.0, 0.0},
2990
{4.8989795, 0.0, 0.0},
2991
{0.0, 0.0, 0.0}};
2992
2993
static const double dmats1[3][3] = \
2994
{{0.0, 0.0, 0.0},
2995
{2.4494897, 0.0, 0.0},
2996
{4.2426407, 0.0, 0.0}};
2997
2998
// Compute reference derivatives.
2999
// Declare pointer to array of derivatives on FIAT element.
3000
double *derivatives = new double[2*num_derivatives];
3001
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3002
{
3003
derivatives[r] = 0.0;
3004
}// end loop over 'r'
3005
3006
// Declare pointer to array of reference derivatives on physical element.
3007
double *derivatives_p = new double[2*num_derivatives];
3008
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3009
{
3010
derivatives_p[r] = 0.0;
3011
}// end loop over 'r'
3012
3013
// Declare derivative matrix (of polynomial basis).
3014
double dmats[3][3] = \
3015
{{1.0, 0.0, 0.0},
3016
{0.0, 1.0, 0.0},
3017
{0.0, 0.0, 1.0}};
3018
3019
// Declare (auxiliary) derivative matrix (of polynomial basis).
3020
double dmats_old[3][3] = \
3021
{{1.0, 0.0, 0.0},
3022
{0.0, 1.0, 0.0},
3023
{0.0, 0.0, 1.0}};
3024
3025
// Loop possible derivatives.
3026
for (unsigned int r = 0; r < num_derivatives; r++)
3027
{
3028
// Resetting dmats values to compute next derivative.
3029
for (unsigned int t = 0; t < 3; t++)
3030
{
3031
for (unsigned int u = 0; u < 3; u++)
3032
{
3033
dmats[t][u] = 0.0;
3034
if (t == u)
3035
{
3036
dmats[t][u] = 1.0;
3037
}
3038
3039
}// end loop over 'u'
3040
}// end loop over 't'
3041
3042
// Looping derivative order to generate dmats.
3043
for (unsigned int s = 0; s < n; s++)
3044
{
3045
// Updating dmats_old with new values and resetting dmats.
3046
for (unsigned int t = 0; t < 3; t++)
3047
{
3048
for (unsigned int u = 0; u < 3; u++)
3049
{
3050
dmats_old[t][u] = dmats[t][u];
3051
dmats[t][u] = 0.0;
3052
}// end loop over 'u'
3053
}// end loop over 't'
3054
3055
// Update dmats using an inner product.
3056
if (combinations[r][s] == 0)
3057
{
3058
for (unsigned int t = 0; t < 3; t++)
3059
{
3060
for (unsigned int u = 0; u < 3; u++)
3061
{
3062
for (unsigned int tu = 0; tu < 3; tu++)
3063
{
3064
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3065
}// end loop over 'tu'
3066
}// end loop over 'u'
3067
}// end loop over 't'
3068
}
3069
3070
if (combinations[r][s] == 1)
3071
{
3072
for (unsigned int t = 0; t < 3; t++)
3073
{
3074
for (unsigned int u = 0; u < 3; u++)
3075
{
3076
for (unsigned int tu = 0; tu < 3; tu++)
3077
{
3078
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3079
}// end loop over 'tu'
3080
}// end loop over 'u'
3081
}// end loop over 't'
3082
}
3083
3084
}// end loop over 's'
3085
for (unsigned int s = 0; s < 3; s++)
3086
{
3087
for (unsigned int t = 0; t < 3; t++)
3088
{
3089
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3090
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
3091
}// end loop over 't'
3092
}// end loop over 's'
3093
3094
// Using contravariant Piola transform to map values back to the physical element.
3095
const double tmp_ref0 = derivatives[r];
3096
const double tmp_ref1 = derivatives[num_derivatives + r];
3097
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
3098
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
3099
}// end loop over 'r'
3100
3101
// Transform derivatives back to physical element
3102
for (unsigned int r = 0; r < num_derivatives; r++)
3103
{
3104
for (unsigned int s = 0; s < num_derivatives; s++)
3105
{
3106
values[r] += transform[r][s]*derivatives_p[s];
3107
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
3108
}// end loop over 's'
3109
}// end loop over 'r'
3110
3111
// Delete pointer to array of derivatives on FIAT element
3112
delete [] derivatives;
3113
3114
// Delete pointer to array of reference derivatives on physical element.
3115
delete [] derivatives_p;
3116
3117
// Delete pointer to array of combinations of derivatives and transform
3118
for (unsigned int r = 0; r < num_derivatives; r++)
3119
{
3120
delete [] combinations[r];
3121
}// end loop over 'r'
3122
delete [] combinations;
3123
for (unsigned int r = 0; r < num_derivatives; r++)
3124
{
3125
delete [] transform[r];
3126
}// end loop over 'r'
3127
delete [] transform;
3128
break;
3129
}
3130
case 3:
3131
{
3132
3133
// Array of basisvalues
3134
double basisvalues[3] = {0.0, 0.0, 0.0};
3135
3136
// Declare helper variables
3137
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3138
3139
// Compute basisvalues
3140
basisvalues[0] = 1.0;
3141
basisvalues[1] = tmp0;
3142
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3143
basisvalues[0] *= std::sqrt(0.5);
3144
basisvalues[2] *= std::sqrt(1.0);
3145
basisvalues[1] *= std::sqrt(3.0);
3146
3147
// Table(s) of coefficients
3148
static const double coefficients0[3] = \
3149
{0.47140452, 0.28867513, 0.83333333};
3150
3151
static const double coefficients1[3] = \
3152
{-0.94280904, 0.0, -0.66666667};
3153
3154
// Tables of derivatives of the polynomial base (transpose).
3155
static const double dmats0[3][3] = \
3156
{{0.0, 0.0, 0.0},
3157
{4.8989795, 0.0, 0.0},
3158
{0.0, 0.0, 0.0}};
3159
3160
static const double dmats1[3][3] = \
3161
{{0.0, 0.0, 0.0},
3162
{2.4494897, 0.0, 0.0},
3163
{4.2426407, 0.0, 0.0}};
3164
3165
// Compute reference derivatives.
3166
// Declare pointer to array of derivatives on FIAT element.
3167
double *derivatives = new double[2*num_derivatives];
3168
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3169
{
3170
derivatives[r] = 0.0;
3171
}// end loop over 'r'
3172
3173
// Declare pointer to array of reference derivatives on physical element.
3174
double *derivatives_p = new double[2*num_derivatives];
3175
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3176
{
3177
derivatives_p[r] = 0.0;
3178
}// end loop over 'r'
3179
3180
// Declare derivative matrix (of polynomial basis).
3181
double dmats[3][3] = \
3182
{{1.0, 0.0, 0.0},
3183
{0.0, 1.0, 0.0},
3184
{0.0, 0.0, 1.0}};
3185
3186
// Declare (auxiliary) derivative matrix (of polynomial basis).
3187
double dmats_old[3][3] = \
3188
{{1.0, 0.0, 0.0},
3189
{0.0, 1.0, 0.0},
3190
{0.0, 0.0, 1.0}};
3191
3192
// Loop possible derivatives.
3193
for (unsigned int r = 0; r < num_derivatives; r++)
3194
{
3195
// Resetting dmats values to compute next derivative.
3196
for (unsigned int t = 0; t < 3; t++)
3197
{
3198
for (unsigned int u = 0; u < 3; u++)
3199
{
3200
dmats[t][u] = 0.0;
3201
if (t == u)
3202
{
3203
dmats[t][u] = 1.0;
3204
}
3205
3206
}// end loop over 'u'
3207
}// end loop over 't'
3208
3209
// Looping derivative order to generate dmats.
3210
for (unsigned int s = 0; s < n; s++)
3211
{
3212
// Updating dmats_old with new values and resetting dmats.
3213
for (unsigned int t = 0; t < 3; t++)
3214
{
3215
for (unsigned int u = 0; u < 3; u++)
3216
{
3217
dmats_old[t][u] = dmats[t][u];
3218
dmats[t][u] = 0.0;
3219
}// end loop over 'u'
3220
}// end loop over 't'
3221
3222
// Update dmats using an inner product.
3223
if (combinations[r][s] == 0)
3224
{
3225
for (unsigned int t = 0; t < 3; t++)
3226
{
3227
for (unsigned int u = 0; u < 3; u++)
3228
{
3229
for (unsigned int tu = 0; tu < 3; tu++)
3230
{
3231
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3232
}// end loop over 'tu'
3233
}// end loop over 'u'
3234
}// end loop over 't'
3235
}
3236
3237
if (combinations[r][s] == 1)
3238
{
3239
for (unsigned int t = 0; t < 3; t++)
3240
{
3241
for (unsigned int u = 0; u < 3; u++)
3242
{
3243
for (unsigned int tu = 0; tu < 3; tu++)
3244
{
3245
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3246
}// end loop over 'tu'
3247
}// end loop over 'u'
3248
}// end loop over 't'
3249
}
3250
3251
}// end loop over 's'
3252
for (unsigned int s = 0; s < 3; s++)
3253
{
3254
for (unsigned int t = 0; t < 3; t++)
3255
{
3256
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3257
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
3258
}// end loop over 't'
3259
}// end loop over 's'
3260
3261
// Using contravariant Piola transform to map values back to the physical element.
3262
const double tmp_ref0 = derivatives[r];
3263
const double tmp_ref1 = derivatives[num_derivatives + r];
3264
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
3265
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
3266
}// end loop over 'r'
3267
3268
// Transform derivatives back to physical element
3269
for (unsigned int r = 0; r < num_derivatives; r++)
3270
{
3271
for (unsigned int s = 0; s < num_derivatives; s++)
3272
{
3273
values[r] += transform[r][s]*derivatives_p[s];
3274
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
3275
}// end loop over 's'
3276
}// end loop over 'r'
3277
3278
// Delete pointer to array of derivatives on FIAT element
3279
delete [] derivatives;
3280
3281
// Delete pointer to array of reference derivatives on physical element.
3282
delete [] derivatives_p;
3283
3284
// Delete pointer to array of combinations of derivatives and transform
3285
for (unsigned int r = 0; r < num_derivatives; r++)
3286
{
3287
delete [] combinations[r];
3288
}// end loop over 'r'
3289
delete [] combinations;
3290
for (unsigned int r = 0; r < num_derivatives; r++)
3291
{
3292
delete [] transform[r];
3293
}// end loop over 'r'
3294
delete [] transform;
3295
break;
3296
}
3297
case 4:
3298
{
3299
3300
// Array of basisvalues
3301
double basisvalues[3] = {0.0, 0.0, 0.0};
3302
3303
// Declare helper variables
3304
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3305
3306
// Compute basisvalues
3307
basisvalues[0] = 1.0;
3308
basisvalues[1] = tmp0;
3309
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3310
basisvalues[0] *= std::sqrt(0.5);
3311
basisvalues[2] *= std::sqrt(1.0);
3312
basisvalues[1] *= std::sqrt(3.0);
3313
3314
// Table(s) of coefficients
3315
static const double coefficients0[3] = \
3316
{-0.47140452, -0.28867513, 0.16666667};
3317
3318
static const double coefficients1[3] = \
3319
{-0.47140452, 0.8660254, 0.16666667};
3320
3321
// Tables of derivatives of the polynomial base (transpose).
3322
static const double dmats0[3][3] = \
3323
{{0.0, 0.0, 0.0},
3324
{4.8989795, 0.0, 0.0},
3325
{0.0, 0.0, 0.0}};
3326
3327
static const double dmats1[3][3] = \
3328
{{0.0, 0.0, 0.0},
3329
{2.4494897, 0.0, 0.0},
3330
{4.2426407, 0.0, 0.0}};
3331
3332
// Compute reference derivatives.
3333
// Declare pointer to array of derivatives on FIAT element.
3334
double *derivatives = new double[2*num_derivatives];
3335
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3336
{
3337
derivatives[r] = 0.0;
3338
}// end loop over 'r'
3339
3340
// Declare pointer to array of reference derivatives on physical element.
3341
double *derivatives_p = new double[2*num_derivatives];
3342
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3343
{
3344
derivatives_p[r] = 0.0;
3345
}// end loop over 'r'
3346
3347
// Declare derivative matrix (of polynomial basis).
3348
double dmats[3][3] = \
3349
{{1.0, 0.0, 0.0},
3350
{0.0, 1.0, 0.0},
3351
{0.0, 0.0, 1.0}};
3352
3353
// Declare (auxiliary) derivative matrix (of polynomial basis).
3354
double dmats_old[3][3] = \
3355
{{1.0, 0.0, 0.0},
3356
{0.0, 1.0, 0.0},
3357
{0.0, 0.0, 1.0}};
3358
3359
// Loop possible derivatives.
3360
for (unsigned int r = 0; r < num_derivatives; r++)
3361
{
3362
// Resetting dmats values to compute next derivative.
3363
for (unsigned int t = 0; t < 3; t++)
3364
{
3365
for (unsigned int u = 0; u < 3; u++)
3366
{
3367
dmats[t][u] = 0.0;
3368
if (t == u)
3369
{
3370
dmats[t][u] = 1.0;
3371
}
3372
3373
}// end loop over 'u'
3374
}// end loop over 't'
3375
3376
// Looping derivative order to generate dmats.
3377
for (unsigned int s = 0; s < n; s++)
3378
{
3379
// Updating dmats_old with new values and resetting dmats.
3380
for (unsigned int t = 0; t < 3; t++)
3381
{
3382
for (unsigned int u = 0; u < 3; u++)
3383
{
3384
dmats_old[t][u] = dmats[t][u];
3385
dmats[t][u] = 0.0;
3386
}// end loop over 'u'
3387
}// end loop over 't'
3388
3389
// Update dmats using an inner product.
3390
if (combinations[r][s] == 0)
3391
{
3392
for (unsigned int t = 0; t < 3; t++)
3393
{
3394
for (unsigned int u = 0; u < 3; u++)
3395
{
3396
for (unsigned int tu = 0; tu < 3; tu++)
3397
{
3398
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3399
}// end loop over 'tu'
3400
}// end loop over 'u'
3401
}// end loop over 't'
3402
}
3403
3404
if (combinations[r][s] == 1)
3405
{
3406
for (unsigned int t = 0; t < 3; t++)
3407
{
3408
for (unsigned int u = 0; u < 3; u++)
3409
{
3410
for (unsigned int tu = 0; tu < 3; tu++)
3411
{
3412
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3413
}// end loop over 'tu'
3414
}// end loop over 'u'
3415
}// end loop over 't'
3416
}
3417
3418
}// end loop over 's'
3419
for (unsigned int s = 0; s < 3; s++)
3420
{
3421
for (unsigned int t = 0; t < 3; t++)
3422
{
3423
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3424
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
3425
}// end loop over 't'
3426
}// end loop over 's'
3427
3428
// Using contravariant Piola transform to map values back to the physical element.
3429
const double tmp_ref0 = derivatives[r];
3430
const double tmp_ref1 = derivatives[num_derivatives + r];
3431
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
3432
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
3433
}// end loop over 'r'
3434
3435
// Transform derivatives back to physical element
3436
for (unsigned int r = 0; r < num_derivatives; r++)
3437
{
3438
for (unsigned int s = 0; s < num_derivatives; s++)
3439
{
3440
values[r] += transform[r][s]*derivatives_p[s];
3441
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
3442
}// end loop over 's'
3443
}// end loop over 'r'
3444
3445
// Delete pointer to array of derivatives on FIAT element
3446
delete [] derivatives;
3447
3448
// Delete pointer to array of reference derivatives on physical element.
3449
delete [] derivatives_p;
3450
3451
// Delete pointer to array of combinations of derivatives and transform
3452
for (unsigned int r = 0; r < num_derivatives; r++)
3453
{
3454
delete [] combinations[r];
3455
}// end loop over 'r'
3456
delete [] combinations;
3457
for (unsigned int r = 0; r < num_derivatives; r++)
3458
{
3459
delete [] transform[r];
3460
}// end loop over 'r'
3461
delete [] transform;
3462
break;
3463
}
3464
case 5:
3465
{
3466
3467
// Array of basisvalues
3468
double basisvalues[3] = {0.0, 0.0, 0.0};
3469
3470
// Declare helper variables
3471
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
3472
3473
// Compute basisvalues
3474
basisvalues[0] = 1.0;
3475
basisvalues[1] = tmp0;
3476
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
3477
basisvalues[0] *= std::sqrt(0.5);
3478
basisvalues[2] *= std::sqrt(1.0);
3479
basisvalues[1] *= std::sqrt(3.0);
3480
3481
// Table(s) of coefficients
3482
static const double coefficients0[3] = \
3483
{0.94280904, 0.57735027, -0.33333333};
3484
3485
static const double coefficients1[3] = \
3486
{-0.47140452, -0.8660254, 0.16666667};
3487
3488
// Tables of derivatives of the polynomial base (transpose).
3489
static const double dmats0[3][3] = \
3490
{{0.0, 0.0, 0.0},
3491
{4.8989795, 0.0, 0.0},
3492
{0.0, 0.0, 0.0}};
3493
3494
static const double dmats1[3][3] = \
3495
{{0.0, 0.0, 0.0},
3496
{2.4494897, 0.0, 0.0},
3497
{4.2426407, 0.0, 0.0}};
3498
3499
// Compute reference derivatives.
3500
// Declare pointer to array of derivatives on FIAT element.
3501
double *derivatives = new double[2*num_derivatives];
3502
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3503
{
3504
derivatives[r] = 0.0;
3505
}// end loop over 'r'
3506
3507
// Declare pointer to array of reference derivatives on physical element.
3508
double *derivatives_p = new double[2*num_derivatives];
3509
for (unsigned int r = 0; r < 2*num_derivatives; r++)
3510
{
3511
derivatives_p[r] = 0.0;
3512
}// end loop over 'r'
3513
3514
// Declare derivative matrix (of polynomial basis).
3515
double dmats[3][3] = \
3516
{{1.0, 0.0, 0.0},
3517
{0.0, 1.0, 0.0},
3518
{0.0, 0.0, 1.0}};
3519
3520
// Declare (auxiliary) derivative matrix (of polynomial basis).
3521
double dmats_old[3][3] = \
3522
{{1.0, 0.0, 0.0},
3523
{0.0, 1.0, 0.0},
3524
{0.0, 0.0, 1.0}};
3525
3526
// Loop possible derivatives.
3527
for (unsigned int r = 0; r < num_derivatives; r++)
3528
{
3529
// Resetting dmats values to compute next derivative.
3530
for (unsigned int t = 0; t < 3; t++)
3531
{
3532
for (unsigned int u = 0; u < 3; u++)
3533
{
3534
dmats[t][u] = 0.0;
3535
if (t == u)
3536
{
3537
dmats[t][u] = 1.0;
3538
}
3539
3540
}// end loop over 'u'
3541
}// end loop over 't'
3542
3543
// Looping derivative order to generate dmats.
3544
for (unsigned int s = 0; s < n; s++)
3545
{
3546
// Updating dmats_old with new values and resetting dmats.
3547
for (unsigned int t = 0; t < 3; t++)
3548
{
3549
for (unsigned int u = 0; u < 3; u++)
3550
{
3551
dmats_old[t][u] = dmats[t][u];
3552
dmats[t][u] = 0.0;
3553
}// end loop over 'u'
3554
}// end loop over 't'
3555
3556
// Update dmats using an inner product.
3557
if (combinations[r][s] == 0)
3558
{
3559
for (unsigned int t = 0; t < 3; t++)
3560
{
3561
for (unsigned int u = 0; u < 3; u++)
3562
{
3563
for (unsigned int tu = 0; tu < 3; tu++)
3564
{
3565
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3566
}// end loop over 'tu'
3567
}// end loop over 'u'
3568
}// end loop over 't'
3569
}
3570
3571
if (combinations[r][s] == 1)
3572
{
3573
for (unsigned int t = 0; t < 3; t++)
3574
{
3575
for (unsigned int u = 0; u < 3; u++)
3576
{
3577
for (unsigned int tu = 0; tu < 3; tu++)
3578
{
3579
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3580
}// end loop over 'tu'
3581
}// end loop over 'u'
3582
}// end loop over 't'
3583
}
3584
3585
}// end loop over 's'
3586
for (unsigned int s = 0; s < 3; s++)
3587
{
3588
for (unsigned int t = 0; t < 3; t++)
3589
{
3590
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3591
derivatives[num_derivatives + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
3592
}// end loop over 't'
3593
}// end loop over 's'
3594
3595
// Using contravariant Piola transform to map values back to the physical element.
3596
const double tmp_ref0 = derivatives[r];
3597
const double tmp_ref1 = derivatives[num_derivatives + r];
3598
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
3599
derivatives_p[num_derivatives + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
3600
}// end loop over 'r'
3601
3602
// Transform derivatives back to physical element
3603
for (unsigned int r = 0; r < num_derivatives; r++)
3604
{
3605
for (unsigned int s = 0; s < num_derivatives; s++)
3606
{
3607
values[r] += transform[r][s]*derivatives_p[s];
3608
values[num_derivatives + r] += transform[r][s]*derivatives_p[num_derivatives + s];
3609
}// end loop over 's'
3610
}// end loop over 'r'
3611
3612
// Delete pointer to array of derivatives on FIAT element
3613
delete [] derivatives;
3614
3615
// Delete pointer to array of reference derivatives on physical element.
3616
delete [] derivatives_p;
3617
3618
// Delete pointer to array of combinations of derivatives and transform
3619
for (unsigned int r = 0; r < num_derivatives; r++)
3620
{
3621
delete [] combinations[r];
3622
}// end loop over 'r'
3623
delete [] combinations;
3624
for (unsigned int r = 0; r < num_derivatives; r++)
3625
{
3626
delete [] transform[r];
3627
}// end loop over 'r'
3628
delete [] transform;
3629
break;
3630
}
3631
case 6:
3632
{
3633
3634
// Array of basisvalues
3635
double basisvalues[1] = {0.0};
3636
3637
// Declare helper variables
3638
3639
// Compute basisvalues
3640
basisvalues[0] = 1.0;
3641
3642
// Table(s) of coefficients
3643
static const double coefficients0[1] = \
3644
{1.0};
3645
3646
// Tables of derivatives of the polynomial base (transpose).
3647
static const double dmats0[1][1] = \
3648
{{0.0}};
3649
3650
static const double dmats1[1][1] = \
3651
{{0.0}};
3652
3653
// Compute reference derivatives.
3654
// Declare pointer to array of derivatives on FIAT element.
3655
double *derivatives = new double[num_derivatives];
3656
for (unsigned int r = 0; r < num_derivatives; r++)
3657
{
3658
derivatives[r] = 0.0;
3659
}// end loop over 'r'
3660
3661
// Declare derivative matrix (of polynomial basis).
3662
double dmats[1][1] = \
3663
{{1.0}};
3664
3665
// Declare (auxiliary) derivative matrix (of polynomial basis).
3666
double dmats_old[1][1] = \
3667
{{1.0}};
3668
3669
// Loop possible derivatives.
3670
for (unsigned int r = 0; r < num_derivatives; r++)
3671
{
3672
// Resetting dmats values to compute next derivative.
3673
for (unsigned int t = 0; t < 1; t++)
3674
{
3675
for (unsigned int u = 0; u < 1; u++)
3676
{
3677
dmats[t][u] = 0.0;
3678
if (t == u)
3679
{
3680
dmats[t][u] = 1.0;
3681
}
3682
3683
}// end loop over 'u'
3684
}// end loop over 't'
3685
3686
// Looping derivative order to generate dmats.
3687
for (unsigned int s = 0; s < n; s++)
3688
{
3689
// Updating dmats_old with new values and resetting dmats.
3690
for (unsigned int t = 0; t < 1; t++)
3691
{
3692
for (unsigned int u = 0; u < 1; u++)
3693
{
3694
dmats_old[t][u] = dmats[t][u];
3695
dmats[t][u] = 0.0;
3696
}// end loop over 'u'
3697
}// end loop over 't'
3698
3699
// Update dmats using an inner product.
3700
if (combinations[r][s] == 0)
3701
{
3702
for (unsigned int t = 0; t < 1; t++)
3703
{
3704
for (unsigned int u = 0; u < 1; u++)
3705
{
3706
for (unsigned int tu = 0; tu < 1; tu++)
3707
{
3708
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3709
}// end loop over 'tu'
3710
}// end loop over 'u'
3711
}// end loop over 't'
3712
}
3713
3714
if (combinations[r][s] == 1)
3715
{
3716
for (unsigned int t = 0; t < 1; t++)
3717
{
3718
for (unsigned int u = 0; u < 1; u++)
3719
{
3720
for (unsigned int tu = 0; tu < 1; tu++)
3721
{
3722
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3723
}// end loop over 'tu'
3724
}// end loop over 'u'
3725
}// end loop over 't'
3726
}
3727
3728
}// end loop over 's'
3729
for (unsigned int s = 0; s < 1; s++)
3730
{
3731
for (unsigned int t = 0; t < 1; t++)
3732
{
3733
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3734
}// end loop over 't'
3735
}// end loop over 's'
3736
}// end loop over 'r'
3737
3738
// Transform derivatives back to physical element
3739
for (unsigned int r = 0; r < num_derivatives; r++)
3740
{
3741
for (unsigned int s = 0; s < num_derivatives; s++)
3742
{
3743
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
3744
}// end loop over 's'
3745
}// end loop over 'r'
3746
3747
// Delete pointer to array of derivatives on FIAT element
3748
delete [] derivatives;
3749
3750
// Delete pointer to array of combinations of derivatives and transform
3751
for (unsigned int r = 0; r < num_derivatives; r++)
3752
{
3753
delete [] combinations[r];
3754
}// end loop over 'r'
3755
delete [] combinations;
3756
for (unsigned int r = 0; r < num_derivatives; r++)
3757
{
3758
delete [] transform[r];
3759
}// end loop over 'r'
3760
delete [] transform;
3761
break;
3762
}
3763
}
3764
3765
}
3766
3767
/// Evaluate order n derivatives of all basis functions at given point x in cell
3768
virtual void evaluate_basis_derivatives_all(std::size_t n,
3769
double* values,
3770
const double* x,
3771
const double* vertex_coordinates,
3772
int cell_orientation) const
3773
{
3774
// Compute number of derivatives.
3775
unsigned int num_derivatives = 1;
3776
for (unsigned int r = 0; r < n; r++)
3777
{
3778
num_derivatives *= 2;
3779
}// end loop over 'r'
3780
3781
// Helper variable to hold values of a single dof.
3782
double *dof_values = new double[3*num_derivatives];
3783
for (unsigned int r = 0; r < 3*num_derivatives; r++)
3784
{
3785
dof_values[r] = 0.0;
3786
}// end loop over 'r'
3787
3788
// Loop dofs and call evaluate_basis_derivatives.
3789
for (unsigned int r = 0; r < 7; r++)
3790
{
3791
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
3792
for (unsigned int s = 0; s < 3*num_derivatives; s++)
3793
{
3794
values[r*3*num_derivatives + s] = dof_values[s];
3795
}// end loop over 's'
3796
}// end loop over 'r'
3797
3798
// Delete pointer.
3799
delete [] dof_values;
3800
}
3801
3802
/// Evaluate linear functional for dof i on the function f
3803
virtual double evaluate_dof(std::size_t i,
3804
const ufc::function& f,
3805
const double* vertex_coordinates,
3806
int cell_orientation,
3807
const ufc::cell& c) const
3808
{
3809
// Declare variables for result of evaluation
3810
double vals[3];
3811
3812
// Declare variable for physical coordinates
3813
double y[2];
3814
3815
double result;
3816
// Compute Jacobian
3817
double J[4];
3818
compute_jacobian_triangle_2d(J, vertex_coordinates);
3819
3820
3821
// Compute Jacobian inverse and determinant
3822
double K[4];
3823
double detJ;
3824
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3825
3826
3827
switch (i)
3828
{
3829
case 0:
3830
{
3831
y[0] = 0.66666667*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
3832
y[1] = 0.66666667*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
3833
f.evaluate(vals, y, c);
3834
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
3835
return result;
3836
break;
3837
}
3838
case 1:
3839
{
3840
y[0] = 0.33333333*vertex_coordinates[2] + 0.66666667*vertex_coordinates[4];
3841
y[1] = 0.33333333*vertex_coordinates[3] + 0.66666667*vertex_coordinates[5];
3842
f.evaluate(vals, y, c);
3843
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
3844
return result;
3845
break;
3846
}
3847
case 2:
3848
{
3849
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[4];
3850
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[5];
3851
f.evaluate(vals, y, c);
3852
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
3853
return result;
3854
break;
3855
}
3856
case 3:
3857
{
3858
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[4];
3859
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[5];
3860
f.evaluate(vals, y, c);
3861
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
3862
return result;
3863
break;
3864
}
3865
case 4:
3866
{
3867
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2];
3868
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3];
3869
f.evaluate(vals, y, c);
3870
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
3871
return result;
3872
break;
3873
}
3874
case 5:
3875
{
3876
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[2];
3877
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[3];
3878
f.evaluate(vals, y, c);
3879
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
3880
return result;
3881
break;
3882
}
3883
case 6:
3884
{
3885
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
3886
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
3887
f.evaluate(vals, y, c);
3888
return vals[2];
3889
break;
3890
}
3891
}
3892
3893
return 0.0;
3894
}
3895
3896
/// Evaluate linear functionals for all dofs on the function f
3897
virtual void evaluate_dofs(double* values,
3898
const ufc::function& f,
3899
const double* vertex_coordinates,
3900
int cell_orientation,
3901
const ufc::cell& c) const
3902
{
3903
// Declare variables for result of evaluation
3904
double vals[3];
3905
3906
// Declare variable for physical coordinates
3907
double y[2];
3908
3909
double result;
3910
// Compute Jacobian
3911
double J[4];
3912
compute_jacobian_triangle_2d(J, vertex_coordinates);
3913
3914
3915
// Compute Jacobian inverse and determinant
3916
double K[4];
3917
double detJ;
3918
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3919
3920
3921
y[0] = 0.66666667*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
3922
y[1] = 0.66666667*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
3923
f.evaluate(vals, y, c);
3924
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
3925
values[0] = result;
3926
y[0] = 0.33333333*vertex_coordinates[2] + 0.66666667*vertex_coordinates[4];
3927
y[1] = 0.33333333*vertex_coordinates[3] + 0.66666667*vertex_coordinates[5];
3928
f.evaluate(vals, y, c);
3929
result = (detJ*(K[0]*vals[0] + K[1]*vals[1])) + (detJ*(K[2]*vals[0] + K[3]*vals[1]));
3930
values[1] = result;
3931
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[4];
3932
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[5];
3933
f.evaluate(vals, y, c);
3934
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
3935
values[2] = result;
3936
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[4];
3937
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[5];
3938
f.evaluate(vals, y, c);
3939
result = (detJ*(K[0]*vals[0] + K[1]*vals[1]));
3940
values[3] = result;
3941
y[0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2];
3942
y[1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3];
3943
f.evaluate(vals, y, c);
3944
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
3945
values[4] = result;
3946
y[0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[2];
3947
y[1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[3];
3948
f.evaluate(vals, y, c);
3949
result = (-1.0)*(detJ*(K[2]*vals[0] + K[3]*vals[1]));
3950
values[5] = result;
3951
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
3952
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
3953
f.evaluate(vals, y, c);
3954
values[6] = vals[2];
3955
}
3956
3957
/// Interpolate vertex values from dof values
3958
virtual void interpolate_vertex_values(double* vertex_values,
3959
const double* dof_values,
3960
const double* vertex_coordinates,
3961
int cell_orientation,
3962
const ufc::cell& c) const
3963
{
3964
// Compute Jacobian
3965
double J[4];
3966
compute_jacobian_triangle_2d(J, vertex_coordinates);
3967
3968
3969
// Compute Jacobian inverse and determinant
3970
double K[4];
3971
double detJ;
3972
compute_jacobian_inverse_triangle_2d(K, detJ, J);
3973
3974
3975
3976
// Evaluate function and change variables
3977
vertex_values[0] = dof_values[2]*(1.0/detJ)*J[0]*2.0 + dof_values[3]*((1.0/detJ)*(J[0]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[1]*(-2.0))) + dof_values[5]*(1.0/detJ)*J[1];
3978
vertex_values[3] = dof_values[0]*(1.0/detJ)*J[0]*2.0 + dof_values[1]*((1.0/detJ)*(J[0]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[0]*(-1.0) + J[1])) + dof_values[5]*((1.0/detJ)*(J[0]*2.0 + J[1]*(-2.0)));
3979
vertex_values[6] = dof_values[0]*((1.0/detJ)*(J[1]*(-1.0))) + dof_values[1]*(1.0/detJ)*J[1]*2.0 + dof_values[2]*((1.0/detJ)*(J[0]*(-1.0) + J[1])) + dof_values[3]*((1.0/detJ)*(J[0]*2.0 + J[1]*(-2.0)));
3980
vertex_values[1] = dof_values[2]*(1.0/detJ)*J[2]*2.0 + dof_values[3]*((1.0/detJ)*(J[2]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[3]*(-2.0))) + dof_values[5]*(1.0/detJ)*J[3];
3981
vertex_values[4] = dof_values[0]*(1.0/detJ)*J[2]*2.0 + dof_values[1]*((1.0/detJ)*(J[2]*(-1.0))) + dof_values[4]*((1.0/detJ)*(J[2]*(-1.0) + J[3])) + dof_values[5]*((1.0/detJ)*(J[2]*2.0 + J[3]*(-2.0)));
3982
vertex_values[7] = dof_values[0]*((1.0/detJ)*(J[3]*(-1.0))) + dof_values[1]*(1.0/detJ)*J[3]*2.0 + dof_values[2]*((1.0/detJ)*(J[2]*(-1.0) + J[3])) + dof_values[3]*((1.0/detJ)*(J[2]*2.0 + J[3]*(-2.0)));
3983
// Evaluate function and change variables
3984
vertex_values[2] = dof_values[6];
3985
vertex_values[5] = dof_values[6];
3986
vertex_values[8] = dof_values[6];
3987
}
3988
3989
/// Map coordinate xhat from reference cell to coordinate x in cell
3990
virtual void map_from_reference_cell(double* x,
3991
const double* xhat,
3992
const ufc::cell& c) const
3993
{
3994
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
3995
}
3996
3997
/// Map from coordinate x in cell to coordinate xhat in reference cell
3998
virtual void map_to_reference_cell(double* xhat,
3999
const double* x,
4000
const ufc::cell& c) const
4001
{
4002
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
4003
}
4004
4005
/// Return the number of sub elements (for a mixed element)
4006
virtual std::size_t num_sub_elements() const
4007
{
4008
return 2;
4009
}
4010
4011
/// Create a new finite element for sub element i (for a mixed element)
4012
virtual ufc::finite_element* create_sub_element(std::size_t i) const
4013
{
4014
switch (i)
4015
{
4016
case 0:
4017
{
4018
return new mixedpoisson_finite_element_0();
4019
break;
4020
}
4021
case 1:
4022
{
4023
return new mixedpoisson_finite_element_1();
4024
break;
4025
}
4026
}
4027
4028
return 0;
4029
}
4030
4031
/// Create a new class instance
4032
virtual ufc::finite_element* create() const
4033
{
4034
return new mixedpoisson_finite_element_2();
4035
}
4036
4037
};
4038
4039
/// This class defines the interface for a local-to-global mapping of
4040
/// degrees of freedom (dofs).
4041
4042
class mixedpoisson_dofmap_0: public ufc::dofmap
4043
{
4044
public:
4045
4046
/// Constructor
4047
mixedpoisson_dofmap_0() : ufc::dofmap()
4048
{
4049
// Do nothing
4050
}
4051
4052
/// Destructor
4053
virtual ~mixedpoisson_dofmap_0()
4054
{
4055
// Do nothing
4056
}
4057
4058
/// Return a string identifying the dofmap
4059
virtual const char* signature() const
4060
{
4061
return "FFC dofmap for FiniteElement('Brezzi-Douglas-Marini', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
4062
}
4063
4064
/// Return true iff mesh entities of topological dimension d are needed
4065
virtual bool needs_mesh_entities(std::size_t d) const
4066
{
4067
switch (d)
4068
{
4069
case 0:
4070
{
4071
return false;
4072
break;
4073
}
4074
case 1:
4075
{
4076
return true;
4077
break;
4078
}
4079
case 2:
4080
{
4081
return false;
4082
break;
4083
}
4084
}
4085
4086
return false;
4087
}
4088
4089
/// Return the topological dimension of the associated cell shape
4090
virtual std::size_t topological_dimension() const
4091
{
4092
return 2;
4093
}
4094
4095
/// Return the geometric dimension of the associated cell shape
4096
virtual std::size_t geometric_dimension() const
4097
{
4098
return 2;
4099
}
4100
4101
/// Return the dimension of the global finite element function space
4102
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4103
num_global_entities) const
4104
{
4105
return 2*num_global_entities[1];
4106
}
4107
4108
/// Return the dimension of the local finite element function space for a cell
4109
virtual std::size_t local_dimension(const ufc::cell& c) const
4110
{
4111
return 6;
4112
}
4113
4114
/// Return the maximum dimension of the local finite element function space
4115
virtual std::size_t max_local_dimension() const
4116
{
4117
return 6;
4118
}
4119
4120
/// Return the number of dofs on each cell facet
4121
virtual std::size_t num_facet_dofs() const
4122
{
4123
return 2;
4124
}
4125
4126
/// Return the number of dofs associated with each cell entity of dimension d
4127
virtual std::size_t num_entity_dofs(std::size_t d) const
4128
{
4129
switch (d)
4130
{
4131
case 0:
4132
{
4133
return 0;
4134
break;
4135
}
4136
case 1:
4137
{
4138
return 2;
4139
break;
4140
}
4141
case 2:
4142
{
4143
return 0;
4144
break;
4145
}
4146
}
4147
4148
return 0;
4149
}
4150
4151
/// Tabulate the local-to-global mapping of dofs on a cell
4152
virtual void tabulate_dofs(std::size_t* dofs,
4153
const std::vector<std::size_t>& num_global_entities,
4154
const ufc::cell& c) const
4155
{
4156
dofs[0] = 2*c.entity_indices[1][0];
4157
dofs[1] = 2*c.entity_indices[1][0] + 1;
4158
dofs[2] = 2*c.entity_indices[1][1];
4159
dofs[3] = 2*c.entity_indices[1][1] + 1;
4160
dofs[4] = 2*c.entity_indices[1][2];
4161
dofs[5] = 2*c.entity_indices[1][2] + 1;
4162
}
4163
4164
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4165
virtual void tabulate_facet_dofs(std::size_t* dofs,
4166
std::size_t facet) const
4167
{
4168
switch (facet)
4169
{
4170
case 0:
4171
{
4172
dofs[0] = 0;
4173
dofs[1] = 1;
4174
break;
4175
}
4176
case 1:
4177
{
4178
dofs[0] = 2;
4179
dofs[1] = 3;
4180
break;
4181
}
4182
case 2:
4183
{
4184
dofs[0] = 4;
4185
dofs[1] = 5;
4186
break;
4187
}
4188
}
4189
4190
}
4191
4192
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4193
virtual void tabulate_entity_dofs(std::size_t* dofs,
4194
std::size_t d, std::size_t i) const
4195
{
4196
if (d > 2)
4197
{
4198
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
4199
}
4200
4201
switch (d)
4202
{
4203
case 0:
4204
{
4205
4206
break;
4207
}
4208
case 1:
4209
{
4210
if (i > 2)
4211
{
4212
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
4213
}
4214
4215
switch (i)
4216
{
4217
case 0:
4218
{
4219
dofs[0] = 0;
4220
dofs[1] = 1;
4221
break;
4222
}
4223
case 1:
4224
{
4225
dofs[0] = 2;
4226
dofs[1] = 3;
4227
break;
4228
}
4229
case 2:
4230
{
4231
dofs[0] = 4;
4232
dofs[1] = 5;
4233
break;
4234
}
4235
}
4236
4237
break;
4238
}
4239
case 2:
4240
{
4241
4242
break;
4243
}
4244
}
4245
4246
}
4247
4248
/// Tabulate the coordinates of all dofs on a cell
4249
virtual void tabulate_coordinates(double** dof_coordinates,
4250
const double* vertex_coordinates) const
4251
{
4252
dof_coordinates[0][0] = 0.66666667*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
4253
dof_coordinates[0][1] = 0.66666667*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
4254
dof_coordinates[1][0] = 0.33333333*vertex_coordinates[2] + 0.66666667*vertex_coordinates[4];
4255
dof_coordinates[1][1] = 0.33333333*vertex_coordinates[3] + 0.66666667*vertex_coordinates[5];
4256
dof_coordinates[2][0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[4];
4257
dof_coordinates[2][1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[5];
4258
dof_coordinates[3][0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[4];
4259
dof_coordinates[3][1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[5];
4260
dof_coordinates[4][0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2];
4261
dof_coordinates[4][1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3];
4262
dof_coordinates[5][0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[2];
4263
dof_coordinates[5][1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[3];
4264
}
4265
4266
/// Return the number of sub dofmaps (for a mixed element)
4267
virtual std::size_t num_sub_dofmaps() const
4268
{
4269
return 0;
4270
}
4271
4272
/// Create a new dofmap for sub dofmap i (for a mixed element)
4273
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4274
{
4275
return 0;
4276
}
4277
4278
/// Create a new class instance
4279
virtual ufc::dofmap* create() const
4280
{
4281
return new mixedpoisson_dofmap_0();
4282
}
4283
4284
};
4285
4286
/// This class defines the interface for a local-to-global mapping of
4287
/// degrees of freedom (dofs).
4288
4289
class mixedpoisson_dofmap_1: public ufc::dofmap
4290
{
4291
public:
4292
4293
/// Constructor
4294
mixedpoisson_dofmap_1() : ufc::dofmap()
4295
{
4296
// Do nothing
4297
}
4298
4299
/// Destructor
4300
virtual ~mixedpoisson_dofmap_1()
4301
{
4302
// Do nothing
4303
}
4304
4305
/// Return a string identifying the dofmap
4306
virtual const char* signature() const
4307
{
4308
return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
4309
}
4310
4311
/// Return true iff mesh entities of topological dimension d are needed
4312
virtual bool needs_mesh_entities(std::size_t d) const
4313
{
4314
switch (d)
4315
{
4316
case 0:
4317
{
4318
return false;
4319
break;
4320
}
4321
case 1:
4322
{
4323
return false;
4324
break;
4325
}
4326
case 2:
4327
{
4328
return true;
4329
break;
4330
}
4331
}
4332
4333
return false;
4334
}
4335
4336
/// Return the topological dimension of the associated cell shape
4337
virtual std::size_t topological_dimension() const
4338
{
4339
return 2;
4340
}
4341
4342
/// Return the geometric dimension of the associated cell shape
4343
virtual std::size_t geometric_dimension() const
4344
{
4345
return 2;
4346
}
4347
4348
/// Return the dimension of the global finite element function space
4349
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4350
num_global_entities) const
4351
{
4352
return num_global_entities[2];
4353
}
4354
4355
/// Return the dimension of the local finite element function space for a cell
4356
virtual std::size_t local_dimension(const ufc::cell& c) const
4357
{
4358
return 1;
4359
}
4360
4361
/// Return the maximum dimension of the local finite element function space
4362
virtual std::size_t max_local_dimension() const
4363
{
4364
return 1;
4365
}
4366
4367
/// Return the number of dofs on each cell facet
4368
virtual std::size_t num_facet_dofs() const
4369
{
4370
return 0;
4371
}
4372
4373
/// Return the number of dofs associated with each cell entity of dimension d
4374
virtual std::size_t num_entity_dofs(std::size_t d) const
4375
{
4376
switch (d)
4377
{
4378
case 0:
4379
{
4380
return 0;
4381
break;
4382
}
4383
case 1:
4384
{
4385
return 0;
4386
break;
4387
}
4388
case 2:
4389
{
4390
return 1;
4391
break;
4392
}
4393
}
4394
4395
return 0;
4396
}
4397
4398
/// Tabulate the local-to-global mapping of dofs on a cell
4399
virtual void tabulate_dofs(std::size_t* dofs,
4400
const std::vector<std::size_t>& num_global_entities,
4401
const ufc::cell& c) const
4402
{
4403
dofs[0] = c.entity_indices[2][0];
4404
}
4405
4406
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4407
virtual void tabulate_facet_dofs(std::size_t* dofs,
4408
std::size_t facet) const
4409
{
4410
switch (facet)
4411
{
4412
case 0:
4413
{
4414
4415
break;
4416
}
4417
case 1:
4418
{
4419
4420
break;
4421
}
4422
case 2:
4423
{
4424
4425
break;
4426
}
4427
}
4428
4429
}
4430
4431
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4432
virtual void tabulate_entity_dofs(std::size_t* dofs,
4433
std::size_t d, std::size_t i) const
4434
{
4435
if (d > 2)
4436
{
4437
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
4438
}
4439
4440
switch (d)
4441
{
4442
case 0:
4443
{
4444
4445
break;
4446
}
4447
case 1:
4448
{
4449
4450
break;
4451
}
4452
case 2:
4453
{
4454
if (i > 0)
4455
{
4456
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
4457
}
4458
4459
dofs[0] = 0;
4460
break;
4461
}
4462
}
4463
4464
}
4465
4466
/// Tabulate the coordinates of all dofs on a cell
4467
virtual void tabulate_coordinates(double** dof_coordinates,
4468
const double* vertex_coordinates) const
4469
{
4470
dof_coordinates[0][0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
4471
dof_coordinates[0][1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
4472
}
4473
4474
/// Return the number of sub dofmaps (for a mixed element)
4475
virtual std::size_t num_sub_dofmaps() const
4476
{
4477
return 0;
4478
}
4479
4480
/// Create a new dofmap for sub dofmap i (for a mixed element)
4481
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4482
{
4483
return 0;
4484
}
4485
4486
/// Create a new class instance
4487
virtual ufc::dofmap* create() const
4488
{
4489
return new mixedpoisson_dofmap_1();
4490
}
4491
4492
};
4493
4494
/// This class defines the interface for a local-to-global mapping of
4495
/// degrees of freedom (dofs).
4496
4497
class mixedpoisson_dofmap_2: public ufc::dofmap
4498
{
4499
public:
4500
4501
/// Constructor
4502
mixedpoisson_dofmap_2() : ufc::dofmap()
4503
{
4504
// Do nothing
4505
}
4506
4507
/// Destructor
4508
virtual ~mixedpoisson_dofmap_2()
4509
{
4510
// Do nothing
4511
}
4512
4513
/// Return a string identifying the dofmap
4514
virtual const char* signature() const
4515
{
4516
return "FFC dofmap for MixedElement(*[FiniteElement('Brezzi-Douglas-Marini', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None), FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)], **{'value_shape': (3,) })";
4517
}
4518
4519
/// Return true iff mesh entities of topological dimension d are needed
4520
virtual bool needs_mesh_entities(std::size_t d) const
4521
{
4522
switch (d)
4523
{
4524
case 0:
4525
{
4526
return false;
4527
break;
4528
}
4529
case 1:
4530
{
4531
return true;
4532
break;
4533
}
4534
case 2:
4535
{
4536
return true;
4537
break;
4538
}
4539
}
4540
4541
return false;
4542
}
4543
4544
/// Return the topological dimension of the associated cell shape
4545
virtual std::size_t topological_dimension() const
4546
{
4547
return 2;
4548
}
4549
4550
/// Return the geometric dimension of the associated cell shape
4551
virtual std::size_t geometric_dimension() const
4552
{
4553
return 2;
4554
}
4555
4556
/// Return the dimension of the global finite element function space
4557
virtual std::size_t global_dimension(const std::vector<std::size_t>&
4558
num_global_entities) const
4559
{
4560
return 2*num_global_entities[1] + num_global_entities[2];
4561
}
4562
4563
/// Return the dimension of the local finite element function space for a cell
4564
virtual std::size_t local_dimension(const ufc::cell& c) const
4565
{
4566
return 7;
4567
}
4568
4569
/// Return the maximum dimension of the local finite element function space
4570
virtual std::size_t max_local_dimension() const
4571
{
4572
return 7;
4573
}
4574
4575
/// Return the number of dofs on each cell facet
4576
virtual std::size_t num_facet_dofs() const
4577
{
4578
return 2;
4579
}
4580
4581
/// Return the number of dofs associated with each cell entity of dimension d
4582
virtual std::size_t num_entity_dofs(std::size_t d) const
4583
{
4584
switch (d)
4585
{
4586
case 0:
4587
{
4588
return 0;
4589
break;
4590
}
4591
case 1:
4592
{
4593
return 2;
4594
break;
4595
}
4596
case 2:
4597
{
4598
return 1;
4599
break;
4600
}
4601
}
4602
4603
return 0;
4604
}
4605
4606
/// Tabulate the local-to-global mapping of dofs on a cell
4607
virtual void tabulate_dofs(std::size_t* dofs,
4608
const std::vector<std::size_t>& num_global_entities,
4609
const ufc::cell& c) const
4610
{
4611
unsigned int offset = 0;
4612
dofs[0] = offset + 2*c.entity_indices[1][0];
4613
dofs[1] = offset + 2*c.entity_indices[1][0] + 1;
4614
dofs[2] = offset + 2*c.entity_indices[1][1];
4615
dofs[3] = offset + 2*c.entity_indices[1][1] + 1;
4616
dofs[4] = offset + 2*c.entity_indices[1][2];
4617
dofs[5] = offset + 2*c.entity_indices[1][2] + 1;
4618
offset += 2*num_global_entities[1];
4619
dofs[6] = offset + c.entity_indices[2][0];
4620
offset += num_global_entities[2];
4621
}
4622
4623
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
4624
virtual void tabulate_facet_dofs(std::size_t* dofs,
4625
std::size_t facet) const
4626
{
4627
switch (facet)
4628
{
4629
case 0:
4630
{
4631
dofs[0] = 0;
4632
dofs[1] = 1;
4633
break;
4634
}
4635
case 1:
4636
{
4637
dofs[0] = 2;
4638
dofs[1] = 3;
4639
break;
4640
}
4641
case 2:
4642
{
4643
dofs[0] = 4;
4644
dofs[1] = 5;
4645
break;
4646
}
4647
}
4648
4649
}
4650
4651
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
4652
virtual void tabulate_entity_dofs(std::size_t* dofs,
4653
std::size_t d, std::size_t i) const
4654
{
4655
if (d > 2)
4656
{
4657
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
4658
}
4659
4660
switch (d)
4661
{
4662
case 0:
4663
{
4664
4665
break;
4666
}
4667
case 1:
4668
{
4669
if (i > 2)
4670
{
4671
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
4672
}
4673
4674
switch (i)
4675
{
4676
case 0:
4677
{
4678
dofs[0] = 0;
4679
dofs[1] = 1;
4680
break;
4681
}
4682
case 1:
4683
{
4684
dofs[0] = 2;
4685
dofs[1] = 3;
4686
break;
4687
}
4688
case 2:
4689
{
4690
dofs[0] = 4;
4691
dofs[1] = 5;
4692
break;
4693
}
4694
}
4695
4696
break;
4697
}
4698
case 2:
4699
{
4700
if (i > 0)
4701
{
4702
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
4703
}
4704
4705
dofs[0] = 6;
4706
break;
4707
}
4708
}
4709
4710
}
4711
4712
/// Tabulate the coordinates of all dofs on a cell
4713
virtual void tabulate_coordinates(double** dof_coordinates,
4714
const double* vertex_coordinates) const
4715
{
4716
dof_coordinates[0][0] = 0.66666667*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
4717
dof_coordinates[0][1] = 0.66666667*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
4718
dof_coordinates[1][0] = 0.33333333*vertex_coordinates[2] + 0.66666667*vertex_coordinates[4];
4719
dof_coordinates[1][1] = 0.33333333*vertex_coordinates[3] + 0.66666667*vertex_coordinates[5];
4720
dof_coordinates[2][0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[4];
4721
dof_coordinates[2][1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[5];
4722
dof_coordinates[3][0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[4];
4723
dof_coordinates[3][1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[5];
4724
dof_coordinates[4][0] = 0.66666667*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2];
4725
dof_coordinates[4][1] = 0.66666667*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3];
4726
dof_coordinates[5][0] = 0.33333333*vertex_coordinates[0] + 0.66666667*vertex_coordinates[2];
4727
dof_coordinates[5][1] = 0.33333333*vertex_coordinates[1] + 0.66666667*vertex_coordinates[3];
4728
dof_coordinates[6][0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[4];
4729
dof_coordinates[6][1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[5];
4730
}
4731
4732
/// Return the number of sub dofmaps (for a mixed element)
4733
virtual std::size_t num_sub_dofmaps() const
4734
{
4735
return 2;
4736
}
4737
4738
/// Create a new dofmap for sub dofmap i (for a mixed element)
4739
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
4740
{
4741
switch (i)
4742
{
4743
case 0:
4744
{
4745
return new mixedpoisson_dofmap_0();
4746
break;
4747
}
4748
case 1:
4749
{
4750
return new mixedpoisson_dofmap_1();
4751
break;
4752
}
4753
}
4754
4755
return 0;
4756
}
4757
4758
/// Create a new class instance
4759
virtual ufc::dofmap* create() const
4760
{
4761
return new mixedpoisson_dofmap_2();
4762
}
4763
4764
};
4765
4766
/// This class defines the interface for the tabulation of the cell
4767
/// tensor corresponding to the local contribution to a form from
4768
/// the integral over a cell.
4769
4770
class mixedpoisson_cell_integral_0_otherwise: public ufc::cell_integral
4771
{
4772
public:
4773
4774
/// Constructor
4775
mixedpoisson_cell_integral_0_otherwise() : ufc::cell_integral()
4776
{
4777
// Do nothing
4778
}
4779
4780
/// Destructor
4781
virtual ~mixedpoisson_cell_integral_0_otherwise()
4782
{
4783
// Do nothing
4784
}
4785
4786
/// Tabulate the tensor for the contribution from a local cell
4787
virtual void tabulate_tensor(double* A,
4788
const double * const * w,
4789
const double* vertex_coordinates,
4790
int cell_orientation) const
4791
{
4792
// Compute Jacobian
4793
double J[4];
4794
compute_jacobian_triangle_2d(J, vertex_coordinates);
4795
4796
// Compute Jacobian inverse and determinant
4797
double K[4];
4798
double detJ;
4799
compute_jacobian_inverse_triangle_2d(K, detJ, J);
4800
4801
// Set scale factor
4802
const double det = std::abs(detJ);
4803
4804
// Cell volume
4805
4806
// Compute circumradius of triangle in 2D
4807
4808
4809
// Facet area
4810
4811
// Array of quadrature weights.
4812
static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
4813
// Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
4814
4815
// Value of basis functions at quadrature points.
4816
static const double FE0_C0[3][7] = \
4817
{{0.33333333, -0.16666667, 1.1666667, -0.33333333, -0.16666667, 0.33333333, 0.0},
4818
{0.33333333, -0.16666667, -0.33333333, 1.1666667, -0.16666667, 0.33333333, 0.0},
4819
{1.3333333, -0.66666667, 0.16666667, 0.16666667, -0.66666667, 1.3333333, 0.0}};
4820
4821
static const double FE0_C0_D01[3][7] = \
4822
{{0.0, 0.0, -3.0, 3.0, 0.0, 0.0, 0.0},
4823
{0.0, 0.0, -3.0, 3.0, 0.0, 0.0, 0.0},
4824
{0.0, 0.0, -3.0, 3.0, 0.0, 0.0, 0.0}};
4825
4826
static const double FE0_C0_D10[3][7] = \
4827
{{2.0, -1.0, -2.0, 1.0, -1.0, 2.0, 0.0},
4828
{2.0, -1.0, -2.0, 1.0, -1.0, 2.0, 0.0},
4829
{2.0, -1.0, -2.0, 1.0, -1.0, 2.0, 0.0}};
4830
4831
static const double FE0_C1[3][7] = \
4832
{{-0.16666667, 0.33333333, 0.16666667, -0.33333333, -1.1666667, 0.33333333, 0.0},
4833
{-0.66666667, 1.3333333, 0.66666667, -1.3333333, -0.16666667, -0.16666667, 0.0},
4834
{-0.16666667, 0.33333333, 0.16666667, -0.33333333, 0.33333333, -1.1666667, 0.0}};
4835
4836
static const double FE0_C1_D01[3][7] = \
4837
{{-1.0, 2.0, 1.0, -2.0, 2.0, -1.0, 0.0},
4838
{-1.0, 2.0, 1.0, -2.0, 2.0, -1.0, 0.0},
4839
{-1.0, 2.0, 1.0, -2.0, 2.0, -1.0, 0.0}};
4840
4841
static const double FE0_C1_D10[3][7] = \
4842
{{0.0, 0.0, 0.0, 0.0, 3.0, -3.0, 0.0},
4843
{0.0, 0.0, 0.0, 0.0, 3.0, -3.0, 0.0},
4844
{0.0, 0.0, 0.0, 0.0, 3.0, -3.0, 0.0}};
4845
4846
static const double FE0_C2[3][7] = \
4847
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0},
4848
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0},
4849
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
4850
4851
// Reset values in the element tensor.
4852
for (unsigned int r = 0; r < 49; r++)
4853
{
4854
A[r] = 0.0;
4855
}// end loop over 'r'
4856
4857
// Compute element tensor using UFL quadrature representation
4858
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
4859
4860
// Loop quadrature points for integral.
4861
// Number of operations to compute element tensor for following IP loop = 13671
4862
for (unsigned int ip = 0; ip < 3; ip++)
4863
{
4864
4865
// Number of operations for primary indices: 4557
4866
for (unsigned int j = 0; j < 7; j++)
4867
{
4868
for (unsigned int k = 0; k < 7; k++)
4869
{
4870
// Number of operations to compute entry: 93
4871
A[j*7 + k] += ((((((K[1]*(1.0/detJ)*J[2]*FE0_C0_D10[ip][j] + K[1]*(1.0/detJ)*J[3]*FE0_C1_D10[ip][j] + K[3]*(1.0/detJ)*J[2]*FE0_C0_D01[ip][j] + K[3]*(1.0/detJ)*J[3]*FE0_C1_D01[ip][j]) + (K[0]*(1.0/detJ)*J[0]*FE0_C0_D10[ip][j] + K[0]*(1.0/detJ)*J[1]*FE0_C1_D10[ip][j] + K[2]*(1.0/detJ)*J[0]*FE0_C0_D01[ip][j] + K[2]*(1.0/detJ)*J[1]*FE0_C1_D01[ip][j])))*FE0_C2[ip][k])*(-1.0) + ((((1.0/detJ)*J[0]*FE0_C0[ip][j] + (1.0/detJ)*J[1]*FE0_C1[ip][j]))*(((1.0/detJ)*J[0]*FE0_C0[ip][k] + (1.0/detJ)*J[1]*FE0_C1[ip][k])) + (((1.0/detJ)*J[2]*FE0_C0[ip][j] + (1.0/detJ)*J[3]*FE0_C1[ip][j]))*(((1.0/detJ)*J[2]*FE0_C0[ip][k] + (1.0/detJ)*J[3]*FE0_C1[ip][k])))) + (((K[1]*(1.0/detJ)*J[2]*FE0_C0_D10[ip][k] + K[1]*(1.0/detJ)*J[3]*FE0_C1_D10[ip][k] + K[3]*(1.0/detJ)*J[2]*FE0_C0_D01[ip][k] + K[3]*(1.0/detJ)*J[3]*FE0_C1_D01[ip][k]) + (K[0]*(1.0/detJ)*J[0]*FE0_C0_D10[ip][k] + K[0]*(1.0/detJ)*J[1]*FE0_C1_D10[ip][k] + K[2]*(1.0/detJ)*J[0]*FE0_C0_D01[ip][k] + K[2]*(1.0/detJ)*J[1]*FE0_C1_D01[ip][k])))*FE0_C2[ip][j])*W3[ip]*det;
4872
}// end loop over 'k'
4873
}// end loop over 'j'
4874
}// end loop over 'ip'
4875
}
4876
4877
};
4878
4879
/// This class defines the interface for the tabulation of the cell
4880
/// tensor corresponding to the local contribution to a form from
4881
/// the integral over a cell.
4882
4883
class mixedpoisson_cell_integral_1_otherwise: public ufc::cell_integral
4884
{
4885
public:
4886
4887
/// Constructor
4888
mixedpoisson_cell_integral_1_otherwise() : ufc::cell_integral()
4889
{
4890
// Do nothing
4891
}
4892
4893
/// Destructor
4894
virtual ~mixedpoisson_cell_integral_1_otherwise()
4895
{
4896
// Do nothing
4897
}
4898
4899
/// Tabulate the tensor for the contribution from a local cell
4900
virtual void tabulate_tensor(double* A,
4901
const double * const * w,
4902
const double* vertex_coordinates,
4903
int cell_orientation) const
4904
{
4905
// Compute Jacobian
4906
double J[4];
4907
compute_jacobian_triangle_2d(J, vertex_coordinates);
4908
4909
// Compute Jacobian inverse and determinant
4910
double K[4];
4911
double detJ;
4912
compute_jacobian_inverse_triangle_2d(K, detJ, J);
4913
4914
// Set scale factor
4915
const double det = std::abs(detJ);
4916
4917
// Cell volume
4918
4919
// Compute circumradius of triangle in 2D
4920
4921
4922
// Facet area
4923
4924
// Array of quadrature weights.
4925
static const double W1 = 0.5;
4926
// Quadrature points on the UFC reference element: (0.33333333, 0.33333333)
4927
4928
// Value of basis functions at quadrature points.
4929
static const double FE0[1][1] = \
4930
{{1.0}};
4931
4932
static const double FE1_C2[1][7] = \
4933
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
4934
4935
// Reset values in the element tensor.
4936
for (unsigned int r = 0; r < 7; r++)
4937
{
4938
A[r] = 0.0;
4939
}// end loop over 'r'
4940
4941
// Compute element tensor using UFL quadrature representation
4942
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
4943
4944
// Loop quadrature points for integral.
4945
// Number of operations to compute element tensor for following IP loop = 30
4946
// Only 1 integration point, omitting IP loop.
4947
4948
// Coefficient declarations.
4949
double F0 = 0.0;
4950
4951
// Total number of operations to compute function values = 2
4952
for (unsigned int r = 0; r < 1; r++)
4953
{
4954
F0 += FE0[0][0]*w[0][0];
4955
}// end loop over 'r'
4956
4957
// Number of operations for primary indices: 28
4958
for (unsigned int j = 0; j < 7; j++)
4959
{
4960
// Number of operations to compute entry: 4
4961
A[j] += FE1_C2[0][j]*F0*W1*det;
4962
}// end loop over 'j'
4963
}
4964
4965
};
4966
4967
/// This class defines the interface for the assembly of the global
4968
/// tensor corresponding to a form with r + n arguments, that is, a
4969
/// mapping
4970
///
4971
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
4972
///
4973
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
4974
/// global tensor A is defined by
4975
///
4976
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
4977
///
4978
/// where each argument Vj represents the application to the
4979
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
4980
/// fixed functions (coefficients).
4981
4982
class mixedpoisson_form_0: public ufc::form
4983
{
4984
public:
4985
4986
/// Constructor
4987
mixedpoisson_form_0() : ufc::form()
4988
{
4989
// Do nothing
4990
}
4991
4992
/// Destructor
4993
virtual ~mixedpoisson_form_0()
4994
{
4995
// Do nothing
4996
}
4997
4998
/// Return a string identifying the form
4999
virtual const char* signature() const
5000
{
5001
return "ba8f7ee8a474704904426d1ae24218877c4fffaff40757e3b3504cd391e7ccb869536395afe73d8e955900506f62cec5a6bb0981dfdd9f562ab4ca394a2cac97";
5002
}
5003
5004
/// Return the rank of the global tensor (r)
5005
virtual std::size_t rank() const
5006
{
5007
return 2;
5008
}
5009
5010
/// Return the number of coefficients (n)
5011
virtual std::size_t num_coefficients() const
5012
{
5013
return 0;
5014
}
5015
5016
/// Return the number of cell domains
5017
virtual std::size_t num_cell_domains() const
5018
{
5019
return 0;
5020
}
5021
5022
/// Return the number of exterior facet domains
5023
virtual std::size_t num_exterior_facet_domains() const
5024
{
5025
return 0;
5026
}
5027
5028
/// Return the number of interior facet domains
5029
virtual std::size_t num_interior_facet_domains() const
5030
{
5031
return 0;
5032
}
5033
5034
/// Return the number of point domains
5035
virtual std::size_t num_point_domains() const
5036
{
5037
return 0;
5038
}
5039
5040
/// Return whether the form has any cell integrals
5041
virtual bool has_cell_integrals() const
5042
{
5043
return true;
5044
}
5045
5046
/// Return whether the form has any exterior facet integrals
5047
virtual bool has_exterior_facet_integrals() const
5048
{
5049
return false;
5050
}
5051
5052
/// Return whether the form has any interior facet integrals
5053
virtual bool has_interior_facet_integrals() const
5054
{
5055
return false;
5056
}
5057
5058
/// Return whether the form has any point integrals
5059
virtual bool has_point_integrals() const
5060
{
5061
return false;
5062
}
5063
5064
/// Create a new finite element for argument function i
5065
virtual ufc::finite_element* create_finite_element(std::size_t i) const
5066
{
5067
switch (i)
5068
{
5069
case 0:
5070
{
5071
return new mixedpoisson_finite_element_2();
5072
break;
5073
}
5074
case 1:
5075
{
5076
return new mixedpoisson_finite_element_2();
5077
break;
5078
}
5079
}
5080
5081
return 0;
5082
}
5083
5084
/// Create a new dofmap for argument function i
5085
virtual ufc::dofmap* create_dofmap(std::size_t i) const
5086
{
5087
switch (i)
5088
{
5089
case 0:
5090
{
5091
return new mixedpoisson_dofmap_2();
5092
break;
5093
}
5094
case 1:
5095
{
5096
return new mixedpoisson_dofmap_2();
5097
break;
5098
}
5099
}
5100
5101
return 0;
5102
}
5103
5104
/// Create a new cell integral on sub domain i
5105
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
5106
{
5107
return 0;
5108
}
5109
5110
/// Create a new exterior facet integral on sub domain i
5111
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
5112
{
5113
return 0;
5114
}
5115
5116
/// Create a new interior facet integral on sub domain i
5117
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
5118
{
5119
return 0;
5120
}
5121
5122
/// Create a new point integral on sub domain i
5123
virtual ufc::point_integral* create_point_integral(std::size_t i) const
5124
{
5125
return 0;
5126
}
5127
5128
/// Create a new cell integral on everywhere else
5129
virtual ufc::cell_integral* create_default_cell_integral() const
5130
{
5131
return new mixedpoisson_cell_integral_0_otherwise();
5132
}
5133
5134
/// Create a new exterior facet integral on everywhere else
5135
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
5136
{
5137
return 0;
5138
}
5139
5140
/// Create a new interior facet integral on everywhere else
5141
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
5142
{
5143
return 0;
5144
}
5145
5146
/// Create a new point integral on everywhere else
5147
virtual ufc::point_integral* create_default_point_integral() const
5148
{
5149
return 0;
5150
}
5151
5152
};
5153
5154
/// This class defines the interface for the assembly of the global
5155
/// tensor corresponding to a form with r + n arguments, that is, a
5156
/// mapping
5157
///
5158
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
5159
///
5160
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
5161
/// global tensor A is defined by
5162
///
5163
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
5164
///
5165
/// where each argument Vj represents the application to the
5166
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
5167
/// fixed functions (coefficients).
5168
5169
class mixedpoisson_form_1: public ufc::form
5170
{
5171
public:
5172
5173
/// Constructor
5174
mixedpoisson_form_1() : ufc::form()
5175
{
5176
// Do nothing
5177
}
5178
5179
/// Destructor
5180
virtual ~mixedpoisson_form_1()
5181
{
5182
// Do nothing
5183
}
5184
5185
/// Return a string identifying the form
5186
virtual const char* signature() const
5187
{
5188
return "ba4bf21fb268523b07e03e4afd5337c82cf9a48e5b60f077e619b7de28c804d51cb9e5832eab0f3ec048a42b64a007c831e68416fd4e56b0f2a4f8a4b5a863ca";
5189
}
5190
5191
/// Return the rank of the global tensor (r)
5192
virtual std::size_t rank() const
5193
{
5194
return 1;
5195
}
5196
5197
/// Return the number of coefficients (n)
5198
virtual std::size_t num_coefficients() const
5199
{
5200
return 1;
5201
}
5202
5203
/// Return the number of cell domains
5204
virtual std::size_t num_cell_domains() const
5205
{
5206
return 0;
5207
}
5208
5209
/// Return the number of exterior facet domains
5210
virtual std::size_t num_exterior_facet_domains() const
5211
{
5212
return 0;
5213
}
5214
5215
/// Return the number of interior facet domains
5216
virtual std::size_t num_interior_facet_domains() const
5217
{
5218
return 0;
5219
}
5220
5221
/// Return the number of point domains
5222
virtual std::size_t num_point_domains() const
5223
{
5224
return 0;
5225
}
5226
5227
/// Return whether the form has any cell integrals
5228
virtual bool has_cell_integrals() const
5229
{
5230
return true;
5231
}
5232
5233
/// Return whether the form has any exterior facet integrals
5234
virtual bool has_exterior_facet_integrals() const
5235
{
5236
return false;
5237
}
5238
5239
/// Return whether the form has any interior facet integrals
5240
virtual bool has_interior_facet_integrals() const
5241
{
5242
return false;
5243
}
5244
5245
/// Return whether the form has any point integrals
5246
virtual bool has_point_integrals() const
5247
{
5248
return false;
5249
}
5250
5251
/// Create a new finite element for argument function i
5252
virtual ufc::finite_element* create_finite_element(std::size_t i) const
5253
{
5254
switch (i)
5255
{
5256
case 0:
5257
{
5258
return new mixedpoisson_finite_element_2();
5259
break;
5260
}
5261
case 1:
5262
{
5263
return new mixedpoisson_finite_element_1();
5264
break;
5265
}
5266
}
5267
5268
return 0;
5269
}
5270
5271
/// Create a new dofmap for argument function i
5272
virtual ufc::dofmap* create_dofmap(std::size_t i) const
5273
{
5274
switch (i)
5275
{
5276
case 0:
5277
{
5278
return new mixedpoisson_dofmap_2();
5279
break;
5280
}
5281
case 1:
5282
{
5283
return new mixedpoisson_dofmap_1();
5284
break;
5285
}
5286
}
5287
5288
return 0;
5289
}
5290
5291
/// Create a new cell integral on sub domain i
5292
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
5293
{
5294
return 0;
5295
}
5296
5297
/// Create a new exterior facet integral on sub domain i
5298
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
5299
{
5300
return 0;
5301
}
5302
5303
/// Create a new interior facet integral on sub domain i
5304
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
5305
{
5306
return 0;
5307
}
5308
5309
/// Create a new point integral on sub domain i
5310
virtual ufc::point_integral* create_point_integral(std::size_t i) const
5311
{
5312
return 0;
5313
}
5314
5315
/// Create a new cell integral on everywhere else
5316
virtual ufc::cell_integral* create_default_cell_integral() const
5317
{
5318
return new mixedpoisson_cell_integral_1_otherwise();
5319
}
5320
5321
/// Create a new exterior facet integral on everywhere else
5322
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
5323
{
5324
return 0;
5325
}
5326
5327
/// Create a new interior facet integral on everywhere else
5328
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
5329
{
5330
return 0;
5331
}
5332
5333
/// Create a new point integral on everywhere else
5334
virtual ufc::point_integral* create_default_point_integral() const
5335
{
5336
return 0;
5337
}
5338
5339
};
5340
5341
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1