v0.15.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Static Private Attributes | List of all members
SetEnrichedIntegration Struct Reference

#include "users_modules/multifield-thermoplasticity-private/src/SolutionMapping.hpp"

Collaboration diagram for SetEnrichedIntegration:
[legend]

Classes

struct  Fe
 

Public Member Functions

 SetEnrichedIntegration ()
 
MoFEMErrorCode operator() (ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
 

Static Private Attributes

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

Detailed Description

Examples
SolutionMapping.hpp.

Definition at line 8 of file SolutionMapping.hpp.

Constructor & Destructor Documentation

◆ SetEnrichedIntegration()

SetEnrichedIntegration::SetEnrichedIntegration ( )
Examples
SolutionMapping.hpp, and thermoplastic.cpp.

Definition at line 643 of file thermoplastic.cpp.

643{}

Member Function Documentation

◆ operator()()

MoFEMErrorCode SetEnrichedIntegration::operator() ( ForcesAndSourcesCore *  fe_raw_ptr,
int  order_row,
int  order_col,
int  order_data 
)
Examples
SolutionMapping.hpp, and thermoplastic.cpp.

Definition at line 646 of file thermoplastic.cpp.

648 {
649 MoFEMFunctionBegin;
650
651 constexpr int numNodes = 3;
652 constexpr int numEdges = 3;
653
654 auto &m_field = fe_raw_ptr->mField;
655 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
656 auto fe_handle = fe_ptr->getFEEntityHandle();
657
658 auto set_base_quadrature = [&]() {
659 MoFEMFunctionBegin;
660 int rule = 2 * (order_data + 1);
661 if (rule < QUAD_2D_TABLE_SIZE) {
662 if (QUAD_2D_TABLE[rule]->dim != 2) {
663 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
664 }
665 if (QUAD_2D_TABLE[rule]->order < rule) {
666 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
667 "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
668 }
669 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
670 fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
671 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
672 &fe_ptr->gaussPts(0, 0), 1);
673 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
674 &fe_ptr->gaussPts(1, 0), 1);
675 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
676 &fe_ptr->gaussPts(2, 0), 1);
677 auto &data = fe_ptr->dataOnElement[H1];
678 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
679 false);
680 double *shape_ptr =
681 &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
682 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
683 1);
684 data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
685 std::copy(
686 Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
687 data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
688
689 } else {
690 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
691 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
692 }
693 MoFEMFunctionReturn(0);
694 };
695
696 CHKERR set_base_quadrature();
697
698 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
699 MoFEMFunctionBegin;
700 fe_ptr->gaussPts.swap(ref_gauss_pts);
701 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
702 auto &data = fe_ptr->dataOnElement[H1];
703 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3);
704 double *shape_ptr =
705 &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
706 CHKERR ShapeMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
707 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
708 MoFEMFunctionReturn(0);
709 };
710
711 auto refine_quadrature = [&]() {
712 MoFEMFunctionBegin;
713
714 moab::Core moab_ref;
715 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
716 EntityHandle nodes[numNodes];
717 for (int nn = 0; nn != numNodes; nn++)
718 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
719 EntityHandle tri;
720 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
721 MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
722 MoFEM::Interface &m_field_ref = mofem_ref_core;
723
724 Range ref_tri(tri, tri);
725 Range edges;
726 CHKERR m_field_ref.get_moab().get_adjacencies(ref_tri, 1, true, edges,
727 moab::Interface::UNION);
728 CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
729 ref_tri, BitRefLevel().set(0), false, VERBOSE);
730
731 Range nodes_at_front;
732 for (int nn = 0; nn != numNodes; nn++) {
733 EntityHandle ent;
734 CHKERR moab_ref.side_element(tri, 0, nn, ent);
735 nodes_at_front.insert(ent);
736 }
737
738 EntityHandle meshset;
739 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
740 for (int ee = 0; ee != numEdges; ee++) {
741 EntityHandle ent;
742 CHKERR moab_ref.side_element(tri, 1, ee, ent);
743 CHKERR moab_ref.add_entities(meshset, &ent, 1);
744 }
745
746 // refine mesh
747 auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
748 Range ref_edges;
749 CHKERR m_field_ref.getInterface<BitRefManager>()
750 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(0),
751 BitRefLevel().set(), MBEDGE, ref_edges);
752
753 // Range tris_to_refine;
754 // CHKERR m_field_ref.getInterface<BitRefManager>()
755 // ->getEntitiesByTypeAndRefLevel(
756 // BitRefLevel().set(0), BitRefLevel().set(), MBTRI,
757 // tris_to_refine);
758 CHKERR m_ref->addVerticesInTheMiddleOfEdges(ref_edges,
759 BitRefLevel().set(1));
760 CHKERR m_ref->addVerticesInTheMiddleOfFaces(ref_tri, BitRefLevel().set(1));
761 CHKERR m_ref->refineTris(ref_tri, BitRefLevel().set(1));
762 CHKERR m_field_ref.getInterface<BitRefManager>()
763 ->updateMeshsetByEntitiesChildren(meshset, BitRefLevel().set(1),
764 meshset, MBEDGE, true);
765 CHKERR m_field_ref.getInterface<BitRefManager>()
766 ->updateMeshsetByEntitiesChildren(meshset, BitRefLevel().set(1),
767 meshset, MBTRI, true);
768
769 // get ref coords
770 Range output_ref_tris;
771 CHKERR m_field_ref.getInterface<BitRefManager>()
772 ->getEntitiesByTypeAndRefLevel(
773 BitRefLevel().set(1), BitRefLevel().set(), MBTRI, output_ref_tris);
774
775#ifndef NDEBUG
776 CHKERR save_range(moab_ref, "ref_tris.vtk", output_ref_tris);
777#endif
778
779 MatrixDouble ref_coords(output_ref_tris.size(), 9, false);
780 int tt = 0;
781 for (Range::iterator tit = output_ref_tris.begin();
782 tit != output_ref_tris.end(); tit++, tt++) {
783 int num_nodes;
784 const EntityHandle *conn;
785 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
786 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
787 }
788
789 auto &data = fe_ptr->dataOnElement[H1];
790 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
791 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
792 MatrixDouble &shape_n = data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
793 int gg = 0;
794 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
795 double *tri_coords = &ref_coords(tt, 0);
796 FTensor::Tensor1<double, 3> t_normal;
797 CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
798 auto det = t_normal.l2();
799 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
800 for (int dd = 0; dd != 2; dd++) {
801 ref_gauss_pts(dd, gg) = shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
802 shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
803 shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
804 }
805 MOFEM_LOG("REMESHING", Sev::noisy)
806 << ref_gauss_pts(0, gg) << ", " << ref_gauss_pts(1, gg);
807 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
808 }
809 }
810
811 int num_nodes;
812 const EntityHandle *conn;
813 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
814 true);
815 std::bitset<numNodes> all_nodes;
816 for (auto nn = 0; nn != numNodes; ++nn) {
817 all_nodes.set(nn);
818 }
819
820 mapRefCoords[all_nodes.to_ulong()].swap(ref_gauss_pts);
821 CHKERR set_gauss_pts(mapRefCoords[all_nodes.to_ulong()]);
822
823 MoFEMFunctionReturn(0);
824 };
825
826 CHKERR refine_quadrature();
827
828 MoFEMFunctionReturn(0);
829}
VERBOSE
@ VERBOSE
Definition definitions.h:222
NOBASE
@ NOBASE
Definition definitions.h:59
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
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
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
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
SetEnrichedIntegration::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition SolutionMapping.hpp:23
save_range
auto save_range
Definition thermo_elastic.cpp:121

Member Data Documentation

◆ mapRefCoords

std::map<long int, MatrixDouble> SetEnrichedIntegration::mapRefCoords
inlinestaticprivate
Examples
SolutionMapping.hpp, and thermoplastic.cpp.

Definition at line 23 of file SolutionMapping.hpp.

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