gmpInterface.h

Go to the documentation of this file.

/************************************************************************
**
** Copyright (C) 2014 by Carlos Augusto Teixera Mendes
** All rights reserved.
**
** This file is part of the "GeMA" software. It's use should respect
** the terms in the license agreement that can be found together
** with this source code.
** It is provided AS IS, with NO WARRANTY OF ANY KIND,
** INCLUDING THE WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR
** A PARTICULAR PURPOSE.
**
************************************************************************/

#ifndef _GEMA_PLUGIN_MECHANICAL_PHYSICS_INTERFACE_H_
#define _GEMA_PLUGIN_MECHANICAL_PHYSICS_INTERFACE_H_


#include "gmpMechanicalInterface.h"
#include <gmTrace.h>
#include <math.h>


template <class T> class GmpInterface : public T
{
public:
  GmpInterface(const char* pluginType, GmSimulationData* simulation, QString id, QString description,
    const GmpFemPhysicsCommonMaterialFactory* matFactory, const GmLogCategory& logger)
    : T(pluginType, simulation, id, description, matFactory, logger)
  {}

protected:
  // Comments on the base class
  virtual double fillBuMatrix(const GmElement* e, const GmShape* shape, const GmVector& ip, const GmMatrix& MX, const GmVector& N, const GmMatrix& J, GmMatrix& Bu)
  {
    S_TRACE();
    assert(e); Q_UNUSED(MX); Q_UNUSED(N);
    int cNodes, inNode, nv;
    // Node Number and intermediate node number for quadratic interface element
    int n = e->numNodes();
    int d = GmpFemPhysicsCommon::nodeDim();
    if (e->type() == GM_INT2DL6)
    {
      n = n - pow(2, d - 1);
    }

    // Compute the number of intermediate nodes and corner nodes

    //int cNodes = pow(2, d - 1);
    if (e->type() == GM_INT3DL6 || e->type() == GM_INT3DQ12)
    {
      inNode = (n - 3.0 * (d - 1)) / 2.0;
    }
    else
    {
      inNode = (n - pow(2.0, d)) / 2.0;
    }

    cNodes = 0.5 * n - inNode;

    double detJ;
    GmVector H;
    GmMatrix R(d,d);

    assert(Bu.n_rows == d && Bu.n_cols == n*d);
    // read the shapes values 
    shape->shapeValues(ip, H);
    // 2D interface element
    if (d == 2)
    {
      // compute determinant of the jacobian matrix
      detJ = shape->scaledJacobianDet(J);
      // fill the 2D rotation matrix
      R(0, 0) = J(0, 0) / detJ;
      R(0, 1) = J(0, 1) / detJ;
      R(1, 0) = -R(0, 1);
      R(1, 1) = R(0, 0);
    }
    // 3D interface element
    else if (d == 3)
    {
      double J1 = 0.0, J2 = 0.0, J3 = 0.0;

      J1 = J(0, 1)*J(1, 2) - J(1, 1)*J(0, 2);
      J2 = J(1, 0)*J(0, 2) - J(0, 0)*J(1, 2);
      J3 = J(0, 0)*J(1, 1) - J(1, 0)*J(0, 1);

      // determinant of the Jacobian matrix
      detJ = sqrt(J1*J1 + J2*J2 + J3*J3);        
      
      
      // fill the 3D rotation matrix
      //           | Vs |      V_xi / |V_xi|
      // _R(3x3) = | Vt |  =  (Vn x Vs) / |(Vn x Vs)|
      //           | Vn |     (V_xi x V_Eta) / |(V_xi x V_Eta)|
      double detJ1, detJ2;

      //Tangential Vs at shear direction 1
      detJ1 = sqrt(J(0, 0)*J(0, 0) + J(0, 1)*J(0, 1) + J(0, 2)*J(0, 2));
      R(0, 0) = J(0, 0) / detJ1;       // Vs(0,0)
      R(0, 1) = J(0, 1) / detJ1;       // Vs(0,1)
      R(0, 2) = J(0, 2) / detJ1;       // Vs(0,2)
      //Normal Vn at normal direction 3   
      detJ2 = sqrt(J1*J1 + J2* J2 + J3* J3);
      R(2, 0) = J1 / detJ2;            // Vn(0,0)
      R(2, 1) = J2 / detJ2;            // Vn(0,1)
      R(2, 2) = J3 / detJ2;            // Vn(0,2)
      //Tangential Vt at shear direction 2 
      R(1, 0) = R(2, 1)*R(0, 2) - R(0, 1)*R(2, 2);  // Vt(0,0)
      R(1, 1) = R(0, 0)*R(2, 2) - R(2, 0)*R(0, 2);  // Vt(0,1)
      R(1, 2) = R(2, 0)*R(0, 1) - R(0, 0)*R(2, 1);  // Vt(0,2)
      //*/
      //
# if 0
      double dipdir, dip2, dip, strike;
      GmVector N2(3);
      // unit vector at normal direction
      N2(0) =  J3 / detJ;           // Un(0)
      N2(1) = -J2 / detJ;           // Un(1)
      N2(2) =  J1 / detJ;           // Un(2)

      dipdir =angleAtanD (N2(0),N2(1));

      if (N2(2) > 0.00)
      {
        dip2 = 90.0 - angleAtanD(abs(N2(2)), sqrt(N2(0)*N2(0) + N2(1)*N2(1)));
      }     
      else
      {
        dip2 = 90.00 + angleAtanD(abs(N2(2)),sqrt(N2(0)*N2(0)+N2(1)*N2(1)));
      }

      dip = -M_PI*(180 + dip2) / 180.0;
      strike = M_PI*(90+dipdir)/180;

      R(0, 0) = sin(strike);            R(0, 1) = cos(strike);            R(0, 2) = 0;
      R(1, 0) = cos(dip)*cos(strike);   R(1, 1) = -cos(dip)*sin(strike);  R(1, 2) = -sin(dip);
      R(2, 0) = -sin(dip)*cos(strike);  R(2, 1) = sin(dip)*sin(strike);   R(2, 2) = -cos(dip);
# endif

  # if 0
      double easting, northing, dip, strike;
      GmVector N2(3);

      // unit vector at normal direction
      N2(0) = J3 / detJ;            // Un(0)
      N2(1) = -J2 / detJ;           // Un(1)
      N2(2) = J1 / detJ;            // Un(2)
      //N(0) = -0.7715;
      //N(1) =  0.6172;
      //N(2) = -0.1543;

      easting = (N2(2) < 0) ? N2(1) : -N2(1);

      northing = (N2(2) > 0) ? N2(0) : -N2(0);

      dip = abs(asin(sqrt(N2(0)*N2(0) + N2(1)*N2(1))));

      strike = acos(northing / sqrt(easting*easting + northing*northing));

      strike = (easting>= 0.0)? strike : (2 * M_PI - strike);
   
      R(0, 0) = sin(strike);            R(0, 1) = cos(strike);            R(0, 2) = 0;
      R(1, 0) = cos(dip)*cos(strike);   R(1, 1) = -cos(dip)*sin(strike);  R(1, 2) = -sin(dip);
      R(2, 0) = -sin(dip)*cos(strike);  R(2, 1) = sin(dip)*sin(strike);   R(2, 2) = -cos(dip);    
      //R = R.t();
 #endif
    }
    else
    {
      gmErrorMsg(logger(), QObject::tr("gmpInterface: Unsupported condition in fillStrainDisplacementMatrix"));
    }

    // Fill matrix B
    // fill the matrix B with zeros.
    Bu.zeros();
    nv = 2 * cNodes;
    // fill the matrix B for corner nodes
    for (int i = 0; i < d; ++i)
    {
      for (int j = 0; j < cNodes; ++j)
      {
        if (d == 2)
        {
          // Bottom nodes 
          Bu(i, j * 2 + i) = -H(j);
          // Top nodes 
          Bu(i, (j + 2) * 2 + i) = H(1 - j);
        }
        else if (d == 3)
        {
          // Bottom nodes 
          Bu(i, j * 3 + i) = -H(j);
          // Top nodes 
          //Bu(i, (j + 4) * 3 + i) = H(j);
          Bu(i, (j + cNodes) * 3 + i) = H(j);
        }
        else
        {
          gmErrorMsg(logger(), QObject::tr("GmpStdCohesive: Unsupported condition to fill B matrix."));
        }
      }
      // fill the matrix B for intermediate nodes
      for (int j = 0; j < inNode; ++j)
      {
        if (d == 2)
        {
          // Bottom nodes 
          Bu(i, (j + 4) * 2 + i) = -H(j + 2);
          // Top nodes 
          Bu(i, (j + 4 + inNode) * 2 + i) = H(1 + inNode - j);
        }
        else if (d == 3)
        {
          // Bottom nodes 
          Bu(i, (nv + j) * 3 + i) = -H(cNodes + j);
          // Top nodes 
          Bu(i, (nv + inNode + j) * 3 + i) = H(cNodes + j);
        }
        else
        {
          gmErrorMsg(logger(), QObject::tr("GmpInterface: Unsupported condition to fill B matrix."));
        }
      }

    }

    // Produto da Bu = R*Bu
    Bu = R*Bu;
   
    return detJ;

  }
    
  //
  double angleAtanD(double y, double x)
  {
    double PI, GRAUS, angle;

    PI = acos(-1.00);
    GRAUS = 180.0 / PI;

    if (x == 0.00)
    {
      if (y > 0.00)
      {
        // I and II
          angle = 90.0;
      }
      else if (y < 0.00)
      {
        // III and IV
          angle = -90.0;
      }
      else if (y == 0.0)
      {
        // Not possible
      }
    }
    else if (x>0.00)
    {
      if (y > 0.00)
      {
        // I(0<ANGLE<90)
          angle = atan(y / x)*GRAUS;
      }      
      else if (y < 0.00)
      {
        // IV(-90<ANGLE<0)
          angle = atan(y / x)*GRAUS;
      }       
      else if (y == 0.00)
      {
        // I and IV
          angle = 0.00;
      }
    }
    else if (x<0.00)
    {
      if (y>0.00)
      {
        // II(90<ANGLE<180)
          angle = 180.00 - abs(atan(y / x)*GRAUS);
      }       
      else if (y < 0.00)
      {
        // III(-180<ANGLE<-90)
          angle = abs(atan(y / x)*GRAUS) - 180.00;
      }     
      else if (y == 0.00)
      {
        // II and III
          angle = 180.00;
      }
    }
    return angle;
  } 
};
// Axisymmetric interface element
template <class T> class GmpAxiInterface : public T
{
public:
  GmpAxiInterface(const char* pluginType, GmSimulationData* simulation, QString id, QString description,
    const GmpFemPhysicsCommonMaterialFactory* matFactory, const GmLogCategory& logger)
    : T(pluginType, simulation, id, description, matFactory, logger)
  {}

  virtual double fillBuMatrix(const GmElement* e, const GmShape* shape, const GmVector& ip, const GmMatrix& MX, const GmVector& N, const GmMatrix& J, GmMatrix& Bu)
  {
    S_TRACE();
    assert(e); Q_UNUSED(MX); Q_UNUSED(N);
    int cNodes, inNode, nv;
    // Node Number and intermediate node number for quadratic interface element
    int n = e->numNodes();
    int d = GmpFemPhysicsCommon::nodeDim();
    assert(d == 2);
    if (e->type() == GM_INT2DL6)
    {
      n = n - pow(2, d - 1);
    }
    // Compute the number of intermediate nodes and corner nodes
    inNode = (n - pow(2.0, d)) / 2.0;
    cNodes = 0.5 * n - inNode;

    double detJ;
    GmVector H;
    GmMatrix R(d, d);

    assert(Bu.n_rows == d && Bu.n_cols == n * d);
    // read the shapes values
    shape->shapeValues(ip, H);
    // 2D interface element
    if (d == 2)
    {
      // compute determinant of the jacobian matrix
      detJ = shape->scaledJacobianDet(J);
      // fill the 2D rotation matrix
      R(0, 0) = J(0, 0) / detJ;
      R(0, 1) = J(0, 1) / detJ;
      R(1, 0) = -R(0, 1);
      R(1, 1) = R(0, 0);
    }
    // 3D interface element
    else
    {
      gmErrorMsg(logger(), QObject::tr("gmpInterface: Axisymmetric interface element must be 2D model"));
    }

    // Fill matrix B
    // fill the matrix B with zeros.
    Bu.zeros();
    nv = 2 * cNodes;
    // fill the matrix B for corner nodes
    for (int i = 0; i < d; ++i)
    {
      for (int j = 0; j < cNodes; ++j)
      {
        if (d == 2)
        {
          // Bottom nodes 
          Bu(i, j * 2 + i) = -H(j);
          // Top nodes 
          Bu(i, (j + 2) * 2 + i) = H(1 - j);
        }
        else
        {
          gmErrorMsg(logger(), QObject::tr("GmpStdCohesive: Unsupported condition to fill B matrix."));
        }
      }
      // fill the matrix B for intermediate nodes
      for (int j = 0; j < inNode; ++j)
      {
        if (d == 2)
        {
          // Bottom nodes 
          Bu(i, (j + 4) * 2 + i) = -H(j + 2);
          // Top nodes 
          Bu(i, (j + 4 + inNode) * 2 + i) = H(1 + inNode - j);
        }
        else
        {
          gmErrorMsg(logger(), QObject::tr("GmpInterface: Unsupported condition to fill B matrix."));
        }
      }

    }

    // Produto da Bu = R*Bu
    Bu = R * Bu;

    return detJ;;
  
  }

protected:
  // Evaluates the axisymmetric factor
  virtual double axisymmetricFactor(const GmElement* e, const GmMatrix& MX, const GmVector& N)
  {
    S_TRACE();
    assert(e);
    // number of nodes
    int n = e->numNodes();
    int d = GmpFemPhysicsCommon::nodeDim();
    // number of nodes for tri-nodded interface elements
    if (e->type() == GM_INT2DL6)
    {
      n = n - pow(2, d - 1);
    }
    assert(MX.n_rows == n / 2 && MX.n_cols == d);

    double rGauss = 0.0;
    for (int j = 0; j < n/2; ++j)
      rGauss += N(j)*MX(j, 0);

    return 2 * arma::datum::pi * rGauss;
  }

  virtual bool isAxisymmetric() { return true; }

};
#endif