Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90029229
stimulation_dislo.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Oct 28, 15:28
Size
23 KB
Mime Type
text/x-c++
Expires
Wed, Oct 30, 15:28 (2 d)
Engine
blob
Format
Raw Data
Handle
21996263
Attached To
rLIBMULTISCALE LibMultiScale
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 "stimulation_dislo.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.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
);
}
inline
Real
round
(
Real
value
)
{
Real
int_value
=
(
Real
)
floor
(
value
);
if
(
value
-
int_value
>=
0.5
)
{
return
(
Real
)
ceil
(
value
);
}
else
{
return
(
Real
)
floor
(
value
);
}
}
/* -------------------------------------------------------------------------- */
StimulationDislo
::
StimulationDislo
(
const
std
::
string
&
name
)
:
LMObject
(
name
)
{
nb_replica
=
0
;
this
->
dir_replica_char
=
'X'
;
this
->
burg_edge
=
this
->
rotation_matrix
;
this
->
normal
=
this
->
burg_edge
+
3
;
this
->
line_dir
=
this
->
burg_edge
+
6
;
for
(
UInt
i
=
0
;
i
<
3
;
++
i
)
{
this
->
pos
[
i
]
=
0.0
;
this
->
burg_edge_quant
[
i
]
=
0
;
this
->
normal_quant
[
i
]
=
0
;
this
->
line_dir_quant
[
i
]
=
0
;
this
->
dipole
[
i
]
=
0
;
}
this
->
normal_quant
[
1
]
=
1
;
this
->
line_dir_quant
[
2
]
=
1
;
this
->
theta
=
90
;
this
->
fix_border_tol
=
0
;
this
->
linear
=
false
;
this
->
aniso
=
false
;
this
->
iso
=
false
;
for
(
UInt
i
=
0
;
i
<
13
;
++
i
)
this
->
elastic_coefs
[
i
]
=
0.0
;
this
->
c11bar
=
0.0
;
this
->
lambda
=
0.0
;
this
->
phi
=
0.0
;
};
/* -------------------------------------------------------------------------- */
StimulationDislo
::~
StimulationDislo
()
{}
/* -------------------------------------------------------------------------- */
void
StimulationDislo
::
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"
);
}
// security check regarding anisotropy
if
((
!
this
->
aniso
)
&&
(
!
this
->
iso
))
{
for
(
UInt
i
=
0
;
i
<
13
;
++
i
)
{
if
(
this
->
elastic_coefs
[
i
]
!=
0.0
)
{
LM_FATAL
(
"In case of (an)isotropy, elastic coeffcients shouldn't be "
"defined."
);
}
}
}
else
if
((
this
->
aniso
)
&&
(
this
->
iso
))
{
LM_FATAL
(
"Both anisotropy and isotropy cannot be defined"
);
}
else
{
Real
sum_elastic_coefs
=
0.0
;
for
(
UInt
i
=
0
;
i
<
13
;
++
i
)
{
sum_elastic_coefs
+=
this
->
elastic_coefs
[
i
];
}
if
(
sum_elastic_coefs
==
0.0
)
{
LM_FATAL
(
"In case of (an)isotropy, elastic coeffcients should be defined."
);
}
if
(
this
->
linear
)
LM_FATAL
(
"In case of (an)isotropy, linear field cannot be applied"
);
Real
c11
=
this
->
elastic_coefs
[
0
];
Real
c12
=
this
->
elastic_coefs
[
1
];
Real
c22
=
this
->
elastic_coefs
[
4
];
Real
c66
=
this
->
elastic_coefs
[
12
];
this
->
c11bar
=
sqrt
(
c11
*
c22
);
if
(
this
->
aniso
)
{
Real
lambda2
=
sqrt
(
c11
/
c22
);
this
->
lambda
=
sqrt
(
lambda2
);
this
->
phi
=
0.5
*
acos
((
c12
*
c12
+
2.0
*
c12
*
c66
-
c11bar
*
c11bar
)
/
(
2.0
*
c11bar
*
c66
));
}
else
{
this
->
lambda
=
1.0
;
this
->
phi
=
M_PI
/
2.0
;
}
}
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
eval_q_coef
(
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
inside
=
X
[
0
]
*
X
[
0
]
+
2.0
*
X
[
0
]
*
X
[
1
]
*
this
->
lambda
*
cos
(
this
->
phi
)
+
X
[
1
]
*
X
[
1
]
*
this
->
lambda
*
this
->
lambda
;
return
sqrt
(
inside
);
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
eval_t_coef
(
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
inside
=
X
[
0
]
*
X
[
0
]
-
2.0
*
X
[
0
]
*
X
[
1
]
*
this
->
lambda
*
cos
(
this
->
phi
)
+
X
[
1
]
*
X
[
1
]
*
this
->
lambda
*
this
->
lambda
;
return
sqrt
(
inside
);
}
/* -------------------------------------------------------------------------- */
inline
Real
StimulationDislo
::
evalArctan
(
Real
X
,
Real
Y
)
{
if
((
X
>=
0.0
)
&&
(
Y
>
0.0
))
{
return
atan
(
Y
/
X
);
}
else
if
((
X
<
0.0
)
&&
(
Y
>=
0.0
))
{
return
atan
(
-
Y
/
-
X
)
+
M_PI
;
}
else
if
((
X
<
0.0
)
&&
(
Y
<
0.0
))
{
return
atan
(
-
Y
/
-
X
)
-
M_PI
;
}
else
if
((
X
>=
0.0
)
&&
(
Y
<
0.0
))
{
return
atan
(
Y
/
X
);
}
else
{
DUMP
(
"X:"
<<
X
<<
", Y: "
<<
Y
,
DBG_MESSAGE
);
DUMP
(
"c11bar:"
<<
this
->
c11bar
<<
", lambda: "
<<
lambda
<<
", phi: "
<<
phi
,
DBG_MESSAGE
);
LM_FATAL
(
"unexpected case"
);
}
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
computeBeDisplacement
(
Real
b_edge
,
Vector
<
_Input
::
Dim
>
&
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
::
computeBeDisplacementAniso
(
Real
b_edge
,
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
first
=
0.0
;
first
+=
evalArctan
((
X
[
0
]
+
X
[
1
]
*
this
->
lambda
*
cos
(
this
->
phi
)),
(
X
[
1
]
*
this
->
lambda
*
sin
(
this
->
phi
)));
first
+=
evalArctan
((
X
[
0
]
-
X
[
1
]
*
this
->
lambda
*
cos
(
this
->
phi
)),
(
X
[
1
]
*
this
->
lambda
*
sin
(
this
->
phi
)));
Real
q
=
eval_q_coef
(
X
);
Real
t
=
eval_t_coef
(
X
);
Real
c12
=
this
->
elastic_coefs
[
1
];
Real
c66
=
this
->
elastic_coefs
[
12
];
Real
logvalue
=
log
(
q
/
t
);
Real
second
=
(
this
->
c11bar
*
this
->
c11bar
-
c12
*
c12
)
/
(
2.0
*
this
->
c11bar
*
c66
*
sin
(
2.0
*
this
->
phi
))
*
logvalue
;
// Real second = 0.0;
Real
res
=
-
1.0
*
b_edge
/
(
4.0
*
M_PI
)
*
(
-
1.0
*
first
+
second
);
return
res
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
computeBeDisplacementIso
(
Real
b_edge
,
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
first
=
0.0
;
Real
second
=
0.0
;
Real
res
=
0.0
;
Real
c12
=
this
->
elastic_coefs
[
1
];
Real
c66
=
this
->
elastic_coefs
[
12
];
first
+=
evalArctan
(
X
[
0
],
X
[
1
]);
first
+=
evalArctan
(
X
[
0
],
X
[
1
]);
second
=
(
this
->
c11bar
*
this
->
c11bar
-
c12
*
c12
)
/
(
4.0
*
this
->
c11bar
*
c66
)
*
((
2.0
*
X
[
0
]
*
X
[
1
])
/
(
X
[
0
]
*
X
[
0
]
+
X
[
1
]
*
X
[
1
]));
res
=
b_edge
/
(
4.0
*
M_PI
)
*
(
first
-
second
);
return
res
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
computeNDisplacement
(
Real
b_edge
,
Vector
<
_Input
::
Dim
>
&
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))
// );
return
-
b_edge
/
(
2.0
*
M_PI
)
*
((
1
-
2.0
*
nu
)
/
(
4.0
*
(
1
-
nu
))
*
log
(
x_2
+
y_2
+
epsilon
)
+
(
x_2
-
y_2
)
/
(
4.0
*
(
1
-
nu
)
*
(
x_2
+
y_2
+
epsilon
)));
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
computeNDisplacementAniso
(
Real
b_edge
,
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
q
=
eval_q_coef
(
X
);
Real
t
=
eval_t_coef
(
X
);
Real
c12
=
this
->
elastic_coefs
[
1
];
Real
first
=
(
this
->
c11bar
-
c12
)
*
cos
(
this
->
phi
)
*
log
(
q
*
t
)
/
sin
(
2.0
*
this
->
phi
);
Real
second
=
(
this
->
c11bar
+
c12
)
*
sin
(
this
->
phi
)
*
atan2
(
(
X
[
1
]
*
X
[
1
]
*
this
->
lambda
*
this
->
lambda
*
sin
(
2.0
*
this
->
phi
)),
(
X
[
0
]
*
X
[
0
]
-
this
->
lambda
*
this
->
lambda
*
X
[
1
]
*
X
[
1
]
*
cos
(
2.0
*
this
->
phi
)))
/
sin
(
2.0
*
this
->
phi
);
Real
res
=
this
->
lambda
*
b_edge
/
4.0
/
M_PI
/
this
->
c11bar
*
(
first
+
second
);
return
res
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
computeNDisplacementIso
(
Real
b_edge
,
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
first
=
0.0
;
Real
second
=
0.0
;
Real
res
=
0.0
;
Real
c12
=
this
->
elastic_coefs
[
1
];
first
=
(
this
->
c11bar
-
c12
)
/
2.
*
log
(
X
[
0
]
*
X
[
0
]
+
X
[
1
]
*
X
[
1
]);
second
=
(
this
->
c11bar
+
c12
)
*
((
X
[
1
]
*
X
[
1
])
/
(
X
[
0
]
*
X
[
0
]
+
X
[
1
]
*
X
[
1
]));
res
=
b_edge
/
4.0
/
M_PI
*
(
first
+
second
);
return
res
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
computeLDisplacement
(
Real
b_screw
,
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
res
=
b_screw
/
(
2.0
*
M_PI
)
*
(
atan2
(
X
[
1
],
X
[
0
]));
return
res
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
inline
Real
StimulationDislo
::
computeLDisplacementAniso
(
Real
b_screw
,
Vector
<
_Input
::
Dim
>
&
X
)
{
Real
c44
=
this
->
elastic_coefs
[
9
];
Real
c45
=
this
->
elastic_coefs
[
10
];
Real
c55
=
this
->
elastic_coefs
[
11
];
Real
res
=
b_screw
/
2.0
/
M_PI
*
atan2
(
sqrt
(
c44
*
c55
-
c45
*
c45
)
*
X
[
1
],
c44
*
X
[
0
]
-
c45
*
X
[
1
]);
return
res
;
}
/* -------------------------------------------------------------------------- */
inline
Real
StimulationDislo
::
computeLinearDisplacement
(
Real
burg
,
Real
*
X
,
Real
xlo
,
Real
xhi
)
{
Real
res
=
burg
/
2.
/
(
xhi
-
xlo
)
*
(
X
[
0
]
-
xhi
);
if
(
X
[
1
]
>
0
)
return
-
1.0
*
res
;
else
return
1.0
*
res
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
_Input
>
void
StimulationDislo
::
fix_boundary
(
typename
_Input
::
Ref
&
atom
,
const
Vector
<
_Input
::
Dim
>
&
xmin
,
const
Vector
<
_Input
::
Dim
>
&
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
)
{
Vector
<
_Input
::
Dim
>
pos
=
atom
.
position
();
pos
[
this
->
dir_replica
]
=
xmax
[
this
->
dir_replica
];
atom
.
position
()
=
pos
;
}
else
if
(
fabs
(
lower_diff
)
<
this
->
fix_border_tol
)
{
Vector
<
_Input
::
Dim
>
pos
=
atom
.
position
();
pos
[
this
->
dir_replica
]
=
xmin
[
this
->
dir_replica
];
atom
.
position
()
=
pos
;
}
}
/* -------------------------------------------------------------------------- */
template
<
typename
Cont
>
void
StimulationDislo
::
stimulate
(
Cont
&
cont
)
{
this
->
stimulate_to_clean
(
cont
,
current_stage
);
LM_TOIMPLEMENT
;
// if (this->fix_border_tol != 0.) {
// for (auto &&at : cont) {
// this->fix_boundary(at, xmin, xmax);
// }
// }
}
/* -------------------------------------------------------------------------- */
template
<
typename
Cont
>
void
StimulationDislo
::
stimulate_to_clean
(
Cont
&
cont
,
UInt
stage
)
{
constexpr
UInt
Dim
=
Cont
::
Dim
;
Cube
&
cube
=
cont
.
getBoundingBox
();
__attribute__
((
unused
))
auto
&&
xmax
=
cube
.
getXmax
<
Cont
::
Dim
>
();
__attribute__
((
unused
))
auto
&&
xmin
=
cube
.
getXmin
<
Cont
::
Dim
>
();
auto
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
);
// identify the direction of propagation (Be)
Real
BeDir0
=
normal_quant
[
1
]
*
line_dir_quant
[
2
]
-
normal_quant
[
2
]
*
line_dir_quant
[
1
];
Real
BeDir1
=
-
(
normal_quant
[
0
]
*
line_dir_quant
[
2
]
-
normal_quant
[
2
]
*
line_dir_quant
[
0
]);
Real
BeDir2
=
normal_quant
[
0
]
*
line_dir_quant
[
1
]
-
normal_quant
[
1
]
*
line_dir_quant
[
0
];
Real
BeMag
=
sqrt
(
BeDir0
*
BeDir0
+
BeDir1
*
BeDir1
+
BeDir2
*
BeDir2
);
BeDir0
/=
BeMag
;
BeDir1
/=
BeMag
;
BeDir2
/=
BeMag
;
// UInt bedir = 4;
// if (fabs(abs(BeDir0) - 1.0) < 1e-4)
// bedir = 0;
// else if (fabs(abs(BeDir1) - 1.0) < 1e-4)
// bedir = 1;
// else if (fabs(abs(BeDir2) - 1.0) < 1e-4)
// bedir = 2;
// else
// LM_FATAL("cannot identify propagation direction");
throw
;
// this was temporarily commented
// Real xloBe = xmin[bedir];
// Real xhiBe = xmax[bedir];
LM_ASSERT
(
Dim
==
3
,
"this stimulator works only in 3D"
);
Cube
c
=
cont
.
getBoundingBox
();
double
dipolelength
=
0
;
for
(
UInt
i
=
0
;
i
<
3
;
i
++
)
{
if
(
this
->
dipole
[
i
]
==
1
)
{
dipolelength
=
c
.
getSize
(
i
);
}
}
for
(
auto
&&
at
:
cont
)
{
auto
X0
=
at
.
position0
();
// position in physical space
auto
X
=
at
.
position
();
// position in physical space
Vector
<
Dim
>
x
=
Vector
<
Dim
>::
Zero
();
// position in the coordinates of the dislo
Vector
<
Dim
>
Disp
=
Vector
<
Dim
>::
Zero
();
Real
disp
[
3
];
for
(
int
rep
=
-
1
*
nb_replica
;
rep
<
nb_replica
+
1
;
++
rep
)
{
X
=
X0
;
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
if
(
this
->
linear
==
false
)
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
];
}
}
if
(
this
->
linear
==
false
)
{
if
(
this
->
aniso
)
{
// compute anisotropy volterra displacement
LM_TOIMPLEMENT
;
// disp[0] = computeBeDisplacementAniso(b_edge, x);
if
(
rep
==
0
)
{
LM_TOIMPLEMENT
;
// disp[1] = computeNDisplacementAniso(b_edge, x);
}
LM_TOIMPLEMENT
;
// disp[2] = computeLDisplacementAniso(b_screw, x);
}
else
if
(
this
->
iso
)
{
// compute isotropy volterra displacement
LM_TOIMPLEMENT
;
// disp[0] = computeBeDisplacementIso(b_edge, x);
if
(
rep
==
0
)
{
LM_TOIMPLEMENT
;
// disp[1] = computeNDisplacementIso(b_edge, x);
}
LM_TOIMPLEMENT
;
// disp[2] = computeLDisplacementAniso(b_screw, x);
}
else
{
// compute isotropy volterra displacement
LM_TOIMPLEMENT
;
// disp[0] = computeBeDisplacement(b_edge, x);
if
(
rep
==
0
)
{
LM_TOIMPLEMENT
;
// disp[1] = computeNDisplacement(b_edge, x);
}
LM_TOIMPLEMENT
;
// disp[2] = computeLDisplacement(b_screw, x);
}
}
else
{
// compute linear displacements
throw
;
// disp[0] = computeLinearDisplacement(b_edge, x, xloBe, xhiBe);
// disp[2] = computeLinearDisplacement(b_screw, x, xloBe, xhiBe);
}
// compensate accumulated images in x-axis
if
((
rep
==
0
)
and
(
this
->
linear
==
false
))
{
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
];
}
}
if
(
this
->
linear
)
break
;
// exit replica for loops
}
at
.
position
()
+=
Disp
;
if
((
fabs
(
at
.
position0
()[
0
]
+
26.3245
)
<
1e-1
)
&&
(
fabs
(
at
.
position0
()[
1
]
-
102.843
)
<
1e-1
)
&&
(
fabs
(
at
.
position0
()[
2
]
+
22.8263
)
<
1e-1
))
{
DUMP
(
"HERE"
,
DBG_MESSAGE
);
}
}
}
/* --------------------------------------------------------------------------
*/
/* 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 */
void
StimulationDislo
::
declareParams
()
{
StimulationInterface
::
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'
);
/* 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
));
/* LMKEYWORD LINEAR
Specify the keyword to employ the linear displacement fields instead of
volterra solution
*/
this
->
parseTag
(
"LINEAR"
,
this
->
linear
,
false
);
/* LMKEYWORD ANISO
Boolean to set anisotropicity:
[1] Elastic Strain Fields and Dislocation Mobility, Volume 31, J.Lothe
Dislocations in anisotropic media
[2] Plane-strain discrete dislocation plasticity incorporating
anisotropic
elasticity IJSS 48(2), January 2011
*/
this
->
parseTag
(
"ANISO"
,
this
->
aniso
,
false
);
/* L2MKEYWORD ELASTIC_COEFS
Specify the thirteen elastic coefficients
of the constitutive matrix of anisotropic media in the following order:
C11, C12, C13, C21, C22, C23, C31, C32, C33, C44, C45, C55, C66
*/
// this->parseVectorKeyword("ELASTIC_COEFS", 13, this->elastic_coefs);
/* LMKEYWORD ISO
Boolean to set isotropicity :
[1] Elastic Strain Fields and Dislocation Mobility, Volume 31, J.Lothe
Dislocations in anisotropic media
[2] Plane-strain discrete dislocation plasticity incorporating
anisotropic
elasticity IJSS 48(2), January 2011
*/
this
->
parseTag
(
"ISO"
,
this
->
iso
,
false
);
}
/* --------------------------------------------------------------------------
*/
DECLARE_STIMULATION_MAKE_CALL
(
StimulationDislo
);
__END_LIBMULTISCALE__
Event Timeline
Log In to Comment