gmpMechanic.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_MECHANIC_H_
#define _GEMA_PLUGIN_MECHANICAL_PHYSICS_MECHANIC_H_


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


template <class T> class GmpMechanic : public T
{
public:
  GmpMechanic(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& ncoord, 
              const GmMatrix& X, const GmVector& N, const GmMatrix& J, GmMatrix& Bu)
  {
    S_TRACE();
    assert(e); Q_UNUSED(X); Q_UNUSED(N);

    int d = GmpFemPhysicsCommon::nodeDim();
    double k, l, m;
    GmMatrix B;
    double detJ = shape->shapeCartesianPartialsFromJacobian(ncoord, J, B, true);

    assert(Bu.n_rows == 2*d && Bu.n_cols == e->numNodes()*d);
    Bu.zeros();
    // 2D solid element
    if (d == 2)
    {
      for (int i = 0; i < e->numNodes(); i++)
      {
        k = 2 * i + 1;
        l = 2 * i;
        Bu(0, l) = B(0, i);
        Bu(1, k) = B(1, i);
        Bu(3, l) = B(1, i);
        Bu(3, k) = B(0, i);
      }
    }
    // 3D solid element
    else if (d == 3)
    {
      for (int i = 0; i < e->numNodes(); i++)
      {
        m = 3 * i + 2;
        k = m - 1;
        l = m - 2;
        Bu(0, l) = B(0, i);
        Bu(1, k) = B(1, i);
        Bu(2, m) = B(2, i);
        Bu(3, l) = B(1, i);
        Bu(3, k) = B(0, i);
        Bu(4, l) = B(2, i);
        Bu(4, m) = B(0, i);
        Bu(5, k) = B(2, i);
        Bu(5, m) = B(1, i);
      }
    }
    return detJ;
  }
};

template <class T> class GmpMechanicPlaneStress : public T
{
public:
  GmpMechanicPlaneStress(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& ncoord, const GmMatrix& X, const GmVector& N, const GmMatrix& J, GmMatrix& Bu)
  {
    S_TRACE();
    assert(e); Q_UNUSED(X); Q_UNUSED(N);
        
    double k, l;
    GmMatrix B;
    double detJ = shape->shapeCartesianPartialsFromJacobian(ncoord, J, B, true);

    assert(Bu.n_rows == 2 * nodeDim() && Bu.n_cols == e->numNodes()*nodeDim());
    Bu.zeros();

    // 2D solid element
    assert(nodeDim() == 2);
    for (int i = 0; i < e->numNodes(); i++)
    {
      k = 2 * i + 1;
      l = 2 * i;
      Bu(0, l) = B(0, i);
      Bu(1, k) = B(1, i);
      Bu(3, l) = B(1, i);
      Bu(3, k) = B(0, i);
    }
    return detJ;
  }
  // Comments on the base class
  virtual bool isPlaneStress() { return true; }
};

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

protected:
  virtual bool isAxisymmetric() { return true; }

  // Comments on the base class
  virtual double fillBuMatrix(const GmElement* e, const GmShape* shape, const GmVector& ncoord, const GmMatrix& X, const GmVector& N, const GmMatrix& J, GmMatrix& Bu)
  {
    S_TRACE();
    assert(e);

    int n = e->numNodes();
    double k, l;
    GmMatrix B;
    double detJ = shape->shapeCartesianPartialsFromJacobian(ncoord, J, B, true);

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

    assert(Bu.n_rows == 2 * nodeDim() && Bu.n_cols == n*nodeDim());
    Bu.zeros();
    // 2D solid element
    assert(nodeDim() == 2);
    for (int i = 0; i < e->numNodes(); i++)
    {
      k = 2 * i + 1;
      l = 2 * i;
      Bu(0, l) = B(0, i);
      Bu(1, k) = B(1, i);
      Bu(2, l) = N(i) / rGauss; 
      Bu(3, l) = B(1, i);
      Bu(3, k) = B(0, i);
    }
    
    return detJ;
  }

  // Comments on the base class
  virtual double axisymmetricFactor(const GmElement* e, const GmMatrix& X, const GmVector& N)
  { 
    int n = e->numNodes();

    double rGauss = 0.0;
    for (int j = 0; j < n; ++j)
      rGauss += N(j)*X(j, 0);
        return 2 * arma::datum::pi * rGauss;
  }

};

#endif