shell_cross_section.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosSolidMechanicsApplication $
//   Created by:          $Author:       Massimo Petracca $
//   Last modified by:    $Co-Author:                     $
//   Date:                $Date:             January 2014 $
//   Revision:            $Revision:                  0.0 $
//
//
 
#if !defined(SHELL_CROSS_SECTION_H_INCLUDED)
#define SHELL_CROSS_SECTION_H_INCLUDED
 
 
// System includes
#include <string>
#include <iostream>
 
// Project includes
#include "includes/define.h"
#include "includes/serializer.h"
#include "includes/constitutive_law.h"
#include "containers/flags.h"
#include "custom_utilities/properties_extensions.hpp"
#include "custom_utilities/solid_mechanics_math_utilities.hpp"
 
namespace Kratos
{
 
class KRATOS_API(SOLID_MECHANICS_APPLICATION) ShellCrossSection : public Flags
{
 
public:
 
    class Ply;
 
    KRATOS_CLASS_POINTER_DEFINITION(ShellCrossSection);
 
    typedef Geometry<Node > GeometryType;
 
    typedef std::vector< Ply > PlyCollection;
 
    typedef std::size_t SizeType;
 
        typedef SolidMechanicsMathUtilities<double> MathUtilsType;
 
 
    enum SectionBehaviorType
    {
        Thick, 
        Thin 
    };
 
 
 
    struct Features
    {
        Flags mOptions;
        double mStrainSize;
        double mSpaceDimension;
        std::vector< ConstitutiveLaw::StrainMeasure > mStrainMeasures;
    };
 
    class Parameters
    {
 
    private:
 
        Flags                mOptions;
 
        Vector*              mpGeneralizedStrainVector;
        Vector*              mpGeneralizedStressVector;
        Matrix*              mpConstitutiveMatrix;
 
        const Vector*        mpShapeFunctionsValues;
        const Matrix*        mpShapeFunctionsDerivatives;
        const ProcessInfo*   mpCurrentProcessInfo;
        const Properties*    mpMaterialProperties;
        const GeometryType*  mpElementGeometry;
 
    public:
 
        Parameters()
            : mpGeneralizedStrainVector(NULL)
            , mpGeneralizedStressVector(NULL)
            , mpConstitutiveMatrix(NULL)
            , mpShapeFunctionsValues(NULL)
            , mpShapeFunctionsDerivatives(NULL)
            , mpCurrentProcessInfo(NULL)
            , mpMaterialProperties(NULL)
            , mpElementGeometry(NULL)
        {}
 
        Parameters (const GeometryType& rElementGeometry,
                    const Properties& rMaterialProperties,
                    const ProcessInfo& rCurrentProcessInfo)
            : mpGeneralizedStrainVector(NULL)
            , mpGeneralizedStressVector(NULL)
            , mpConstitutiveMatrix(NULL)
            , mpShapeFunctionsValues(NULL)
            , mpShapeFunctionsDerivatives(NULL)
            , mpCurrentProcessInfo(&rCurrentProcessInfo)
            , mpMaterialProperties(&rMaterialProperties)
            , mpElementGeometry(&rElementGeometry)
        {}
 
        Parameters (const Parameters & rNewParameters)
            : mOptions(rNewParameters.mOptions)
            , mpGeneralizedStrainVector(rNewParameters.mpGeneralizedStrainVector)
            , mpGeneralizedStressVector(rNewParameters.mpGeneralizedStressVector)
            , mpConstitutiveMatrix(rNewParameters.mpConstitutiveMatrix)
            , mpShapeFunctionsValues(rNewParameters.mpShapeFunctionsValues)
            , mpShapeFunctionsDerivatives(rNewParameters.mpShapeFunctionsDerivatives)
            , mpCurrentProcessInfo(rNewParameters.mpCurrentProcessInfo)
            , mpMaterialProperties(rNewParameters.mpMaterialProperties)
            , mpElementGeometry(rNewParameters.mpElementGeometry)
        {}
 
    public:
 
        bool CheckShapeFunctions ()
        {
            if(!mpShapeFunctionsValues)
                KRATOS_THROW_ERROR(std::invalid_argument,"ShapeFunctionsValues NOT SET","");
 
            if(!mpShapeFunctionsDerivatives)
                KRATOS_THROW_ERROR(std::invalid_argument,"ShapeFunctionsDerivatives NOT SET","");
 
            return 1;
        }
 
        bool CheckInfoMaterialGeometry ()
        {
            if(!mpCurrentProcessInfo)
                KRATOS_THROW_ERROR(std::invalid_argument,"CurrentProcessInfo NOT SET","");
 
            if(!mpMaterialProperties)
                KRATOS_THROW_ERROR(std::invalid_argument,"MaterialProperties NOT SET","");
 
            if(!mpElementGeometry)
                KRATOS_THROW_ERROR(std::invalid_argument,"ElementGeometry NOT SET","");
 
            return 1;
        }
 
        bool CheckMechanicalVariables ()
        {
            if(!mpGeneralizedStrainVector)
                KRATOS_THROW_ERROR(std::invalid_argument,"GenralizedStrainVector NOT SET","");
 
            if(!mpGeneralizedStressVector)
                KRATOS_THROW_ERROR(std::invalid_argument,"GenralizedStressVector NOT SET","");
 
            if(!mpConstitutiveMatrix)
                KRATOS_THROW_ERROR(std::invalid_argument,"ConstitutiveMatrix NOT SET","");
 
            return 1;
        }
 
        void Set                             (Flags ThisFlag)                           {mOptions.Set(ThisFlag);};
        void Reset                           (Flags ThisFlag)                           {mOptions.Reset(ThisFlag);};
 
        void SetOptions                      (const Flags&  rOptions)                   {mOptions=rOptions;};
 
        void SetGeneralizedStrainVector      (Vector& rGeneralizedStrainVector)         {mpGeneralizedStrainVector=&rGeneralizedStrainVector;};
        void SetGeneralizedStressVector      (Vector& rGeneralizedStressVector)         {mpGeneralizedStressVector=&rGeneralizedStressVector;};
        void SetConstitutiveMatrix           (Matrix& rConstitutiveMatrix)              {mpConstitutiveMatrix =&rConstitutiveMatrix;};
 
        void SetShapeFunctionsValues         (const Vector& rShapeFunctionsValues)      {mpShapeFunctionsValues=&rShapeFunctionsValues;};
        void SetShapeFunctionsDerivatives    (const Matrix& rShapeFunctionsDerivatives) {mpShapeFunctionsDerivatives=&rShapeFunctionsDerivatives;};
        void SetProcessInfo                  (const ProcessInfo& rProcessInfo)          {mpCurrentProcessInfo =&rProcessInfo;};
        void SetMaterialProperties           (const Properties&  rMaterialProperties)   {mpMaterialProperties =&rMaterialProperties;};
        void SetElementGeometry              (const GeometryType& rElementGeometry)     {mpElementGeometry =&rElementGeometry;};
 
        Flags& GetOptions () {return mOptions;};
 
        Vector& GetGeneralizedStrainVector         () {return *mpGeneralizedStrainVector;};
        Vector& GetGeneralizedStressVector         () {return *mpGeneralizedStressVector;};
        Matrix& GetConstitutiveMatrix              () {return *mpConstitutiveMatrix;};
 
        const Vector& GetShapeFunctionsValues      () {return *mpShapeFunctionsValues;};
        const Matrix& GetShapeFunctionsDerivatives () {return *mpShapeFunctionsDerivatives;};
        const ProcessInfo&  GetProcessInfo         () {return *mpCurrentProcessInfo;};
        const Properties&   GetMaterialProperties  () {return *mpMaterialProperties;};
        const GeometryType& GetElementGeometry     () {return *mpElementGeometry;};
    };
 
    class IntegrationPoint
    {
 
    private:
 
        double mWeight;
        double mLocation;
        ConstitutiveLaw::Pointer mConstitutiveLaw;
 
    public:
 
        IntegrationPoint()
            : mWeight(0.0)
            , mLocation(0.0)
            , mConstitutiveLaw(ConstitutiveLaw::Pointer())
        {}
 
        IntegrationPoint(double location, double weight, const ConstitutiveLaw::Pointer pMaterial)
            : mWeight(weight)
            , mLocation(location)
            , mConstitutiveLaw(pMaterial)
        {}
 
        IntegrationPoint(const IntegrationPoint& other)
            : mWeight(other.mWeight)
            , mLocation(other.mLocation)
            , mConstitutiveLaw(other.mConstitutiveLaw != NULL ? other.mConstitutiveLaw->Clone() : ConstitutiveLaw::Pointer())
        {}
 
        IntegrationPoint & operator = (const IntegrationPoint & other) {
            if(this != &other) {
                mWeight = other.mWeight;
                mLocation = other.mLocation;
                mConstitutiveLaw = other.mConstitutiveLaw != NULL ? other.mConstitutiveLaw->Clone() : ConstitutiveLaw::Pointer();
            }
            return *this;
        }
 
                virtual ~IntegrationPoint(){};
 
    public:
 
        inline double GetWeight()const { return mWeight; }
        inline void SetWeight(double w) { mWeight = w; }
 
        inline double GetLocation()const { return mLocation; }
        inline void SetLocation(double l) { mLocation = l; }
 
        inline const ConstitutiveLaw::Pointer& GetConstitutiveLaw()const { return mConstitutiveLaw; }
        inline void SetConstitutiveLaw(const ConstitutiveLaw::Pointer& pLaw) { mConstitutiveLaw = pLaw; }
 
    private:
 
        friend class Serializer;
 
        virtual void save(Serializer& rSerializer) const
        {
            rSerializer.save("W", mWeight);
            rSerializer.save("L", mLocation);
            rSerializer.save("CLaw", mConstitutiveLaw);
        }
 
        virtual void load(Serializer& rSerializer)
        {
            rSerializer.load("W", mWeight);
            rSerializer.load("L", mLocation);
            rSerializer.load("CLaw", mConstitutiveLaw);
        }
    };
 
    class Ply
    {
 
    public:
 
        typedef std::vector< IntegrationPoint > IntegrationPointCollection;
 
    private:
 
        double mThickness;
        double mLocation;
        double mOrientationAngle;
        IntegrationPointCollection mIntegrationPoints;
        Properties::Pointer mpProperties;
 
    public:
 
        Ply()
            : mThickness(0.0)
            , mLocation(0.0)
            , mOrientationAngle(0.0)
            , mIntegrationPoints()
            , mpProperties(Properties::Pointer())
        {}
 
        Ply(double thickness, double location, double orientationAngle, int numPoints, const Properties::Pointer & pProperties)
            : mThickness(thickness)
            , mLocation(location)
            , mIntegrationPoints()
            , mpProperties(pProperties)
        {
            this->SetOrientationAngle(orientationAngle);
            this->SetUpIntegrationPoints(numPoints);
        }
 
        Ply(const Ply& other)
            : mThickness(other.mThickness)
            , mLocation(other.mLocation)
            , mOrientationAngle(other.mOrientationAngle)
            , mIntegrationPoints(other.mIntegrationPoints)
            , mpProperties(other.mpProperties)
        {}
 
        Ply & operator = (const Ply & other) {
            if(this != &other) {
                mThickness = other.mThickness;
                mLocation = other.mLocation;
                mOrientationAngle = other.mOrientationAngle;
                mIntegrationPoints = other.mIntegrationPoints;
                mpProperties = other.mpProperties;
            }
            return *this;
        }
 
                virtual ~Ply(){};
 
    public:
 
        inline double GetThickness()const { return mThickness; }
        inline void SetThickness(double thickness) { mThickness = thickness; }
 
        inline double GetLocation()const { return mLocation; }
        inline void SetLocation(double location) {
            if(location != mLocation) {
                for(IntegrationPointCollection::iterator it = mIntegrationPoints.begin(); it != mIntegrationPoints.end(); ++it)
                    (*it).SetLocation((*it).GetLocation() + location - mLocation); // remove the last location and add the new one (this avoids to re-setup the integration points.
                mLocation = location; // update the current location
            }
        }
 
        inline double GetOrientationAngle()const { return mOrientationAngle; }
        inline void SetOrientationAngle(double degrees) {
            mOrientationAngle = std::fmod(degrees, 360.0);
            if(mOrientationAngle < 0.0)
                mOrientationAngle += 360.0;
        }
 
        inline const IntegrationPointCollection& GetIntegrationPoints()const { return mIntegrationPoints; }
        inline IntegrationPointCollection& GetIntegrationPoints() { return mIntegrationPoints; }
 
        inline const Properties::Pointer & GetPropertiesPointer()const { return mpProperties; }
 
        inline const Properties & GetProperties()const { return *mpProperties; }
 
        inline double CalculateMassPerUnitArea()const { return mpProperties->GetValue(DENSITY) * mThickness; }
 
        inline IntegrationPointCollection::size_type NumberOfIntegrationPoints()const { return mIntegrationPoints.size(); }
 
        inline void SetConstitutiveLawAt(IntegrationPointCollection::size_type integrationPointID, const ConstitutiveLaw::Pointer& pNewConstitutiveLaw)
        {
            if(integrationPointID < mIntegrationPoints.size())
                mIntegrationPoints[integrationPointID].SetConstitutiveLaw(pNewConstitutiveLaw);
        }
 
    private:
 
        void SetUpIntegrationPoints(int n)
        {
            KRATOS_TRY
 
                            const ConstitutiveLaw::Pointer & pMaterial = GetProperties()[CONSTITUTIVE_LAW];
            if(pMaterial == NULL)
                          KRATOS_THROW_ERROR(std::logic_error, "A Ply needs a constitutive law to be set. Missing constitutive law in property : ", GetProperties().Id());
 
            // make sure the number is greater than 0 and odd
            if(n < 0) n = -n;
            if(n == 0) n = 5;
            if(n % 2 == 0) n += 1;
 
            // generate the weights (composite simpson rule)
            Vector ip_w(n, 1.0);
            if(n >= 3) {
                          for(int i = 1; i < n-1; i++) {
                            double iw = (i % 2 == 0) ? 2.0 : 4.0;
                            ip_w(i) = iw;
                          }
                          ip_w /= sum( ip_w );
            }
 
            // generate locations (direction: top(+thickness/2) to bottom(-thickness/2)
            Vector ip_loc(n, 0.0);
            if(n >= 3) {
                          double loc_start = mLocation + 0.5 * mThickness;
                          double loc_incr = mThickness / double(n-1);
                          for(int i = 0; i < n; i++) {
                            ip_loc(i) = loc_start;
                            loc_start -= loc_incr;
                          }
            }
 
            // generate the integration points
            mIntegrationPoints.clear();
            mIntegrationPoints.resize(n);
            for(int i = 0; i < n; i++) {
                          IntegrationPoint& intp = mIntegrationPoints[i];
                          intp.SetWeight(ip_w(i) * mThickness);
                          intp.SetLocation(ip_loc(i));
                          intp.SetConstitutiveLaw(pMaterial->Clone());
            }
 
            KRATOS_CATCH("")
                            }
 
         private:
 
          friend class Serializer;
 
          virtual void save(Serializer& rSerializer) const
          {
            rSerializer.save("T", mThickness);
            rSerializer.save("L", mLocation);
            rSerializer.save("O", mOrientationAngle);
            rSerializer.save("IntP", mIntegrationPoints);
            rSerializer.save("Prop", mpProperties);
          }
 
          virtual void load(Serializer& rSerializer)
          {
            rSerializer.load("T", mThickness);
            rSerializer.load("L", mLocation);
            rSerializer.load("O", mOrientationAngle);
            rSerializer.load("IntP", mIntegrationPoints);
            rSerializer.load("Prop", mpProperties);
          }
 
    };
 
 protected:
 
  struct ElementVariables
  {
    double DeterminantF;
    double DeterminantF0;
 
    Vector StrainVector_2D;
    Vector StressVector_2D;
    Matrix ConstitutiveMatrix_2D;
    Matrix DeformationGradientF_2D;
    Matrix DeformationGradientF0_2D;
 
    Vector StrainVector_3D;
    Vector StressVector_3D;
    Matrix ConstitutiveMatrix_3D;
    Matrix DeformationGradientF_3D;
    Matrix DeformationGradientF0_3D;
 
    double GYZ;
    double GXZ;
 
    Matrix H;
    Matrix L;
    Matrix LT;
    Vector CondensedStressVector;
  };
 
 
 public:
 
 
  ShellCrossSection();
 
  ShellCrossSection(const ShellCrossSection & other);
 
  ~ShellCrossSection() override;
 
 
 
  ShellCrossSection & operator = (const ShellCrossSection & other);
 
 
 
  void BeginStack();
 
  void AddPly(double thickness, double orientationAngle, int numPoints, const Properties::Pointer & pProperties);
 
  void EndStack();
 
  virtual std::string GetInfo()const;
 
  virtual ShellCrossSection::Pointer Clone()const;
 
  virtual bool Has(const Variable<double>& rThisVariable);
 
  virtual bool Has(const Variable<Vector>& rThisVariable);
 
  virtual bool Has(const Variable<Matrix>& rThisVariable);
 
  virtual bool Has(const Variable<array_1d<double, 3 > >& rThisVariable);
 
  virtual bool Has(const Variable<array_1d<double, 6 > >& rThisVariable);
 
  virtual double& GetValue(const Variable<double>& rThisVariable, double& rValue);
 
  virtual Vector& GetValue(const Variable<Vector>& rThisVariable, Vector& rValue);
 
  virtual Matrix& GetValue(const Variable<Matrix>& rThisVariable, Matrix& rValue);
 
  virtual array_1d<double, 3 > & GetValue(const Variable<array_1d<double, 3 > >& rVariable,
                                          array_1d<double, 3 > & rValue);
 
  virtual array_1d<double, 6 > & GetValue(const Variable<array_1d<double, 6 > >& rVariable,
                                          array_1d<double, 6 > & rValue);
 
  virtual void SetValue(const Variable<double>& rVariable,
                        const double& rValue,
                        const ProcessInfo& rCurrentProcessInfo);
 
  virtual void SetValue(const Variable<Vector >& rVariable,
                        const Vector& rValue,
                        const ProcessInfo& rCurrentProcessInfo);
 
  virtual void SetValue(const Variable<Matrix >& rVariable,
                        const Matrix& rValue,
                        const ProcessInfo& rCurrentProcessInfo);
 
  virtual void SetValue(const Variable<array_1d<double, 3 > >& rVariable,
                        const array_1d<double, 3 > & rValue,
                        const ProcessInfo& rCurrentProcessInfo);
 
  virtual void SetValue(const Variable<array_1d<double, 6 > >& rVariable,
                        const array_1d<double, 6 > & rValue,
                        const ProcessInfo& rCurrentProcessInfo);
 
  virtual bool ValidateInput(const Properties& rMaterialProperties);
 
  virtual void InitializeCrossSection(const Properties& rMaterialProperties,
                                      const GeometryType& rElementGeometry,
                                      const Vector& rShapeFunctionsValues);
 
  virtual void InitializeSolutionStep(const Properties& rMaterialProperties,
                                      const GeometryType& rElementGeometry,
                                      const Vector& rShapeFunctionsValues,
                                      const ProcessInfo& rCurrentProcessInfo);
 
  virtual void FinalizeSolutionStep(const Properties& rMaterialProperties,
                                    const GeometryType& rElementGeometry,
                                    const Vector& rShapeFunctionsValues,
                                    const ProcessInfo& rCurrentProcessInfo);
 
  virtual void InitializeNonLinearIteration(const Properties& rMaterialProperties,
                                            const GeometryType& rElementGeometry,
                                            const Vector& rShapeFunctionsValues,
                                            const ProcessInfo& rCurrentProcessInfo);
 
  virtual void FinalizeNonLinearIteration(const Properties& rMaterialProperties,
                                          const GeometryType& rElementGeometry,
                                          const Vector& rShapeFunctionsValues,
                                          const ProcessInfo& rCurrentProcessInfo);
 
  virtual void CalculateSectionResponse(Parameters& rValues, const ConstitutiveLaw::StressMeasure& rStressMeasure);
 
  virtual void FinalizeSectionResponse(Parameters& rValues, const ConstitutiveLaw::StressMeasure& rStressMeasure);
 
  virtual void ResetCrossSection(const Properties& rMaterialProperties,
                                 const GeometryType& rElementGeometry,
                                 const Vector& rShapeFunctionsValues);
 
  virtual int Check(const Properties& rMaterialProperties,
                    const GeometryType& rElementGeometry,
                    const ProcessInfo& rCurrentProcessInfo);
 
  inline void GetRotationMatrixForGeneralizedStrains(double radians, Matrix & T)
  {
    double c = std::cos(radians);
    double s = std::sin(radians);
 
    SizeType strain_size = GetStrainSize();
 
    if(T.size1() != strain_size || T.size2() != strain_size)
      T.resize(strain_size, strain_size, false);
    noalias( T ) = ZeroMatrix(strain_size, strain_size);
 
    T(0, 0) = c * c;            T(0, 1) =   s * s;              T(0, 2) = - s * c;
    T(1, 0) = s * s;            T(1, 1) =   c * c;              T(1, 2) =   s * c;
    T(2, 0) = 2.0 * s * c;      T(2, 1) = - 2.0 * s * c;        T(2, 2) = c * c - s * s;
 
    Matrix H = MathUtilsType::project(T, MathUtilsType::range(0, 3), MathUtilsType::range(0, 3));
    MathUtilsType::project(T, MathUtilsType::range(3, 6), MathUtilsType::range(3, 6), H);
 
    if(strain_size == 8)
    {
      T(6, 6) =   c;        T(6, 7) = s;
      T(7, 6) = - s;        T(7, 7) = c;
    }
  }
 
  inline void GetRotationMatrixForCondensedStrains(double radians, Matrix & T)
  {
    SizeType strain_size = GetCondensedStrainSize();
 
    if(T.size1() != strain_size || T.size2() != strain_size)
      T.resize(strain_size, strain_size, false);
    noalias( T ) = ZeroMatrix(strain_size, strain_size);
 
    T(0, 0) = 1.0; // condensed strain E.zz is always at index 0
 
    if(strain_size == 3) // if section is thin the condensed strains are (in order): E.zz E.yz E.xz
    {
      double c = std::cos(radians);
      double s = std::sin(radians);
 
      T(1, 1) =   c;        T(1, 2) = s;
      T(2, 1) = - s;        T(2, 2) = c;
    }
  }
 
  inline void GetRotationMatrixForGeneralizedStresses(double radians, Matrix & T)
  {
    double c = std::cos(radians);
    double s = std::sin(radians);
 
    SizeType strain_size = GetStrainSize();
 
    if(T.size1() != strain_size || T.size2() != strain_size)
      T.resize(strain_size, strain_size, false);
    noalias( T ) = ZeroMatrix(strain_size, strain_size);
 
    T(0, 0) = c * c;        T(0, 1) =   s * s;      T(0, 2) = - 2.0 * s * c;
    T(1, 0) = s * s;        T(1, 1) =   c * c;      T(1, 2) =   2.0 * s * c;
    T(2, 0) = s * c;        T(2, 1) = - s * c;      T(2, 2) = c * c - s * s;
 
    Matrix H = MathUtilsType::project(T, MathUtilsType::range(0, 3), MathUtilsType::range(0, 3));
    MathUtilsType::project(T, MathUtilsType::range(3, 6), MathUtilsType::range(3, 6), H);
 
    if(strain_size == 8)
    {
      T(6, 6) =   c;        T(6, 7) = s;
      T(7, 6) = - s;        T(7, 7) = c;
    }
  }
 
  inline void GetRotationMatrixForCondensedStresses(double radians, Matrix & T)
  {
    SizeType strain_size = GetCondensedStrainSize();
 
    if(T.size1() != strain_size || T.size2() != strain_size)
      T.resize(strain_size, strain_size, false);
    noalias( T ) = ZeroMatrix(strain_size, strain_size);
 
    T(0, 0) = 1.0; // condensed stresse S.zz is always at index 0
 
    if(strain_size == 3) // if section is thin the condensed stresses are (in order): S.zz S.yz S.xz
    {
      double c = std::cos(radians);
      double s = std::sin(radians);
 
      T(1, 1) =   c;        T(1, 2) = s;
      T(2, 1) = - s;        T(2, 2) = c;
    }
  }
 
 
 public:
 
 
  inline const double GetThickness()const
  {
    return mThickness;
  }
 
  inline const double GetOffset()const
  {
    return mOffset;
  }
 
  inline void SetOffset(double offset)
  {
    if((mOffset != offset) && (!mEditingStack))
    {
      for(PlyCollection::iterator it = mStack.begin(); it != mStack.end(); ++it)
        (*it).SetLocation((*it).GetLocation() + offset - mOffset);
      mOffset = offset;
    }
  }
 
  inline PlyCollection::size_type NumberOfPlies()const
  {
    return mStack.size();
  }
 
  inline SizeType NumberOfIntegrationPointsAt(SizeType ply_id)const
  {
    if(ply_id < mStack.size())
      return mStack[ply_id].NumberOfIntegrationPoints();
    return 0;
  }
 
  inline void SetConstitutiveLawAt(SizeType ply_id, SizeType point_id, const ConstitutiveLaw::Pointer& pNewConstitutiveLaw)
  {
    if(ply_id < mStack.size())
      mStack[ply_id].SetConstitutiveLawAt(point_id, pNewConstitutiveLaw);
  }
 
  inline double CalculateMassPerUnitArea()const
  {
    double vol(0.0);
    for(PlyCollection::const_iterator it = mStack.begin(); it != mStack.end(); ++it)
      vol += (*it).CalculateMassPerUnitArea();
    return vol;
  }
 
  inline double CalculateAvarageDensity()const
  {
    return CalculateMassPerUnitArea() / mThickness;
  }
 
  inline double GetOrientationAngle()const
  {
    return mOrientation;
  }
 
  inline void SetOrientationAngle(double radians)
  {
    mOrientation = radians;
  }
 
  inline SectionBehaviorType GetSectionBehavior()const
  {
    return mBehavior;
  }
 
  inline void SetSectionBehavior(SectionBehaviorType behavior)
  {
    mBehavior = behavior;
  }
 
  inline SizeType GetStrainSize()
  {
    return (mBehavior == Thick) ? 8 : 6;
  }
 
  inline SizeType GetCondensedStrainSize()
  {
    return (mBehavior == Thick) ? 1 : 3;
  }
 
  inline double GetDrillingStiffness()const
  {
    return mDrillingPenalty;
  }
 
 
 private:
 
 
  void InitializeParameters(Parameters& rValues, ConstitutiveLaw::Parameters& rMaterialValues, ElementVariables& rVariables);
 
  void UpdateIntegrationPointParameters(IntegrationPoint& rPoint, ConstitutiveLaw::Parameters& rMaterialValues, ElementVariables& rVariables);
 
  void CalculateIntegrationPointResponse(IntegrationPoint& rPoint,
                                         ConstitutiveLaw::Parameters& rMaterialValues,
                                         Parameters& rValues,
                                         ElementVariables& rVariables,
                                         const ConstitutiveLaw::StressMeasure& rStressMeasure);
 
  void PrivateCopy(const ShellCrossSection & other);
 
 
 public:
 
 
 
 private:
 
 
  double mThickness;
  double mOffset;
  PlyCollection mStack;
  bool mEditingStack;
  bool mHasDrillingPenalty;
  double mDrillingPenalty;
  double mOrientation;
  SectionBehaviorType mBehavior;
  bool mInitialized;
  bool mNeedsOOPCondensation;
  Vector mOOP_CondensedStrains;
  Vector mOOP_CondensedStrains_converged;
 
 
 
  friend class Serializer;
 
  void save(Serializer& rSerializer) const override
  {
    KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Flags );
    rSerializer.save("th", mThickness);
    rSerializer.save("offs", mOffset);
    rSerializer.save("stack", mStack);
    rSerializer.save("edit", mEditingStack);
    rSerializer.save("dr", mHasDrillingPenalty);
    rSerializer.save("bdr", mDrillingPenalty);
    rSerializer.save("or", mOrientation);
 
    rSerializer.save("behav", (int)mBehavior);
 
    rSerializer.save("init", mInitialized);
    rSerializer.save("hasOOP", mNeedsOOPCondensation);
    rSerializer.save("OOP_eps", mOOP_CondensedStrains_converged);
  }
 
  void load(Serializer& rSerializer) override
  {
    KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Flags );
    rSerializer.load("th", mThickness);
    rSerializer.load("offs", mOffset);
    rSerializer.load("stack", mStack);
    rSerializer.load("edit", mEditingStack);
    rSerializer.load("dr", mHasDrillingPenalty);
    rSerializer.load("bdr", mDrillingPenalty);
    rSerializer.load("or", mOrientation);
 
    int temp;
    rSerializer.load("behav", temp);
    mBehavior = (SectionBehaviorType)temp;
 
    rSerializer.load("init", mInitialized);
    rSerializer.load("hasOOP", mNeedsOOPCondensation);
    rSerializer.load("OOP_eps", mOOP_CondensedStrains_converged);
  }
 
 
 public:
 
  DECLARE_ADD_THIS_TYPE_TO_PROPERTIES
  DECLARE_GET_THIS_TYPE_FROM_PROPERTIES
 
};
 
 
inline std::istream & operator >> (std::istream & rIStream, ShellCrossSection & rThis);
 
inline std::ostream & operator << (std::ostream & rOStream, const ShellCrossSection & rThis)
{
  return rOStream << rThis.GetInfo();
}
 
 
}
 
 
#endif // SHELL_CROSS_SECTION_H_INCLUDED