ElasticTie::SetIntegrationOnTiePrismFaces Struct Reference

#include "tutorials/vec-11_elastic_tie_mesh/src/ElasticTie.hpp"

Collaboration diagram for ElasticTie::SetIntegrationOnTiePrismFaces:

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

Public Member Functions
	SetIntegrationOnTiePrismFaces ()=delete
	SetIntegrationOnTiePrismFaces (int tie_ref_level=1)
	SetIntegrationOnTiePrismFaces (FunRule fun_rule, int tie_ref_level=1)
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; }
int	tieRefLevel = 1

Static Private Member Functions
static MoFEMErrorCode	setBaseQuadrature (FlatPrismElementForcesAndSourcesCore &fe, const int rule)
static MoFEMErrorCode	projectPointToMasterPlane (const double *master_coords, const double *point, double *projected_point)
static MoFEMErrorCode	projectSlavePointsToMasterPlane (const double *master_coords, const MatrixDouble &slave_global_coords, MatrixDouble &master_global_coords)
static MoFEMErrorCode	refineQuadrature (FlatPrismElementForcesAndSourcesCore &fe, const int refinement_levels, MatrixDouble &ref_gauss_pts)

Static Private Attributes
static std::map< std::pair< int, int >, MatrixDouble >	mapRefCoords

Detailed Description

Definition at line 237 of file ElasticTie.hpp.

Member Typedef Documentation

FunRule

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

Definition at line 239 of file ElasticTie.hpp.

Constructor & Destructor Documentation

SetIntegrationOnTiePrismFaces() [1/3]

ElasticTie::SetIntegrationOnTiePrismFaces::SetIntegrationOnTiePrismFaces ( )

delete

SetIntegrationOnTiePrismFaces() [2/3]

ElasticTie::SetIntegrationOnTiePrismFaces::SetIntegrationOnTiePrismFaces ( int  tie_ref_level = 1 )

inline

Definition at line 245 of file ElasticTie.hpp.

      : tieRefLevel(tie_ref_level) {}

SetIntegrationOnTiePrismFaces() [3/3]

ElasticTie::SetIntegrationOnTiePrismFaces::SetIntegrationOnTiePrismFaces ( FunRule  fun_rule, int  tie_ref_level = 1  )

inline

Definition at line 248 of file ElasticTie.hpp.

      : funRule(fun_rule), tieRefLevel(tie_ref_level) {}

Member Function Documentation

operator()()

MoFEMErrorCode ElasticTie::SetIntegrationOnTiePrismFaces::operator() ( ForcesAndSourcesCore *  fe_raw_ptr, int  order_row, int  order_col, int  order_data  )

inline

Definition at line 251 of file ElasticTie.hpp.

                                                           {
    MoFEMFunctionBegin;
 
    auto fe_ptr = static_cast<FlatPrismElementForcesAndSourcesCore *>(fe_raw_ptr);
    const int rule = funRule(std::max(order_data, std::max(order_row, order_col)));
 
    if (tieRefLevel < 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "-tie_ref_level has to be >= 0");
    }
 
    if (rule < 0 || rule >= QUAD_2D_TABLE_SIZE) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Triangle quadrature rule %d is not available for TIE prism face "
              "integration",
              rule);
    }
 
    if (QUAD_2D_TABLE[rule]->dim != 2) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Expected 2D quadrature for TIE prism face integration");
    }
 
    if (QUAD_2D_TABLE[rule]->order < rule) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong quadrature order %d < %d for TIE prism face integration",
              QUAD_2D_TABLE[rule]->order, rule);
    }
 
    CHKERR setBaseQuadrature(*fe_ptr, rule);
 
    if (tieRefLevel > 0) {
      const auto cache_key = std::make_pair(rule, tieRefLevel);
      auto it = mapRefCoords.find(cache_key);
      if (it == mapRefCoords.end()) {
        MatrixDouble ref_gauss_pts;
        CHKERR refineQuadrature(*fe_ptr, tieRefLevel, ref_gauss_pts);
        mapRefCoords[cache_key].swap(ref_gauss_pts);
        it = mapRefCoords.find(cache_key);
      }
      fe_ptr->gaussPts = it->second;
    }
 
    const MatrixDouble slave_gauss_pts = fe_ptr->gaussPts;
    const auto nb_gauss_pts = slave_gauss_pts.size2();
 
    auto &moab = fe_ptr->mField.get_moab();
    auto get_prism_face_coords = [&](const int side, std::array<double, 9> &coords) {
      MoFEMFunctionBegin;
      const EntityHandle prism = fe_ptr->getFEEntityHandle();
      const EntityHandle *prism_conn = nullptr;
      int num_prism_nodes = 0;
      CHKERR moab.get_connectivity(prism, prism_conn, num_prism_nodes, true);
      if (num_prism_nodes != 6) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "TIE prism is expected to have 6 nodes");
      }
 
      std::array<EntityHandle, 3> face_conn;
      switch (side) {
      case SLAVE_FACE_SIDE:
        face_conn = {prism_conn[0], prism_conn[1], prism_conn[2]};
        break;
      case MASTER_FACE_SIDE:
        face_conn = {prism_conn[3], prism_conn[4], prism_conn[5]};
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Unsupported TIE prism triangle side %d", side);
      }
 
      CHKERR moab.get_coords(face_conn.data(), 3, coords.data());
      MoFEMFunctionReturn(0);
    };
 
    std::array<double, 9> slave_coords; // nodal coords
    std::array<double, 9> master_coords; // nodal coords
    CHKERR get_prism_face_coords(SLAVE_FACE_SIDE, slave_coords);
    CHKERR get_prism_face_coords(MASTER_FACE_SIDE, master_coords);
 
    MatrixDouble slave_global_coords(nb_gauss_pts, 3, false); // gp coords
    MatrixDouble master_global_coords(nb_gauss_pts, 3, false); //gp coords
    MatrixDouble slave_shape(nb_gauss_pts, 3, false);
    CHKERR ShapeMBTRI(&slave_shape(0, 0), &slave_gauss_pts(0, 0),
                      &slave_gauss_pts(1, 0), nb_gauss_pts);
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      for (int dd = 0; dd != 3; ++dd) {
        slave_global_coords(gg, dd) =
            slave_shape(gg, 0) * slave_coords[0 + dd] +
            slave_shape(gg, 1) * slave_coords[3 + dd] +
            slave_shape(gg, 2) * slave_coords[6 + dd];
      }
    }
 
    CHKERR projectSlavePointsToMasterPlane(master_coords.data(),
                                           slave_global_coords,
                                           master_global_coords);
 
    MatrixDouble master_local_coords(nb_gauss_pts, 2, false);
    CHKERR Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
        master_coords.data(), &master_global_coords(0, 0),
        master_global_coords.size1(), &master_local_coords(0, 0));
 
    MatrixDouble combined_gauss_pts(3, 2 * nb_gauss_pts, false);
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      const bool valid_master_point =
          isInsideReferenceTriangle(master_local_coords(gg, 0),
                                    master_local_coords(gg, 1));
      const double paired_weight = valid_master_point ? slave_gauss_pts(2, gg) : 0.;
 
      combined_gauss_pts(0, gg) = slave_gauss_pts(0, gg);
      combined_gauss_pts(1, gg) = slave_gauss_pts(1, gg);
      combined_gauss_pts(2, gg) = paired_weight;
 
      combined_gauss_pts(0, gg + nb_gauss_pts) = master_local_coords(gg, 0);
      combined_gauss_pts(1, gg + nb_gauss_pts) = master_local_coords(gg, 1);
      combined_gauss_pts(2, gg + nb_gauss_pts) = paired_weight;
    }
    fe_ptr->gaussPts.swap(combined_gauss_pts);
 
    MoFEMFunctionReturn(0);
  }

projectPointToMasterPlane()

static MoFEMErrorCode ElasticTie::SetIntegrationOnTiePrismFaces::projectPointToMasterPlane ( const double *  master_coords, const double *  point, double *  projected_point  )

inlinestaticprivate

Definition at line 390 of file ElasticTie.hpp.

                                                                           {
    MoFEMFunctionBegin;
 
    // For the infinite master plane, the closest point is the orthogonal
    // projection along the plane normal.
    const double ax = master_coords[0];
    const double ay = master_coords[1];
    const double az = master_coords[2];
    const double bx = master_coords[3];
    const double by = master_coords[4];
    const double bz = master_coords[5];
    const double cx = master_coords[6];
    const double cy = master_coords[7];
    const double cz = master_coords[8];
 
    const double abx = bx - ax;
    const double aby = by - ay;
    const double abz = bz - az;
    const double acx = cx - ax;
    const double acy = cy - ay;
    const double acz = cz - az;
 
    const double nx = aby * acz - abz * acy;
    const double ny = abz * acx - abx * acz;
    const double nz = abx * acy - aby * acx;
    const double n_norm_sq = nx * nx + ny * ny + nz * nz;
    if (n_norm_sq <= std::numeric_limits<double>::epsilon()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Degenerated master triangle plane");
    }
 
    const double px = point[0] - ax;
    const double py = point[1] - ay;
    const double pz = point[2] - az;
    const double signed_distance = (px * nx + py * ny + pz * nz) / n_norm_sq;
    projected_point[0] = point[0] - signed_distance * nx;
    projected_point[1] = point[1] - signed_distance * ny;
    projected_point[2] = point[2] - signed_distance * nz;
 
    MoFEMFunctionReturn(0);
  }

projectSlavePointsToMasterPlane()

static MoFEMErrorCode ElasticTie::SetIntegrationOnTiePrismFaces::projectSlavePointsToMasterPlane ( const double *  master_coords, const MatrixDouble &  slave_global_coords, MatrixDouble &  master_global_coords  )

inlinestaticprivate

Definition at line 434 of file ElasticTie.hpp.

                                          {
    MoFEMFunctionBegin;
 
    // The weak form is integrated on the slave surface. Each slave Gauss point
    // is paired with its orthogonal projection on the infinite plane defined
    // by the master triangle.
    master_global_coords.resize(slave_global_coords.size1(), 3, false);
    for (size_t gg = 0; gg != slave_global_coords.size1(); ++gg) {
      CHKERR projectPointToMasterPlane(master_coords,
                                       &slave_global_coords(gg, 0),
                                       &master_global_coords(gg, 0));
    }
 
    MoFEMFunctionReturn(0);
  }

refineQuadrature()

static MoFEMErrorCode ElasticTie::SetIntegrationOnTiePrismFaces::refineQuadrature ( FlatPrismElementForcesAndSourcesCore &  fe, const int  refinement_levels, MatrixDouble &  ref_gauss_pts  )

inlinestaticprivate

Definition at line 452 of file ElasticTie.hpp.

                                   {
    MoFEMFunctionBegin;
 
    constexpr int num_nodes = 3;
    moab::Core moab_ref;
    double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
    EntityHandle nodes[num_nodes];
    for (int nn = 0; nn != num_nodes; ++nn)
      CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
 
    EntityHandle tri;
    CHKERR moab_ref.create_element(MBTRI, nodes, num_nodes, tri);
    MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
    MoFEM::Interface &m_field_ref = mofem_ref_core;
 
    {
      Range tris(tri, tri);
      Range edges;
      // Force MOAB to materialise the triangle edges before the bit-ref
      // database is seeded. The working reference-refinement snippets all do
      // this, and the side-number lookup below verifies that MOAB sees a
      // consistent local edge numbering on the temporary reference triangle.
      CHKERR m_field_ref.get_moab().get_adjacencies(
          tris, 1, true, edges, moab::Interface::UNION);
      for (auto edge : edges) {
        int side_number = -1, sense = 0, offset = 0;
        CHKERR moab_ref.side_number(tri, edge, side_number, sense, offset);
      }
      CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
          tris, BitRefLevel().set(0), false, VERBOSE);
    }
 
    auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
    for (int ll = 0; ll != refinement_levels; ++ll) {
      Range tris;
      CHKERR m_field_ref.getInterface<BitRefManager>()
          ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                         BitRefLevel().set(), MBTRI, tris);
 
      Range edges;
      CHKERR moab_ref.get_adjacencies(tris, 1, true, edges,
                                      moab::Interface::UNION);
 
      const BitRefLevel child_bit = BitRefLevel().set(ll + 1);
      CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges, child_bit);
      CHKERR m_ref->refineTris(tris, child_bit, QUIET, false);
    }
 
    Range tris;
    CHKERR m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(refinement_levels), BitRefLevel().set(), MBTRI, tris);
 
    MatrixDouble ref_coords(tris.size(), 9, false);
    int tt = 0;
    for (auto tit = tris.begin(); tit != tris.end(); ++tit, ++tt) {
      int num_conn;
      const EntityHandle *conn;
      CHKERR moab_ref.get_connectivity(*tit, conn, num_conn, false);
      CHKERR moab_ref.get_coords(conn, num_conn, &ref_coords(tt, 0));
    }
 
    const size_t nb_gauss_pts = fe.gaussPts.size2();
    MatrixDouble shape_n(nb_gauss_pts, 3, false);
    CHKERR ShapeMBTRI(&shape_n(0, 0), &fe.gaussPts(0, 0), &fe.gaussPts(1, 0),
                      nb_gauss_pts);
    ref_gauss_pts.resize(3, nb_gauss_pts * ref_coords.size1(), false);
    int gg = 0;
    for (size_t rr = 0; rr != ref_coords.size1(); ++rr) {
      double *tri_coords = &ref_coords(rr, 0);
      FTensor::Tensor1<double, 3> t_normal;
      CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
      const double det = t_normal.l2();
      for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
        for (int dd = 0; dd != 2; ++dd) {
          ref_gauss_pts(dd, gg) =
              shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
              shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
              shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
        }
        ref_gauss_pts(2, gg) = fe.gaussPts(2, ggg) * det;
      }
    }
 
    MoFEMFunctionReturn(0);
  }

setBaseQuadrature()

static MoFEMErrorCode ElasticTie::SetIntegrationOnTiePrismFaces::setBaseQuadrature ( FlatPrismElementForcesAndSourcesCore &  fe, const int  rule  )

inlinestaticprivate

Definition at line 376 of file ElasticTie.hpp.

                                                                {
    MoFEMFunctionBegin;
    const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
    fe.gaussPts.resize(3, nb_gauss_pts, false);
    cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
                &fe.gaussPts(0, 0), 1);
    cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
                &fe.gaussPts(1, 0), 1);
    cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
                &fe.gaussPts(2, 0), 1);
    MoFEMFunctionReturn(0);
  }

Member Data Documentation

funRule

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

Definition at line 240 of file ElasticTie.hpp.

{ return 2 * p + 1; };

mapRefCoords

std::map<std::pair<int, int>, MatrixDouble> ElasticTie::SetIntegrationOnTiePrismFaces::mapRefCoords

inlinestaticprivate

Definition at line 540 of file ElasticTie.hpp.

tieRefLevel

int ElasticTie::SetIntegrationOnTiePrismFaces::tieRefLevel = 1

Definition at line 241 of file ElasticTie.hpp.

---

The documentation for this struct was generated from the following file:

• tutorials/vec-11_elastic_tie_mesh/src/ElasticTie.hpp