1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
// Copyright (c) 2020 Chris Richardson
// FEniCS Project
// SPDX-License-Identifier: MIT
#include "finite-element.h"
#include "brezzi-douglas-marini.h"
#include "crouzeix-raviart.h"
#include "lagrange.h"
#include "nedelec.h"
#include "polyset.h"
#include "raviart-thomas.h"
#include "regge.h"
#include <iostream>
#include <numeric>
#define str_macro(X) #X
#define str(X) str_macro(X)
using namespace basix;
//-----------------------------------------------------------------------------
basix::FiniteElement basix::create_element(std::string family,
std::string cell, int degree)
{
if (family == "Lagrange" or family == "P" or family == "Q")
return create_lagrange(cell::str_to_type(cell), degree, family);
else if (family == "Discontinuous Lagrange")
return create_dlagrange(cell::str_to_type(cell), degree, family);
else if (family == "Brezzi-Douglas-Marini")
return create_bdm(cell::str_to_type(cell), degree, family);
else if (family == "Raviart-Thomas")
return create_rt(cell::str_to_type(cell), degree, family);
else if (family == "Nedelec 1st kind H(curl)")
return create_nedelec(cell::str_to_type(cell), degree, family);
else if (family == "Nedelec 2nd kind H(curl)")
return create_nedelec2(cell::str_to_type(cell), degree, family);
else if (family == "Regge")
return create_regge(cell::str_to_type(cell), degree, family);
else if (family == "Crouzeix-Raviart")
return cr::create(cell::str_to_type(cell), degree);
else
throw std::runtime_error("Family not found: \"" + family + "\"");
}
//-----------------------------------------------------------------------------
Eigen::MatrixXd
basix::compute_expansion_coefficients(const Eigen::MatrixXd& coeffs,
const Eigen::MatrixXd& dual,
bool condition_check)
{
#ifndef NDEBUG
std::cout << "Initial coeffs = \n[" << coeffs << "]\n";
std::cout << "Dual matrix = \n[" << dual << "]\n";
#endif
const Eigen::MatrixXd A = coeffs * dual.transpose();
if (condition_check)
{
Eigen::JacobiSVD svd(A);
const int size = svd.singularValues().size();
const double kappa
= svd.singularValues()(0) / svd.singularValues()(size - 1);
if (kappa > 1e6)
{
throw std::runtime_error("Poorly conditioned B.D^T when computing "
"expansion coefficients");
}
}
Eigen::MatrixXd new_coeffs = A.colPivHouseholderQr().solve(coeffs);
#ifndef NDEBUG
std::cout << "New coeffs = \n[" << new_coeffs << "]\n";
#endif
return new_coeffs;
}
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
FiniteElement::FiniteElement(
std::string name, cell::type cell_type, int degree,
const std::vector<int>& value_shape, const Eigen::ArrayXXd& coeffs,
const std::vector<std::vector<int>>& entity_dofs,
const std::vector<Eigen::MatrixXd>& base_permutations,
const Eigen::ArrayXXd& points, const Eigen::MatrixXd interpolation_matrix,
const std::string mapping_name)
: _cell_type(cell_type), _family_name(name), _degree(degree),
_value_shape(value_shape), _mapping_name(mapping_name), _coeffs(coeffs),
_entity_dofs(entity_dofs), _base_permutations(base_permutations),
_points(points), _interpolation_matrix(interpolation_matrix)
{
// Check that entity dofs add up to total number of dofs
int sum = 0;
for (const std::vector<int>& q : entity_dofs)
sum = std::accumulate(q.begin(), q.end(), sum);
if (sum != _coeffs.rows())
{
throw std::runtime_error(
"Number of entity dofs does not match total number of dofs");
}
}
//-----------------------------------------------------------------------------
cell::type FiniteElement::cell_type() const { return _cell_type; }
//-----------------------------------------------------------------------------
int FiniteElement::degree() const { return _degree; }
//-----------------------------------------------------------------------------
int FiniteElement::value_size() const
{
int value_size = 1;
for (const int& d : _value_shape)
value_size *= d;
return value_size;
}
//-----------------------------------------------------------------------------
const std::vector<int>& FiniteElement::value_shape() const
{
return _value_shape;
}
//-----------------------------------------------------------------------------
int FiniteElement::dim() const { return _coeffs.rows(); }
//-----------------------------------------------------------------------------
const std::string& FiniteElement::family_name() const { return _family_name; }
//-----------------------------------------------------------------------------
const std::string& FiniteElement::mapping_name() const { return _mapping_name; }
//-----------------------------------------------------------------------------
const Eigen::MatrixXd& FiniteElement::interpolation_matrix() const
{
return _interpolation_matrix;
}
//-----------------------------------------------------------------------------
const std::vector<std::vector<int>>& FiniteElement::entity_dofs() const
{
return _entity_dofs;
}
//-----------------------------------------------------------------------------
std::vector<Eigen::ArrayXXd>
FiniteElement::tabulate(int nd, const Eigen::ArrayXXd& x) const
{
const int tdim = cell::topological_dimension(_cell_type);
if (x.cols() != tdim)
throw std::runtime_error("Point dim does not match element dim.");
std::vector<Eigen::ArrayXXd> basis
= polyset::tabulate(_cell_type, _degree, nd, x);
const int psize = polyset::dim(_cell_type, _degree);
const int ndofs = _coeffs.rows();
const int vs = value_size();
std::vector<Eigen::ArrayXXd> dresult(basis.size(),
Eigen::ArrayXXd(x.rows(), ndofs * vs));
for (std::size_t p = 0; p < dresult.size(); ++p)
{
for (int j = 0; j < vs; ++j)
{
dresult[p].block(0, ndofs * j, x.rows(), ndofs)
= basis[p].matrix()
* _coeffs.block(0, psize * j, _coeffs.rows(), psize).transpose();
}
}
return dresult;
}
//-----------------------------------------------------------------------------
std::vector<Eigen::MatrixXd> FiniteElement::base_permutations() const
{
return _base_permutations;
}
//-----------------------------------------------------------------------------
const Eigen::ArrayXXd& FiniteElement::points() const { return _points; }
//-----------------------------------------------------------------------------
const std::string& basix::version()
{
static const std::string version_str = str(BASIX_VERSION);
return version_str;
}
//-----------------------------------------------------------------------------
|
