ASDShellQ4.h

// -*-c++-*-
//----------------------------------------------------------------------------
//  XC program; finite element analysis code
//  for structural analysis and design.
//
//  Copyright (C)  Luis C. Pérez Tato
//
//  This program derives from OpenSees <http://opensees.berkeley.edu>
//  developed by the  «Pacific earthquake engineering research center».
//
//  Except for the restrictions that may arise from the copyright
//  of the original program (see copyright_opensees.txt)
//  XC is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or 
//  (at your option) any later version.
//
//  This software is distributed in the hope that it will be useful, but 
//  WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details. 
//
//
// You should have received a copy of the GNU General Public License 
// along with this program.
// If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------------
/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */
                                                                        
// $Revision: 1.10 $
// $Date: 2020/05/18 22:51:21 $

// Original implementation: Massimo Petracca (ASDEA)
//
// A 4-node general shell element based on 
// the AGQ6-I formulation for the membrane part, 
// and the MITC4 formulation for the shear-deformable bending part.
// 
// It supports both linear and corotational kinematics. Warped geometries
// can be modelled since this element is not assumed flat.
//
//  References:
// 
// 1 Chen, Xiao-Ming, et al. "Membrane elements insensitive to distortion 
//   using the quadrilateral area coordinate method." 
//   Computers & Structures 82.1 (2004): 35-54.
//   http://www.paper.edu.cn/scholar/showpdf/MUT2ANwINTT0Ax5h
//
// 2 Dvorkin, Eduardo N., and Klaus-Jurgen Bathe. 
//   "A continuum mechanics based four-node shell element for 
//   general non-linear analysis." Engineering computations (1984).
//   https://www.researchgate.net/profile/Eduardo_Dvorkin/publication/235313212_A_Continuum_mechanics_based_four-node_shell_element_for_general_nonlinear_analysis/links/00b7d52611d8813ffe000000.pdf
//
// 3 Bathe, Klaus-Jurgen, and Eduardo N. Dvorkin. 
//   "A four-node plate bending element based on Mindlin/Reissner plate theory 
//   and a mixed interpolation." 
//   International Journal for Numerical Methods in Engineering 21.2 (1985): 367-383.
//   http://www.simytec.com/docs/Short_communicaion_%20four_node_plate.pdf
//
// 4 Hughes, Thomas JR, and F. Brezzi. 
//   "On drilling degrees of freedom." 
//   Computer methods in applied mechanics and engineering 72.1 (1989): 105-121.
//   https://www.sciencedirect.com/science/article/pii/0045782589901242
//
// Notes:
//
// - The AGQ6-I formulation greatly enhances the membrane behavior of the
//   standard 4-node iso-parametric element. However it does not pass the constant-strain
//   patch test as most of the enhanced elements based on incompatible-modes.
//   In this implementation the enhancing BQ matrix is corrected to make the element
//   pass the constant strain patch test. At each gauss point the BQ = BQ - 1/V*BQ_mean
//
// - The drilling DOF is treated using the Hughes-Brezzi formulation. The drilling
//   stiffness is taken equal to the initial (elastic) in-plane section shear modulus.
//   Taking the real shear modulus (and not a scaled version of it, as suggested in many
//   articles) is mandatory to obtain a correct response for warped shells, i.e. when
//   the drilling rotation affects the bending rotation. However using the full shear modulus
//   deteriorates the element membrane behavior, making the element stiffer. For example,
//   when using softening material models, this drilling stiffness makes the element lock
//   failing in representing crack propagation correctly.
//   Here we solved this problem using a reduced integration for the drilling part. However
//   a pure reduced integration for the drilling DOFs, leads to 3 spurious zero-energy modes.
//   Therefore we also include a fully integrated contribution of the drilling DOFs
//   scaling the drilling stiffness by 1.0e-4, so just to suppress these spurious zero energy modes.
// 

#ifndef ASDShellQ4_h
#define ASDShellQ4_h

#include "utility/matrix/Vector.h"
#include "utility/matrix/Matrix.h"
#include "utility/matrix/ID.h"
#include "domain/mesh/element/plane/QuadBase4N.h"
#include "domain/mesh/element/utils/physical_properties/SectionFDPhysicalProperties.h"
#include "domain/mesh/element/utils/damping/DampingVector.h"

namespace XC {

class SectionForceDeformation;
class ASDShellQ4Transformation;
class ASDShellQ4LocalCoordinateSystem;

class ASDShellQ4: public QuadBase4N<SectionFDPhysicalProperties>
  {
  public:
    enum DrillingDOFMode
      {
        DrillingDOF_Elastic = 0,
        DrillingDOF_NonLinear = 1,
      };
    class NLDrillingData
      {
      public:
        std::vector<Vector> strain_comm = { Vector(8), Vector(8), Vector(8), Vector(8) };
        std::vector<Vector> stress_comm = { Vector(8), Vector(8), Vector(8), Vector(8) };
        std::vector<double> damage = { 0.0, 0.0, 0.0, 0.0 };
        std::vector<double> damage_comm = { 0.0, 0.0, 0.0, 0.0 };
      };
    class EASData
      {
      public:
        Vector Q = Vector(4);
        Vector Q_converged = Vector(4);
        Vector U = Vector(24);
        Vector U_converged = Vector(24);
        Vector Q_residual = Vector(4);
        Matrix KQQ_inv = Matrix(4, 4);
        Matrix KQU = Matrix(4, 24); // L = G'*C*B
        Matrix KUQ = Matrix(24, 4); // L^T = B'*C'*G
      };
    
  private:
    // coordinate transformation
    ASDShellQ4Transformation *m_transformation= nullptr;

    // vectors for applying load.
    Vector m_load;

    // drilling strain for the indipendent rotation field (Hughes-Brezzi)
    mutable std::vector<double> m_drill_strain= { 0.0, 0.0, 0.0, 0.0 };
    double m_drill_stiffness = 0.0;

    // section orientation with respect to the local coordinate system
    Vector m_local_x;
    double m_angle = 0.0;

    // initialization flag
    bool m_initialized= false;
    // damping
    DampingVector m_damping;
    
    // members for non-linear treatment of AGQI internal DOFs:
    // it has 24 displacement DOFs
    // and 4 internal DOFs for membrane enhancement
    EASData *m_eas= nullptr;
    
    // drilling strain for the independent rotation field (Hughes-Brezzi)
    DrillingDOFMode m_drill_mode = DrillingDOF_Elastic;
    double m_drill_stab = 0.01;
    NLDrillingData* m_nldrill= nullptr;
    
  private:

    void free_mem(void);
    void free_crd_transf(void);
    void free_eas(void);
    void free_nldrill(void);
    void alloc(bool corotational);
    void alloc(const ASDShellQ4Transformation *);
    void alloc_eas(void);
    void copy_eas(const EASData &);
    void alloc_nldrill(void);
    void copy_nldrill(const NLDrillingData &);
    // internal method to compute everything using switches...
    int calculateAll(Matrix& LHS, Vector& RHS, int options) const;

    void AGQIinitialize(void);
    void AGQIupdate(const Vector& UL) const;
    void AGQIbeginGaussLoop(const ASDShellQ4LocalCoordinateSystem& reference_cs) const;
  protected:
    int sendData(Communicator &);
    int recvData(const Communicator &);
    
  public:

    // life cycle
    ASDShellQ4(int tag= 0, bool corotational= true);
    ASDShellQ4(int tag,const SectionForceDeformation *ptr_mat);
    ASDShellQ4(int tag,
        int node1,
        int node2,
        int node3,
        int node4,
        SectionForceDeformation* section,
    const Vector &local_x,
    bool corotational= true,
        bool use_eas = true,
        DrillingDOFMode drill_mode = DrillingDOF_Elastic,
        double drilling_stab = 0.01,
    Damping *theDamping= nullptr);
    ASDShellQ4(const ASDShellQ4 &);
    ASDShellQ4 &operator=(const ASDShellQ4 &);
    virtual ~ASDShellQ4();
    Element *getCopy(void) const;

    // domain
    int resetNodalCoordinates(void);
    void setDomain(Domain* theDomain);
    int setDamping(Domain *theDomain, Damping *theDamping);

    // print
    void Print(std::ostream &, int flag);

    int getNumDOF(void) const;

    virtual const SectionForceDeformation *getSectionPtr(const size_t &) const;
    
    // methods dealing with committed state and update
    int commitState(void);
    int revertToLastCommit(void);
    int revertToStart(void);
    int update(void);

    // methods to return the current linearized stiffness,
    // damping and mass matrices
    const Matrix &getTangentStiff(void) const;
    const Matrix &getInitialStiff(void) const;
    const Matrix &getMass(void) const;

    double getCharacteristicLength(void) const;
    
    // methods for applying loads
    void alive(void);
    void zeroLoad(void);
    int addLoad(ElementalLoad* theLoad, double loadFactor);
    int addInertiaLoadToUnbalance(const Vector& accel);

    // methods for obtaining resisting force (force includes elemental loads)
    const Vector &getResistingForce(void) const;
    const Vector &getResistingForceIncInertia(void) const;

    // public methods for element output
    int sendSelf(Communicator &);
    int recvSelf(const Communicator &);

    Response *setResponse(const std::vector<std::string> &argv, Information &eleInformation);
    int getResponse(int responseID, Information &eleInformation);

    int setParameter(const std::vector<std::string> &argv, Parameter &param);
  };
} // end of XC namespace
#endif // ASDShellQ4_h