testMls.cpp

Demonstrates the use of Moving Least Squares (Mls) to compute deformation gradients and push stress tensors to a spatial grid.

This test sets up a simulation involving a particle configuration and grid data, computes the Piola-Kirchhoff stress using the LDAD method, and then applies the Moving Least Squares (MLS) method to compute the deformation gradient and push the Piola stress to Cauchy stress on a grid.

Key steps in this example:

• Load atomic configuration and potential.
• Load reference grid data.
• Compute LDAD Piola stress using MethodLdadConstant.
• Construct an Mls object and compute deformation gradients.
• Push stress field to obtain Cauchy stress tensors.
• Write deformation gradient, pushed grid, and Cauchy stress to output files.

Output files:

• mls.DeformationGradient
• mlsPushed.grid
• mlsPushed.stress

Requirements:

• Input file config.data with atomic data in BoxConfiguration format.
• Grid files reference.grid.
• A valid KIM model (here: "LJ_Smoothed_Bernardes_1958_Ar__MO_764178710049_001").
• Expected Cauchy stress field mlsPushedReference.stress for unit testing.

1. Open the configuration file (*.data in MDStressLab format) and read the number of particles

       int referenceAndFinal= true;
       std::string configFileName= "config.data";
       std::ifstream file(configFileName);
       if(!file) MY_ERROR("ERROR: config.data could not be opened for reading!");
       int numberOfParticles;
       file >> numberOfParticles;
       if (numberOfParticles < 0) MY_ERROR("Error: Negative number of particles.\n");

2. Initialize the BoxConfiguration and read the reference and current atomic coordinates from the configuration file. The reference configuration is an fcc Ar crystal in the relaxed state. The deformed/current configuration is the reference crystal strained in the \(y\)-direction.

       BoxConfiguration body{numberOfParticles,referenceAndFinal};
       body.read(configFileName,referenceAndFinal);

3. Initialize the Kim model

       std::string modelname= "LJ_Smoothed_Bernardes_1958_Ar__MO_764178710049_001";
       Kim kim(modelname);

4. Load reference grid coordinates from files.

       int ngrid = 125;
       Grid<Reference> gridFromFile_ref(ngrid);
       gridFromFile_ref.read("reference.grid");

5. Initialize LDAD vectors and construct LDAD method MethodLdadConstant. The LDAD vectors (columns) are the three lattice vectors of the conventional unit cell.

       Matrix3d ldadVectors;
       ldadVectors << 5.29216036151419, 0.0, 0.0, 
                      0.0, 5.29216036151419, 0.0,
                      0.0, 0.0, 5.29216036151419;
       MethodLdadConstant ldad(ldadVectors);

6. Construct the Piola stress using the above method

       //TODO The bond function should be accepted as a reference
       Stress<MethodLdadConstant,Piola> ldad_constant_stress_ref("ldad",ldad,&gridFromFile_ref);

7. Construct an Mls object to compute deformation gradients and push the stress to Cauchy form.

       double MlsRadius = 15.87648108454257;
       Mls testMls(body,&gridFromFile_ref,MlsRadius,"mls");

8. Write deformation gradients, pushed grid points, and Cauchy stresses to output files.

       std::vector<Matrix3d> cauchyPushedField;
       testMls.pushToCauchy(ldad_constant_stress_ref.field,cauchyPushedField);
       testMls.writeDeformationGradient();
       testMls.writeGridPushed();
       testMls.writePushedCauchyStress(cauchyPushedField);

9. Compare the pushed fields to expected data.

       compareStress("mlsPushed");

Full code:

/*
 * main.cpp
 *
 *  Created on: Nov 3, 2019
 *      Author: Nikhil Admal
 */
 
#include <string>
#include <iostream>
#include <tuple>
#include <fstream>
#include "MethodLdad.h"
#include "BoxConfiguration.h"
#include "calculateStress.h"
#include "Mls.h"
#include "Grid.h"
#include "typedef.h"
#include "../compareStress.cpp"
 
 
 
int main()
{
    int referenceAndFinal= true;
    std::string configFileName= "config.data";
    std::ifstream file(configFileName);
    if(!file) MY_ERROR("ERROR: config.data could not be opened for reading!");
    int numberOfParticles;
    file >> numberOfParticles;
    if (numberOfParticles < 0) MY_ERROR("Error: Negative number of particles.\n");
    BoxConfiguration body{numberOfParticles,referenceAndFinal};
    body.read(configFileName,referenceAndFinal);
    std::string modelname= "LJ_Smoothed_Bernardes_1958_Ar__MO_764178710049_001";
    Kim kim(modelname);
    int ngrid = 125;
    Grid<Reference> gridFromFile_ref(ngrid);
    gridFromFile_ref.read("reference.grid");
    Matrix3d ldadVectors;
    ldadVectors << 5.29216036151419, 0.0, 0.0, 
                   0.0, 5.29216036151419, 0.0,
                   0.0, 0.0, 5.29216036151419;
    MethodLdadConstant ldad(ldadVectors);
    //TODO The bond function should be accepted as a reference
    Stress<MethodLdadConstant,Piola> ldad_constant_stress_ref("ldad",ldad,&gridFromFile_ref);
    calculateStress(body,kim,std::tie(ldad_constant_stress_ref));
    double MlsRadius = 15.87648108454257;
    Mls testMls(body,&gridFromFile_ref,MlsRadius,"mls");
    std::vector<Matrix3d> cauchyPushedField;
    testMls.pushToCauchy(ldad_constant_stress_ref.field,cauchyPushedField);
    testMls.writeDeformationGradient();
    testMls.writeGridPushed();
    testMls.writePushedCauchyStress(cauchyPushedField);
    compareStress("mlsPushed");
    return 0;
}