v0.14.0
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)
 
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 332 of file EshelbianPlasticity.cpp.

Constructor & Destructor Documentation

◆ SetIntegrationAtFrontVolume()

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

Definition at line 334 of file EshelbianPlasticity.cpp.

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

Member Function Documentation

◆ operator()()

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

Definition at line 338 of file EshelbianPlasticity.cpp.

339  {
340  MoFEMFunctionBegin;
341 
342  constexpr bool debug = false;
343 
344  constexpr int numNodes = 4;
345  constexpr int numEdges = 6;
346  constexpr int refinementLevels = 4;
347 
348  auto &m_field = fe_raw_ptr->mField;
349  auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
350  auto fe_handle = fe_ptr->getFEEntityHandle();
351 
352  auto set_base_quadrature = [&]() {
353  MoFEMFunctionBegin;
354  int rule = 2 * order_data + 1;
355  if (rule < QUAD_3D_TABLE_SIZE) {
356  if (QUAD_3D_TABLE[rule]->dim != 3) {
357  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
358  "wrong dimension");
359  }
360  if (QUAD_3D_TABLE[rule]->order < rule) {
361  SETERRQ2(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
362  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
363  }
364  const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
365  auto &gauss_pts = fe_ptr->gaussPts;
366  gauss_pts.resize(4, nb_gauss_pts, false);
367  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
368  &gauss_pts(0, 0), 1);
369  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
370  &gauss_pts(1, 0), 1);
371  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
372  &gauss_pts(2, 0), 1);
373  cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
374  &gauss_pts(3, 0), 1);
375  auto &data = fe_ptr->dataOnElement[H1];
376  data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
377  false);
378  double *shape_ptr =
379  &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
380  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
381  1);
382  } else {
383  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
384  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
385  }
386  MoFEMFunctionReturn(0);
387  };
388 
389  CHKERR set_base_quadrature();
390 
391  if (frontNodes->size() && EshelbianCore::setSingularity) {
392 
393  auto get_singular_nodes = [&]() {
394  int num_nodes;
395  const EntityHandle *conn;
396  CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
397  true);
398  std::bitset<numNodes> singular_nodes;
399  for (auto nn = 0; nn != numNodes; ++nn) {
400  if (frontNodes->find(conn[nn]) != frontNodes->end()) {
401  singular_nodes.set(nn);
402  } else {
403  singular_nodes.reset(nn);
404  }
405  }
406  return singular_nodes;
407  };
408 
409  auto get_singular_edges = [&]() {
410  std::bitset<numEdges> singular_edges;
411  for (int ee = 0; ee != numEdges; ee++) {
412  EntityHandle edge;
413  CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
414  if (frontEdges->find(edge) != frontEdges->end()) {
415  singular_edges.set(ee);
416  } else {
417  singular_edges.reset(ee);
418  }
419  }
420  return singular_edges;
421  };
422 
423  auto set_gauss_pts = [&](auto &ref_gauss_pts) {
424  MoFEMFunctionBegin;
425  fe_ptr->gaussPts.swap(ref_gauss_pts);
426  const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
427  auto &data = fe_ptr->dataOnElement[H1];
428  data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
429  double *shape_ptr =
430  &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
431  CHKERR ShapeMBTET(shape_ptr, &fe_ptr->gaussPts(0, 0),
432  &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
433  nb_gauss_pts);
434  MoFEMFunctionReturn(0);
435  };
436 
437  auto singular_nodes = get_singular_nodes();
438  if (singular_nodes.count()) {
439  auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
440  if (it_map_ref_coords != mapRefCoords.end()) {
441  CHKERR set_gauss_pts(it_map_ref_coords->second);
442  MoFEMFunctionReturnHot(0);
443  } else {
444 
445  auto refine_quadrature = [&]() {
446  MoFEMFunctionBegin;
447 
448  const int max_level = refinementLevels;
449  EntityHandle tet;
450 
451  moab::Core moab_ref;
452  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
453  EntityHandle nodes[4];
454  for (int nn = 0; nn != 4; nn++)
455  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
456  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
457  MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
458  MoFEM::Interface &m_field_ref = mofem_ref_core;
459  {
460  Range tets(tet, tet);
461  Range edges;
462  CHKERR m_field_ref.get_moab().get_adjacencies(
463  tets, 1, true, edges, moab::Interface::UNION);
464  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
465  tets, BitRefLevel().set(0), false, VERBOSE);
466  }
467 
468  Range nodes_at_front;
469  for (int nn = 0; nn != numNodes; nn++) {
470  if (singular_nodes[nn]) {
471  EntityHandle ent;
472  CHKERR moab_ref.side_element(tet, 0, nn, ent);
473  nodes_at_front.insert(ent);
474  }
475  }
476 
477  auto singular_edges = get_singular_edges();
478 
479  EntityHandle meshset;
480  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
481  for (int ee = 0; ee != numEdges; ee++) {
482  if (singular_edges[ee]) {
483  EntityHandle ent;
484  CHKERR moab_ref.side_element(tet, 1, ee, ent);
485  CHKERR moab_ref.add_entities(meshset, &ent, 1);
486  }
487  }
488 
489  // refine mesh
490  auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
491  for (int ll = 0; ll != max_level; ll++) {
492  Range edges;
493  CHKERR m_field_ref.getInterface<BitRefManager>()
494  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
495  BitRefLevel().set(), MBEDGE,
496  edges);
497  Range ref_edges;
498  CHKERR moab_ref.get_adjacencies(
499  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
500  ref_edges = intersect(ref_edges, edges);
501  Range ents;
502  CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
503  ref_edges = intersect(ref_edges, ents);
504  Range tets;
505  CHKERR m_field_ref.getInterface<BitRefManager>()
506  ->getEntitiesByTypeAndRefLevel(
507  BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
508  CHKERR m_ref->addVerticesInTheMiddleOfEdges(
509  ref_edges, BitRefLevel().set(ll + 1));
510  CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
511  CHKERR m_field_ref.getInterface<BitRefManager>()
512  ->updateMeshsetByEntitiesChildren(meshset,
513  BitRefLevel().set(ll + 1),
514  meshset, MBEDGE, true);
515  }
516 
517  // get ref coords
518  Range tets;
519  CHKERR m_field_ref.getInterface<BitRefManager>()
520  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
521  BitRefLevel().set(), MBTET,
522  tets);
523 
524  if (debug) {
525  CHKERR save_range(moab_ref, "ref_tets.vtk", tets);
526  }
527 
528  MatrixDouble ref_coords(tets.size(), 12, false);
529  int tt = 0;
530  for (Range::iterator tit = tets.begin(); tit != tets.end();
531  tit++, tt++) {
532  int num_nodes;
533  const EntityHandle *conn;
534  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
535  CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
536  }
537 
538  auto &data = fe_ptr->dataOnElement[H1];
539  const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
540  MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
541  MatrixDouble &shape_n =
542  data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
543  int gg = 0;
544  for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
545  double *tet_coords = &ref_coords(tt, 0);
546  double det = Tools::tetVolume(tet_coords);
547  det *= 6;
548  for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
549  for (int dd = 0; dd != 3; dd++) {
550  ref_gauss_pts(dd, gg) =
551  shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
552  shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
553  shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
554  shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
555  }
556  ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
557  }
558  }
559 
560  mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
561  CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
562 
563  MoFEMFunctionReturn(0);
564  };
565 
566  CHKERR refine_quadrature();
567  }
568  }
569  }
570 
571  MoFEMFunctionReturn(0);
572  }

Member Data Documentation

◆ frontEdges

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

Definition at line 583 of file EshelbianPlasticity.cpp.

◆ frontNodes

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

Definition at line 582 of file EshelbianPlasticity.cpp.

◆ mapRefCoords

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

Definition at line 585 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
EntityHandle
EshelbianPlasticity::SetIntegrationAtFrontVolume::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition: EshelbianPlasticity.cpp:585
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
QUAD_::order
int order
Definition: quad.h:28
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
QUAD_::npoints
int npoints
Definition: quad.h:29
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
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
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
EshelbianPlasticity::SetIntegrationAtFrontVolume::frontNodes
boost::shared_ptr< Range > frontNodes
Definition: EshelbianPlasticity.cpp:582
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
EshelbianPlasticity::SetIntegrationAtFrontVolume::frontEdges
boost::shared_ptr< Range > frontEdges
Definition: EshelbianPlasticity.cpp:583
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
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow