oILAB
Loading...
Searching...
No Matches
testGenerateGBs.cpp

Given a tilt-axis of a 3D lattice, this example demonstrates the construction of multiple tilt GBs of varying misorientations and inclinations, and the calculation of grain boundaries' disconnection modes

  1. Define types
    using VectorDimI = LatticeCore<3>::VectorDimI;
    using IntScalarType = LatticeCore<3>::IntScalarType;
  2. Instantiate a lattice \(\mathcal A\)
    const auto A(
    TextFileParser("bicrystal_3d.txt").readMatrix<double, 3, 3>("A", true));
    Lattice<3> lattice(A);
    std::cout << "Lattice A = " << std::endl;
    std::cout << lattice.latticeBasis << std::endl;
  3. Specify the Cartesian coordinates of an axis. This should be parallel to a reciprocal lattice vector of \(\mathcal A\)
    const auto axis(TextFileParser("bicrystal_3d.txt")
    .readMatrix<double, 3, 1>("axis", true));
    ReciprocalLatticeVector<3> rv(
    lattice.reciprocalLatticeDirection(axis).reciprocalLatticeVector());
    std::cout << "Cartesian coordinates of axis = " << std::endl;
    std::cout << rv.cartesian().transpose() << std::endl;
  4. Generate all rotations \(\mathbf R\) about the given axis that result in a coincidence relation between \(\mathcal A\) and \(\mathbf R\mathcal A\), and form the corresponding bicrystals
    const auto &coincidentLattices =
    lattice.generateCoincidentLattices(rv, 150, 150);
  5. Loop over the generated bicrystals, i.e., loop over misorientation angles
    for (const auto &rotation : coincidentLattices) {
    // Loop over misorientation angles
    std::cout << "###################################################"
    << std::endl;
  6. SNF output for each misorientation
    double theta =
    acos((rotation.trace() - 1.0) / 2.0) * 180 / std::numbers::pi;
    std::cout << "Misorientation angle = " << std::setprecision(20) << theta
    << "; ";
    Lattice<3> latticeB(lattice.latticeBasis, rotation);
    // BiCrystal<3>
    // bc(lattice,Lattice<3>(lattice.latticeBasis,rotation),false);
    BiCrystal<3> bc(lattice, latticeB, false);
    std::cout << "Sigma = " << std::setprecision(20) << bc.sigma << std::endl;
    std::cout << std::endl;
    std::cout << "Lattice B = " << std::endl;
    std::cout << std::setprecision(20) << rotation * lattice.latticeBasis
    << std::endl;
    std::cout << std::endl;
    std::cout << "Parallel CSL basis Cp= " << std::endl;
    std::cout << std::setprecision(20) << bc.csl.latticeBasis << std::endl;
    std::cout << std::endl;
    std::cout << "Parallel DSCL basis Dp = " << std::endl;
    std::cout << std::setprecision(20) << bc.dscl.latticeBasis << std::endl;
    std::cout << std::endl;
  7. Output the invariance property of the CSL in the following steps: 1) compute reduced basis vectors \(\mathbf d_1\) and \(\mathbf d_2\) for the DSCL, 2) compute the CSL shifts \(\mathbf s_1\) and \(\mathbf s_2\) when lattice \(\mathcal A\) is displaced by \(\mathbf d_1\) and \(\mathbf d_2\), respectively.
    auto reducedDsclBasis = RLLL(bc.dscl.latticeBasis, 0.75);
    auto U_Dscl = reducedDsclBasis.unimodularMatrix();
    LatticeVector<3> d1(bc.dscl);
    d1 << U_Dscl.col(0).template cast<IntScalarType>();
    std::cout << "Reduced DSCL basis vectors:" << std::endl;
    std::cout << "d1 = ";
    std::cout << std::setprecision(20) << d1.cartesian().transpose()
    << std::endl;
    std::cout << "Integer coordinates of d1:";
    std::cout << std::setprecision(20) << d1.transpose() << std::endl;
    std::cout << std::endl;
    LatticeVector<3> d2(bc.dscl);
    d2 << U_Dscl.col(1).template cast<IntScalarType>();
    std::cout << "d2 = ";
    std::cout << std::setprecision(20) << d2.cartesian().transpose()
    << std::endl;
    std::cout << "Integer coordinates of d2:";
    std::cout << std::setprecision(20) << d2.transpose() << std::endl;
    std::cout << std::endl;
    LatticeVector<3> d3(bc.dscl);
    d3 << U_Dscl.col(2).template cast<IntScalarType>();
    std::cout << "d3 = ";
    std::cout << std::setprecision(20) << d3.cartesian().transpose()
    << std::endl;
    std::cout << "Integer coordinates of d3:";
    std::cout << std::setprecision(20) << d3.transpose() << std::endl;
    std::cout << std::endl;
    // shift vectors corresponding to d1, d2, and d3
    LatticeVector<3> s1(bc.dscl), s2(bc.dscl), s3(bc.dscl);
    s1 << bc.LambdaA * d1;
    s2 << bc.LambdaA * d2;
    s3 << bc.LambdaA * d3;
    Lattice<3> reducedCsl(RLLL(bc.csl.latticeBasis, 0.75).reducedBasis());
    std::cout << "Reduced shift vectors: " << std::endl;
    Eigen::Vector3d s1_coordinates_in_reduced_csl =
    reducedCsl.latticeBasis.inverse() * s1.cartesian();
    Eigen::Vector3d s2_coordinates_in_reduced_csl =
    reducedCsl.latticeBasis.inverse() * s2.cartesian();
    Eigen::Vector3d s3_coordinates_in_reduced_csl =
    reducedCsl.latticeBasis.inverse() * s3.cartesian();
    Eigen::Vector3d s1_coordinates_modulo =
    s1_coordinates_in_reduced_csl.array() -
    s1_coordinates_in_reduced_csl.array().round();
    Eigen::Vector3d s2_coordinates_modulo =
    s2_coordinates_in_reduced_csl.array() -
    s2_coordinates_in_reduced_csl.array().round();
    Eigen::Vector3d s3_coordinates_modulo =
    s3_coordinates_in_reduced_csl.array() -
    s3_coordinates_in_reduced_csl.array().round();
    std::cout << "s1 = ";
    std::cout << std::setprecision(20)
    << (reducedCsl.latticeBasis * s1_coordinates_modulo).transpose()
    << std::endl;
    std::cout << "s2 = ";
    std::cout << std::setprecision(20)
    << (reducedCsl.latticeBasis * s2_coordinates_modulo).transpose()
    << std::endl;
    std::cout << "s3 = ";
    std::cout << std::setprecision(20)
    << (reducedCsl.latticeBasis * s3_coordinates_modulo).transpose()
    << std::endl;
    std::cout << std::endl;
    1. Generate grain boundaries of varying inclinations
      auto gbSet(bc.generateGrainBoundaries(
      bc.A.latticeDirection(rv.cartesian()), 60));
    2. Loop over inclination angles
      int gbCount = 0;
      ReciprocalLatticeVector<3> refnA(bc.A);
      std::cout << "GBs of varying inclination (measured with respect to the "
      "first grain boundary)"
      << std::endl;
      std::cout << "-----------------------------------------------------------"
      "------------------"
      << std::endl;
      std::cout << "-- CAUTION: The integer coordinates of GB normals are "
      "w.r.t the reciprocal --"
      << std::endl;
      std::cout << "-- basis of the primitive unit cell. "
      " --"
      << std::endl;
      std::cout << "-----------------------------------------------------------"
      "------------------"
      << std::endl;
      for (const auto &gb : gbSet) {
      1. For each GB, compute its period and the Burgers vector of the glide disconnection
        try {
        LatticeVector<3> glide(
        rv.cross(gb.second.nA.reciprocalLatticeVector()).latticeVector());
        LatticeVector<3> burgersVector(
        bc.getLatticeDirectionInD(glide).latticeVector());
        LatticeVector<3> periodVector(
        bc.getLatticeDirectionInC(glide).latticeVector());
        ReciprocalLatticeVector<3> rvInA =
        gb.second.nA.reciprocalLatticeVector();
        ReciprocalLatticeVector<3> rvInCsl =
        bc.getReciprocalLatticeDirectionInC(rvInA)
        .reciprocalLatticeVector();
      2. Note the reference GB with respect to which inclination angles are measured
        if (gbCount == 0)
        refnA = gb.second.nA.reciprocalLatticeVector();
      3. Output properties of GBs whose period is less than 100 Angstrom
        double epsilon = 1e-8;
        if (gbCount == 0 || periodVector.cartesian().norm() < 100) {
        double cosAngle = refnA.cartesian().normalized().dot(
        gb.second.nA.cartesian().normalized());
        if (cosAngle - 1 > -epsilon)
        cosAngle = 1.0;
        if (cosAngle + 1 < epsilon)
        cosAngle = -1.0;
        std::cout << gbCount + 1
        << ") Inclination = " << std::setprecision(20)
        << acos(cosAngle) * 180 / std::numbers::pi << std::endl;
        Eigen::Vector3d nAglobalCoords = gb.second.nA.cartesian();
        Eigen::Vector3d nBglobalCoords =
        rotation.transpose() * gb.second.nB.cartesian();
        std::cout << "nA = " << std::setprecision(20)
        << nAglobalCoords.transpose() << std::endl;
        std::cout << "nB = " << std::setprecision(20)
        << nBglobalCoords.transpose() << std::endl;
        /* Change the above two lines as follows if you need the GB normals
        * in the coordinate system of A */
        /*
        std::cout << "nA = " << gb.second.nA << std::endl;
        std::cout << "nB = " << gb.second.nB << std::endl;
        */
        nAglobalCoords = nAglobalCoords.array().round().cwiseAbs();
        nBglobalCoords = nBglobalCoords.array().round().cwiseAbs();
        std::sort(nAglobalCoords.data(), nAglobalCoords.data() + 3);
        std::sort(nBglobalCoords.data(), nBglobalCoords.data() + 3);
        if ((nAglobalCoords - nBglobalCoords).norm() < FLT_EPSILON)
        std::cout << "STGB" << std::endl;
        std::cout << "GB period = " << std::setprecision(20)
        << periodVector.cartesian().norm() << std::endl;
        std::cout << "CSL plane distance (Height)= "
        << std::setprecision(20)
        << 1.0 / rvInCsl.cartesian().norm() << std::endl;
        std::cout << "Glide disconnection Burgers vector = "
        << std::setprecision(20)
        << burgersVector.cartesian().transpose()
        << "; norm = " << burgersVector.cartesian().norm()
        << std::endl;
        std::cout << "Step height of glide disconnection = "
        << std::setprecision(20)
        << gb.second.stepHeight(burgersVector) << std::endl;
        std::cout << "Step height of disconnection with Burgers vector d1= "
        << std::setprecision(20) << gb.second.stepHeight(d1)
        << std::endl;
        std::cout << "Step height of disconnection with Burgers vector d2= "
        << std::setprecision(20) << gb.second.stepHeight(d2)
        << std::endl;
        std::cout << "Step height of disconnection with Burgers vector d3= "
        << std::setprecision(20) << gb.second.stepHeight(d3)
        << std::endl;
        std::cout << "-----------------------------------------------------"
        "------------------------"
        << std::endl;
        gbCount++;
        }

Full code:

#include "../../include/IO/TextFileParser.h"
#include "../../include/Lattices/BiCrystal.h"
#include "../../include/Lattices/LatticeModule.h"
#include <numbers>
using namespace oILAB;
int main() {
using VectorDimI = LatticeCore<3>::VectorDimI;
using IntScalarType = LatticeCore<3>::IntScalarType;
const auto A(
TextFileParser("bicrystal_3d.txt").readMatrix<double, 3, 3>("A", true));
Lattice<3> lattice(A);
std::cout << "Lattice A = " << std::endl;
std::cout << lattice.latticeBasis << std::endl;
const auto axis(TextFileParser("bicrystal_3d.txt")
.readMatrix<double, 3, 1>("axis", true));
ReciprocalLatticeVector<3> rv(
lattice.reciprocalLatticeDirection(axis).reciprocalLatticeVector());
std::cout << "Cartesian coordinates of axis = " << std::endl;
std::cout << rv.cartesian().transpose() << std::endl;
const auto &coincidentLattices =
lattice.generateCoincidentLattices(rv, 150, 150);
for (const auto &rotation : coincidentLattices) {
// Loop over misorientation angles
std::cout << "###################################################"
<< std::endl;
try {
double theta =
acos((rotation.trace() - 1.0) / 2.0) * 180 / std::numbers::pi;
std::cout << "Misorientation angle = " << std::setprecision(20) << theta
<< "; ";
Lattice<3> latticeB(lattice.latticeBasis, rotation);
// BiCrystal<3>
// bc(lattice,Lattice<3>(lattice.latticeBasis,rotation),false);
BiCrystal<3> bc(lattice, latticeB, false);
std::cout << "Sigma = " << std::setprecision(20) << bc.sigma << std::endl;
std::cout << std::endl;
std::cout << "Lattice B = " << std::endl;
std::cout << std::setprecision(20) << rotation * lattice.latticeBasis
<< std::endl;
std::cout << std::endl;
std::cout << "Parallel CSL basis Cp= " << std::endl;
std::cout << std::setprecision(20) << bc.csl.latticeBasis << std::endl;
std::cout << std::endl;
std::cout << "Parallel DSCL basis Dp = " << std::endl;
std::cout << std::setprecision(20) << bc.dscl.latticeBasis << std::endl;
std::cout << std::endl;
auto reducedDsclBasis = RLLL(bc.dscl.latticeBasis, 0.75);
auto U_Dscl = reducedDsclBasis.unimodularMatrix();
LatticeVector<3> d1(bc.dscl);
d1 << U_Dscl.col(0).template cast<IntScalarType>();
std::cout << "Reduced DSCL basis vectors:" << std::endl;
std::cout << "d1 = ";
std::cout << std::setprecision(20) << d1.cartesian().transpose()
<< std::endl;
std::cout << "Integer coordinates of d1:";
std::cout << std::setprecision(20) << d1.transpose() << std::endl;
std::cout << std::endl;
LatticeVector<3> d2(bc.dscl);
d2 << U_Dscl.col(1).template cast<IntScalarType>();
std::cout << "d2 = ";
std::cout << std::setprecision(20) << d2.cartesian().transpose()
<< std::endl;
std::cout << "Integer coordinates of d2:";
std::cout << std::setprecision(20) << d2.transpose() << std::endl;
std::cout << std::endl;
LatticeVector<3> d3(bc.dscl);
d3 << U_Dscl.col(2).template cast<IntScalarType>();
std::cout << "d3 = ";
std::cout << std::setprecision(20) << d3.cartesian().transpose()
<< std::endl;
std::cout << "Integer coordinates of d3:";
std::cout << std::setprecision(20) << d3.transpose() << std::endl;
std::cout << std::endl;
// shift vectors corresponding to d1, d2, and d3
LatticeVector<3> s1(bc.dscl), s2(bc.dscl), s3(bc.dscl);
s1 << bc.LambdaA * d1;
s2 << bc.LambdaA * d2;
s3 << bc.LambdaA * d3;
Lattice<3> reducedCsl(RLLL(bc.csl.latticeBasis, 0.75).reducedBasis());
std::cout << "Reduced shift vectors: " << std::endl;
Eigen::Vector3d s1_coordinates_in_reduced_csl =
reducedCsl.latticeBasis.inverse() * s1.cartesian();
Eigen::Vector3d s2_coordinates_in_reduced_csl =
reducedCsl.latticeBasis.inverse() * s2.cartesian();
Eigen::Vector3d s3_coordinates_in_reduced_csl =
reducedCsl.latticeBasis.inverse() * s3.cartesian();
Eigen::Vector3d s1_coordinates_modulo =
s1_coordinates_in_reduced_csl.array() -
s1_coordinates_in_reduced_csl.array().round();
Eigen::Vector3d s2_coordinates_modulo =
s2_coordinates_in_reduced_csl.array() -
s2_coordinates_in_reduced_csl.array().round();
Eigen::Vector3d s3_coordinates_modulo =
s3_coordinates_in_reduced_csl.array() -
s3_coordinates_in_reduced_csl.array().round();
std::cout << "s1 = ";
std::cout << std::setprecision(20)
<< (reducedCsl.latticeBasis * s1_coordinates_modulo).transpose()
<< std::endl;
std::cout << "s2 = ";
std::cout << std::setprecision(20)
<< (reducedCsl.latticeBasis * s2_coordinates_modulo).transpose()
<< std::endl;
std::cout << "s3 = ";
std::cout << std::setprecision(20)
<< (reducedCsl.latticeBasis * s3_coordinates_modulo).transpose()
<< std::endl;
std::cout << std::endl;
auto gbSet(bc.generateGrainBoundaries(
bc.A.latticeDirection(rv.cartesian()), 60));
int gbCount = 0;
ReciprocalLatticeVector<3> refnA(bc.A);
std::cout << "GBs of varying inclination (measured with respect to the "
"first grain boundary)"
<< std::endl;
std::cout << "-----------------------------------------------------------"
"------------------"
<< std::endl;
std::cout << "-- CAUTION: The integer coordinates of GB normals are "
"w.r.t the reciprocal --"
<< std::endl;
std::cout << "-- basis of the primitive unit cell. "
" --"
<< std::endl;
std::cout << "-----------------------------------------------------------"
"------------------"
<< std::endl;
for (const auto &gb : gbSet) {
// Loop over inclinations
try {
LatticeVector<3> glide(
rv.cross(gb.second.nA.reciprocalLatticeVector()).latticeVector());
LatticeVector<3> burgersVector(
bc.getLatticeDirectionInD(glide).latticeVector());
LatticeVector<3> periodVector(
bc.getLatticeDirectionInC(glide).latticeVector());
ReciprocalLatticeVector<3> rvInA =
gb.second.nA.reciprocalLatticeVector();
ReciprocalLatticeVector<3> rvInCsl =
bc.getReciprocalLatticeDirectionInC(rvInA)
.reciprocalLatticeVector();
if (gbCount == 0)
refnA = gb.second.nA.reciprocalLatticeVector();
double epsilon = 1e-8;
if (gbCount == 0 || periodVector.cartesian().norm() < 100) {
double cosAngle = refnA.cartesian().normalized().dot(
gb.second.nA.cartesian().normalized());
if (cosAngle - 1 > -epsilon)
cosAngle = 1.0;
if (cosAngle + 1 < epsilon)
cosAngle = -1.0;
std::cout << gbCount + 1
<< ") Inclination = " << std::setprecision(20)
<< acos(cosAngle) * 180 / std::numbers::pi << std::endl;
Eigen::Vector3d nAglobalCoords = gb.second.nA.cartesian();
Eigen::Vector3d nBglobalCoords =
rotation.transpose() * gb.second.nB.cartesian();
std::cout << "nA = " << std::setprecision(20)
<< nAglobalCoords.transpose() << std::endl;
std::cout << "nB = " << std::setprecision(20)
<< nBglobalCoords.transpose() << std::endl;
/* Change the above two lines as follows if you need the GB normals
* in the coordinate system of A */
/*
std::cout << "nA = " << gb.second.nA << std::endl;
std::cout << "nB = " << gb.second.nB << std::endl;
*/
nAglobalCoords = nAglobalCoords.array().round().cwiseAbs();
nBglobalCoords = nBglobalCoords.array().round().cwiseAbs();
std::sort(nAglobalCoords.data(), nAglobalCoords.data() + 3);
std::sort(nBglobalCoords.data(), nBglobalCoords.data() + 3);
if ((nAglobalCoords - nBglobalCoords).norm() < FLT_EPSILON)
std::cout << "STGB" << std::endl;
std::cout << "GB period = " << std::setprecision(20)
<< periodVector.cartesian().norm() << std::endl;
std::cout << "CSL plane distance (Height)= "
<< std::setprecision(20)
<< 1.0 / rvInCsl.cartesian().norm() << std::endl;
std::cout << "Glide disconnection Burgers vector = "
<< std::setprecision(20)
<< burgersVector.cartesian().transpose()
<< "; norm = " << burgersVector.cartesian().norm()
<< std::endl;
std::cout << "Step height of glide disconnection = "
<< std::setprecision(20)
<< gb.second.stepHeight(burgersVector) << std::endl;
std::cout << "Step height of disconnection with Burgers vector d1= "
<< std::setprecision(20) << gb.second.stepHeight(d1)
<< std::endl;
std::cout << "Step height of disconnection with Burgers vector d2= "
<< std::setprecision(20) << gb.second.stepHeight(d2)
<< std::endl;
std::cout << "Step height of disconnection with Burgers vector d3= "
<< std::setprecision(20) << gb.second.stepHeight(d3)
<< std::endl;
std::cout << "-----------------------------------------------------"
"------------------------"
<< std::endl;
gbCount++;
}
} catch (std::runtime_error &e) {
std::cout << e.what() << std::endl;
std::cout << "Moving onto the next inclination" << std::endl;
}
}
} catch (std::runtime_error &e) {
std::cout << e.what() << std::endl;
std::cout << "Moving on the the next misorientation" << std::endl;
}
}
return 0;
}
oILAB::BiCrystal::generateGrainBoundaries
std::enable_if< dm==2||dm==3, std::map< IntScalarType, Gb< dm > > >::type generateGrainBoundaries(const LatticeDirection< dim > &d, int div=30) const
Given a tilt axis , that belongs to lattices or , this function generate a set of tilt GBs....
Definition BiCrystal.cpp:466
oILAB::BiCrystal::LambdaA
const MatrixDimI LambdaA
Shift tensor describes the shift in the CSL when lattice is shifted by a DSCL vector.
Definition BiCrystal.h:119
oILAB::BiCrystal::getLatticeDirectionInD
LatticeDirection< dim > getLatticeDirectionInD(const LatticeVector< dim > &v) const
Definition BiCrystal.cpp:336
oILAB::BiCrystal::getLatticeDirectionInC
LatticeDirection< dim > getLatticeDirectionInC(const LatticeVector< dim > &v) const
Definition BiCrystal.cpp:306
oILAB::BiCrystal::A
const Lattice< dim > & A
Definition BiCrystal.h:52
oILAB::BiCrystal::csl
const Lattice< dim > csl
CSL lattice .
Definition BiCrystal.h:102
oILAB::BiCrystal::getReciprocalLatticeDirectionInC
ReciprocalLatticeDirection< dim > getReciprocalLatticeDirectionInC(const ReciprocalLatticeVector< dim > &v) const
Definition BiCrystal.cpp:400
oILAB::BiCrystal::dscl
const Lattice< dim > dscl
DCSL lattice .
Definition BiCrystal.h:106
oILAB::BiCrystal::sigma
const int sigma
if , else
Definition BiCrystal.h:78
oILAB::Lattice
Lattice class.
Definition Lattice.h:31
oILAB::Lattice::latticeBasis
const MatrixDimD latticeBasis
Definition Lattice.h:41
oILAB::LatticeVector
LatticeVector class.
Definition LatticeVector.h:18
oILAB::LatticeVector::cartesian
VectorDimD cartesian() const
Definition LatticeVector.cpp:132
oILAB::RLLL::reducedBasis
const MatrixType & reducedBasis() const
Definition RLLL.cpp:185
oILAB::ReciprocalLatticeVector::cross
std::enable_if< dm==2, LatticeDirection< dim > >::type cross(const ReciprocalLatticeVector< dim > &other) const
Definition ReciprocalLatticeVector.cpp:170
oILAB::TextFileParser
Definition TextFileParser.h:125
main
int main(int argc, char **argv)
Definition main.cpp:12
oILAB::LatticeCore::VectorDimI
Eigen::Matrix< IntScalarType, dim, 1 > VectorDimI
Definition LatticeCore.h:23
oILAB::LatticeCore::IntScalarType
long long int IntScalarType
Definition LatticeCore.h:22