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

Classes

struct  Fe
 

Public Types

using FunRule = boost::function<int(int)>
 

Public Member Functions

 SetIntegrationAtFrontFace (boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges)
 
 SetIntegrationAtFrontFace (boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges, FunRule fun_rule)
 
MoFEMErrorCode operator() (ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
 

Public Attributes

FunRule funRule = [](int p) { return 2 * (p + 1); }
 

Private Attributes

boost::shared_ptr< Range > frontNodes
 
boost::shared_ptr< Range > frontEdges
 

Static Private Attributes

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

Detailed Description

Definition at line 616 of file EshelbianPlasticity.cpp.

Member Typedef Documentation

◆ FunRule

using EshelbianPlasticity::SetIntegrationAtFrontFace::FunRule = boost::function<int(int)>
Examples
EshelbianPlasticity.cpp.

Definition at line 618 of file EshelbianPlasticity.cpp.

Constructor & Destructor Documentation

◆ SetIntegrationAtFrontFace() [1/2]

EshelbianPlasticity::SetIntegrationAtFrontFace::SetIntegrationAtFrontFace ( boost::shared_ptr< Range > front_nodes,
boost::shared_ptr< Range > front_edges )
inline
Examples
EshelbianPlasticity.cpp.

Definition at line 621 of file EshelbianPlasticity.cpp.

623 : frontNodes(front_nodes), frontEdges(front_edges){};

◆ SetIntegrationAtFrontFace() [2/2]

EshelbianPlasticity::SetIntegrationAtFrontFace::SetIntegrationAtFrontFace ( boost::shared_ptr< Range > front_nodes,
boost::shared_ptr< Range > front_edges,
FunRule fun_rule )
inline

Definition at line 625 of file EshelbianPlasticity.cpp.

628 : frontNodes(front_nodes), frontEdges(front_edges), funRule(fun_rule){};

Member Function Documentation

◆ operator()()

MoFEMErrorCode EshelbianPlasticity::SetIntegrationAtFrontFace::operator() ( ForcesAndSourcesCore * fe_raw_ptr,
int order_row,
int order_col,
int order_data )
inline
Examples
EshelbianPlasticity.cpp.

Definition at line 630 of file EshelbianPlasticity.cpp.

631 {
632 MoFEMFunctionBegin;
633
634 constexpr bool debug = false;
635
636 constexpr int numNodes = 3;
637 constexpr int numEdges = 3;
638 constexpr int refinementLevels = 4;
639
640 auto &m_field = fe_raw_ptr->mField;
641 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
642 auto fe_handle = fe_ptr->getFEEntityHandle();
643
644 auto set_base_quadrature = [&]() {
645 MoFEMFunctionBegin;
646 int rule = funRule(order_data);
647 if (rule < QUAD_2D_TABLE_SIZE) {
648 if (QUAD_2D_TABLE[rule]->dim != 2) {
649 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
650 }
651 if (QUAD_2D_TABLE[rule]->order < rule) {
652 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
653 "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
654 }
655 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
656 fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
657 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
658 &fe_ptr->gaussPts(0, 0), 1);
659 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
660 &fe_ptr->gaussPts(1, 0), 1);
661 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
662 &fe_ptr->gaussPts(2, 0), 1);
663 auto &data = fe_ptr->dataOnElement[H1];
664 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
665 false);
666 double *shape_ptr =
667 &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
668 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
669 1);
670 data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
671 std::copy(
672 Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
673 data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
674
675 } else {
676 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
677 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
678 }
679 MoFEMFunctionReturn(0);
680 };
681
682 CHKERR set_base_quadrature();
683
684 if (frontNodes->size() && EshelbianCore::setSingularity) {
685
686 auto get_singular_nodes = [&]() {
687 int num_nodes;
688 const EntityHandle *conn;
689 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
690 true);
691 std::bitset<numNodes> singular_nodes;
692 for (auto nn = 0; nn != numNodes; ++nn) {
693 if (frontNodes->find(conn[nn]) != frontNodes->end()) {
694 singular_nodes.set(nn);
695 } else {
696 singular_nodes.reset(nn);
697 }
698 }
699 return singular_nodes;
700 };
701
702 auto get_singular_edges = [&]() {
703 std::bitset<numEdges> singular_edges;
704 for (int ee = 0; ee != numEdges; ee++) {
705 EntityHandle edge;
706 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
707 if (frontEdges->find(edge) != frontEdges->end()) {
708 singular_edges.set(ee);
709 } else {
710 singular_edges.reset(ee);
711 }
712 }
713 return singular_edges;
714 };
715
716 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
717 MoFEMFunctionBegin;
718 fe_ptr->gaussPts.swap(ref_gauss_pts);
719 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
720 auto &data = fe_ptr->dataOnElement[H1];
721 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
722 double *shape_ptr =
723 &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
724 CHKERR ShapeMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
725 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
726 MoFEMFunctionReturn(0);
727 };
728
729 auto singular_nodes = get_singular_nodes();
730 if (singular_nodes.count()) {
731 auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
732 if (it_map_ref_coords != mapRefCoords.end()) {
733 CHKERR set_gauss_pts(it_map_ref_coords->second);
734 MoFEMFunctionReturnHot(0);
735 } else {
736
737 auto refine_quadrature = [&]() {
738 MoFEMFunctionBegin;
739
740 const int max_level = refinementLevels;
741
742 moab::Core moab_ref;
743 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
744 EntityHandle nodes[numNodes];
745 for (int nn = 0; nn != numNodes; nn++)
746 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
747 EntityHandle tri;
748 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
749 MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
750 MoFEM::Interface &m_field_ref = mofem_ref_core;
751 {
752 Range tris(tri, tri);
753 Range edges;
754 CHKERR m_field_ref.get_moab().get_adjacencies(
755 tris, 1, true, edges, moab::Interface::UNION);
756 CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
757 tris, BitRefLevel().set(0), false, VERBOSE);
758 }
759
760 Range nodes_at_front;
761 for (int nn = 0; nn != numNodes; nn++) {
762 if (singular_nodes[nn]) {
763 EntityHandle ent;
764 CHKERR moab_ref.side_element(tri, 0, nn, ent);
765 nodes_at_front.insert(ent);
766 }
767 }
768
769 auto singular_edges = get_singular_edges();
770
771 EntityHandle meshset;
772 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
773 for (int ee = 0; ee != numEdges; ee++) {
774 if (singular_edges[ee]) {
775 EntityHandle ent;
776 CHKERR moab_ref.side_element(tri, 1, ee, ent);
777 CHKERR moab_ref.add_entities(meshset, &ent, 1);
778 }
779 }
780
781 // refine mesh
782 auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
783 for (int ll = 0; ll != max_level; ll++) {
784 Range edges;
785 CHKERR m_field_ref.getInterface<BitRefManager>()
786 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
787 BitRefLevel().set(), MBEDGE,
788 edges);
789 Range ref_edges;
790 CHKERR moab_ref.get_adjacencies(
791 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
792 ref_edges = intersect(ref_edges, edges);
793 Range ents;
794 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
795 ref_edges = intersect(ref_edges, ents);
796 Range tris;
797 CHKERR m_field_ref.getInterface<BitRefManager>()
798 ->getEntitiesByTypeAndRefLevel(
799 BitRefLevel().set(ll), BitRefLevel().set(), MBTRI, tris);
800 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
801 ref_edges, BitRefLevel().set(ll + 1));
802 CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
803 CHKERR m_field_ref.getInterface<BitRefManager>()
804 ->updateMeshsetByEntitiesChildren(meshset,
805 BitRefLevel().set(ll + 1),
806 meshset, MBEDGE, true);
807 }
808
809 // get ref coords
810 Range tris;
811 CHKERR m_field_ref.getInterface<BitRefManager>()
812 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
813 BitRefLevel().set(), MBTRI,
814 tris);
815
816 if (debug) {
817 CHKERR save_range(moab_ref, "ref_tris.vtk", tris);
818 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "debug");
819 }
820
821 MatrixDouble ref_coords(tris.size(), 9, false);
822 int tt = 0;
823 for (Range::iterator tit = tris.begin(); tit != tris.end();
824 tit++, tt++) {
825 int num_nodes;
826 const EntityHandle *conn;
827 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
828 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
829 }
830
831 auto &data = fe_ptr->dataOnElement[H1];
832 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
833 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
834 MatrixDouble &shape_n =
835 data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
836 int gg = 0;
837 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
838 double *tri_coords = &ref_coords(tt, 0);
839 FTensor::Tensor1<double, 3> t_normal;
840 CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
841 auto det = t_normal.l2();
842 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
843 for (int dd = 0; dd != 2; dd++) {
844 ref_gauss_pts(dd, gg) =
845 shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
846 shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
847 shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
848 }
849 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
850 }
851 }
852
853 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
854 CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
855
856 MoFEMFunctionReturn(0);
857 };
858
859 CHKERR refine_quadrature();
860 }
861 }
862 }
863
864 MoFEMFunctionReturn(0);
865 }
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
Definition dg_projection.cpp:18
debug
static const bool debug
Definition dm_create_subdm.cpp:12
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition fem_tools.c:182
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_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition quad.h:174
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition quad.h:186
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition quad.h:175
EshelbianPlasticity::SetIntegrationAtFrontFace::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition EshelbianPlasticity.cpp:878
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::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition Tools.hpp:104
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition Tools.cpp:353
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
QUAD_::npoints
int npoints
Definition quad.h:29
save_range
auto save_range
Definition thermo_elastic.cpp:121

Member Data Documentation

◆ frontEdges

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontFace::frontEdges
private
Examples
EshelbianPlasticity.cpp.

Definition at line 876 of file EshelbianPlasticity.cpp.

◆ frontNodes

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontFace::frontNodes
private
Examples
EshelbianPlasticity.cpp.

Definition at line 875 of file EshelbianPlasticity.cpp.

◆ funRule

FunRule EshelbianPlasticity::SetIntegrationAtFrontFace::funRule = [](int p) { return 2 * (p + 1); }
Examples
EshelbianPlasticity.cpp.

Definition at line 619 of file EshelbianPlasticity.cpp.

619{ return 2 * (p + 1); };

◆ mapRefCoords

std::map<long int, MatrixDouble> EshelbianPlasticity::SetIntegrationAtFrontFace::mapRefCoords
inlinestaticprivate
Examples
EshelbianPlasticity.cpp.

Definition at line 878 of file EshelbianPlasticity.cpp.

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