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
|
// Copyright (c) 2020 Chris Richardson
// FEniCS Project
// SPDX-License-Identifier: MIT
#include "brezzi-douglas-marini.h"
#include "dof-permutations.h"
#include "lagrange.h"
#include "moments.h"
#include "nedelec.h"
#include "polyset.h"
#include "quadrature.h"
#include <Eigen/Dense>
#include <numeric>
#include <vector>
using namespace basix;
//----------------------------------------------------------------------------
FiniteElement basix::create_bdm(cell::type celltype, int degree,
const std::string& name)
{
if (celltype != cell::type::triangle and celltype != cell::type::tetrahedron)
throw std::runtime_error("Unsupported cell type");
const int tdim = cell::topological_dimension(celltype);
const cell::type facettype
= (tdim == 2) ? cell::type::interval : cell::type::triangle;
// The number of order (degree) scalar polynomials
const int npoly = polyset::dim(celltype, degree);
const int ndofs = npoly * tdim;
// Create coefficients for order (degree-1) vector polynomials
Eigen::MatrixXd wcoeffs = Eigen::MatrixXd::Identity(ndofs, ndofs);
// Dual space
Eigen::MatrixXd dual = Eigen::MatrixXd::Zero(ndofs, ndofs);
// quadrature degree
int quad_deg = 5 * degree;
// Add rows to dualmat for integral moments on facets
const int facet_count = tdim + 1;
const int facet_dofs = polyset::dim(facettype, degree);
dual.block(0, 0, facet_count * facet_dofs, ndofs)
= moments::make_normal_integral_moments(
create_dlagrange(facettype, degree), celltype, tdim, degree,
quad_deg);
const int internal_dofs = ndofs - facet_count * facet_dofs;
// Add rows to dualmat for integral moments on interior
if (degree > 1)
{
// Interior integral moment
dual.block(facet_count * facet_dofs, 0, internal_dofs, ndofs)
= moments::make_dot_integral_moments(
create_nedelec(celltype, degree - 1), celltype, tdim, degree,
quad_deg);
}
const std::vector<std::vector<std::vector<int>>> topology
= cell::topology(celltype);
int perm_count = 0;
for (int i = 1; i < tdim; ++i)
perm_count += topology[i].size() * i;
std::vector<Eigen::MatrixXd> base_permutations(
perm_count, Eigen::MatrixXd::Identity(ndofs, ndofs));
if (tdim == 2)
{
Eigen::ArrayXi edge_ref = dofperms::interval_reflection(degree + 1);
Eigen::ArrayXXd edge_dir
= dofperms::interval_reflection_tangent_directions(degree + 1);
for (int edge = 0; edge < facet_count; ++edge)
{
const int start = edge_ref.size() * edge;
for (int i = 0; i < edge_ref.size(); ++i)
{
base_permutations[edge](start + i, start + i) = 0;
base_permutations[edge](start + i, start + edge_ref[i]) = 1;
}
Eigen::MatrixXd directions = Eigen::MatrixXd::Identity(ndofs, ndofs);
directions.block(edge_dir.rows() * edge, edge_dir.cols() * edge,
edge_dir.rows(), edge_dir.cols())
= edge_dir;
base_permutations[edge] *= directions;
}
}
else if (tdim == 3)
{
Eigen::ArrayXi face_ref = dofperms::triangle_reflection(degree + 1);
Eigen::ArrayXi face_rot = dofperms::triangle_rotation(degree + 1);
for (int face = 0; face < facet_count; ++face)
{
const int start = face_ref.size() * face;
for (int i = 0; i < face_rot.size(); ++i)
{
base_permutations[6 + 2 * face](start + i, start + i) = 0;
base_permutations[6 + 2 * face](start + i, start + face_rot[i]) = 1;
base_permutations[6 + 2 * face + 1](start + i, start + i) = 0;
base_permutations[6 + 2 * face + 1](start + i, start + face_ref[i])
= -1;
}
}
}
// BDM has facet_dofs dofs on each facet, and ndofs-facet_count*facet_dofs in
// the interior
std::vector<std::vector<int>> entity_dofs(topology.size());
for (int i = 0; i < tdim - 1; ++i)
entity_dofs[i].resize(topology[i].size(), 0);
entity_dofs[tdim - 1].resize(topology[tdim - 1].size(), facet_dofs);
entity_dofs[tdim] = {internal_dofs};
Eigen::MatrixXd coeffs = compute_expansion_coefficients(wcoeffs, dual);
return FiniteElement(name, celltype, degree, {tdim}, coeffs, entity_dofs,
base_permutations, {}, {}, "contravariant piola");
}
//-----------------------------------------------------------------------------
|
