Fem2DElementTriangle.h

// This file is a part of the OpenSurgSim project.
// Copyright 2013, SimQuest Solutions Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef SURGSIM_PHYSICS_FEM2DELEMENTTRIANGLE_H
#define SURGSIM_PHYSICS_FEM2DELEMENTTRIANGLE_H

#include <array>

#include "SurgSim/Physics/Fem.h"
#include "SurgSim/Physics/FemElement.h"

namespace SurgSim
{

namespace Physics
{
SURGSIM_STATIC_REGISTRATION(Fem2DElementTriangle);

class Fem2DElementTriangle : public FemElement
{
    typedef Eigen::Matrix<double, 3, 3> Matrix33Type;

    typedef Eigen::Matrix<double, 3, 6> Matrix36Type;
    typedef Eigen::Matrix<double, 6, 6> Matrix66Type;

    typedef Eigen::Matrix<double, 3, 9> Matrix39Type;
    typedef Eigen::Matrix<double, 9, 9> Matrix99Type;

public:
    Fem2DElementTriangle();

    explicit Fem2DElementTriangle(std::array<size_t, 3> nodeIds);

    explicit Fem2DElementTriangle(std::shared_ptr<FemElementStructs::FemElementParameter> elementData);

    SURGSIM_CLASSNAME(SurgSim::Physics::Fem2DElementTriangle)

    
    void setThickness(double thickness);

    double getThickness() const;

    void initialize(const SurgSim::Math::OdeState& state) override;

    double getVolume(const SurgSim::Math::OdeState& state) const override;

    SurgSim::Math::Vector computeCartesianCoordinate(const SurgSim::Math::OdeState& state,
            const SurgSim::Math::Vector& naturalCoordinate) const override;

    SurgSim::Math::Vector computeNaturalCoordinate(const SurgSim::Math::OdeState& state,
            const SurgSim::Math::Vector& cartesianCoordinate) const override;

protected:
    void initializeMembers();

    SurgSim::Math::Matrix33d computeRotation(const SurgSim::Math::OdeState& state);

    virtual void computeLocalStiffness(const SurgSim::Math::OdeState& state,
                                       Eigen::Matrix<double, 18, 18>* localStiffnessMatrix);

    void computeStiffness(const SurgSim::Math::OdeState& state, SurgSim::Math::Matrix* stiffnessMatrix);

    virtual void computeLocalMass(const SurgSim::Math::OdeState& state,
                                  Eigen::Matrix<double, 18, 18>* localMassMatrix);

    void computeMass(const SurgSim::Math::OdeState& state, SurgSim::Math::Matrix* massMatrix);

    void doUpdateFMDK(const Math::OdeState& state, int options) override;

    Eigen::Matrix<double, 18, 1> m_x0;

    SurgSim::Math::Matrix33d m_initialRotation;

    Eigen::Matrix<double, 18, 18> m_MLocal;
    Eigen::Matrix<double, 18, 18> m_KLocal;

    double m_restArea;

    double m_thickness;

    void computeShapeFunctionsParameters(const SurgSim::Math::OdeState& restState);

    SurgSim::Math::Matrix33d m_membraneShapeFunctionsParameters;

    // Thin-plate (bending/twisting) specific data structure
    // DOF simulated: (z, thetaX, thetaY)
    // "A Study Of Three-Node Triangular Plate Bending Elements", Jean-Louis Batoz
    // Numerical Methods in Engineering, vol 15, 1771-1812 (1980)
    // Indices are as follow:
    // 0 1 2 denotes triangle's points ABC:
    // 4 (mid-edge 12) 5 (mid-edge 20) 6 (mid-edge 01) denotes mid-edge points
    // Data structures having only mid-edge information are 0 based (0->4 (mid-egde 12) ; 1->5 ; 2->6)

    SurgSim::Math::Vector3d m_xij;     
    SurgSim::Math::Vector3d m_yij;     
    SurgSim::Math::Vector3d m_lij_sqr; 
    SurgSim::Math::Vector3d m_ak;      
    SurgSim::Math::Vector3d m_bk;      
    SurgSim::Math::Vector3d m_ck;      
    SurgSim::Math::Vector3d m_dk;      
    SurgSim::Math::Vector3d m_ek;      


    SurgSim::Math::Vector3d m_Pk;      
    SurgSim::Math::Vector3d m_qk;      
    SurgSim::Math::Vector3d m_tk;      
    SurgSim::Math::Vector3d m_rk;      


    SurgSim::Math::Matrix m_integral_dT_d;  
    SurgSim::Math::Matrix m_integralHyiHyj; 
    SurgSim::Math::Matrix m_integralHxiHxj; 


    std::array<double, 9> batozDhxDxi(double xi, double eta) const;
    std::array<double, 9> batozDhxDeta(double xi, double eta) const;
    std::array<double, 9> batozDhyDxi(double xi, double eta) const;
    std::array<double, 9> batozDhyDeta(double xi, double eta) const;
    Matrix39Type batozStrainDisplacement(double xi, double eta) const;

private:
    void computeLocalMembraneMass(const SurgSim::Math::OdeState& state,
                                  Eigen::Matrix<double, 18, 18>* localMassMatrix);

    void computeLocalPlateMass(const SurgSim::Math::OdeState& state,
                               Eigen::Matrix<double, 18, 18>* localMassMatrix);

    void computeIntegral_dTd();

    void computeIntegral_HxHxT();

    void computeIntegral_HyHyT();
};

} // namespace Physics

} // namespace SurgSim

#endif // SURGSIM_PHYSICS_FEM2DELEMENTTRIANGLE_H