243 {
245
248 const int nb_gauss_pts = pts.size2();
249
251 (unsigned int)nb_gauss_pts)
253 "Base functions or nodes has wrong number of integration points "
254 "for base %s",
257
259 auto &ptr = data.getBBAlphaIndicesSharedPtr(
field_name);
260 if (!ptr)
262 return *ptr;
263 };
264
267 if (!ptr)
269 return *ptr;
270 };
271
273 auto &ptr = data.getBBDiffNSharedPtr(
field_name);
274 if (!ptr)
276 return *ptr;
277 };
278
279 auto get_alpha_by_name_ptr =
280 [](auto &data,
281 const std::string &
field_name) -> boost::shared_ptr<MatrixInt> & {
282 return data.getBBAlphaIndicesSharedPtr(
field_name);
283 };
284
285 auto get_base_by_name_ptr =
286 [](auto &data,
287 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
289 };
290
291 auto get_diff_base_by_name_ptr =
292 [](auto &data,
293 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
295 };
296
297 auto get_alpha_by_order_ptr =
298 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
299 return data.getBBAlphaIndicesByOrderSharedPtr(o);
300 };
301
302 auto get_base_by_order_ptr =
303 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
304 return data.getBBNByOrderSharedPtr(o);
305 };
306
307 auto get_diff_base_by_order_ptr =
308 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
309 return data.getBBDiffNByOrderSharedPtr(o);
310 };
311
313 auto &vertex_alpha = get_alpha(vert_ent_data);
314 vertex_alpha.resize(4, 4, false);
315 vertex_alpha.clear();
316 for (
int n = 0;
n != 4; ++
n)
318
319 auto &vert_get_n = get_base(vert_ent_data);
320 auto &vert_get_diff_n = get_diff_base(vert_ent_data);
321 vert_get_n.resize(nb_gauss_pts, 4, false);
322 vert_get_diff_n.resize(nb_gauss_pts, 12, false);
324 1,
lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
326 &vert_get_diff_n(0, 0));
327 for (
int n = 0;
n != 4; ++
n) {
328 const double f = boost::math::factorial<double>(
330 for (
int g = 0;
g != nb_gauss_pts; ++
g) {
331 vert_get_n(
g,
n) *=
f;
332 for (
int d = 0;
d != 3; ++
d)
333 vert_get_diff_n(
g, 3 *
n + d) *=
f;
334 }
335 }
336
337
341 "Wrong size of ent data");
342
343 constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
344 {0, 3}, {1, 3}, {2, 3}};
345 for (int ee = 0; ee != 6; ++ee) {
347 const int sense = ent_data.getSense();
348 if (sense == 0)
350 "Sense of the edge unknown");
351 const int order = ent_data.getOrder();
353
354 if (nb_dofs) {
355 if (get_alpha_by_order_ptr(ent_data,
order)) {
357 get_alpha_by_order_ptr(ent_data,
order);
359 get_base_by_order_ptr(ent_data,
order);
360 get_diff_base_by_name_ptr(ent_data,
field_name) =
361 get_diff_base_by_order_ptr(ent_data,
order);
362 } else {
363 auto &get_n = get_base(ent_data);
364 auto &get_diff_n = get_diff_base(ent_data);
365 get_n.resize(nb_gauss_pts, nb_dofs, false);
366 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
367
369 edge_alpha.resize(nb_dofs, 4, false);
371 &edge_alpha(0, 0));
372 if (sense == -1) {
373 for (
int i = 0;
i != edge_alpha.size1(); ++
i) {
374 int a = edge_alpha(
i, edges_nodes[ee][0]);
375 edge_alpha(
i, edges_nodes[ee][0]) =
376 edge_alpha(
i, edges_nodes[ee][1]);
377 edge_alpha(
i, edges_nodes[ee][1]) =
a;
378 }
379 }
381 order,
lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
383 &get_diff_n(0, 0));
384
385 get_alpha_by_order_ptr(ent_data,
order) =
387 get_base_by_order_ptr(ent_data,
order) =
389 get_diff_base_by_order_ptr(ent_data,
order) =
390 get_diff_base_by_name_ptr(ent_data,
field_name);
391 }
392 }
393 }
394 } else {
395 for (int ee = 0; ee != 6; ++ee) {
397 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
398 auto &get_n = get_base(ent_data);
399 auto &get_diff_n = get_diff_base(ent_data);
400 get_n.resize(nb_gauss_pts, 0, false);
401 get_diff_n.resize(nb_gauss_pts, 0, false);
402 }
403 }
404
405
409 "Wrong size of ent data");
414
415 for (int ff = 0; ff != 4; ++ff) {
417 const int order = ent_data.getOrder();
419
420 if (nb_dofs) {
421 if (get_alpha_by_order_ptr(ent_data,
order)) {
423 get_alpha_by_order_ptr(ent_data,
order);
425 get_base_by_order_ptr(ent_data,
order);
426 get_diff_base_by_name_ptr(ent_data,
field_name) =
427 get_diff_base_by_order_ptr(ent_data,
order);
428 } else {
429
430 auto &get_n = get_base(ent_data);
431 auto &get_diff_n = get_diff_base(ent_data);
432 get_n.resize(nb_gauss_pts, nb_dofs, false);
433 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
434
435 auto &face_alpha = get_alpha(ent_data);
436 face_alpha.resize(nb_dofs, 4, false);
437
439 &face_alpha(0, 0));
440 senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(),
false);
442 constexpr int tri_nodes[4][3] = {
443 {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
444 for (
int d = 0;
d != nb_dofs; ++
d)
445 for (
int n = 0;
n != 3; ++
n)
447 face_alpha(d, tri_nodes[ff][
n]);
450 order,
lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
452 &get_diff_n(0, 0));
453
454 get_alpha_by_order_ptr(ent_data,
order) =
456 get_base_by_order_ptr(ent_data,
order) =
458 get_diff_base_by_order_ptr(ent_data,
order) =
459 get_diff_base_by_name_ptr(ent_data,
field_name);
460 }
461 }
462 }
463 } else {
464 for (int ff = 0; ff != 4; ++ff) {
466 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
467 auto &get_n = get_base(ent_data);
468 auto &get_diff_n = get_diff_base(ent_data);
469 get_n.resize(nb_gauss_pts, 0, false);
470 get_diff_n.resize(nb_gauss_pts, 0, false);
471 }
472 }
473
477 "Wrong size ent of ent data");
478
480 const int order = ent_data.getOrder();
482 if (get_alpha_by_order_ptr(ent_data,
order)) {
484 get_alpha_by_order_ptr(ent_data,
order);
486 get_base_by_order_ptr(ent_data,
order);
487 get_diff_base_by_name_ptr(ent_data,
field_name) =
488 get_diff_base_by_order_ptr(ent_data,
order);
489 } else {
490
491 auto &get_n = get_base(ent_data);
492 auto &get_diff_n = get_diff_base(ent_data);
493 get_n.resize(nb_gauss_pts, nb_dofs, false);
494 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
495 if (nb_dofs) {
496 auto &tet_alpha = get_alpha(ent_data);
497 tet_alpha.resize(nb_dofs, 4, false);
498
501 order,
lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
503 &get_diff_n(0, 0));
504
505 get_alpha_by_order_ptr(ent_data,
order) =
507 get_base_by_order_ptr(ent_data,
order) =
509 get_diff_base_by_order_ptr(ent_data,
order) =
510 get_diff_base_by_name_ptr(ent_data,
field_name);
511 }
512 }
513 } else {
515 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
516 auto &get_n = get_base(ent_data);
517 auto &get_diff_n = get_diff_base(ent_data);
518 get_n.resize(nb_gauss_pts, 0, false);
519 get_diff_n.resize(nb_gauss_pts, 0, false);
520 }
521
523}
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 generateIndicesTriTet(const int N[], int *alpha[])
static MoFEMErrorCode baseFunctionsTet(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 generateIndicesEdgeTet(const int N[], int *alpha[])
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
const std::string fieldName
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Public Member Functions inherited from MoFEM::BaseFunction