v0.15.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Private Attributes | Static Private Attributes | List of all members
EshelbianPlasticity::SetIntegrationAtFrontVolume Struct Reference
Collaboration diagram for EshelbianPlasticity::SetIntegrationAtFrontVolume:
[legend]

Classes

struct  Fe
 

Public Member Functions

 SetIntegrationAtFrontVolume (boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges, boost::shared_ptr< CGGUserPolynomialBase::CachePhi > cache_phi=nullptr)
 
MoFEMErrorCode operator() (ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
 

Private Attributes

boost::shared_ptr< Range > frontNodes
 
boost::shared_ptr< Range > frontEdges
 
boost::shared_ptr< CGGUserPolynomialBase::CachePhi > cachePhi
 

Static Private Attributes

static std::map< long int, MatrixDouble > mapRefCoords
 

Detailed Description

Definition at line 363 of file EshelbianPlasticity.cpp.

Constructor & Destructor Documentation

◆ SetIntegrationAtFrontVolume()

EshelbianPlasticity::SetIntegrationAtFrontVolume::SetIntegrationAtFrontVolume ( boost::shared_ptr< Range >  front_nodes,
boost::shared_ptr< Range >  front_edges,
boost::shared_ptr< CGGUserPolynomialBase::CachePhi >  cache_phi = nullptr 
)
inline

Definition at line 365 of file EshelbianPlasticity.cpp.

369 : frontNodes(front_nodes), frontEdges(front_edges), cachePhi(cache_phi){};
EshelbianPlasticity::SetIntegrationAtFrontVolume::cachePhi
boost::shared_ptr< CGGUserPolynomialBase::CachePhi > cachePhi
Definition EshelbianPlasticity.cpp:625

Member Function Documentation

◆ operator()()

MoFEMErrorCode EshelbianPlasticity::SetIntegrationAtFrontVolume::operator() ( ForcesAndSourcesCore *  fe_raw_ptr,
int  order_row,
int  order_col,
int  order_data 
)
inline
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 371 of file EshelbianPlasticity.cpp.

372 {
373 MoFEMFunctionBegin;
374
375 constexpr bool debug = false;
376
377 constexpr int numNodes = 4;
378 constexpr int numEdges = 6;
379 constexpr int refinementLevels = 6;
380
381 auto &m_field = fe_raw_ptr->mField;
382 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
383 auto fe_handle = fe_ptr->getFEEntityHandle();
384
385 auto set_base_quadrature = [&]() {
386 MoFEMFunctionBegin;
387 int rule = 2 * order_data + 1;
388 if (rule < QUAD_3D_TABLE_SIZE) {
389 if (QUAD_3D_TABLE[rule]->dim != 3) {
390 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
391 "wrong dimension");
392 }
393 if (QUAD_3D_TABLE[rule]->order < rule) {
394 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
395 "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
396 }
397 const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
398 auto &gauss_pts = fe_ptr->gaussPts;
399 gauss_pts.resize(4, nb_gauss_pts, false);
400 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
401 &gauss_pts(0, 0), 1);
402 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
403 &gauss_pts(1, 0), 1);
404 cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
405 &gauss_pts(2, 0), 1);
406 cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
407 &gauss_pts(3, 0), 1);
408 auto &data = fe_ptr->dataOnElement[H1];
409 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
410 false);
411 double *shape_ptr =
412 &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
413 cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
414 1);
415 } else {
416 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
417 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
418 }
419 MoFEMFunctionReturn(0);
420 };
421
422 CHKERR set_base_quadrature();
423
424 if (frontNodes->size() && EshelbianCore::setSingularity) {
425
426 auto get_singular_nodes = [&]() {
427 int num_nodes;
428 const EntityHandle *conn;
429 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
430 true);
431 std::bitset<numNodes> singular_nodes;
432 for (auto nn = 0; nn != numNodes; ++nn) {
433 if (frontNodes->find(conn[nn]) != frontNodes->end()) {
434 singular_nodes.set(nn);
435 } else {
436 singular_nodes.reset(nn);
437 }
438 }
439 return singular_nodes;
440 };
441
442 auto get_singular_edges = [&]() {
443 std::bitset<numEdges> singular_edges;
444 for (int ee = 0; ee != numEdges; ee++) {
445 EntityHandle edge;
446 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
447 if (frontEdges->find(edge) != frontEdges->end()) {
448 singular_edges.set(ee);
449 } else {
450 singular_edges.reset(ee);
451 }
452 }
453 return singular_edges;
454 };
455
456 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
457 MoFEMFunctionBegin;
458 fe_ptr->gaussPts.swap(ref_gauss_pts);
459 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
460 auto &data = fe_ptr->dataOnElement[H1];
461 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
462 double *shape_ptr =
463 &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
464 CHKERR ShapeMBTET(shape_ptr, &fe_ptr->gaussPts(0, 0),
465 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
466 nb_gauss_pts);
467 MoFEMFunctionReturn(0);
468 };
469
470 auto singular_nodes = get_singular_nodes();
471 if (singular_nodes.count()) {
472 auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
473 if (it_map_ref_coords != mapRefCoords.end()) {
474 CHKERR set_gauss_pts(it_map_ref_coords->second);
475 MoFEMFunctionReturnHot(0);
476 } else {
477
478 auto refine_quadrature = [&]() {
479 MoFEMFunctionBegin;
480
481 const int max_level = refinementLevels;
482 EntityHandle tet;
483
484 moab::Core moab_ref;
485 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
486 EntityHandle nodes[4];
487 for (int nn = 0; nn != 4; nn++)
488 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
489 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
490 MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
491 MoFEM::Interface &m_field_ref = mofem_ref_core;
492 {
493 Range tets(tet, tet);
494 Range edges;
495 CHKERR m_field_ref.get_moab().get_adjacencies(
496 tets, 1, true, edges, moab::Interface::UNION);
497 CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
498 tets, BitRefLevel().set(0), false, VERBOSE);
499 }
500
501 Range nodes_at_front;
502 for (int nn = 0; nn != numNodes; nn++) {
503 if (singular_nodes[nn]) {
504 EntityHandle ent;
505 CHKERR moab_ref.side_element(tet, 0, nn, ent);
506 nodes_at_front.insert(ent);
507 }
508 }
509
510 auto singular_edges = get_singular_edges();
511
512 EntityHandle meshset;
513 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
514 for (int ee = 0; ee != numEdges; ee++) {
515 if (singular_edges[ee]) {
516 EntityHandle ent;
517 CHKERR moab_ref.side_element(tet, 1, ee, ent);
518 CHKERR moab_ref.add_entities(meshset, &ent, 1);
519 }
520 }
521
522 // refine mesh
523 auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
524 for (int ll = 0; ll != max_level; ll++) {
525 Range edges;
526 CHKERR m_field_ref.getInterface<BitRefManager>()
527 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
528 BitRefLevel().set(), MBEDGE,
529 edges);
530 Range ref_edges;
531 CHKERR moab_ref.get_adjacencies(
532 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
533 ref_edges = intersect(ref_edges, edges);
534 Range ents;
535 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
536 ref_edges = intersect(ref_edges, ents);
537 Range tets;
538 CHKERR m_field_ref.getInterface<BitRefManager>()
539 ->getEntitiesByTypeAndRefLevel(
540 BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
541 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
542 ref_edges, BitRefLevel().set(ll + 1));
543 CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
544 CHKERR m_field_ref.getInterface<BitRefManager>()
545 ->updateMeshsetByEntitiesChildren(meshset,
546 BitRefLevel().set(ll + 1),
547 meshset, MBEDGE, true);
548 }
549
550 // get ref coords
551 Range tets;
552 CHKERR m_field_ref.getInterface<BitRefManager>()
553 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
554 BitRefLevel().set(), MBTET,
555 tets);
556
557 if (debug) {
558 CHKERR save_range(moab_ref, "ref_tets.vtk", tets);
559 }
560
561 MatrixDouble ref_coords(tets.size(), 12, false);
562 int tt = 0;
563 for (Range::iterator tit = tets.begin(); tit != tets.end();
564 tit++, tt++) {
565 int num_nodes;
566 const EntityHandle *conn;
567 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
568 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
569 }
570
571 auto &data = fe_ptr->dataOnElement[H1];
572 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
573 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
574 MatrixDouble &shape_n =
575 data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
576 int gg = 0;
577 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
578 double *tet_coords = &ref_coords(tt, 0);
579 double det = Tools::tetVolume(tet_coords);
580 det *= 6;
581 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
582 for (int dd = 0; dd != 3; dd++) {
583 ref_gauss_pts(dd, gg) =
584 shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
585 shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
586 shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
587 shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
588 }
589 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
590 }
591 }
592
593 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
594 CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
595
596 // clear cache bubble
597 cachePhi->get<0>() = 0;
598 cachePhi->get<1>() = 0;
599 // tet base cache
600 TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
601 TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
602
603 MoFEMFunctionReturn(0);
604 };
605
606 CHKERR refine_quadrature();
607 }
608 }
609 }
610
611 MoFEMFunctionReturn(0);
612 }
VERBOSE
@ VERBOSE
Definition definitions.h:222
NOBASE
@ NOBASE
Definition definitions.h:59
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
H1
@ H1
continuous field
Definition definitions.h:85
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
order
constexpr int order
Order displacement.
Definition dg_projection.cpp:38
debug
static const bool debug
Definition dm_create_subdm.cpp:12
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition fem_tools.c:306
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition quad.h:186
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition quad.h:187
EshelbianCore::setSingularity
static PetscBool setSingularity
Definition EshelbianCore.hpp:28
EshelbianPlasticity::SetIntegrationAtFrontVolume::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition EshelbianPlasticity.cpp:627
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::MeshRefinement
Mesh refinement interface.
Definition MeshRefinement.hpp:26
MoFEM::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition Tools.cpp:30
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
QUAD_::order
int order
Definition quad.h:28
QUAD_::npoints
int npoints
Definition quad.h:29
save_range
auto save_range
Definition thermo_elastic.cpp:121

Member Data Documentation

◆ cachePhi

boost::shared_ptr<CGGUserPolynomialBase::CachePhi> EshelbianPlasticity::SetIntegrationAtFrontVolume::cachePhi
private
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 625 of file EshelbianPlasticity.cpp.

◆ frontEdges

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontVolume::frontEdges
private
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 623 of file EshelbianPlasticity.cpp.

◆ frontNodes

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontVolume::frontNodes
private
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 622 of file EshelbianPlasticity.cpp.

◆ mapRefCoords

std::map<long int, MatrixDouble> EshelbianPlasticity::SetIntegrationAtFrontVolume::mapRefCoords
inlinestaticprivate
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 627 of file EshelbianPlasticity.cpp.

The documentation for this struct was generated from the following file:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow