testThermalCrossSlipTopology.cpp

Go to the documentation of this file.

#include <DislocationDynamics/DislocationDynamicsBase.h>
#include <MicrostructureGenerator.h>
#include <DislocationNetwork.h>
#include <DefectiveCrystal/MicrostructureContainer.h>
#include <DefectiveCrystal/DefectiveCrystalBase.h>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <chrono>
#include <iostream>
#include <set>
#include <array>
#include <sstream>
#include <DislocationCrossSlip.h>
namespace model
{
struct CrossSlipTestAccess
{
    template<class N> static void queue(N& network,const std::deque<std::shared_ptr<typename N::LoopNodeType>>& nodes,int slip)
    {
        network.crossSlipMaker.csNodes.clear();
        network.crossSlipMaker.csNodes.emplace_back(nodes,slip);
        network.crossSlipMaker.csNodes.back().thermal=model::ThermalCrossSlipEvent{slip,300.0,0.5};
    }
};
}
 
 
void require(bool condition, const char* message)
{
    if(!condition)
    {
        throw std::runtime_error(message);
    }
}
 
int main(int argc, char** argv)
{
    try
    {
        require(argc == 2, "Expected repository root");
        const std::filesystem::path root(argv[1]);
        const auto folder = std::filesystem::temp_directory_path() /
            ("thermal-cross-slip-" + std::to_string(std::chrono::steady_clock::now().time_since_epoch().count()));
        const auto input = folder / "inputFiles";
        std::filesystem::create_directories(input);
        std::filesystem::create_directories(folder / "F");
        std::filesystem::copy_file(root / "Library/Materials/W.txt", input / "W.txt");
 
        std::filesystem::copy_file(root / "Library/Meshes/bicrystal_48.msh", input / "mesh.msh");
        std::filesystem::copy_file(root / "Library/DislocationDynamics/DD.txt", input / "DD.txt");
        {
            std::ifstream inputDD(input/"DD.txt");
            std::string contents((std::istreambuf_iterator<char>(inputDD)),std::istreambuf_iterator<char>());
            const auto pos=contents.find("crossSlipModel=0;");
            require(pos!=std::string::npos,"Missing cross-slip setting");
            contents.replace(pos,17,"crossSlipModel=1;");
            std::ofstream outputDD(input/"DD.txt");
            outputDD<<contents;
        }
        std::filesystem::create_directories(folder / "config");
        model::SimplicialMesh<3> mesh((input/"mesh.msh").string(),Eigen::Matrix3d::Identity(),Eigen::Vector3d::Zero(),{-1});
        std::array<int,3> seeds{-1,-1,-1};
        for(const auto& region : mesh.regions())
        {
            for(const auto& face : region.second->faces())
            {
                const auto& pair=face.second->periodicFacePair;
                if(pair.second)
                {
                    for(int axis=0;axis<3;++axis)
                    {
                        if(pair.first.normalized().cross(Eigen::Vector3d::Unit(axis)).norm()<1.e-8)
                        {
                            seeds[axis]=face.first;
                        }
                    }
                }
            }
        }
        require(seeds[0]>=0 && seeds[1]>=0,"Missing periodic seeds");
        model::StaticID<model::Polyhedron<3,3-1>>::force_count(0);
        const double length=3671.9477120460206;
        for(int refinement : {1,2,12})
        {
        for(int scenario=0;scenario<6;++scenario)
        {
            const bool full=scenario==1;
            const bool periodic=scenario!=2;
            {
                std::ofstream f(input/"DefectiveCrystal.txt");
                f<<"physics=DislocationDynamics;\nuseFEM=0;\nNsteps=0;\nmaxResolveSteps=0;\ndtMax=1;\noutputFrequency=1;\noutputBinary=0;\nstartAtTimeStep=0;\nperiodicImageSize="<<(full ? "0 0 0" : "0")<<";\nEwaldLengthFactor=0;\n";
                std::ofstream p(input/"polycrystal.txt");
                p<<std::setprecision(17)<<"materialFile=W.txt;\nmeshFile=mesh.msh;\nabsoluteTemperature=300;\nC2G1=1 0 0 0 1 0 0 0 1;\nC2G2=0.8 0 0.6 0 1 0 -0.6 0 0.8;\nF="<<length<<" 0 0 0 "<<length<<" 0 0 0 "<<length<<";\nX0=0 0 0;\nperiodicFaceIDs=";
                if(periodic)
                {
                    p<<(full ? -1 : seeds[1]);
                }
                p<<";\n";
            }
            model::DefectiveCrystalBase<3> base(folder.string());
            model::MicrostructureContainer<3> container(base);
            using Network=model::DislocationNetwork<3>;
            using V=Eigen::Vector3d;
            Network network(container);
            const auto grain=base.poly.grain(1);
            const auto slip=grain->slipSystems().begin()->second;
            int other=-1;
            for(size_t k=0;k<grain->slipSystems().size();++k)
            {
                const auto candidate=grain->slipSystems().at(k);
                if((candidate->s-slip->s).squaredNorm()==0 && candidate->unitNormal.cross(slip->unitNormal).norm()>0.1)
                {
                    other=k;
                    break;
                }
            }
            require(other>=0,"No cross-slip system");
            const auto css=grain->slipSystems().at(other);
            const V b=slip->s.cartesian();
            const V center=scenario==4 ? V::Zero().eval() : V(length*V(-0.5,periodic ? 0.49 : 0.0,0));
            const auto plane=base.dislocationDynamics().glidePlaneFactory.getFromKey(model::GlidePlaneKey<3>(center,slip->n));
            const long spacing=std::max(1L,static_cast<long>(20.0/css->n.planeSpacing()));
            const long height=std::lround(double(css->n.closestPlaneIndexOfPoint(center))/spacing)*spacing;
            model::PlanePlaneIntersection<3> intersection(plane->P,plane->unitNormal,height*css->n.planeSpacing()*css->unitNormal,css->unitNormal);
            V origin=intersection.P+(center-intersection.P).dot(intersection.d)*intersection.d;
            V along=intersection.d;
            if(scenario==4 && along.x()<0.0)
            {
                along=-along;
            }
            V side=plane->unitNormal.cross(along).normalized();
            if(scenario==4)
            {
                origin=-0.15*length*along;
                if(side.x()>0.0)
                {
                    side=-side;
                }
            }
            if(scenario==5)
            {
                origin+=(0.5*length+2.0-origin.y())/along.y()*along;
            }
            const V a=origin-0.15*length*along;
            const V middle=origin;
            const V end=origin+0.15*length*along;
            std::vector<V> polygon{a,middle,end,end+100.0*side,a+100.0*side};
            if(scenario==5)
            {
                polygon[1]-=4.0/side.y()*side;
            }
            if(scenario==3)
            {
                polygon[0]+=0.2*side;
                polygon[1]+=0.2*side;
                polygon[2]+=2.0*css->n.planeSpacing()/css->unitNormal.dot(side)*side;
            }
            if(scenario!=3 && scenario!=5)
            {
                polygon.clear();
                for(int k=0;k<=refinement;++k)
                {
                    polygon.push_back(a+(end-a)*(double(k)/refinement));
                }
                polygon.push_back(end+100.0*side);
                polygon.push_back(a+100.0*side);
            }
            model::MicrostructureGenerator generator(base);
            require(generator.insertJunctionLoop(polygon,base.dislocationDynamics().periodicGlidePlaneFactory.get(plane->key),b,plane->unitNormal,plane->P,1,model::DislocationLoopIO<3>::GLISSILELOOP),"Parent insertion failed");
            network.initializeConfiguration(generator.config(),std::ofstream(),std::ofstream());
 
            const auto parent=network.loops().begin()->second.lock();
            std::deque<std::shared_ptr<Network::LoopNodeType>> branch;
            for(const auto& point : std::vector<V>{polygon.front(),polygon[polygon.size()-3]})
            {
                for(const auto* link : parent->linkSequence())
                {
                    if((link->source->get_P()-point).norm()<1.e-6 && !link->source->isPeriodicSubdivision())
                    {
                        branch.push_back(link->source);
                        break;
                    }
                }
            }
            require(branch.size()==2,"Missing branch anchors");
            if(scenario==3 || scenario==5)
            {
                for(const auto* link : parent->linkSequence())
                {
                    if((link->source->get_P()-polygon[1]).norm()<1.e-6 && !link->source->isPeriodicSubdivision())
                    {
                        branch.insert(branch.begin()+1,link->source);
                        break;
                    }
                }
            }
            // A second loop shares the final endpoint. Its plane either permits
            // alignment, or locks the displaced endpoint against that motion.
            const auto endpoint=branch.back()->networkNode;
            const V normal=scenario==3 ? css->unitNormal : plane->unitNormal;
            const auto extraPlane=base.dislocationDynamics().glidePlaneFactory.getFromKey(model::GlidePlaneKey<3>(endpoint->get_P(),grain->reciprocalLatticeDirection(normal)));
            const auto extra=network.loops().create(b,extraPlane);
            V u=normal.cross(along).normalized();
            if(scenario==4 && u.x()>0.0)
            {
                u=-u;
            }
            const auto e1=network.networkNodes().create(V(endpoint->get_P()-30.0*along));
            const auto e2=network.networkNodes().create(V(endpoint->get_P()-30.0*along+20.0*u));
            std::vector<std::shared_ptr<Network::LoopNodeType>> ring;
            for(const auto& node : std::vector{endpoint,e1,e2})
            {
                ring.push_back(network.loopNodes().create(extra,node,node->get_P(),extra->periodicGlidePlane->getPatch(V::Zero()),std::make_pair(nullptr,nullptr)));
            }
            network.insertLoop(extra,ring);
            network.updateGeometry();
            const auto snapshot=[&]()
            {
                std::ostringstream out;
                out<<std::setprecision(17);
                for(const auto& entry : network.networkNodes())
                {
                    const auto node=entry.second.lock();
                    out<<node->sID<<' '<<node->get_P().transpose()<<' '<<node->get_V().transpose()<<'\n';
                }
                for(const auto& entry : network.loops())
                {
                    out<<entry.first<<'\n';
                    for(const auto* link : entry.second.lock()->linkSequence())
                    {
                        out<<link->source->sID<<' '<<link->source->networkNode->sID<<' '<<link->source->get_P().transpose()<<' '<<link->sink->sID<<'\n';
                    }
                }
                return out.str();
            };
            const auto before=snapshot();
            const auto endpointID=endpoint->sID;
            const auto middleNode=branch[1]->networkNode;
            const V oldMiddle=middleNode->get_P();
            model::CrossSlipTestAccess::queue(network,branch,other);
            network.crossSlipMaker.execute();
            if(scenario==3)
            {
                require(snapshot()==before,"Rejected alignment changed network");
                continue;
            }
            if(scenario==5)
            {
                require((middleNode->get_P()-oldMiddle).norm()>0.5*length,"Alignment did not cross the periodic face");
                require(network.networkNodes().get(middleNode->sID)==middleNode,"Alignment lost network identity at seam");
            }
            branch.clear();
            ring.clear();
            require(network.loops().size()==3,"Cross-slip seed was discarded");
            require(network.networkNodes().get(endpointID)==endpoint && endpoint->loopNodes().size()==2,"Shared endpoint was detached");
            std::shared_ptr<Network::LoopType> seed;
            for(const auto& entry : network.loops())
            {
                const auto loop=entry.second.lock();
                if(loop!=parent && loop!=extra)
                {
                    seed=loop;
                }
            }
            require(seed && seed->loopLinks().size()>=4 && seed->slipSystem(),"Invalid zero-area seed");
            double perimeter=0.0;
            for(const auto* link : seed->linkSequence())
            {
                perimeter+=(link->sink->get_P()-link->source->get_P()).norm();
            }
            require(std::abs(perimeter-600.0)<1.e-6,"Nucleus length depends on segmentation");
            require(seed->currentAreaVector().norm()<1.e-6,"Seed added slipped area before motion");
            size_t seams=0,cancelled=0;
            for(const auto* link : seed->linkSequence())
            {
                seams+=link->source->isPeriodicSubdivision();
                if(link->networkLink() && link->networkLink()->loopLinks().size()==2)
                {
                    require(link->networkLink()->hasZeroBurgers(),"Parent branch Burgers vector not cancelled");
                    ++cancelled;
                }
                require((link->source->get_P()+link->source->periodicPlanePatch()->shift-link->source->networkNode->get_P()).norm()<1.e-6,"Incorrect periodic mapping");
            }
            require(cancelled>=2,"No cancellation on original branch");
            require((seams>0)==(periodic && scenario!=4),"Periodic crossing was not preserved");
            const auto expected=network.glideSolver->getNodeVelocities();
            size_t k=0;
            for(const auto& entry : network.networkNodes())
            {
                const auto node=entry.second.lock();
                const V stored=node->get_V();
                node->set_V(expected.segment<3>(3*k),false);
                require((stored-node->get_V()).norm()<1.e-7,"Stale nodal velocity after cross-slip");
                ++k;
            }
            // Drive only the free side of the seed on its new plane. This
            // checks the ordinary motion/boundary walk after seed creation.
            for(const auto& entry : network.networkNodes())
            {
                entry.second.lock()->set_V(V::Zero(),false);
            }
            for(const auto* link : seed->linkSequence())
            {
                const auto node=link->source->networkNode;
                if(node->loopNodes().size()==1)
                {
                    node->set_V(css->unitNormal.cross(along).normalized(),false);
                }
            }
            network.moveNodes(0.1);
            require(seed->currentAreaVector().norm()>1.e-6,"Seed failed to expand under motion");
            require(endpoint->loopNodes().size()==2,"Motion lost shared endpoint");
        }
        }
        std::filesystem::remove_all(folder);
        std::cout<<"PASS: thermal nucleus length, periodic mapping, shared endpoint, fresh velocities and atomic rejection\n";
        return 0;
    }
    catch(const std::exception& error)
    {
        std::cerr<<error.what()<<'\n';
        return 1;
    }
}