DislocationLoop.cpp

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_DislocationLoop_cpp_
#define model_DislocationLoop_cpp_
 
 
#include <DislocationDynamics/DislocationDynamicsBase.h>
#include <algorithm>
#include <DislocationLoop.h>
#include <DislocationDynamics/DislocationFieldBase.h>
#include <SutherlandHodgman.h>
#include <DislocationDynamics/DislocationLoopPatches.h>
 
namespace model
{
    template <int dim>
    DislocationLoop<dim>::DislocationLoop(LoopNetworkType* const net,
                                                 const VectorDim& B,
                                                 const std::shared_ptr<GlidePlaneType>& glidePlane_in) :
    /* init */ Loop<DislocationLoop>(net,glidePlane_in->grain.rationalLatticeDirection(B))
    /* init */,glidePlane(glidePlane_in)
    /* init */,periodicGlidePlane(this->network().ddBase.dislocationDynamics().periodicGlidePlaneFactory.getFromKey(glidePlane->key))
    /* init */,_patches(periodicGlidePlane)
    /* init */,grain(glidePlane->grain)
    /* init */,loopType(this->flow().dot(glidePlane->n)==0? DislocationLoopIO<dim>::GLISSILELOOP : DislocationLoopIO<dim>::SESSILELOOP)
    /* init */,nA(VectorDim::Zero())
    /* init */,nAR(VectorDim::Zero())
    /* init */,_slippedArea(0.0)
    /* init */,_slippedAreaRate(0.0)
    /* init */,_rightHandedUnitNormal(VectorDim::Zero())
    /* init */,_rightHandedUnitNormal_old(VectorDim::Zero())
    /* init */,_rightHandedNormal(new ReciprocalLatticeDirectionType(grain))
    /* init */,_slipSystem(nullptr)
 
    {
        VerboseDislocationLoop(1,"Constructing DislocationLoop "<<this->tag()<<std::endl;);
    }
 
    template <int dim>
    void DislocationLoop<dim>::crossSlipBranches(std::deque<CrossSlipBranch<dim>>& csNodes,const double dt) const
    {
        std::deque<std::deque<std::pair<const LoopLinkType*,int>>> csBranches;
        if(this->network().crossSlipModel)
        {
            const auto atJunction=[this](const auto& node)
            {
                for(const auto* occurrence : node->networkNode->loopNodes())
                {
                    if(occurrence->loop().get()!=this ||
                       (occurrence!=node.get() && !occurrence->isPeriodicSubdivision()))
                    {
                        return true;
                    }
                }
                return false;
            };
            std::deque<std::pair<const LoopLinkType*,int>> currentBranch; // pair<link,cross-slip slipSystem ID>
            auto sequence=this->linkSequence();
            // Start on a physical segment so an initial seam connector cannot
            // be lost before currentBranch exists. This makes cyclic branch
            // merging independent of the loop's arbitrary starting link.
            const auto start=std::find_if(sequence.begin(),sequence.end(),[](const auto* link)
            {
                return link->hasNetworkLink();
            });
            std::rotate(sequence.begin(),start,sequence.end());
            for(const auto& link : sequence)
            {
                if(link->hasNetworkLink())
                {
                    if(this->network().crossSlipModel->isThermal() && !currentBranch.empty())
                    {
                        if(atJunction(link->source))
                        {
                            csBranches.push_back(currentBranch);
                            currentBranch.clear();
                        }
                    }
                    const auto isCSLink(this->network().crossSlipModel->isCrossSlipLink(*link->networkLink()));
                    
                    if(isCSLink.first)
                    {// a cross-slip segment
                        if(currentBranch.empty())
                        {// start of new branch
                            currentBranch.emplace_back(link,isCSLink.second.second);
                        }
                        else
                        {// existing branch
                            if(currentBranch.back().second==isCSLink.second.second)
                            {
                                currentBranch.emplace_back(link,isCSLink.second.second);
                            }
                            else
                            {// close, store, and push to currentBranch,
                                csBranches.push_back(currentBranch);
                                currentBranch.clear();
                                currentBranch.emplace_back(link,isCSLink.second.second);
                            }
                        }
                    }
                    else
                    {// not a cross-slip segment
                        if(!currentBranch.empty())
                        {// close and store branch if not empty
                            csBranches.push_back(currentBranch);
                            currentBranch.clear();
                        }
                    }
                }
                else
                {
                    if(!currentBranch.empty())
                    {// close and store branch if not empty
                        currentBranch.emplace_back(link,currentBranch.back().second);
                    }
                }
            }
            if(!currentBranch.empty())
            {// close and store branch if not empty
                csBranches.push_back(currentBranch);
                currentBranch.clear();
            }
            
            if(csBranches.size()>1)
            {// Inserted two or more branches. Merge last and first branch if possible
                if(   csBranches[0].front().first->prev==csBranches.back().back().first
                   && csBranches[0].front().second==csBranches.back().back().second
                   && (!this->network().crossSlipModel->isThermal() ||
                       !atJunction(csBranches[0].front().first->source)))
                {
                    for(typename std::deque<std::pair<const LoopLinkType*,int>>::reverse_iterator rIter = csBranches.back().rbegin();
                        rIter != csBranches.back().rend(); ++rIter)
                    {
                        csBranches[0].push_front(*rIter);
                    }
                    csBranches.pop_back();
                }
            }
            
            for(const auto& branch : csBranches)
            {
                if(branch.size())
                {
                    std::optional<ThermalCrossSlipEvent> event;
                    if(this->network().crossSlipModel->isThermal())
                    {
                        std::vector<const NetworkLinkType*> links;
                        for(const auto& pair : branch)
                        {
                            if(pair.first->hasNetworkLink())
                            {
                                links.push_back(pair.first->networkLink().get());
                            }
                        }
                        event=this->network().crossSlipModel->sampleBranch(links,dt>=0.0? dt : this->network().ddBase.simulationParameters.dt);
                        if(!event)
                        {
                            continue;
                        }
                    }
                    csNodes.emplace_back(std::deque<std::shared_ptr<LoopNodeType>>(),event? event->target : branch.back().second);
                    csNodes.back().thermal=event;
                    for(const auto& pair : branch)
                    {
                        if(!pair.first->source->isPeriodicSubdivision())
                        {// Keep physical endpoints; only periodic subdivisions are omitted.
                            csNodes.back().first.emplace_back(pair.first->source);
                        }
                    }
                    
                    if(csNodes.back().first.size())
                    {
                        if(    branch.back().first->sink!=csNodes.back().first.back()
                           && !branch.back().first->sink->isPeriodicSubdivision())
                        {
                            csNodes.back().first.emplace_back(branch.back().first->sink);
                        }
                    }
                }
            }
        }
    }
 
    template <int dim>
    DislocationLoop<dim>::~DislocationLoop()
    {
        VerboseDislocationLoop(1,"Destroying DislocationLoop "<<this->sID<<std::endl;);
    }
 
    template <int dim>
    const DislocationLoopPatches<dim>& DislocationLoop<dim>::patches() const
    {
        return _patches;
    }
 
    template <int dim>
    std::shared_ptr<typename TypeTraits<DislocationLoop<dim>>::LoopType> DislocationLoop<dim>::clone() const
    {
        return std::shared_ptr<typename TypeTraits<DislocationLoop<dim>>::LoopType>(new DislocationLoop(this->p_network(),burgers(),glidePlane));
    }
 
    template <int dim>
    void DislocationLoop<dim>::initFromFile(const std::string& fileName)
    {
        verboseDislocationLoop=TextFileParser(fileName).readScalar<int>("verboseDislocationLoop",false);
    }
 
    template <int dim>
    typename DislocationLoop<dim>::VectorDim DislocationLoop<dim>::currentAreaVector() const
    {
        VectorDim area=VectorDim::Zero();
        if(!this->loopLinks().empty())
        {
            const VectorDim origin=(*this->loopLinks().begin())->source->get_P();
            for(const auto* link : this->loopLinks())
            {
                area+=0.5*(link->source->get_P()-origin).cross(link->sink->get_P()-origin);
            }
        }
        return area;
    }
 
    template <int dim>
    double DislocationLoop<dim>::currentAreaRate() const
    {
        const VectorDim area=currentAreaVector();
        if(area.norm()<=FLT_EPSILON)
        {
            return 0.0;
        }
        VectorDim rate=VectorDim::Zero();
        for(const auto* link : this->loopLinks())
        {
            rate+=0.5*(link->source->networkNode->get_V()+link->sink->networkNode->get_V()).cross(
                link->sink->get_P()-link->source->get_P());
        }
        return area.normalized().dot(rate);
    }
 
    template <int dim>
    const double& DislocationLoop<dim>::slippedArea() const
    {
        return _slippedArea;
    }
 
    template <int dim>
    const double& DislocationLoop<dim>::slippedAreaRate() const
    {
        return _slippedAreaRate;
    }
 
    template <int dim>
    const typename DislocationLoop<dim>::VectorDim& DislocationLoop<dim>::rightHandedUnitNormal() const
    {
        return _rightHandedUnitNormal;
    }
 
    template <int dim>
    const typename DislocationLoop<dim>::VectorDim& DislocationLoop<dim>::rightHandedUnitNormalOld() const
    {
        return _rightHandedUnitNormal_old;
    }
 
    template <int dim>
    const typename DislocationLoop<dim>::ReciprocalLatticeDirectionType& DislocationLoop<dim>::rightHandedNormal() const
    {
        return *_rightHandedNormal;
    }
 
    template <int dim>
    bool DislocationLoop<dim>::isVirtualBoundaryLoop() const
    {
        return loopType==DislocationLoopIO<dim>::VIRTUALLOOP;
    }
 
    template <int dim>
    double DislocationLoop<dim>::planarSolidAngle(const VectorDim& x,
                                                         const VectorDim& planePoint,
                                                         const VectorDim& rhN,
                                                         const std::vector<std::pair<VectorDim,VectorDim>>& polygonSegments)
    {
        double temp(0.0);
        const double posNorm((x-planePoint).norm());
        const double dotProd((x-planePoint).dot(rhN));
        if(std::fabs(dotProd)>FLT_EPSILON*posNorm)
        {// x is outside the plane of the loop
            const VectorDim s(sgn(dotProd)*rhN); // s points along +n for points above, and along -n for points below
            for(const auto& pair : polygonSegments)
            {
                VectorDim e1(pair.first-x);
                const double e1Norm(e1.norm());
                if(e1Norm>FLT_EPSILON)
                {
                    e1/=e1Norm;
                    VectorDim Y1(pair.second-x);
                    const double Y1norm(Y1.norm());
                    if(Y1norm>FLT_EPSILON)
                    {
                        Y1/=Y1norm;
                        VectorDim e3(e1.cross(Y1));
                        const double e3Norm(e3.norm());
                        if(e3Norm>FLT_EPSILON)
                        {// e1 and Y1 are not align. If they are the projection on the unit sphere is a point and therefore there is no contribution to solid angle
                            e3/=e3Norm; // normalize e3
                            const VectorDim e2(e3.cross(e1));
                            const double ydy(e1.dot(Y1));
                            const double w=sqrt((1.0-ydy)/(1.0+ydy));
                            const double oneA2=sqrt(1.0+DislocationFieldBase<dim>::a2);
                            const double s3(s.dot(e3));
                            const double s3A2=sqrt(std::pow(s3,2)+DislocationFieldBase<dim>::a2);
                            temp+=2.0*s3/oneA2/s3A2*atan(s3A2*w/(oneA2-s.dot(e1)-s.dot(e2)*w));
                        }
                    }
                }
            }
        }
        return temp;
    }
 
    template <int dim>
    template <typename T>
    int DislocationLoop<dim>::sgn(const T& val)
    {
        return (val > T(0)) - (val < T(0));
    }
 
    template <int dim>
    double DislocationLoop<dim>::solidAngle(const VectorDim& x) const
    {
        double temp(0.0);
        if(glidePlane)
        {
            if(_slippedArea>FLT_EPSILON)
            {// a right-handed normal for the loop can be determined
                temp+=_patches.solidAngle(x);
            }
        }
        else
        {
            const auto linkSeq(this->linkSequence());
            assert(linkSeq.size()==4);
            
            std::vector<std::pair<VectorDim,VectorDim>> triangle0;
            triangle0.emplace_back(linkSeq[0]->source->get_P(),linkSeq[0]->sink->get_P());
            triangle0.emplace_back(linkSeq[1]->source->get_P(),linkSeq[1]->sink->get_P());
            triangle0.emplace_back(linkSeq[1]->sink->get_P(),linkSeq[0]->source->get_P());
            const VectorDim planePoint0(triangle0[0].first);
            VectorDim rhN0(VectorDim::Zero());
            for(const auto& pair : triangle0)
            {
                rhN0+= 0.5*(pair.first-planePoint0).cross(pair.second-pair.first);
            }
            const double rhN0norm(rhN0.norm());
            if(rhN0norm>FLT_EPSILON)
            {
                temp+=planarSolidAngle(x,planePoint0,rhN0/rhN0norm,triangle0);
            }
            
            std::vector<std::pair<VectorDim,VectorDim>> triangle1;
            triangle1.emplace_back(linkSeq[2]->source->get_P(),linkSeq[2]->sink->get_P());
            triangle1.emplace_back(linkSeq[3]->source->get_P(),linkSeq[3]->sink->get_P());
            triangle1.emplace_back(linkSeq[3]->sink->get_P(),linkSeq[2]->source->get_P());
            const VectorDim planePoint1(triangle1[0].first);
            VectorDim rhN1(VectorDim::Zero());
            for(const auto& pair : triangle1)
            {
                rhN1+= 0.5*(pair.first-planePoint1).cross(pair.second-pair.first);
            }
            const double rhN1norm(rhN1.norm());
            if(rhN1norm>FLT_EPSILON)
            {
                temp+=planarSolidAngle(x,planePoint1,rhN1/rhN1norm,triangle1);
            }
        }
        return temp;
    }
 
    template <int dim>
    std::tuple<double,double,double,double> DislocationLoop<dim>::loopLength() const
    {
        double freeLength=0.0;
        double boundaryLength=0.0;
        double junctionLength=0.0;
        for(const auto& link : this->loopLinks())
        {
            if(link->networkLink())
            {
                if(link->networkLink()->isBoundarySegment())
                {
                    boundaryLength+=(link->sink->get_P()-link->source->get_P()).norm();
                }
                else
                {
                    if(link->networkLink()->loopLinks().size()==1)
                    {
                        freeLength+=(link->sink->get_P()-link->source->get_P()).norm();
                    }
                    else
                    {
                        junctionLength+=(link->sink->get_P()-link->source->get_P()).norm();
                    }
                }
            }
        }
        return std::make_tuple(freeLength,junctionLength,boundaryLength,slippedArea());
    }
 
    template <int dim>
    std::shared_ptr<SlipSystem> DislocationLoop<dim>::searchSlipSystem() const
    {
        for(const auto& ss : grain.slipSystems())
        {
            if(ss.second->isSameAs(this->flow(),rightHandedNormal()))
            {
                return ss.second;
            }
        }
        return std::shared_ptr<SlipSystem>(nullptr);
    }
 
    template <int dim>
    typename DislocationLoop<dim>::MatrixDim DislocationLoop<dim>::averagePlasticDistortion() const
    {
        return -burgers()*nA.transpose()/this->network().ddBase.mesh.volume();
    }
 
    template <int dim>
    typename DislocationLoop<dim>::MatrixDim DislocationLoop<dim>::averagePlasticDistortionRate() const
    {
        return -burgers()*nAR.transpose()/this->network().ddBase.mesh.volume();
    }
 
    template <int dim>
    typename DislocationLoop<dim>::VectorDim DislocationLoop<dim>::burgers() const
    {
        return this->flow().cartesian();
    }
 
    template <int dim>
    void DislocationLoop<dim>::updateSlipSystem()
    {
        VerboseDislocationLoop(3,"DislocationLoop "<<this->tag()<<std::endl;);
        
        const SlipSystem* const oldSlipSystem(_slipSystem.get());
        if(glidePlane)
        {// a glide plane exists
            if(_slipSystem)
            {// a current slip system exists
                if(_slipSystem->isSameAs(this->flow(),rightHandedNormal()))
                {// currenst slip system still valid, don't do anything
                }
                else
                {// currenst slip system not valid
                    _slipSystem=searchSlipSystem();
                }
            }
            else
            {// a current slip system does not exist
                _slipSystem=searchSlipSystem();
            }
        }
        else
        {// no glide plane
            _slipSystem=nullptr;
        }
        
        if(_slipSystem.get()!=oldSlipSystem)
        {
            for(const auto& link : this->loopLinks())
            {
                if(link->networkLink())
                {
                    link->networkLink()->updateSlipSystem();
                }
            }
        }
        if(_slipSystem)
        {
            VerboseDislocationLoop(3,"_slipSystem.s= "<<_slipSystem->s.cartesian().transpose()<<std::endl;);
            VerboseDislocationLoop(3,"_slipSystem.unitNormal= "<<_slipSystem->unitNormal.transpose()<<std::endl;);
        }
        else
        {
            VerboseDislocationLoop(3,"slipSystem NOT FOUND"<<std::endl;);
        }
    }
 
    template <int dim>
    const std::shared_ptr<SlipSystem>&  DislocationLoop<dim>::slipSystem() const
    {
        return _slipSystem;
    }
 
    template <int dim>
    bool DislocationLoop<dim>::flippedInsideOut(const int& N,const double& areaThreshold) const
    {
        
//        for(int n=0;n<N;++n)
//        {
//            const double du((1.0*n)/N);
//        
//        
//        VectorDim newAreaVector(VectorDim::Zero());
//        const VectorDim P0((*this->loopLinks().begin())->source->get_P()-du*(*this->loopLinks().begin())->source->get_deltaP());
//        for(const auto& loopLink : this->loopLinks())
//        {
//            newAreaVector+= 0.5*(loopLink->source->get_P()-du*loopLink->source->get_deltaP()-P0).cross(loopLink->sink->get_P()-du*loopLink->sink->get_deltaP()-loopLink->source->get_P()+du*loopLink->source->get_deltaP());
//            }
//            const double newArea(newAreaVector.norm());
//            if(newArea<areaThreshold && newAreaVector.dot(_rightHandedUnitNormal)>=0.0)
//            {// area fell close to zero and rightHandedNormal did not change
//
//                CHECK THIS IMPLEMENTATION ESPECIALLY  CONDITION ON THE NORMAL.
//                ALSO, CONSIDER JUNCTIONS< SO MAYBE LOOP MUST BE ISOLATED?
//                
//                return true;
//            }
//        }
        return false;
    }
 
 
    template <int dim>
    void DislocationLoop<dim>::updateGeometry()
    {
        VerboseDislocationLoop(2,"DislocationLoop "<<this->sID<<" updateGeometry"<<std::endl;);
        _rightHandedUnitNormal_old=_rightHandedUnitNormal;
        nA.setZero();
        //    nAR.setZero();
        if(glidePlane)
        {// remove numerical errors by projecting along glidePlane->unitNormal
            if(this->loopLinks().size())
            {
                const VectorDim P0((*this->loopLinks().begin())->source->get_P());
                for(const auto& loopLink : this->loopLinks())
                {
                    nA+= 0.5*(loopLink->source->get_P()-P0).cross(loopLink->sink->get_P()-loopLink->source->get_P());
                    }
            }
            _slippedArea=nA.norm();
            _rightHandedUnitNormal= _slippedArea>FLT_EPSILON? (nA/_slippedArea).eval() : VectorDim::Zero();
            const double nnDot(_rightHandedUnitNormal.dot(glidePlane->unitNormal));
//            _rightHandedNormal= nnDot>=0.0? glidePlane->n : ReciprocalLatticeDirection<dim>(glidePlane->n*(-1));
            _rightHandedNormal.reset(nnDot>=0.0? new ReciprocalLatticeDirectionType(glidePlane->n) : new ReciprocalLatticeDirection<dim>(glidePlane->n*(-1)));
            _rightHandedUnitNormal=rightHandedNormal().cartesian().normalized();
            VerboseDislocationLoop(3,"_rightHandedUnitNormal= "<<_rightHandedUnitNormal.transpose()<<std::endl;);
        }
        else
        {
            if(this->loopLinks().size())
            {
                const VectorDim P0((*this->loopLinks().begin())->source->get_P());
                for(const auto& loopLink : this->loopLinks())
                {
                    nA+= 0.5*(loopLink->source->get_P()-P0).cross(loopLink->sink->get_P()-loopLink->source->get_P());
                }
                Plane<dim> tempPlane(P0,nA);
                bool isPlanar(true);
                for(const auto& loopLink : this->loopLinks())
                {
                    isPlanar= (isPlanar && tempPlane.contains(loopLink->source->get_P()));
                }
                if(!isPlanar)
                {
                    nA.setZero();
                }
                _slippedArea=nA.norm();
                _rightHandedUnitNormal= _slippedArea>FLT_EPSILON? (nA/_slippedArea).eval() : VectorDim::Zero();
//                _rightHandedNormal= grain.reciprocalLatticeDirection(_rightHandedUnitNormal);
                _rightHandedNormal.reset(new ReciprocalLatticeDirectionType(grain.reciprocalLatticeDirection(_rightHandedUnitNormal)));
                VerboseDislocationLoop(3,"non-glide _rightHandedUnitNormal= "<<_rightHandedUnitNormal.transpose()<<std::endl;);
            }
        }
        updateSlipSystem();
        
        if(periodicGlidePlane)
        {
            std::vector<VectorDim> linkShifts;
            std::vector<VectorDim> globalNodePos;
            
            const auto linkSeq(this->linkSequence());
            for(const auto& link : linkSeq)
            {// Collect patches of the loop
                if(link->periodicPlanePatch())
                {
                    linkShifts.push_back(link->periodicPlanePatch()->shift);
                }
                // Physical-boundary nodes define the trimmed contour, including
                // its corners. Only periodic seam subdivisions may be omitted.
                if(!link->source->isPeriodicSubdivision())
                {
                    globalNodePos.push_back(link->source->get_P());
                }
            }
            _patches.update(linkShifts,globalNodePos);
        }
    }
 
    template <int dim>
    void DislocationLoop<dim>::updateRates()
    {
        VerboseDislocationLoop(2,"DislocationLoop "<<this->sID<<" updateRates"<<std::endl;);
        nAR.setZero();
        if(glidePlane)
        {// remove numerical errors by projecting along glidePlane->unitNormal
            if(this->loopLinks().size())
            {
                for(const auto& loopLink : this->loopLinks())
                {
                    const VectorDim linkChord(loopLink->sink->get_P()-loopLink->source->get_P());
                    nAR+=0.5*(loopLink->source->networkNode->get_V()+loopLink->sink->networkNode->get_V()).cross(linkChord);
                }
            }
            _slippedAreaRate=nA.dot(nAR)/_slippedArea;
        }
    }
 
    template <int dim>
    std::vector<MeshedDislocationLoop> DislocationLoop<dim>::meshed(const double& meshSize,const double& localMeshSize) const
    {
        std::vector<MeshedDislocationLoop> temp;
        for(const auto& gPatch : patches().globalPatches())
        {
            temp.emplace_back(burgers(),*gPatch.first->glidePlane,gPatch.second,this->network(),meshSize,localMeshSize);
        }
        return temp;
    }
 
    template <int dim>
    int DislocationLoop<dim>::verboseDislocationLoop=0;
 
    template class DislocationLoop<3>;
 
}
#endif