SphericalInclusionsGenerator.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_SphericalInclusionsGenerator_cpp_
#define model_SphericalInclusionsGenerator_cpp_
 
 
#include <chrono>
#include <random>
#include <cmath>
#include <list>
#include <assert.h>
#include <Eigen/LU>
#include <Eigen/Cholesky>
#include <limits>
 
//#include <Simplex.h>
#include <SimplicialMesh.h>
#include <Polycrystal.h>
#include <PolycrystallineMaterialBase.h>
#include <LatticeModule.h>
//#include <PlaneMeshIntersection.h>
#include <DislocationDynamics/DislocationNodeIO.h>
#include <DislocationDynamics/DislocationLoopIO.h>
#include <DislocationDynamics/DislocationLoopLinkIO.h>
#include <DislocationDynamics/DislocationLoopNodeIO.h>
#include <DefectiveCrystal/DefectiveCrystalConfigIO.h>
#include <DefectiveCrystal/DefectiveCrystalAuxIO.h>
 
#include <DislocationDynamics/DislocationLinkingNumber.h>
#include <TextFileParser.h>
#include <DislocationDynamics/DislocationInjector.h>
#include <MeshBoundarySegment.h>
//#include <ConfinedDislocationObject.h>
#include <GlidePlaneModule.h>
#include <MeshModule.h>
#include <Plane.h>
#include <MicrostructureGenerator.h>
#include <SphericalInclusionsGenerator.h>
#include <PlaneLineIntersection.h>
 
namespace model
{
 
    SphericalInclusionsGenerator::SphericalInclusionsGenerator(const SphericalInclusionDensitySpecification& spec,MicrostructureGenerator& mg)
    {
        std::cout<<magentaBoldColor<<"Generating Spherical Inclusions Density"<<defaultColor<<std::endl;
        //    const double targetInclusionDensity(this->parser.readScalar<double>("targetDensity",true));
        if(spec.targetDensity>0)
        {
            //        const double inclusionsDiameterLognormalDistribution_M(this->parser.readScalar<double>("diameterLognormalDistribution_M",true) );
            //        const double inclusionsDiameterLognormalDistribution_S(this->parser.readScalar<double>("diameterLognormalDistribution_S",true));
            //        const double inclusionsDiameterLognormalDistribution_A(this->parser.readScalar<double>("diameterLognormalDistribution_A",true));
            //        const Eigen::Matrix<double,1,dim*dim> inclusionsTransformationEigenDistortion(this->parser.readMatrix<double>("crystalEigenDistortion",1,dim*dim,true));
            //        const Eigen::Matrix<double,1,dim> patternVector(this->parser.readMatrix<double>("patternVector",1,dim,true)/mg.ddBase.poly.b_SI);
            //        const double velocityReductionFactor(this->parser.readScalar<double>("velocityReductionFactor",true));
            //        const int phaseIDs(this->parser.readScalar<int>("phaseID",true));
            
            
            std::lognormal_distribution<double> distribution(log(spec.diameterLognormalDistribution_M/spec.diameterLognormalDistribution_A),spec.diameterLognormalDistribution_S);
            
            const VectorDimD patternVector(spec.patternVector_SI/mg.ddBase.poly.b_SI);
            const double patternHeigth(patternVector.norm());
            const bool applyPattern(patternHeigth>0.0);
            const VectorDimD patternDir(applyPattern? (patternVector/patternHeigth).eval() : VectorDimD::Zero());
            
            std::mt19937 generator;
            double density(0.0);
            while(density<spec.targetDensity)
            {
                
                const double diameter = distribution(generator)*spec.diameterLognormalDistribution_A/mg.ddBase.poly.b_SI;
                const double radius(0.5*diameter);
                std::pair<LatticeVector<dim>,int> pointPair=mg.ddBase.poly.randomLatticePointInMesh();
                VectorDimD P=pointPair.first.cartesian();
                const int& grainID(pointPair.second);
                
                if(applyPattern)
                {
                    const VectorDimD globalVector(mg.ddBase.poly.grain(grainID)->C2G*patternVector);
                    const VectorDimD globalDir(mg.ddBase.poly.grain(grainID)->C2G*patternDir);
                    const long long pointHeigth=std::round(P.dot(globalDir)/patternHeigth);
                    const VectorDimD O(pointHeigth*globalVector);
                    P-=(P-O).dot(globalDir)*globalDir;
                }
                if(generateSingle(mg,P,radius,spec.crystalEigenDistortion,spec.velocityReductionFactor,spec.phaseID,spec.allowOutside,spec.allowOverlap))
                {
                    density+=1.0/mg.ddBase.mesh.volume()/std::pow(mg.ddBase.poly.b_SI,3);
                    std::cout<<"inclusion density="<<density<<std::endl;
                }
                
                
            }
            
            
            //                            const double patternHeigth(currentPattern.norm());
            //                            const bool applyPattern(patternHeigth>0.0);
            
            
            //            const double inclusionsDiameterLognormalDistribution_S(isInclusionsUsed()? TextFileParser("./inputFiles/initialMicrostructure.txt").readArray<double>("inclusionsDiameterLognormalDistribution_S",true) : std::vector<double>())
            //            const double inclusionsDiameterLognormalDistribution_A(isInclusionsUsed()? TextFileParser("./inputFiles/initialMicrostructure.txt").readArray<double>("inclusionsDiameterLognormalDistribution_A",true) : std::vector<double>())
            //            const double inclusionsTransformationStrains(isInclusionsUsed()? TextFileParser("./inputFiles/initialMicrostructure.txt").readMatrix<double>("inclusionsTransformationStrains",targetInclusionDensities.size(),MicrostructureGenerator::dim*MicrostructureGenerator::dim,true) : Eigen::Matrix<double,Eigen::Dynamic,MicrostructureGenerator::dim*MicrostructureGenerator::dim>::Zero(1,MicrostructureGenerator::dim*MicrostructureGenerator::dim))
            //            const double patternVectors(isInclusionsUsed()? TextFileParser("./inputFiles/initialMicrostructure.txt").readMatrix<double>("patternVectors",targetInclusionDensities.size(),MicrostructureGenerator::dim,true) : Eigen::Matrix<double,Eigen::Dynamic,MicrostructureGenerator::dim>::Zero(1,MicrostructureGenerator::dim))
            
        }
    }
 
    SphericalInclusionsGenerator::SphericalInclusionsGenerator(const SphericalInclusionIndividualSpecification& spec,MicrostructureGenerator& mg)
    {
        //    const std::vector<double> inclusionRadii_SI(this->parser.readArray<double>("inclusionRadii_SI",true));
        
        std::cout<<magentaBoldColor<<"Generating Individual Sperical Inclusions"<<defaultColor<<std::endl;
        if(spec.radii_SI.size())
        {
            
            for(size_t k=0;k<spec.radii_SI.size();++k)
            {
                if(generateSingle(mg,spec.centers.row(k),spec.radii_SI[k]/mg.ddBase.poly.b_SI,spec.crystalEigenDistortions.row(k),spec.velocityReductionFactors[k],spec.phaseIDs[k],spec.allowOutside,spec.allowOverlap))
                {
                    std::cout<<"generated inclusion"<<std::endl;
                }
            }
        }
        
    }
 
    bool SphericalInclusionsGenerator::generateSingle(MicrostructureGenerator& mg,const VectorDimD& C,const double& R, const Eigen::Matrix<double,1,dim*dim>& eTrow, const double& vrc,const int&type, const bool& allowOutside, const bool& allowOverlap)
    {
        
        const auto searchPair = mg.ddBase.mesh.search(C);
        if(!searchPair.first)
        {
            std::cout << "Spherical inclusion center outside mesh; cannot determine crystal frame; not generated\n";
            return false;
        }
        const auto &C2G = mg.ddBase.poly.grain(searchPair.second->region->regionID)->C2G;
        const Eigen::Matrix<double,dim,dim> crystalET(
            Eigen::Map<const Eigen::Matrix<double,dim,dim>>(eTrow.data()).transpose());
        const Eigen::Matrix<double,dim,dim> eT = C2G * crystalET * C2G.transpose();
        
        
        bool isInside(allowOutside);
        if(!isInside)
        {
            isInside=(   mg.ddBase.mesh.search(C+VectorDimD::UnitX()*R).first
                      && mg.ddBase.mesh.search(C-VectorDimD::UnitX()*R).first
                      && mg.ddBase.mesh.search(C+VectorDimD::UnitY()*R).first
                      && mg.ddBase.mesh.search(C-VectorDimD::UnitY()*R).first
                      && mg.ddBase.mesh.search(C+VectorDimD::UnitZ()*R).first
                      && mg.ddBase.mesh.search(C-VectorDimD::UnitZ()*R).first);
        }
        
        if(isInside)
        {
            if(allowOverlap)
            {
                mg.insertInclusion(C,R,eT,vrc,type);
                return true;
            }
            else
            {
                bool isOutsidePrecipitates(true);
                for(const auto& incl : mg.config().sphericalInclusions())
                {
                    isOutsidePrecipitates=(isOutsidePrecipitates && ((C-incl.C).norm()>R+incl.a));
                }
                if(isOutsidePrecipitates)
                {
                    mg.insertInclusion(C,R,eT,vrc,type);
                    return true;
                }
                else
                {
                    return false;
                }
            }
        }
        else
        {
            std::cout<<"Inclusion outside mesh"<<std::endl;
            return false;
        }
    }
 
}
#endif