Page MenuHomec4science
Files F91977630

stimulation_dislo.cc
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 16, 08:03

stimulation_dislo.cc
View Options

/**
* @file stimulation_dislo.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief This stmulation command imposes the Volterra displacement fields of strainght screw edge or mixed dislocations
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "stimulation_dislo.hh"
#include "lib_md.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "filter.hh"
#include "ref_point_data.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
inline Real magnitude(int * vec, UInt nb) {
Real mag = 0.0;
for (UInt i = 0; i < nb; ++i) {
mag += Real(vec[i])* Real(vec[i]);
}
return sqrt(mag);
}
/* -------------------------------------------------------------------------- */
template <typename _Input>
void StimulationDislo<_Input>::init() {
// first step: normalize all the vectors to build a coordinate system
Real mag_n = magnitude(this->normal_quant, 3);
Real mag_l = magnitude(this->line_dir_quant, 3);
this->burg_edge_quant[0] = this->normal_quant[1]*this->line_dir_quant[2] - this->normal_quant[2]*this->line_dir_quant[1];
this->burg_edge_quant[1] = this->normal_quant[2]*this->line_dir_quant[0] - this->normal_quant[0]*this->line_dir_quant[2];
this->burg_edge_quant[2] = this->normal_quant[0]*this->line_dir_quant[1] - this->normal_quant[1]*this->line_dir_quant[0];
for (UInt i = 0; i < 3; ++i ) {
this->normal[i] = Real(this->normal_quant[i])/mag_n;
this->line_dir[i] = Real(this->line_dir_quant[i])/mag_l;
}
this->burg_edge[0] = this->normal[1]*this->line_dir[2] - this->normal[2]*this->line_dir[1];
this->burg_edge[1] = this->normal[2]*this->line_dir[0] - this->normal[0]*this->line_dir[2];
this->burg_edge[2] = this->normal[0]*this->line_dir[1] - this->normal[1]*this->line_dir[0];
if(this->dir_replica_char=='X' || this->dir_replica_char=='x'){
dir_replica = 0;
}
else if(this->dir_replica_char=='Y' || this->dir_replica_char=='y'){
dir_replica = 1;
}
else if(this->dir_replica_char=='Z' || this->dir_replica_char=='z'){
dir_replica = 2;
}
else{
LM_FATAL("replica direction can only be X, Y or Z");
}
}
/* -------------------------------------------------------------------------- */
template <typename _Input>
inline Real StimulationDislo<_Input>::computeBeDisplacement(Real b_edge, Real * X){
Real x_2 = X[0]*X[0];
Real y_2 = X[1]*X[1];
Real epsilon = 1e-15;
Real res = b_edge/(2.0*M_PI)*(atan2(X[1],X[0]) +
X[0]*X[1]/(2.0*(1-nu)*(x_2 + y_2+epsilon))
);
return res;
}
/* -------------------------------------------------------------------------- */
template <typename _Input>
inline Real StimulationDislo<_Input>::computeNDisplacement(Real b_edge, Real * X){
Real x_2 = X[0]*X[0];
Real y_2 = X[1]*X[1];
Real b_edge2 = b_edge*b_edge;
Real epsilon = 1e-15;
return -b_edge/(2.0*M_PI)*((1-2.0*nu)/(4.0*(1-nu))*(log((x_2 + y_2+epsilon)/b_edge2))
+ (x_2 - y_2)/
(4.0*(1-nu)*(x_2 + y_2+epsilon))
);
}
/* -------------------------------------------------------------------------- */
template <typename _Input>
inline Real StimulationDislo<_Input>::computeLDisplacement(Real b_screw, Real * X){
Real res = b_screw/(2.0*M_PI)*(atan2(X[1],X[0]));
return res;
}
/* -------------------------------------------------------------------------- */
template <typename _Input>
void StimulationDislo<_Input>::fix_boundary(typename _Input::Ref & atom,
const Quantity<Length,3> & xmin,
const Quantity<Length,3> & xmax) {
Real upper_diff = atom.position(this->dir_replica) - xmax[this->dir_replica];
Real lower_diff = atom.position(this->dir_replica) - xmin[this->dir_replica];
if (fabs(upper_diff) < this->fix_border_tol) {
atom.position(this->dir_replica) = xmax[this->dir_replica];
} else if (fabs(lower_diff) < this->fix_border_tol) {
atom.position(this->dir_replica) = xmin[this->dir_replica];
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void StimulationDislo<Cont>::stimulate(Cont & cont, UInt stage) {
if (this->fix_border_tol == 0.) {
this->stimulate_fix<false>(cont, stage);
} else {
this->stimulate_fix<true>(cont, stage);
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
template <bool do_the_fix>
void StimulationDislo<Cont>::stimulate_fix(Cont & cont, UInt stage)
{
#ifndef LM_OPTIMIZED
static const UInt Dim = Cont::Dim;
#endif// LM_OPTIMIZED
Cube & cube = cont.getBoundingBox();
const Quantity<Length,3> & xmax = cube.getXmax();
const Quantity<Length,3> & xmin = cube.getXmin();
Real * domainSize = cube.getSize();
Real b_edge;
Real b_screw;
b_edge = burgers*sin(theta*M_PI/180.0);
b_screw=-burgers*cos(theta*M_PI/180.0);
//b_edge = burgers*cos((theta-90.0)*M_PI/180.0);
//b_screw= burgers*sin((theta-90.0)*M_PI/180.0);
LM_ASSERT(Dim == 3,"this stimulator works only in 3D");
typedef typename Cont::Ref RefDof;
typename Cont::iterator it = cont.getIterator();
Cube c = cont.getBoundingBox();
double dipolelength = 0;
for (UInt i = 0; i < 3; i++) {
if(this->dipole[i] == 1) {
dipolelength = c.getSize(i);
}
}
for (RefDof at = it.getFirst() ; !it.end() ; at = it.getNext()) {
Real X0[3] = {0.,0.,0.}; // position in physical space
Real X[3] = {0.,0.,0.}; // position in physical space
Real x[3] = {0.,0.,0.}; // position in the coordinates of the dislo
Real Disp[3] = {0.,0.,0.};
Real disp[3];
at.getPositions0(X0);
for (int rep = -1*nb_replica ; rep < nb_replica+1 ; ++rep) {
X[0] = X0[0];
X[1] = X0[1];
X[2] = X0[2];
for (UInt i = 0; i < 3; i++) {
if(this->dipole[i] == 1) {
if(X[i]<0) {
X[i] = -1.0*(X[i] + (dipolelength/2.0)*0.5);
}
else {
X[i] = 1.0*(X[i] - (dipolelength/2.0)*0.5);
}
}
}
for (UInt i = 0; i < 3; ++i) {
X[i] -= pos[i];
disp[i] = 0.0;
}
// calculate coordinates of images
X[dir_replica] += domainSize[dir_replica]*Real(rep);
// rotate point
for (UInt i = 0; i < 3; ++i) {
x[i] = 0.0;
for (UInt j = 0; j < 3; ++j) {
x[i] += this->rotation_matrix[3*i+j]*X[j];
}
}
// compute displacement
disp[0] = computeBeDisplacement(b_edge, x);
if(rep==0) disp[1] = computeNDisplacement(b_edge, x);
disp[2] = computeLDisplacement(b_screw, x);
// compensate accumulated images in x-axis
if (rep == 0) {
if (x[1] > 0) {
disp[0] -= Real(nb_replica)*b_edge/2;
disp[2] -= Real(nb_replica)*b_screw/2;
} else {
disp[0] += Real(nb_replica)*b_edge/2;
disp[2] += Real(nb_replica)*b_screw/2;
}
}
// rotate back
for (UInt i = 0; i < 3; ++i) {
for (UInt j = 0; j < 3; ++j) {
Disp[i] += this->rotation_matrix[3*j+i] * disp[j];
}
}
}
at.position(0) += Disp[0];
at.position(1) += Disp[1];
at.position(2) += Disp[2];
if (do_the_fix){
this->fix_boundary(at, xmin, xmax);
}
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC DISLO
This stmulation command imposes the Volterra displacement\\
fields of straight screw, edge, or combination of two dislocations\\\\
Displacement field of Screw dislocation:
\begin{equation}
u_{z}(x, y) = \frac{b}{2\pi} \cdot atan2(y,x)
\end{equation}
Displacement fields of Edge dislocation:
\begin{equation}
u_{x}(x, y) = \frac{b}{2\pi} \cdot
\left(
atan2(y,x) + \frac{xy}{2 (1-\nu) (x^2 + y^2)} \right)
\end{equation}
\begin{equation}
u_{y}(x, y) = -\frac{b}{2\pi} \cdot
\left(
\frac{1-2\nu}{4(1 - \nu)}\ln(x^2 + y^2)
+ \frac{x^2 - y^2}{4(1-\nu)(x^2 + y^2)}
\right)
\end{equation}
are introduced.
*/
/* LMEXAMPLE STIMULATION myDislo DISLO INPUT md STAGE PRE_DUMP BURGERS burgers POS 0 0 0 REPLICA 10 0 POISSON 0.3468 THETA 180 ONESHOT 0 */
/* LMHERITANCE action_interface */
template <typename _Input>
void StimulationDislo<_Input>::declareParams(){
Stimulation<_Input>::declareParams();
/* LMKEYWORD BURGERS
Specify the magnitude of Burger's vector
*/
this->parseKeyword("BURGERS", this->burgers);
/* LMKEYWORD POS
Specify the coordinates of a point on the dislocation line
*/
this->parseVectorKeyword("POS", 3, this->pos);
/* LMKEYWORD SLIP_PLANE
Specifies the normal vector to the slip plane
*/
this->parseVectorKeyword("SLIP_PLANE", 3, this->normal_quant);
/* LMKEYWORD LINE_DIR
Specifies the normal vector to the slip plane
*/
this->parseVectorKeyword("LINE_DIR", 3, this->line_dir_quant);
/* LMKEYWORD POISSON
Specify Poisson's ratio of the material
*/
this->parseKeyword("POISSON", this->nu);
// Resolution and width are temporarily taken out, because they don't do anything,
// plus, it's not clear how to spread a screw dislo, since it doesn't have a slip plane
///* LM KEYW ORD RESOLUTION
// Specify the number of increments in the spread of dislocation density in Z direction
//*/
//parseKey word("RESOLUTION",resolution);
//
///* LM KEY esWORD WIDTH
// Specify the width of the core dislocation in Z direction
//*/
//parseKe yword("WIDTH",width);
/* LMKEYWORD NB_REPLICA
Specify the number of replicas in burgers edge components
*/
this->parseKeyword("NB_REPLICA", this->nb_replica);
/* LMKEYWORD DIR_REPLICA
Specify the direction of replication (X, Y or Z)
*/
this->parseKeyword("DIR_REPLICA", this->dir_replica_char,'X');
// temporarily taken out, because we can probably get along without it
///* LMKEY WORD SIGN
// Specify the direction of dislocation by changing signs of x and y coordinates,//
// For example,
// 1 1 means dislocation with Z and Y coordinates,
// -1 1 means dislocation with -Z and Y coordinates,
// 1 -1 means dislocation with Z and -Y coordinates, and
// -1 -1 means dislocation with -Z and -Y coordinates
//*/
//parseVe ctorKeyword("SIGN",2,sign);
/* LMKEYWORD THETA
Specify the angle of burgers vector with dislocation line in unit of degree\\
For example, 90 and 270 impose pure edge burgers vector,\\
0 and 180 impose pure screw burgers vector, and\\
other angles impose combination of two kinds of burgers vector
*/
this->parseKeyword("THETA", this->theta);
/* LMKEYWORD FIX_BORDER_TOL
specify a tolerance whithin which border points are forces to be exacly
on the border:
------------- -------------
| | | |
| | | |
/ | => | |
\ | => | |
| | | |
| | | |
------------- -------------
*/
this->parseKeyword("FIX_BORDER_TOL", this->fix_border_tol,0.);
/* LMKEYWORD DIPOLE
Specify the direction where the dislocation dipole creation in X, Y or Z.
As results, two dislocation are inserted with full periodicity in all directions.
%% For example, DIPOLE 0 1 0 creates:
%% ------------- -------------
%% | | | |
%% | | | T |
%% | _|_ | => | |
%% | | | _|_ |
%% | | | |
%% ------------- -------------
%% where T and _|_ are same burger vectors with opposite signs.
*/
this->parseVectorKeyword("DIPOLE", 3, this->dipole, VEC_DEFAULTS(0,0,0));
//this->parseKeyword("DIPOLE", this->dipole);
}
/* -------------------------------------------------------------------------- */
DECLARE_STIMULATION(StimulationDislo,LIST_ATOM_MODEL)
DECLARE_STIMULATION(StimulationDislo,LIST_CONTINUUM_MODEL)
DECLARE_STIMULATION(StimulationDislo,LIST_DD_MODEL)
DECLARE_STIMULATION_REFPOINT(StimulationDislo)
DECLARE_STIMULATION_GENERIC_MESH(StimulationDislo)
__END_LIBMULTISCALE__

Event Timeline

c4science · Help