DislocationLinkingNumber.h

Go to the documentation of this file.

/* This file is part of MoDELib, the Mechanics Of Defects Evolution Library.
 *
 *
 * MoDELib is distributed without any warranty under the
 * GNU General Public License (GPL) v2 <http://www.gnu.org/licenses/>.
 */
 
#ifndef model_DislocationLinkingNumber_H_
#define model_DislocationLinkingNumber_H_
 
#include <numbers>
#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <cfloat>
#include <math.h>       /* asin */
 
#include <SegmentSegmentDistance.h>
 
namespace model
{
    template<int dim> class DislocationNetwork;
    
    
    template<int dim>
    struct LinkingNumber
    {
        
        static constexpr double tol=FLT_EPSILON;
        
        
        typedef Eigen::Matrix<double,dim,1> VectorDim;
        
        //        /**********************************************************************/
        //        static double F(const double& t1,
        //                        const double& t2,
        //                        const double& cosBeta,
        //                        const double& sinBeta2,
        //                        const double& a0)
        //        {
        //
        //            const double num=t1*t2+a0*a0*cosBeta;
        //            const double den=a0*sqrt(t1*t1+t2*t2-2.0*t1*t2*cosBeta+a0*a0*sinBeta2);
        //            //            return -atan2(den+tol,num)/(4.0*std::numbers::pi);
        //            //            return fabs(den)>tol? -atan2(den,num)/(4.0*std::numbers::pi) : -atan2(0.0,num)/(4.0*std::numbers::pi);
        //            return fabs(den)>tol? -atan2(den,num)/(4.0*std::numbers::pi) : 0.0;
        //
        //        }
        
        /**********************************************************************/
 
        static VectorDim unitCross(const VectorDim& v1,
                                   const VectorDim& v2);
        
        /**********************************************************************/
 
        static double myAsin(const double& x);
        
        /**********************************************************************/
 
        static double segmentPairLN(const VectorDim& R1,
                                    const VectorDim& R2,
                                    const VectorDim& R3,
                                    const VectorDim& R4);
        
        /**********************************************************************/
 
        static double loopPairLN(const std::vector<VectorDim>& loop1,
                                 const std::vector<VectorDim>& loop2);
        
        /**********************************************************************/
 
        static double loopPairHelicity(const std::vector<VectorDim>& loop1,
                                       const VectorDim& b1,
                                       const std::vector<VectorDim>& loop2,
                                       const VectorDim& b2);
        
    };
    
    template<int dim>
    class DislocationLinkingNumber
    {
        typedef typename DislocationNetwork<dim>::LoopType LoopType;
        typedef typename DislocationNetwork<dim>::LoopLinkType LoopLink;
        typedef typename DislocationNetwork<dim>::VectorDim VectorDim;
        
        
        const DislocationNetwork<dim>& DN;
        
        static constexpr double tol=FLT_EPSILON;
        
        
        
        /**********************************************************************/
 
        static double segmentPairLN(const LoopLink& link1,
                                    const LoopLink& link2);
        
        
        
        
        
        //        /**********************************************************************/
        //        static double segmentPairLN(const VectorDim& R1,
        //                                      const VectorDim& R2,
        //                                      const VectorDim& R3,
        //                                      const VectorDim& R4)
        //        {
        //
        //            const VectorDim R13(R3-R1);
        //            const VectorDim R12(R2-R1);
        //            const VectorDim R34(R4-R3);
        //
        //            const double R13norm(R13.norm());
        //            const double R12norm(R12.norm());
        //            const double R34norm(R34.norm());
        //
        //            if(R13norm>tol && R12norm>tol && R34norm>tol)
        //            {
        //                const double prod(R34.cross(R12).dot(R13)/R13norm/R12norm/R34norm);
        //                if(fabs(prod)>tol)
        //                {
        //                    const VectorDim R14(R4-R1);
        //                    const VectorDim R24(R4-R2);
        //                    const VectorDim R23(R3-R2);
        //
        //                    const VectorDim n1(unitCross(R13,R14));
        //                    const VectorDim n2(unitCross(R14,R24));
        //                    const VectorDim n3(unitCross(R24,R23));
        //                    const VectorDim n4(unitCross(R23,R13));
        //
        //
        //                    const double Os= myAsin(n1.dot(n2))
        //                    /*            */+myAsin(n2.dot(n3))
        //                    /*            */+myAsin(n3.dot(n4))
        //                    /*            */+myAsin(n4.dot(n1));
        //
        //                    return Os/(4.0*std::numbers::pi)*prod/fabs(prod);
        //
        //                }
        //                else
        //                {
        //                    return 0.0;
        //                }
        //            }
        //            else
        //            {
        //                return 0.0;
        //            }
        //        }
        
        //        /**********************************************************************/
        //        static double segmentPairLN_b(const LoopLink& link1,
        //                                      const LoopLink& link2)
        //        {
        //            const VectorDim R1(link1.source->get_P());
        //            const VectorDim R2(link2.source->get_P());
        //
        //            const VectorDim S1(link1.sink->get_P()-link1.source->get_P());
        //            const VectorDim S2(link2.sink->get_P()-link2.source->get_P());
        //            const double s1(S1.norm());
        //            const double s2(S2.norm());
        //            const VectorDim e1(S1.normalized());
        //            const VectorDim e2(S2.normalized());
        //            const VectorDim R12(R2-R1);
        //
        //            const double cosBeta(e1.dot(e2));
        //            const double sinBeta2(1.0-cosBeta*cosBeta);
        //            const double a1=R12.dot(e2*cosBeta-e1)/(sinBeta2+tol);
        //            const double a2=R12.dot(e2-e1*cosBeta)/(sinBeta2+tol);
        //            const double a0=R12.dot(e1.cross(e2))/(sinBeta2+tol);
        //
        //
        //            return F(a1+s1,a2+s2,cosBeta,sinBeta2,a0)
        //            /* */ -F(a1+s1,a2   ,cosBeta,sinBeta2,a0)
        //            /* */ -F(a1   ,a2+s2,cosBeta,sinBeta2,a0)
        //            /* */ +F(a1   ,a2   ,cosBeta,sinBeta2,a0);
        //        }
        
        
        /**********************************************************************/
 
        static std::pair<double,bool> linkingNumber(const LoopType& loop1,
                                                    const LoopType& loop2);
        
    public:
        
        /**********************************************************************/
 
        DislocationLinkingNumber(const DislocationNetwork<dim>& DN_in);
        
        
        /**********************************************************************/
        template <class T>
        friend T& operator << (T& os, const DislocationLinkingNumber<dim>& ln)
        {
            ln.write(os);
            return os;
        }
 
        std::ostream& write(std::ostream& os) const;
        
    };
    
}
#endif