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();
    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++)
      {
        int k = 2 * i + 1;
        int 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)
    {
      constexpr int bstride = 3;
      constexpr int bustride = 6;
      const double* bptr = B.memptr();
      double * buptr = Bu.memptr();
      for (int i = 0, n = e->numNodes(); i < n; i++)
      {
        int m = 3 * i + 2;
        int k = m - 1;
        int l = m - 2;
 
        const double* bi = bptr + i*bstride;
        double b0 = bi[0], b1 = bi[1], b2 = bi[2];
        /*
        Bu(0, l) = B(0, i);
        Bu(3, l) = B(1, i);
        Bu(4, l) = B(2, i);
        */
        double* bu = buptr + bustride * l;
        bu[0] = b0;
        bu[3] = b1;
        bu[4] = b2;
 
        /*
        Bu(1, k) = B(1, i);
        Bu(3, k) = B(0, i);
        Bu(5, k) = B(2, i);
        */
 
        bu = buptr + bustride * k;
        bu[1] = b1;
        bu[3] = b0;
        bu[5] = b2;
 
        /*
        Bu(2, m) = B(2, i);
        Bu(4, m) = B(0, i);
        Bu(5, m) = B(1, i);
        */
 
        bu = buptr + bustride * m;
        bu[2] = b2;
        bu[4] = b0;
        bu[5] = b1;
      }
    }
    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) const
  { 
    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