v0.14.0
Classes | Public Member Functions | Private Attributes | Static Private Attributes | List of all members
EshelbianPlasticity::SetIntegrationAtFrontFace Struct Reference
Collaboration diagram for EshelbianPlasticity::SetIntegrationAtFrontFace:
[legend]

Classes

struct  Fe
 

Public Member Functions

 SetIntegrationAtFrontFace (boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges)
 
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
 

Static Private Attributes

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

Detailed Description

Definition at line 588 of file EshelbianPlasticity.cpp.

Constructor & Destructor Documentation

◆ SetIntegrationAtFrontFace()

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

Definition at line 590 of file EshelbianPlasticity.cpp.

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

Member Function Documentation

◆ operator()()

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

Definition at line 594 of file EshelbianPlasticity.cpp.

595  {
596  MoFEMFunctionBegin;
597 
598  constexpr bool debug = false;
599 
600  constexpr int numNodes = 3;
601  constexpr int numEdges = 3;
602  constexpr int refinementLevels = 4;
603 
604  auto &m_field = fe_raw_ptr->mField;
605  auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
606  auto fe_handle = fe_ptr->getFEEntityHandle();
607 
608  auto set_base_quadrature = [&]() {
609  MoFEMFunctionBegin;
610  int rule = 2 * (order_data + 1);
611  if (rule < QUAD_2D_TABLE_SIZE) {
612  if (QUAD_2D_TABLE[rule]->dim != 2) {
613  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
614  }
615  if (QUAD_2D_TABLE[rule]->order < rule) {
616  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
617  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
618  }
619  const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
620  fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
621  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
622  &fe_ptr->gaussPts(0, 0), 1);
623  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
624  &fe_ptr->gaussPts(1, 0), 1);
625  cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
626  &fe_ptr->gaussPts(2, 0), 1);
627  auto &data = fe_ptr->dataOnElement[H1];
628  data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
629  false);
630  double *shape_ptr =
631  &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
632  cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
633  1);
634  data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
635  std::copy(
636  Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
637  data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
638 
639  } else {
640  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
641  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
642  }
643  MoFEMFunctionReturn(0);
644  };
645 
646  CHKERR set_base_quadrature();
647 
648  if (frontNodes->size() && EshelbianCore::setSingularity) {
649 
650  auto get_singular_nodes = [&]() {
651  int num_nodes;
652  const EntityHandle *conn;
653  CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
654  true);
655  std::bitset<numNodes> singular_nodes;
656  for (auto nn = 0; nn != numNodes; ++nn) {
657  if (frontNodes->find(conn[nn]) != frontNodes->end()) {
658  singular_nodes.set(nn);
659  } else {
660  singular_nodes.reset(nn);
661  }
662  }
663  return singular_nodes;
664  };
665 
666  auto get_singular_edges = [&]() {
667  std::bitset<numEdges> singular_edges;
668  for (int ee = 0; ee != numEdges; ee++) {
669  EntityHandle edge;
670  CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
671  if (frontEdges->find(edge) != frontEdges->end()) {
672  singular_edges.set(ee);
673  } else {
674  singular_edges.reset(ee);
675  }
676  }
677  return singular_edges;
678  };
679 
680  auto set_gauss_pts = [&](auto &ref_gauss_pts) {
681  MoFEMFunctionBegin;
682  fe_ptr->gaussPts.swap(ref_gauss_pts);
683  const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
684  auto &data = fe_ptr->dataOnElement[H1];
685  data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
686  double *shape_ptr =
687  &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
688  CHKERR ShapeMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
689  &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
690  MoFEMFunctionReturn(0);
691  };
692 
693  auto singular_nodes = get_singular_nodes();
694  if (singular_nodes.count()) {
695  auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
696  if (it_map_ref_coords != mapRefCoords.end()) {
697  CHKERR set_gauss_pts(it_map_ref_coords->second);
698  MoFEMFunctionReturnHot(0);
699  } else {
700 
701  auto refine_quadrature = [&]() {
702  MoFEMFunctionBegin;
703 
704  const int max_level = refinementLevels;
705 
706  moab::Core moab_ref;
707  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
708  EntityHandle nodes[numNodes];
709  for (int nn = 0; nn != numNodes; nn++)
710  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
711  EntityHandle tri;
712  CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
713  MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
714  MoFEM::Interface &m_field_ref = mofem_ref_core;
715  {
716  Range tris(tri, tri);
717  Range edges;
718  CHKERR m_field_ref.get_moab().get_adjacencies(
719  tris, 1, true, edges, moab::Interface::UNION);
720  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
721  tris, BitRefLevel().set(0), false, VERBOSE);
722  }
723 
724  Range nodes_at_front;
725  for (int nn = 0; nn != numNodes; nn++) {
726  if (singular_nodes[nn]) {
727  EntityHandle ent;
728  CHKERR moab_ref.side_element(tri, 0, nn, ent);
729  nodes_at_front.insert(ent);
730  }
731  }
732 
733  auto singular_edges = get_singular_edges();
734 
735  EntityHandle meshset;
736  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
737  for (int ee = 0; ee != numEdges; ee++) {
738  if (singular_edges[ee]) {
739  EntityHandle ent;
740  CHKERR moab_ref.side_element(tri, 1, ee, ent);
741  CHKERR moab_ref.add_entities(meshset, &ent, 1);
742  }
743  }
744 
745  // refine mesh
746  auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
747  for (int ll = 0; ll != max_level; ll++) {
748  Range edges;
749  CHKERR m_field_ref.getInterface<BitRefManager>()
750  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
751  BitRefLevel().set(), MBEDGE,
752  edges);
753  Range ref_edges;
754  CHKERR moab_ref.get_adjacencies(
755  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
756  ref_edges = intersect(ref_edges, edges);
757  Range ents;
758  CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
759  ref_edges = intersect(ref_edges, ents);
760  Range tris;
761  CHKERR m_field_ref.getInterface<BitRefManager>()
762  ->getEntitiesByTypeAndRefLevel(
763  BitRefLevel().set(ll), BitRefLevel().set(), MBTRI, tris);
764  CHKERR m_ref->addVerticesInTheMiddleOfEdges(
765  ref_edges, BitRefLevel().set(ll + 1));
766  CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
767  CHKERR m_field_ref.getInterface<BitRefManager>()
768  ->updateMeshsetByEntitiesChildren(meshset,
769  BitRefLevel().set(ll + 1),
770  meshset, MBEDGE, true);
771  }
772 
773  // get ref coords
774  Range tris;
775  CHKERR m_field_ref.getInterface<BitRefManager>()
776  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
777  BitRefLevel().set(), MBTRI,
778  tris);
779 
780  if (debug) {
781  CHKERR save_range(moab_ref, "ref_tris.vtk", tris);
782  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "debug");
783  }
784 
785  MatrixDouble ref_coords(tris.size(), 9, false);
786  int tt = 0;
787  for (Range::iterator tit = tris.begin(); tit != tris.end();
788  tit++, tt++) {
789  int num_nodes;
790  const EntityHandle *conn;
791  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
792  CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
793  }
794 
795  auto &data = fe_ptr->dataOnElement[H1];
796  const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
797  MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
798  MatrixDouble &shape_n =
799  data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
800  int gg = 0;
801  for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
802  double *tri_coords = &ref_coords(tt, 0);
803  FTensor::Tensor1<double, 3> t_normal;
804  CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
805  auto det = t_normal.l2();
806  for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
807  for (int dd = 0; dd != 2; dd++) {
808  ref_gauss_pts(dd, gg) =
809  shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
810  shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
811  shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
812  }
813  ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
814  }
815  }
816 
817  mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
818  CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
819 
820  MoFEMFunctionReturn(0);
821  };
822 
823  CHKERR refine_quadrature();
824  }
825  }
826  }
827 
828  MoFEMFunctionReturn(0);
829  }

Member Data Documentation

◆ frontEdges

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontFace::frontEdges
private

Definition at line 840 of file EshelbianPlasticity.cpp.

◆ frontNodes

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontFace::frontNodes
private

Definition at line 839 of file EshelbianPlasticity.cpp.

◆ mapRefCoords

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

Definition at line 842 of file EshelbianPlasticity.cpp.

The documentation for this struct was generated from the following file:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
H1
@ H1
continuous field
Definition: definitions.h:85
save_range
static auto save_range(moab::Interface &moab, const std::string name, const Range r, std::vector< Tag > tags={})
Definition: EshelbianPlasticity.cpp:128
FTensor::Tensor1< double, 3 >
EntityHandle
EshelbianPlasticity::SetIntegrationAtFrontFace::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition: EshelbianPlasticity.cpp:842
NOBASE
@ NOBASE
Definition: definitions.h:59
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
QUAD_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EshelbianCore::setSingularity
static PetscBool setSingularity
Definition: EshelbianCore.hpp:19
VERBOSE
@ VERBOSE
Definition: definitions.h:222
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MoFEM::CoreTmp
Definition: Core.hpp:36
EshelbianPlasticity::SetIntegrationAtFrontFace::frontEdges
boost::shared_ptr< Range > frontEdges
Definition: EshelbianPlasticity.cpp:840
QUAD_::npoints
int npoints
Definition: quad.h:29
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
Range
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::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
EshelbianPlasticity::SetIntegrationAtFrontFace::frontNodes
boost::shared_ptr< Range > frontNodes
Definition: EshelbianPlasticity.cpp:839
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
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
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow