33 {
37 int nb_gauss_pts = pts.size2();
38
40 auto &ptr = data.getBBAlphaIndicesSharedPtr(
field_name);
41 if (!ptr)
43 return *ptr;
44 };
45
48 if (!ptr)
50 return *ptr;
51 };
52
54 auto &ptr = data.getBBDiffNSharedPtr(
field_name);
55 if (!ptr)
57 return *ptr;
58 };
59
60 auto get_alpha_by_name_ptr =
61 [](auto &data,
62 const std::string &
field_name) -> boost::shared_ptr<MatrixInt> & {
63 return data.getBBAlphaIndicesSharedPtr(
field_name);
64 };
65
66 auto get_base_by_name_ptr =
67 [](auto &data,
68 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
70 };
71
72 auto get_diff_base_by_name_ptr =
73 [](auto &data,
74 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
76 };
77
78 auto get_alpha_by_order_ptr =
79 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
80 return data.getBBAlphaIndicesByOrderSharedPtr(o);
81 };
82
83 auto get_base_by_order_ptr =
84 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
85 return data.getBBNByOrderSharedPtr(o);
86 };
87
88 auto get_diff_base_by_order_ptr =
89 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
90 return data.getBBDiffNByOrderSharedPtr(o);
91 };
92
93
94
95
97 auto &n_ptr = ent_data.getBBNSharedPtr(
field_name);
98 if (!n_ptr)
99 return;
100 for (auto &ptr : ent_data.getBBNByOrderArray()) {
101 if (ptr && ptr.get() == n_ptr.get()) {
102 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
103 ent_data.getBBDiffNSharedPtr(
field_name).reset();
104 n_ptr.reset();
105 return;
106 }
107 }
108 };
109
110
111
112 auto set_empty = [&](auto &ent_data) {
113 make_own(ent_data);
114 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
115 get_base(ent_data).resize(nb_gauss_pts, 0, false);
116 get_diff_base(ent_data).resize(nb_gauss_pts, 0, false);
117 };
118
119
120
121
122
123
124 auto check_uniform_order = [&]() {
126 auto &node_order = data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder();
127 if (node_order.size() != 3)
129 "Wrong number of vertex orders for BB base: %d",
130 static_cast<int>(node_order.size()));
132 for (
int n = 0;
n != 3; ++
n)
134 auto max_order_on = [&](
const EntityType
type) {
136 for (
auto &ent_data : data.dataOnEntities[
type])
138 };
139 max_order_on(MBEDGE);
140 max_order_on(MBTRI);
143 "Bernstein-Bezier H1 base requires approximation order >= 1");
144 if (
static_cast<size_t>(
order) >=
147 "Bernstein-Bezier order %d exceeds the supported maximum %d",
149 static_cast<int>(
151 1);
152 for (
int n = 0;
n != 3; ++
n)
153 if (node_order[
n] < 1)
155 "Bernstein-Bezier H1 base requires vertex order between 1 "
156 "and the element order %d; vertex %d has order %d",
158 auto check_type = [&](
const EntityType
type,
const int nb_dofs_at_order) {
162 int side = 0;
163 for (
auto &ent_data : data.dataOnEntities[
type]) {
164 if (ent_data.getOrder() !=
order)
166 "Bernstein-Bezier base in H1 space requires uniform "
167 "approximation order on the element: %s %d has order %d, "
168 "element order is %d. BB functions are not hierarchical, "
169 "so mixed degrees do not span a complete polynomial space. "
170 "Set one order for the whole H1 field (L2 fields are not "
171 "restricted).",
172 moab::CN::EntityTypeName(
type), side, ent_data.getOrder(),
174 ++side;
175 }
177 };
181 };
182 CHKERR check_uniform_order();
183
185 auto &vert_get_diff_n = get_diff_base(data.
dataOnEntities[MBVERTEX][0]);
186 vert_get_n.resize(nb_gauss_pts, 3, false);
187 vert_get_diff_n.resize(nb_gauss_pts, 6, false);
188
190 (unsigned int)nb_gauss_pts)
192 "Base functions or nodes has wrong number of integration points "
193 "for base %s",
196
198 vertex_alpha.resize(3, 3, false);
199 vertex_alpha.clear();
200 for (
int n = 0;
n != 3; ++
n)
202
204 1,
lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
206 &vert_get_diff_n(0, 0));
207 for (
int n = 0;
n != 3; ++
n) {
208 const int f = boost::math::factorial<double>(
210 for (
int g = 0;
g != nb_gauss_pts; ++
g) {
211 vert_get_n(
g,
n) *=
f;
212 for (
int d = 0;
d != 2; ++
d)
213 vert_get_diff_n(
g, 2 *
n + d) *=
f;
214 }
215 }
216
217
221 "Wrong size of ent data");
222
223 constexpr int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
224 for (int ee = 0; ee != 3; ++ee) {
226
227 if (ent_data.getSense() == 0)
229 "Sense of the edge unknown");
230
231 const int sense = ent_data.getSense();
232 const int order = ent_data.getOrder();
233 const int nb_dofs =
NBEDGE_H1(ent_data.getOrder());
234
235 if (nb_dofs) {
236 if (get_alpha_by_order_ptr(ent_data,
order)) {
238 get_alpha_by_order_ptr(ent_data,
order);
240 get_base_by_order_ptr(ent_data,
order);
241 get_diff_base_by_name_ptr(ent_data,
field_name) =
242 get_diff_base_by_order_ptr(ent_data,
order);
243 } else {
244 make_own(ent_data);
245 auto &get_n = get_base(ent_data);
246 auto &get_diff_n = get_diff_base(ent_data);
247 get_n.resize(nb_gauss_pts, nb_dofs, false);
248 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
249
250 auto &edge_alpha = get_alpha(ent_data);
251 edge_alpha.resize(nb_dofs, 3, false);
253 &edge_alpha(0, 0));
254 if (sense == -1) {
255 for (
int i = 0;
i != edge_alpha.size1(); ++
i) {
256 int a = edge_alpha(
i, edges_nodes[ee][0]);
257 edge_alpha(
i, edges_nodes[ee][0]) =
258 edge_alpha(
i, edges_nodes[ee][1]);
259 edge_alpha(
i, edges_nodes[ee][1]) =
a;
260 }
261 }
263 order,
lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
265 &get_diff_n(0, 0));
266
267 get_alpha_by_order_ptr(ent_data,
order) =
269 get_base_by_order_ptr(ent_data,
order) =
271 get_diff_base_by_order_ptr(ent_data,
order) =
272 get_diff_base_by_name_ptr(ent_data,
field_name);
273 }
274 } else {
275 set_empty(ent_data);
276 }
277 }
278 } else {
279 for (int ee = 0; ee != 3; ++ee)
281 }
282
286 "Wrong size ent of ent data");
287
289 const int order = ent_data.getOrder();
291 if (get_alpha_by_order_ptr(ent_data,
order)) {
293 get_alpha_by_order_ptr(ent_data,
order);
295 get_base_by_order_ptr(ent_data,
order);
296 get_diff_base_by_name_ptr(ent_data,
field_name) =
297 get_diff_base_by_order_ptr(ent_data,
order);
298 } else {
299
300 make_own(ent_data);
301 auto &get_n = get_base(ent_data);
302 auto &get_diff_n = get_diff_base(ent_data);
303 get_n.resize(nb_gauss_pts, nb_dofs, false);
304 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
305 if (nb_dofs) {
306 auto &tri_alpha = get_alpha(ent_data);
307 tri_alpha.resize(nb_dofs, 3, false);
308
311 order,
lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
313 &get_diff_n(0, 0));
314
315 get_alpha_by_order_ptr(ent_data,
order) =
317 get_base_by_order_ptr(ent_data,
order) =
319 get_diff_base_by_order_ptr(ent_data,
order) =
320 get_diff_base_by_name_ptr(ent_data,
field_name);
321 } else {
322 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
323 }
324 }
325 } else {
327 }
328
330}
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
UBlasMatrix< double > MatrixDouble
UBlasMatrix< int > MatrixInt
constexpr auto field_name
static MoFEMErrorCode baseFunctionsTri(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])
const std::string fieldName
static constexpr size_t MaxBernsteinBezierOrder
Public Member Functions inherited from MoFEM::BaseFunction