!/*****************************************************************************/
! *
! *  Elmer, A Finite Element Software for Multiphysical Problems
! *
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
! * 
! *  This library is free software; you can redistribute it and/or
! *  modify it under the terms of the GNU Lesser General Public
! *  License as published by the Free Software Foundation; either
! *  version 2.1 of the License, or (at your option) any later version.
! *
! *  This library 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
! *  Lesser General Public License for more details.
! * 
! *  You should have received a copy of the GNU Lesser General Public
! *  License along with this library (in file ../LGPL-2.1); if not, write 
! *  to the Free Software Foundation, Inc., 51 Franklin Street, 
! *  Fifth Floor, Boston, MA  02110-1301  USA
! *
! *****************************************************************************/
!
!/******************************************************************************
!
!/******************************************************************************
! *
! *  A module for solving the one-dimensional Timoshenko beam equations
! *
! *                  m*d^2u/dt^2 - N'= F  on (0,L)
! *  IM*d^2(theta)/dt^2 - M' - t x N = G  on (0,L).
! *
! *  with the stress resultants N and M representing forces and moments 
! *  experienced by the cross section of a beam and t giving the tangent
! *  vector to the beam axis. By using the constitutive equation, the 
! *  displacement u:[0,L] -> R^3 and the so-called rotation theta:[0,L] -> R^3 
! *  are approximated by a finite element discretization. The fields
! *  m and IM are the beam mass per unit length and the mass moment of
! *  inertia of the cross section. In addition to Dirichlet constraints,
! *  the force and the moment can be specified as nodal loads.
! *
! *  One may assume the quadratic approximation of the deflection and 
! *  the linear approximation of the rotation and eliminate the quadratic
! *  bubble part of the deflection by static condensation. This gives
! *  an approximation over the lowest-order nodal mesh. This approximation
! *  may be generated by giving the element definition "Element = p:1 b:1". 
! *  Higher-order versions are not yet available.
! *
! *  Authors: Mika Malinen
! *  Email:   mika.malinen@csc.fi
! *  Web:     http://www.csc.fi/elmer
! *  Address: CSC - IT Center for Science Ltd.
! *           Keilaranta 14
! *           02101 Espoo, Finland 
! *
! *  Original Date: May 27, 2019
! *
! *****************************************************************************/

!------------------------------------------------------------------------------
SUBROUTINE TimoshenkoSolver_Init0(Model, Solver, dt, Transient)
!------------------------------------------------------------------------------
  USE DefUtils
  IMPLICIT NONE
!------------------------------------------------------------------------------
  TYPE(Model_t) :: Model
  TYPE(Solver_t) :: Solver
  REAL(KIND=dp) :: dt
  LOGICAL :: Transient
!------------------------------------------------------------------------------
  TYPE(ValueList_t), POINTER :: SolverPars, Simulation
  LOGICAL :: Found, MeshDisplacementActive
!------------------------------------------------------------------------------
  SolverPars => GetSolverParams()

  CALL ListAddInteger(SolverPars, 'Variable DOFs', 6)
  CALL ListAddNewString(SolverPars, 'Variable', 'Deflection[U:3 Theta:3]')
  CALL ListAddNewString(SolverPars, 'Element', 'p:1 b:1')

  CALL ListAddNewLogical(SolverPars, 'Bubbles in Global System', .FALSE.)
  CALL ListAddNewLogical(SolverPars, 'Use Global Mass Matrix',.TRUE.)
  IF (Transient) THEN
    CALL ListAddInteger(SolverPars, 'Time derivative order', 2)
    CALL ListAddNewString(SolverPars, 'Timestepping Method', 'Bossak')
  END IF

  MeshDisplacementActive = GetLogical(SolverPars, 'Displace Mesh', Found)
  IF (MeshDisplacementActive) THEN
    Simulation => GetSimulation()
    CALL ListAddLogical(Simulation, 'Initialize Dirichlet Conditions', .FALSE.) 
  END IF
  
  CALL ListAddNewLogical(SolverPars,'Beam Solver',.TRUE.)
!------------------------------------------------------------------------------
END SUBROUTINE TimoshenkoSolver_Init0
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
SUBROUTINE TimoshenkoSolver(Model, Solver, dt, TransientSimulation)
!------------------------------------------------------------------------------
  USE DefUtils
  USE SolidMechanicsUtils

  IMPLICIT NONE
!------------------------------------------------------------------------------
  TYPE(Model_t) :: Model
  TYPE(Solver_t) :: Solver
  REAL(KIND=dp) :: dt
  LOGICAL :: TransientSimulation
!------------------------------------------------------------------------------
! Local variables
!------------------------------------------------------------------------------
  TYPE(Element_t), POINTER :: Element
  TYPE(Mesh_t), POINTER :: Mesh
  LOGICAL :: Found
  INTEGER :: K, Active, n, nb, nd
  INTEGER :: iter, maxiter
  REAL(KIND=dp) :: Norm
  LOGICAL :: HarmonicAssembly, MassAssembly, MeshDisplacementActive
  TYPE(ValueList_t), POINTER :: Params
!------------------------------------------------------------------------------

  CALL DefaultStart()

  IF (.NOT. ListCheckPresentAnyMaterial(Model, 'Principal Direction 2') .AND. &
      .NOT. ListCheckPresentAnyMaterial(Model, 'Director')) THEN
    CALL Warn('TimoshenkoSolver', &
        'Principal axes unspecified, assuming a circular cross section')
  END IF
  
  Params => GetSolverParams()
  
  maxiter = ListGetInteger(Params, &
      'Nonlinear System Max Iterations', Found, minv=1)
  IF (.NOT. Found ) maxiter = 1

  HarmonicAssembly = EigenOrHarmonicAnalysis() &
      .OR. ListGetLogical( Params,'Harmonic Mode',Found ) 
  MassAssembly = TransientSimulation .OR. HarmonicAssembly

  MeshDisplacementActive = GetLogical(Params, 'Displace Mesh', Found)
  IF (MeshDisplacementActive) THEN
    Mesh => GetMesh()
    CALL Info('TimoshenkoSolver', 'Returning the mesh to its reference position', Level=4)     
    CALL DisplaceMesh(Mesh, Solver % Variable % Values, -1, Solver % Variable % Perm, &
        6, .FALSE., 3)      
  END IF

  !--------------------------
  ! Nonlinear iteration loop:
  !--------------------------
  DO iter=1,maxiter
    !-----------------
    ! System assembly:
    !-----------------
    CALL DefaultInitialize()
    Active = GetNOFActive()
    DO K=1,Active
      Element => GetActiveElement(K)

      IF ( .NOT. (GetElementFamily(Element) == 2) ) CYCLE

      n  = GetElementNOFNodes()
      nd = GetElementNOFDOFs()
      nb = GetElementNOFBDOFs()

      CALL BeamStiffnessMatrix(Element, n, nd+nb, nb, TransientSimulation, &
          MassAssembly=MassAssembly, HarmonicAssembly=HarmonicAssembly)      
    END DO

    CALL DefaultFinishBulkAssembly()

    CALL DefaultFinishBoundaryAssembly()
    CALL DefaultFinishAssembly()
    CALL DefaultDirichletBCs()

    !-----------------------
    ! Call a linear solver:
    !-----------------------
    Norm = DefaultSolve()
    IF ( DefaultConverged() ) EXIT    

  END DO

  CALL DefaultFinish()

  IF (MeshDisplacementActive) THEN
    CALL Info('TimoshenkoSolver', 'Displacing the mesh with computed displacement field', Level=4)
    CALL DisplaceMesh(Mesh, Solver % Variable % Values, 1, Solver % Variable % Perm, &
        6, .FALSE., 3)
  END IF
!------------------------------------------------------------------------------
END SUBROUTINE TimoshenkoSolver
!------------------------------------------------------------------------------
