EshelbianPlasticity.cpp File Reference

Implementation of automatic differentiation. More...

#include <MoFEM.hpp>
#include <CGGTonsorialBubbleBase.hpp>
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
#include <cholesky.hpp>
#include <EshelbianAux.hpp>
#include <EshelbianContact.hpp>
#include <phg-quadrule/quad.h>
#include <queue>
#include "impl/EshelbianMonitor.cpp"
#include "impl/SetUpSchurImpl.cpp"
#include "impl/TSElasticPostStep.cpp"
#include <impl/EshelbianContact.cpp>

Go to the source code of this file.

Classes
struct	EshelbianPlasticity::VolUserDataOperatorStabAssembly
struct	EshelbianPlasticity::FaceUserDataOperatorStabAssembly
struct	EshelbianPlasticity::SetIntegrationAtFrontVolume
struct	EshelbianPlasticity::SetIntegrationAtFrontVolume::Fe
struct	EshelbianPlasticity::SetIntegrationAtFrontFace
struct	EshelbianPlasticity::SetIntegrationAtFrontFace::Fe
struct	EshelbianPlasticity::VolRule
	Set integration rule on element. More...
struct	EshelbianPlasticity::FaceRule
struct	EshelbianPlasticity::EshelbianMonitor
struct	EshelbianPlasticity::SetUpSchurImpl
struct	EshelbianPlasticity::SetUpSchurImpl::P_MultiGridData
struct	EshelbianPlasticity::TSElasticPostStep
struct	EshelbianPlasticity::solve_elastic_setup

Namespaces
namespace	EshelbianPlasticity

Macros
#define	SINGULARITY

Functions
static auto	send_type (MoFEM::Interface &m_field, Range r, const EntityType type)
static auto	get_range_from_block (MoFEM::Interface &m_field, const std::string block_name, int dim)
static auto	get_range_from_block_map (MoFEM::Interface &m_field, const std::string block_name, int dim)
static auto	get_block_meshset (MoFEM::Interface &m_field, const int ms_id, const unsigned int cubit_bc_type)
static auto	save_range (moab::Interface &moab, const std::string name, const Range r, std::vector< Tag > tags={})
static auto	filter_true_skin (MoFEM::Interface &m_field, Range &&skin)
static auto	filter_owners (MoFEM::Interface &m_field, Range skin)
static auto	get_skin (MoFEM::Interface &m_field, Range body_ents)
static auto	get_crack_front_edges (MoFEM::Interface &m_field, Range crack_faces)
static auto	get_two_sides_of_crack_surface (MoFEM::Interface &m_field, Range crack_faces)

Detailed Description

Implementation of automatic differentiation.

Definition in file EshelbianPlasticity.cpp.

Macro Definition Documentation

SINGULARITY

#define SINGULARITY

Definition at line 8 of file EshelbianPlasticity.cpp.

Function Documentation

filter_owners()

static auto filter_owners ( MoFEM::Interface & m_field, Range skin )

static

Examples

EshelbianPlasticity.cpp.

Definition at line 177 of file EshelbianPlasticity.cpp.

                                                               {
  Range owner_ents;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  CHK_MOAB_THROW(pcomm->filter_pstatus(skin, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
                                       &owner_ents),
                 "filter_pstatus");
  return owner_ents;
};

filter_true_skin()

static auto filter_true_skin ( MoFEM::Interface & m_field, Range && skin )

static

Examples

EshelbianPlasticity.cpp, adolc_plasticity.cpp, plastic.cpp, seepage.cpp, tensor_divergence_operator.cpp, and thermo_elastic.cpp.

Definition at line 166 of file EshelbianPlasticity.cpp.

                                                                    {
  Range boundary_ents;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  CHK_MOAB_THROW(pcomm->filter_pstatus(skin,
                                       PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                       PSTATUS_NOT, -1, &boundary_ents),
                 "filter_pstatus");
  return boundary_ents;
};

get_block_meshset()

static auto get_block_meshset ( MoFEM::Interface & m_field, const int ms_id, const unsigned int cubit_bc_type )

static

Examples

EshelbianPlasticity.cpp.

Definition at line 144 of file EshelbianPlasticity.cpp.

                                                                {
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  EntityHandle meshset;
  CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
  return meshset;
};

get_crack_front_edges()

static auto get_crack_front_edges ( MoFEM::Interface & m_field, Range crack_faces )

static

Examples

EshelbianPlasticity.cpp.

Definition at line 194 of file EshelbianPlasticity.cpp.

                                                     {
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  auto &moab = m_field.get_moab();
  Range crack_skin_without_bdy;
  if (pcomm->rank() == 0) {
    Range crack_edges;
    CHKERR moab.get_adjacencies(crack_faces, 1, true, crack_edges,
                                moab::Interface::UNION);
    auto crack_skin = get_skin(m_field, crack_faces);
    Range body_ents;
    CHK_MOAB_THROW(
        m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents),
        "get_entities_by_dimension");
    auto body_skin = get_skin(m_field, body_ents);
    Range body_skin_edges;
    CHK_MOAB_THROW(moab.get_adjacencies(body_skin, 1, true, body_skin_edges,
                                        moab::Interface::UNION),
                   "get_adjacencies");
    crack_skin_without_bdy = subtract(crack_skin, body_skin_edges);
    auto front_edges_map = get_range_from_block_map(m_field, "FRONT", 1);
    for (auto &m : front_edges_map) {
      auto add_front = subtract(m.second, crack_edges);
      auto i = intersect(m.second, crack_edges);
      if (i.empty()) {
        crack_skin_without_bdy.merge(add_front);
      } else {
        auto i_skin = get_skin(m_field, i);
        Range adj_i_skin;
        CHKERR moab.get_adjacencies(i_skin, 1, true, adj_i_skin,
                                    moab::Interface::UNION);
        adj_i_skin = subtract(intersect(adj_i_skin, m.second), crack_edges);
        crack_skin_without_bdy.merge(adj_i_skin);
      }
    }
  }
  return send_type(m_field, crack_skin_without_bdy, MBEDGE);
}

get_range_from_block()

static auto get_range_from_block ( MoFEM::Interface & m_field, const std::string block_name, int dim )

static

Definition at line 100 of file EshelbianPlasticity.cpp.

                                                                      {
  Range r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r.merge(faces);
  }
 
  return r;
};

get_range_from_block_map()

static auto get_range_from_block_map ( MoFEM::Interface & m_field, const std::string block_name, int dim )

static

Examples

EshelbianPlasticity.cpp.

Definition at line 122 of file EshelbianPlasticity.cpp.

                                                                          {
  std::map<std::string, Range> r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r[bc->getName()] = faces;
  }
 
  return r;
}

get_skin()

static auto get_skin ( MoFEM::Interface & m_field, Range body_ents )

static

Definition at line 187 of file EshelbianPlasticity.cpp.

                                                               {
  Skinner skin(&m_field.get_moab());
  Range skin_ents;
  CHK_MOAB_THROW(skin.find_skin(0, body_ents, false, skin_ents), "find_skin");
  return skin_ents;
};

get_two_sides_of_crack_surface()

static auto get_two_sides_of_crack_surface ( MoFEM::Interface & m_field, Range crack_faces )

static

Examples

EshelbianPlasticity.cpp.

Definition at line 234 of file EshelbianPlasticity.cpp.

                                                              {
 
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
 
  MOFEM_LOG("EP", Sev::noisy) << "get_two_sides_of_crack_surface";
 
  if (!pcomm->rank()) {
 
    auto impl = [&](auto &saids) {
      MoFEMFunctionBegin;
 
      auto &moab = m_field.get_moab();
 
      auto get_adj = [&](auto e, auto dim) {
        Range adj;
        CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(
                           e, dim, true, adj, moab::Interface::UNION),
                       "get adj");
        return adj;
      };
 
      auto get_conn = [&](auto e) {
        Range conn;
        CHK_MOAB_THROW(m_field.get_moab().get_connectivity(e, conn, true),
                       "get connectivity");
        return conn;
      };
 
      constexpr bool debug = false;
      Range body_ents;
      CHKERR m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                          body_ents);
      auto body_skin = get_skin(m_field, body_ents);
      auto body_skin_edges = get_adj(body_skin, 1);
 
      auto crack_skin =
          subtract(get_skin(m_field, crack_faces), body_skin_edges);
      auto crack_skin_conn = get_conn(crack_skin);
      auto crack_skin_conn_edges = get_adj(crack_skin_conn, 1);
      auto crack_edges = get_adj(crack_faces, 1);
      crack_edges = subtract(crack_edges, crack_skin);
      auto all_tets = get_adj(crack_edges, 3);
      crack_edges = subtract(crack_edges, crack_skin_conn_edges);
      auto crack_conn = get_conn(crack_edges);
      all_tets.merge(get_adj(crack_conn, 3));
 
      if (debug) {
        CHKERR save_range(m_field.get_moab(), "crack_faces.vtk", crack_faces);
        CHKERR save_range(m_field.get_moab(), "all_crack_tets.vtk", all_tets);
        CHKERR save_range(m_field.get_moab(), "crack_edges_all.vtk",
                          crack_edges);
      }
 
      if (crack_faces.size()) {
        auto grow = [&](auto r) {
          auto crack_faces_conn = get_conn(crack_faces);
          Range v;
          auto size_r = 0;
          while (size_r != r.size() && r.size() > 0) {
            size_r = r.size();
            CHKERR moab.get_connectivity(r, v, true);
            v = subtract(v, crack_faces_conn);
            if (v.size()) {
              CHKERR moab.get_adjacencies(v, SPACE_DIM, true, r,
                                          moab::Interface::UNION);
              r = intersect(r, all_tets);
            }
            if (r.empty()) {
              break;
            }
          }
          return r;
        };
 
        Range all_tets_ord = all_tets;
        while (all_tets.size()) {
          Range faces = get_adj(unite(saids.first, saids.second), 2);
          faces = subtract(crack_faces, faces);
          if (faces.size()) {
            Range tets;
            auto fit = faces.begin();
            for (; fit != faces.end(); ++fit) {
              tets = intersect(get_adj(Range(*fit, *fit), 3), all_tets);
              if (tets.size() == 2) {
                break;
              }
            }
            if (tets.empty()) {
              break;
            } else {
              saids.first.insert(tets[0]);
              saids.first = grow(saids.first);
              all_tets = subtract(all_tets, saids.first);
              if (tets.size() == 2) {
                saids.second.insert(tets[1]);
                saids.second = grow(saids.second);
                all_tets = subtract(all_tets, saids.second);
              }
            }
          } else {
            break;
          }
        }
 
        saids.first = subtract(all_tets_ord, saids.second);
        saids.second = subtract(all_tets_ord, saids.first);
      }
 
      MoFEMFunctionReturn(0);
    };
 
    std::pair<Range, Range> saids;
    if (crack_faces.size())
      CHK_THROW_MESSAGE(impl(saids), "get crack both sides");
    return saids;
  }
 
  MOFEM_LOG("EP", Sev::noisy) << "get_two_sides_of_crack_surface <- done";
 
  return std::pair<Range, Range>();
}

save_range()

static auto save_range ( moab::Interface & moab, const std::string name, const Range r, std::vector< Tag > tags = {} )

static

Definition at line 152 of file EshelbianPlasticity.cpp.

                                                            {}) {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  if (r.size()) {
    CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1,
                           tags.data(), tags.size());
  } else {
    MOFEM_LOG("SELF", Sev::warning) << "Empty range for " << name;
  }
  MoFEMFunctionReturn(0);
};

send_type()

static auto send_type ( MoFEM::Interface & m_field, Range r, const EntityType type )

static

Examples

EshelbianPlasticity.cpp.

Definition at line 54 of file EshelbianPlasticity.cpp.

                                             {
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
 
  auto dim = CN::Dimension(type);
 
  std::vector<int> sendcounts(pcomm->size());
  std::vector<int> displs(pcomm->size());
  std::vector<int> sendbuf(r.size());
  if (pcomm->rank() == 0) {
    for (auto p = 1; p != pcomm->size(); p++) {
      auto part_ents = m_field.getInterface<CommInterface>()
                           ->getPartEntities(m_field.get_moab(), p)
                           .subset_by_dimension(SPACE_DIM);
      Range faces;
      CHKERR m_field.get_moab().get_adjacencies(part_ents, dim, true, faces,
                                                moab::Interface::UNION);
      faces = intersect(faces, r);
      sendcounts[p] = faces.size();
      displs[p] = sendbuf.size();
      for (auto f : faces) {
        auto id = id_from_handle(f);
        sendbuf.push_back(id);
      }
    }
  }
 
  int recv_data;
  MPI_Scatter(sendcounts.data(), 1, MPI_INT, &recv_data, 1, MPI_INT, 0,
              pcomm->comm());
  std::vector<int> recvbuf(recv_data);
  MPI_Scatterv(sendbuf.data(), sendcounts.data(), displs.data(), MPI_INT,
               recvbuf.data(), recv_data, MPI_INT, 0, pcomm->comm());
 
  if (pcomm->rank() > 0) {
    Range r;
    for (auto &f : recvbuf) {
      r.insert(ent_form_type_and_id(type, f));
    }
    return r;
  }
 
  return r;
}