Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69632011
template_cadd.hh
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
Tue, Jul 2, 18:28
Size
30 KB
Mime Type
text/x-c++
Expires
Thu, Jul 4, 18:28 (2 d)
Engine
blob
Format
Raw Data
Handle
18734233
Attached To
rLIBMULTISCALE LibMultiScale
template_cadd.hh
View Options
/**
* @file template_cadd.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Wed Jan 15 17:00:43 2014
*
* @brief Simple CADD-like coupler
*
* @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/>.
*
*/
#ifndef __LIBMULTISCALE_TEMPLATE_CADD_HH__
#define __LIBMULTISCALE_TEMPLATE_CADD_HH__
#include "bridging.hh"
#include "compute_compatibility.hh"
#include "coupling_dd_md.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_parsable.hh"
#include "lm_parser.hh"
#include <string>
#include "container_mesh.hh"
#include "ref_point_data.hh"
//#include "twod_template.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
// struct JaehyunDDNeighbor {
// Real idx; // 0
// Vector<3> X; // 1,2,3
// Vector<3> burg; // 4,5,6
// Vector<3> normal; // 7,8,9
// };
// struct JaehyunDDPoint {
// Real idx; // 0
// Real num_neigh; // 1
// Vector<3> X; // 2,3,4
// Real idx0; // 5
// std::vector<JaehyunDDNeighbor> neighs; // 6
// Vector<3> tmp_X; // 27,28,29
// Vector<3> tmp2_X; // 30,31,32
// };
// inline std::ostream &operator<<(std::ostream &os, JaehyunDDNeighbor &ref) {
// os << "\t\tnbr_idx:" << ref.idx << std::endl;
// os << "\t\tnbr_x:" << ref.X << std::endl;
// os << "\t\tbx:" << ref.burg << std::endl;
// os << "\t\tnx:" << ref.normal << std::endl;
// return os;
// }
// inline std::ostream &operator<<(std::ostream &os, JaehyunDDPoint &ref) {
// os << "JaehyunDD point:" << std::endl;
// os << "\tidx:" << ref.idx << std::endl;
// os << "\tnum_neigh:" << ref.num_neigh << std::endl;
// os << "\tcrr_x:" << ref.X << std::endl;
// os << "\tidx0:" << ref.idx0 << std::endl;
// for (UInt i = 0; i < ref.num_neigh; ++i) {
// os << ref.neighs[i];
// }
// return os;
// }
/* -------------------------------------------------------------------------- */
class
CADD_template
:
public
CouplingDDMD
{
public
:
DECLARE_COUPLER
(
CADD_template
,
domDD
,
domMD
);
template
<
typename
DomainDD
,
typename
DomainA
>
void
coupling
(
CouplingStage
stage
,
DomainDD
&
domDD
,
DomainA
&
domA
);
template
<
typename
DomainDD
,
typename
DomainA
>
void
init
(
DomainDD
&
domDD
,
DomainA
&
domA
);
protected
:
template
<
typename
DomainA
>
ContainerGenericDDMesh
<
3
>
sendDDSegments2MD
(
DomainA
&
domA
);
template
<
typename
RefAtom
>
void
functionifAniso
(
UInt
num_dislo
,
Vector
<
3
>
&
pos0_pbc
,
Vector
<
3
>
&
pos0_pbc0
,
std
::
vector
<
std
::
vector
<
Real
>>
&
dislocations
,
RefAtom
&
at
);
template
<
typename
RefAtom
>
void
functionifVolterra
(
UInt
num_dislo
,
Vector
<
3
>
&
pos0_pbc
,
Vector
<
3
>
&
pos0_pbc0
,
std
::
vector
<
std
::
vector
<
Real
>>
&
dislocations
,
RefAtom
&
at
);
template
<
typename
RefAtom
>
void
functionElse
(
ContainerGenericDDMesh
<
3
>
&
dd_points
,
Vector
<
3
>
&
pos0_pbc
,
RefAtom
&
at
);
inline
void
setGeometries
(
int
pad_id
,
int
atom_id
,
int
geom_id
,
int
boundary_id
);
inline
virtual
void
clearAll
();
inline
Real
compute_dd_theta
(
Real
x3
[],
Real
x4
[]);
template
<
typename
Vec
>
inline
bool
intersectTwoD
(
Vec
&
p1
,
Vec
&
p2
,
Vec
&
p3
,
Vec
&
p4
,
Vector
<
2
>
&
pts
,
UInt
slip_case
);
template
<
typename
Vec
>
inline
bool
intersectThreeD
(
Vec
&
p1
,
Vec
&
p2
,
Vec
&
p3
,
Vec
&
p4
,
Vec
&
pts
,
UInt
slip_case
);
inline
bool
intersectTwoD_infinite
(
Real
p1
[],
Real
p2
[],
Real
p3
[],
Real
p4
[],
Real
*
ptsx
,
Real
*
ptsz
);
template
<
typename
Vec
>
inline
bool
intersect
(
Vec
&
x1
,
Vec
&
x2
,
Vec
&
x3
,
Vec
&
x4
,
UInt
slip_case
);
inline
Real
ccw
(
Real
Ax
,
Real
Ay
,
Real
Bx
,
Real
By
,
Real
Cx
,
Real
Cy
);
inline
Real
determinant
(
Real
a
,
Real
b
,
Real
c
,
Real
d
,
Real
e
,
Real
f
,
Real
g
,
Real
h
,
Real
i
);
inline
Real
radian2degree
(
Real
rad
);
inline
Real
angle_on_plane
(
Real
pos
[],
UInt
axis
);
// inline void RotationMatrixTemplateToRealLattices(Real theta, Real theta0,
// std::vector<Real>
// &resultR);
// inline void MatrixMultiplication(std::vector<Real> &mt, Real vec[],
// Real result[]);
// inline void CrossProduct(Real vec1[], Real vec2[], Real result[]);
inline
Real
distance_between_line_and_point
(
Real
x0
,
Real
z0
,
Real
pts1
[],
Real
pts2
[]);
inline
void
intersect_of_point_on_line_3d
(
Real
p0
[],
Real
p1
[],
Real
p2
[],
Real
result
[]);
template
<
typename
Vec
>
inline
Real
distance_between_line_and_point3D
(
Vec
&
pos0
,
Vec
&
pts1
,
Vec
&
pts2
);
template
<
typename
Vec
>
inline
Real
AreaTriangle
(
Vec
&
p1
,
Vec
&
p2
,
Vec
&
p3
,
UInt
slip_case
);
template
<
typename
Vec
>
inline
bool
is_it_closed
(
Vec
cp
,
Vec
&
cp0
,
Real
&
ci
,
Real
&
ci0
);
inline
UInt
next_index
(
UInt
next_idx
,
UInt
num_dd_pts0
,
std
::
vector
<
UInt
>
indices
);
template
<
typename
Vec
>
inline
UInt
find_one_direction
(
Vec
&
nb1
,
Vec
&
nb2
,
Real
b1mag
,
Real
b2mag
,
UInt
idx1
,
UInt
idx2
,
UInt
*
next_i
);
inline
Real
DistanceBetweenTwoNodes
(
Real
x1
,
Real
y1
,
Real
z1
,
Real
x2
,
Real
y2
,
Real
z2
);
void
mergeDuplicateNodes
();
/**
Performs the following tasks on processors MD and DD
MD: 1. detect the MD dislocation segments (DXA or user-defined compute)
2. send to DD processor(s)
DD: 1. (at every step) fix the DD nodes located within MD domain
3. receive the detected dislocation segments from MD processors
4. move DD nodes to the closest detected segments
5. register the remaining detected nodes as new DD nodes
6. check topology - remesh
*/
//! Impose the position of detected segments to coupled DDD segments
template
<
typename
DomainDD
,
typename
DomainMD
>
void
applyDetection
(
DomainDD
&
domDD
,
DomainMD
&
domA
);
// first coupling step
template
<
typename
DomainDD
,
typename
DomainMD
>
void
applyTemplateToPad
(
DomainDD
&
domDD
,
DomainMD
&
domA
);
// second coupling step
//! lower level auxiliary coupling methods
void
detect_core
();
// void detect_core_gui();
ContainerGenericDDMesh
<
3
>
detect_cores_by_compute
();
template
<
typename
RefAtom
>
void
applyTemplateToPoint
(
RefAtom
&
at
,
const
Real
pos_core
[
3
]);
template
<
typename
DomainA
>
void
close_dislocation_loop
(
ContainerGenericDDMesh
<
3
>
&
dd_points
,
ContainerGenericDDMesh
<
3
>
&
dd_points_global
,
DomainA
&
domA
);
bool
coarse_dislocation_loop
(
ContainerGenericDDMesh
<
3
>
&
dd_points
,
ContainerGenericDDMesh
<
3
>
&
dd_points0
);
void
find_indices_for_neighbors
(
ContainerGenericDDMesh
<
3
>
&
dd_points
);
void
number_of_dislocation_loops
(
UInt
num_dd_pts
,
UInt
*
num_loops
,
ContainerGenericDDMesh
<
3
>
&
dd_points
);
// void set_loop_number(Real node[], UInt num_dd_pts, Real **dd_points, UInt
// num_loops, UInt *i_index);
UInt
startNewDDPoint
(
UInt
num_dd_pts
,
std
::
vector
<
UInt
>
indices
);
UInt
findIndex_nbrnode
(
Real
node
[],
UInt
num_dd_pts
,
Real
**
dd_points
);
// void position_of_atom(Real pos0[], int *position, Real *template_disp,
// Real *barnett_disp, UInt num_dd_pts, Real **dd_points, Real maxSeg,
// Real *f);
void
center_of_circumcircle
(
Real
pts1
[],
Real
pts2
[],
Real
pts3
[],
Real
*
cx
,
Real
*
cy
,
Real
*
cz
);
void
coordinate_for_2d_template
(
Real
distance
,
Real
distance_from_dd_core
[],
Real
global_coord
[],
Real
template_coords
[]);
void
compute_angle_ratio
(
Real
angle
,
UInt
*
two_indices
,
Real
*
two_angles
);
void
computeNextPositionBasedOnCompute
(
std
::
vector
<
Real
>
computed_x
,
std
::
vector
<
Real
>
computed_z
,
Real
x
,
Real
z
,
Real
mindist
,
Real
maxdist
,
Real
*
next_position
);
void
newNodesRegistration
(
std
::
vector
<
Real
>
NextDDXs
,
std
::
vector
<
Real
>
NextDDYs
,
std
::
vector
<
Real
>
NextDDZs
);
template
<
typename
RefDDNode
>
void
moveDDconstNode2MDcore
(
std
::
vector
<
Real
>
&
computed_x
,
std
::
vector
<
Real
>
&
computed_y
,
std
::
vector
<
Real
>
&
computed_z
,
RefDDNode
node
);
template
<
typename
RefDDNode
,
typename
IteratorDDNodes
>
void
letAsFreeCombinedNode
(
RefDDNode
node
,
IteratorDDNodes
nodes
);
void
findPointInDDLine
(
std
::
vector
<
Real
>
ddx
,
std
::
vector
<
Real
>
ddy
,
std
::
vector
<
Real
>
ddz
,
Real
pts0
[],
Real
pts1
[],
Real
pts2
[],
Real
pts3
[],
std
::
vector
<
Real
>
&
NextDDXs
,
std
::
vector
<
Real
>
&
NextDDYs
,
std
::
vector
<
Real
>
&
NextDDZs
);
int
evaluate_edge
(
Real
x
,
Real
y
,
Real
z
,
Real
pts1
[],
Real
pts2
[],
Real
pts3
[],
Real
pts4
[],
Real
cent
[]);
bool
in_side_ellipse
(
Real
pos
[],
Real
rmax
,
Real
rmin
,
Real
*
factor
);
bool
is_it_infinite_connection
(
UInt
index0
,
UInt
*
index2
);
template
<
typename
RefAtom
>
void
applyDisplacement
(
RefAtom
at
,
Vector
<
3
>
&
total_disp
);
template
<
typename
DomainA
>
void
SetPadAtoms
(
DomainA
&
domA
);
void
SetP0Atoms
();
void
barnett_disp_atom_all_segments
(
Real
pos0
[],
Real
barnett_disp
[],
Real
**
dd_points
,
UInt
num_dd_pts
,
Real
burgMag
,
Real
pois
,
UInt
loop_number
);
void
compute_closure_point
(
Real
**
dd_points
,
UInt
num_dd_pts
,
UInt
loop_number
,
Real
A
[],
Real
B
[],
Real
b
[],
Real
C
[]);
template
<
typename
IteratorDDNodes
>
void
compute_pad_dd
(
UInt
*
Npad_dd
,
UInt
*
Npad_dd_column
,
IteratorDDNodes
nodes
,
Real
md_xrange
[],
Real
md_yrange
[],
Real
md_zrange
[]);
void
constructElems
(
std
::
vector
<
Real
>
nextDDXs
,
std
::
vector
<
Real
>
nextDDYs
,
std
::
vector
<
Real
>
nextDDZs
,
std
::
vector
<
RefGenericElem
<
3
>>
&
els
,
std
::
vector
<
Real
>
&
normal_ptr
,
std
::
vector
<
Real
>
&
burgers_ptr
);
template
<
typename
Vec
>
bool
PointInSquare
(
Vec
&
pts1
,
Vec
&
pts2
,
Vec
&
pts3
,
Vec
&
pts4
,
Vec
&
pts
,
UInt
slip_case
)
{
LM_TOIMPLEMENT
;
}
bool
is_it_pad_atom
(
Real
pos0
[],
Real
pos
[],
Geometry
*
pad
);
template
<
typename
Vec
>
bool
ComputeTemplateDisplacements
(
Vec
&
template_disp
,
Vec
&
current_node
,
Vec
&
next_node
,
Vec
&
burgers
,
Vec
&
normal
,
Vec
&
pos0
)
{
LM_TOIMPLEMENT
;
}
void
ComputeAnalyticDisplacements
(
Vector
<
3
>
&
barnett_disp
,
Vector
<
3
>
&
B
,
Vector
<
3
>
&
A
,
Vector
<
3
>
&
burg
,
Vector
<
3
>
&
pos0
);
void
ComputeBoundaryPointsInCoreRegion
(
Real
pts1
[],
Real
pts2
[],
Real
currnet_node
[],
Real
center
[]);
UInt
FindNeighborIndex
(
UInt
next_i
,
Real
**
dd_points
,
UInt
num_dd_pts
);
template
<
typename
DomainA
>
UInt
CollectingLocalDDpoints
(
ContainerGenericDDMesh
<
3
>
&
dd_points0_local
,
DomainA
&
domA
);
Real
*
CollectingDDpoints
();
template
<
typename
RefAtom
>
inline
void
shear_strain
(
RefAtom
at
,
Vector
<
3
>
&
disp
);
template
<
typename
RefAtom
>
inline
void
tilt_strain
(
RefAtom
at
,
Vector
<
3
>
&
disp
);
inline
void
ComputeOffsetParallelCurve
(
ContainerGenericDDMesh
<
3
>
&
dd_points
);
inline
void
computeOffset
(
Vector
<
3
>
&
current_node
,
Real
slope
,
Vector
<
3
>
&
pts1
,
Vector
<
3
>
&
pts2
,
UInt
slip_case
);
inline
void
FindOffsetPoints
(
Real
node_idx
,
Real
*
parallel_curve
,
UInt
num_dd_pts
,
Real
*
pts1
,
Real
*
pts2
);
inline
bool
near_md_boundary
(
Real
*
A
,
Real
*
B
,
Real
multiple
);
inline
void
findCorrectOffsets
(
Real
ang
,
Vector
<
3
>
&
pts1
,
Vector
<
3
>
&
pts2
,
UInt
slip_case
);
inline
void
CorrectIntersetingPoints
(
UInt
num_pts
,
Real
**
dd_points
);
//, Real *parallel_curve);
// inline void ComputeOffsetParallelCurve(UInt num_dd_pts, Real **dd_points,
// Real *parallel_curve);
// inline void computeOffset(Real *current_node, Real slope, Real *pts1, Real
// *pts2);
// inline void FindOffsetPoints(Real node_idx, Real *parallel_curve, UInt
// num_dd_pts, Real *normal);
inline
void
RelativePositionToPlane
(
Real
*
node
,
Real
*
normal
,
Real
*
pos
,
Real
*
position
);
// inline void findCorrectOffsets(Real ang, Real *pts1, Real *pts2);
// inline void CorrectIntersetingPoints(UInt num_pts, Real **dd_points, Real
// *parallel_curve);
inline
void
find_another_point_on_same_plane
(
Real
*
dd_node1
,
Real
*
burg
,
Real
*
normal1
,
Real
*
normal2
,
UInt
num_dd_pts
,
Real
**
dd_points
,
Real
*
dX1
,
Real
*
dX2
);
// inline Real BeDisplacement(Real b_edge, Real nu, Real X[]);
// inline Real NDisplacement(Real b_edge, Real nu, Real X[]);
// inline Real LDisplacement(Real b_screw, Real nu, Real X[]);
// inline void volterra_strain(Real pos[], Real disp[], Real pos0[]);
inline
Real
eval_q_coef
(
Real
*
X
);
inline
Real
eval_t_coef
(
Real
*
X
);
inline
Real
evalArctan
(
Real
X
,
Real
Y
);
inline
Real
computeBeDisplacementAniso
(
Real
b_edge
,
Real
*
X
);
inline
Real
computeNDisplacementAniso
(
Real
b_edge
,
Real
*
X
);
inline
Real
computeLDisplacementAniso
(
Real
b_screw
,
Real
*
X
);
inline
Real
computeBeDisplacement
(
Real
b_edge
,
Real
*
X
);
inline
Real
computeNDisplacement
(
Real
b_edge
,
Real
*
X
);
inline
Real
computeLDisplacement
(
Real
b_screw
,
Real
*
X
);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
//! detected core positions
std
::
vector
<
std
::
vector
<
Real
>>
md_cores
;
int
num_cores
;
int
detect_freq
;
Real
cluster_radius
;
//! The pad zone and DD nodes in atomistic domain
OutputContainer
pad_atoms
;
//! geometry ids for the pad zone and atomistic domain
LMID
pad_geom
;
// GeomID top_moving_geom;
// GeomID bot_moving_geom;
// GeomID slip_moving_geom;
// FilterID pad_filter;
int
normal_quant
[
3
];
Real
center
[
3
];
Real
angle_between_b_and_x
;
int
close_direction
[
3
];
int
num_close_loop_in_infinity
;
bool
close_in_infinity
;
bool
multiple_loops
;
bool
current_position_pad
;
bool
coarsen_dd
;
bool
use_specific_template
;
Real
specific_template_angle
;
Real
min_area
;
int
detect_simple_direction
[
3
];
bool
detect_simple_bool
;
UInt
simple_direction
;
LMID
md_geom
;
LMID
dd_geom
;
Vector
<
3
>
md_bbox_min
,
md_bbox_max
;
std
::
shared_ptr
<
Cube
>
md_cube
;
std
::
shared_ptr
<
Cube
>
dd_cube
;
// Real top_moving_xrange[2], top_moving_yrange[2], top_moving_zrange[2];
// Real bot_moving_xrange[2], bot_moving_yrange[2], bot_moving_zrange[2];
// Real slip_moving_xrange[2],slip_moving_yrange[2],slip_moving_zrange[2];
Real
md_box_lengths
[
3
];
Real
md_top_z
,
md_bot_z
;
LMID
detect_geom
;
bool
detect_dxa
;
bool
detect_partial
;
bool
dd_points_out
;
bool
dd_boundary
;
UInt
num_templates
;
std
::
vector
<
std
::
string
>
templatefilename_list
;
std
::
string
directory_dd_points
;
//! template related data
// std::vector<TwoDTemplate<RefAtom, IteratorAtomsSubset> *> core_templates;
// TwoDTemplate<RefAtom, IteratorAtomsSubset> core_template;
//! list of computes as added input
std
::
vector
<
LMID
>
compute_list
;
// std::vector<FilterID> truedisp;
// std::vector<FilterID> barnett_compute;
//! computed values filtered out
std
::
vector
<
ComputeCompatibility
*>
out_computes
;
Real
epsilon
[
2
];
Real
core_rcut
;
// core template correction radius
Real
rmin
;
// minor radius of the region where core template applied
Real
factor
;
Real
mu0
;
Real
shear_stresses
[
4
];
// # sigma_xy, sigma_yx, sigma_zy, sigma_yz
Real
tilt_stresses
[
3
];
// # xy, xz, yz
// UInt nb_slips;
// Real len_inter_slips;
// UInt node_index_list[8];
// UInt num_pairs_in_infinity;
LMID
domdd
;
// int volterra_replica;
// Real volterra_angle;
Real
distance_ddline_from_pad
;
Real
pad_thickness
;
// Real release_Xhi;
// Real release_Xlo;
// Real release_Zhi;
// Real release_Zlo;
Real
pois0
;
// Real maxSeg0;
// Real minSeg0;
Real
burgMag0
;
Quantity
<
Length
>
position_tolerance
;
// Real maxDDY0;
// Real minDDY0;
bool
aniso
;
bool
volterra
;
Real
elastic_coefs
[
13
];
Real
lambda
;
Real
phi
;
Real
c11bar
;
};
/* -------------------------------------------------------------------------- */
template
<>
inline
UInt
Parser
::
parse
(
dislo_orientation
&
orientation
,
std
::
stringstream
&
line
,
UInt
)
{
int
val
;
UInt
nb
=
parse
(
val
,
line
);
orientation
=
static_cast
<
dislo_orientation
>
(
val
);
return
nb
;
}
/* -------------------------------------------------------------------------- */
inline
void
CADD_template
::
setGeometries
(
int
pad_id
,
int
atom_id
,
int
geom_id
,
[[
gnu
::
unused
]]
int
boundary_id
)
{
this
->
pad_geom
=
pad_id
;
this
->
md_geom
=
atom_id
;
this
->
detect_geom
=
geom_id
;
}
/* -------------------------------------------------------------------------- */
inline
void
CADD_template
::
clearAll
()
{
// this->md_cores.clear();
// this->pad_atoms.empty();
// this->core_templates.clear();
this
->
out_computes
.
clear
();
}
/* -------------------------------------------------------------------------- */
inline
Real
CADD_template
::
compute_dd_theta
(
Real
x3
[],
Real
x4
[])
{
Real
norm1
,
bx1
,
bz1
,
bx2
,
bz2
,
norm2
,
inner
,
angle_with_ONEZERO
,
angdist
;
bx1
=
0.0
;
bz1
=
0.0
;
bx2
=
1.0
;
bz2
=
0.0
;
norm1
=
1.0
;
norm2
=
sqrt
(
pow
((
x4
[
0
]
-
x3
[
0
]),
2
)
+
pow
((
x4
[
1
]
-
x3
[
1
]),
2
));
inner
=
(
bx2
-
bx1
)
*
(
x4
[
0
]
-
x3
[
0
])
+
(
bz2
-
bz1
)
*
(
x4
[
1
]
-
x3
[
1
]);
angle_with_ONEZERO
=
acos
(
inner
/
norm1
/
norm2
)
*
180.0
/
M_PI
;
angdist
=
fabs
(
90.0
-
abs
(
angle_with_ONEZERO
));
return
(
180.0
-
angdist
);
}
/* -------------------------------------------------------------------------- */
template
<
typename
Vec
>
inline
bool
CADD_template
::
intersectTwoD
(
Vec
&
p1
,
Vec
&
p2
,
Vec
&
p3
,
Vec
&
p4
,
Vector
<
2
>
&
pts
,
UInt
slip_case
)
{
if
(
intersect
(
p1
,
p2
,
p3
,
p4
,
slip_case
))
{
UInt
id0
=
0
;
UInt
id1
=
1
;
if
(
slip_case
==
0
)
{
id0
=
1
;
id1
=
2
;
}
else
if
(
slip_case
==
1
)
{
id0
=
0
;
id1
=
2
;
}
else
if
(
slip_case
==
2
)
{
id0
=
0
;
id1
=
1
;
}
Real
det
=
(
p1
[
id0
]
-
p2
[
id0
])
*
(
p3
[
id1
]
-
p4
[
id1
])
-
(
p1
[
id1
]
-
p2
[
id1
])
*
(
p3
[
id0
]
-
p4
[
id0
]);
Real
nominator1
=
(
p1
[
id0
]
*
p2
[
id1
]
-
p1
[
id1
]
*
p2
[
id0
])
*
(
p3
[
id0
]
-
p4
[
id0
])
-
(
p1
[
id0
]
-
p2
[
id0
])
*
(
p3
[
id0
]
*
p4
[
id1
]
-
p3
[
id1
]
*
p4
[
id0
]);
Real
nominator2
=
(
p1
[
id0
]
*
p2
[
id1
]
-
p1
[
id1
]
*
p2
[
id0
])
*
(
p3
[
id1
]
-
p4
[
id1
])
-
(
p1
[
id1
]
-
p2
[
id1
])
*
(
p3
[
id0
]
*
p4
[
id1
]
-
p3
[
id1
]
*
p4
[
id0
]);
Real
px
=
nominator1
/
det
;
Real
py
=
nominator2
/
det
;
Real
ratioid0
;
if
(
fabs
(
p2
[
id0
]
-
p1
[
id0
])
<
1e-5
)
ratioid0
=
0.0
;
else
ratioid0
=
roundf
((
px
-
p1
[
id0
])
*
10.0
)
/
10.0
/
(
p2
[
id0
]
-
p1
[
id0
]);
Real
ratioid1
;
if
(
fabs
(
p2
[
id1
]
-
p1
[
id1
])
<
1e-5
)
ratioid1
=
0.0
;
else
ratioid1
=
roundf
((
py
-
p1
[
id1
])
*
10.0
)
/
10.0
/
(
p2
[
id1
]
-
p1
[
id1
]);
if
((
ratioid0
<=
1.0
)
&&
(
0.0
<=
ratioid0
)
&&
(
ratioid1
<=
1.0
)
&&
(
0.0
<=
ratioid1
))
{
pts
[
0
]
=
px
;
pts
[
1
]
=
py
;
return
true
;
}
else
{
pts
[
0
]
=
-
1.0
;
pts
[
1
]
=
-
1.0
;
return
false
;
}
}
pts
[
0
]
=
-
1.0
;
pts
[
1
]
=
-
1.0
;
return
false
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
Vec
>
inline
bool
CADD_template
::
intersectThreeD
(
Vec
&
p1
,
Vec
&
p2
,
Vec
&
p3
,
Vec
&
p4
,
Vec
&
pts
,
UInt
slip_case
)
{
// different slip plane
if
((
fabs
(
p1
[
1
]
-
p3
[
1
])
>
1.0
)
||
(
fabs
(
p2
[
1
]
-
p3
[
1
])
>
1.0
)
||
(
fabs
(
p1
[
1
]
-
p4
[
1
])
>
1.0
)
||
(
fabs
(
p2
[
1
]
-
p4
[
1
])
>
1.0
))
{
return
false
;
}
else
{
Vector
<
2
>
pts2D
=
{
pts
[
0
],
pts
[
1
]};
if
(
intersectTwoD
(
p1
,
p2
,
p3
,
p4
,
pts2D
,
slip_case
))
{
pts
[
1
]
=
p1
[
1
];
return
true
;
}
else
{
return
false
;
}
}
return
false
;
}
/* -------------------------------------------------------------------------- */
inline
bool
CADD_template
::
intersectTwoD_infinite
(
Real
p1
[],
Real
p2
[],
Real
p3
[],
Real
p4
[],
Real
*
ptsx
,
Real
*
ptsz
)
{
Real
x1
=
p1
[
0
];
Real
y1
=
p1
[
2
];
Real
x2
=
p2
[
0
];
Real
y2
=
p2
[
2
];
Real
x3
=
p3
[
0
];
Real
y3
=
p3
[
2
];
Real
x4
=
p4
[
0
];
Real
y4
=
p4
[
2
];
Real
denom
=
(
x1
-
x2
)
*
(
y3
-
y4
)
-
(
y1
-
y2
)
*
(
x3
-
x4
);
if
(
fabs
(
denom
)
<=
1e-3
)
{
*
ptsx
=
-
1.0
;
*
ptsz
=
-
1.0
;
return
false
;
}
else
{
*
ptsx
=
((
x1
*
y2
-
y1
*
x2
)
*
(
x3
-
x4
)
-
(
x1
-
x2
)
*
(
x3
*
y4
-
y3
*
x4
))
/
denom
;
*
ptsz
=
((
x1
*
y2
-
y1
*
x2
)
*
(
y3
-
y4
)
-
(
y1
-
y2
)
*
(
x3
*
y4
-
y3
*
x4
))
/
denom
;
return
true
;
}
}
/* -------------------------------------------------------------------------- */
template
<
typename
Vec
>
inline
bool
CADD_template
::
intersect
(
Vec
&
x1
,
Vec
&
x2
,
Vec
&
x3
,
Vec
&
x4
,
UInt
slip_case
)
{
if
(
slip_case
==
0
)
return
ccw
(
x1
[
1
],
x1
[
2
],
x3
[
1
],
x3
[
2
],
x4
[
1
],
x4
[
2
])
!=
ccw
(
x2
[
1
],
x2
[
2
],
x3
[
1
],
x3
[
2
],
x4
[
1
],
x4
[
2
])
&&
ccw
(
x1
[
1
],
x1
[
2
],
x2
[
1
],
x2
[
2
],
x3
[
1
],
x3
[
2
])
!=
ccw
(
x1
[
1
],
x1
[
2
],
x2
[
1
],
x2
[
2
],
x4
[
1
],
x4
[
2
]);
else
if
(
slip_case
==
1
)
return
ccw
(
x1
[
0
],
x1
[
2
],
x3
[
0
],
x3
[
2
],
x4
[
0
],
x4
[
2
])
!=
ccw
(
x2
[
0
],
x2
[
2
],
x3
[
0
],
x3
[
2
],
x4
[
0
],
x4
[
2
])
&&
ccw
(
x1
[
0
],
x1
[
2
],
x2
[
0
],
x2
[
2
],
x3
[
0
],
x3
[
2
])
!=
ccw
(
x1
[
0
],
x1
[
2
],
x2
[
0
],
x2
[
2
],
x4
[
0
],
x4
[
2
]);
else
if
(
slip_case
==
2
)
return
ccw
(
x1
[
0
],
x1
[
1
],
x3
[
0
],
x3
[
1
],
x4
[
0
],
x4
[
1
])
!=
ccw
(
x2
[
0
],
x2
[
1
],
x3
[
0
],
x3
[
1
],
x4
[
0
],
x4
[
1
])
&&
ccw
(
x1
[
0
],
x1
[
1
],
x2
[
0
],
x2
[
1
],
x3
[
0
],
x3
[
1
])
!=
ccw
(
x1
[
0
],
x1
[
1
],
x2
[
0
],
x2
[
1
],
x4
[
0
],
x4
[
1
]);
return
false
;
}
/* -------------------------------------------------------------------------- */
inline
Real
CADD_template
::
ccw
(
Real
Ax
,
Real
Ay
,
Real
Bx
,
Real
By
,
Real
Cx
,
Real
Cy
)
{
return
(
Cy
-
Ay
)
*
(
Bx
-
Ax
)
>
(
By
-
Ay
)
*
(
Cx
-
Ax
);
}
/* -------------------------------------------------------------------------- */
inline
Real
CADD_template
::
determinant
(
Real
a
,
Real
b
,
Real
c
,
Real
d
,
Real
e
,
Real
f
,
Real
g
,
Real
h
,
Real
i
)
{
return
a
*
(
e
*
i
-
f
*
h
)
-
b
*
(
d
*
i
-
f
*
g
)
+
c
*
(
d
*
h
-
e
*
g
);
}
/* -------------------------------------------------------------------------- */
inline
Real
CADD_template
::
radian2degree
(
Real
rad
)
{
// switch radian to degree having a range (0,360)
Real
ang
=
rad
*
180.0
/
M_PI
;
if
(
ang
<=
0.0
)
return
360.0
+
ang
;
else
return
ang
;
}
/* -------------------------------------------------------------------------- */
inline
Real
CADD_template
::
angle_on_plane
(
Real
pos
[],
UInt
axis
)
{
if
(
axis
==
0
)
{
return
atan
(
pos
[
1
]
/
pos
[
2
]);
}
else
if
(
axis
==
1
)
{
return
atan
(
pos
[
2
]
/
pos
[
0
]);
}
else
if
(
axis
==
2
)
{
return
atan
(
pos
[
0
]
/
pos
[
1
]);
}
else
{
LM_FATAL
(
"wrong axis"
);
}
}
// inline void CADD_template::RotationMatrixTemplateToRealLattices(Real theta,
// Real theta0,
// std::vector<Real>
// &resultR) {
// // |R[0] R[1] R[2]|
// // |R[3] R[4] R[5]|
// // |R[6] R[7] R[8]|
// resultR[0] = cos((theta - theta0) * M_PI / 180.0);
// resultR[2] = -sin((theta - theta0) * M_PI / 180.0);
// resultR[6] = sin((theta - theta0) * M_PI / 180.0);
// resultR[8] = cos((theta - theta0) * M_PI / 180.0);
// resultR[1] = this->normal_quant[0];
// resultR[4] = this->normal_quant[1];
// resultR[7] = this->normal_quant[2];
// resultR[3] = 0.0;
// resultR[5] = 0.0;
// }
// inline void CADD_template::MatrixMultiplication(std::vector<Real> &mt, Real
// vec[],
// Real result[]) {
// result[0] = mt[0] * vec[0] + mt[1] * vec[1] + mt[2] * vec[2];
// result[1] = mt[3] * vec[0] + mt[4] * vec[1] + mt[5] * vec[2];
// result[2] = mt[6] * vec[0] + mt[7] * vec[1] + mt[8] * vec[2];
// }
// inline void CADD_template::CrossProduct(Real vec1[], Real vec2[], Real
// result[]) {
// result[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
// result[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
// result[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
// }
/* -------------------------------------------------------------------------- */
inline
Real
CADD_template
::
distance_between_line_and_point
(
Real
x0
,
Real
z0
,
Real
pts1
[],
Real
pts2
[])
{
Real
x2
=
pts2
[
0
];
Real
x1
=
pts1
[
0
];
Real
z2
=
pts2
[
2
];
Real
z1
=
pts1
[
2
];
Real
dist
=
fabs
((
x1
-
x2
)
*
(
z1
-
z0
)
-
(
x1
-
x0
)
*
(
z2
-
z1
))
/
sqrt
((
x2
-
x1
)
*
(
x2
-
x1
)
+
(
z2
-
z1
)
*
(
z2
-
z1
));
return
dist
;
}
/* -------------------------------------------------------------------------- */
inline
void
CADD_template
::
intersect_of_point_on_line_3d
(
Real
p0
[],
Real
p1
[],
Real
p2
[],
Real
result
[])
{
Real
nom
=
(
p1
[
0
]
-
p0
[
0
])
*
(
p2
[
0
]
-
p1
[
0
])
+
(
p1
[
1
]
-
p0
[
1
])
*
(
p2
[
1
]
-
p1
[
1
])
+
(
p1
[
2
]
-
p0
[
2
])
*
(
p2
[
2
]
-
p1
[
2
]);
Real
denom
=
(
p2
[
0
]
-
p1
[
0
])
*
(
p2
[
0
]
-
p1
[
0
])
+
(
p2
[
1
]
-
p1
[
1
])
*
(
p2
[
1
]
-
p1
[
1
])
+
(
p2
[
2
]
-
p1
[
2
])
*
(
p2
[
2
]
-
p1
[
2
]);
Real
t
=
-
1.0
*
nom
/
(
denom
+
1e-8
);
result
[
0
]
=
p1
[
0
]
+
(
p2
[
0
]
-
p1
[
0
])
*
t
;
result
[
1
]
=
p1
[
1
]
+
(
p2
[
1
]
-
p1
[
1
])
*
t
;
result
[
2
]
=
p1
[
2
]
+
(
p2
[
2
]
-
p1
[
2
])
*
t
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
Vec
>
inline
Real
CADD_template
::
distance_between_line_and_point3D
(
Vec
&
pos0
,
Vec
&
pts1
,
Vec
&
pts2
)
{
Real
x0
=
pos0
[
0
];
Real
y0
=
pos0
[
1
];
Real
z0
=
pos0
[
2
];
Real
u1
=
x0
-
pts1
[
0
];
Real
u2
=
y0
-
pts1
[
1
];
Real
u3
=
z0
-
pts1
[
2
];
Real
v1
=
x0
-
pts2
[
0
];
Real
v2
=
y0
-
pts2
[
1
];
Real
v3
=
z0
-
pts2
[
2
];
Real
cross
[
3
]
=
{
u2
*
v3
-
u3
*
v2
,
u3
*
v1
-
u1
*
v3
,
u1
*
v2
-
u2
*
v1
};
Real
norm_cross
=
sqrt
(
cross
[
0
]
*
cross
[
0
]
+
cross
[
1
]
*
cross
[
1
]
+
cross
[
2
]
*
cross
[
2
]);
Real
r1
=
pts2
[
0
]
-
pts1
[
0
];
Real
r2
=
pts2
[
1
]
-
pts1
[
1
];
Real
r3
=
pts2
[
2
]
-
pts1
[
2
];
Real
denom
=
sqrt
(
r1
*
r1
+
r2
*
r2
+
r3
*
r3
);
Real
dist
=
norm_cross
/
denom
;
return
dist
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
Vec
>
inline
Real
CADD_template
::
AreaTriangle
(
Vec
&
p1
,
Vec
&
p2
,
Vec
&
p3
,
UInt
slip_case
)
{
Real
nom
=
0.0
;
if
(
slip_case
==
0
)
nom
=
p1
[
1
]
*
(
p2
[
2
]
-
p3
[
2
])
+
p2
[
1
]
*
(
p3
[
2
]
-
p1
[
2
])
+
p3
[
1
]
*
(
p1
[
2
]
-
p2
[
2
]);
else
if
(
slip_case
==
1
)
nom
=
p1
[
0
]
*
(
p2
[
2
]
-
p3
[
2
])
+
p2
[
0
]
*
(
p3
[
2
]
-
p1
[
2
])
+
p3
[
0
]
*
(
p1
[
2
]
-
p2
[
2
]);
else
if
(
slip_case
==
2
)
nom
=
p1
[
0
]
*
(
p2
[
1
]
-
p3
[
1
])
+
p2
[
0
]
*
(
p3
[
1
]
-
p1
[
1
])
+
p3
[
0
]
*
(
p1
[
1
]
-
p2
[
1
]);
Real
area
=
fabs
(
nom
/
2.0
)
+
1e-8
;
return
area
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
Vec
>
inline
bool
CADD_template
::
is_it_closed
(
Vec
cp
,
Vec
&
cp0
,
Real
&
ci
,
Real
&
ci0
)
{
Real
dist
=
(
cp
-
cp0
).
norm
();
if
((
dist
<
1e-8
)
&&
(
fabs
(
ci0
-
ci
)
<
1e-5
))
{
return
true
;
}
else
{
return
false
;
}
}
/* -------------------------------------------------------------------------- */
inline
UInt
CADD_template
::
next_index
(
UInt
next_idx
,
UInt
num_dd_pts0
,
std
::
vector
<
UInt
>
indices
)
{
for
(
UInt
i
=
0
;
i
<
num_dd_pts0
;
++
i
)
{
if
(
!
(
std
::
find
(
indices
.
begin
(),
indices
.
end
(),
i
)
!=
indices
.
end
()))
{
return
i
;
}
}
return
next_idx
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
Vec
>
inline
UInt
CADD_template
::
find_one_direction
(
Vec
&
nb1
[[
gnu
::
unused
]],
Vec
&
nb2
[[
gnu
::
unused
]],
Real
b1mag
,
Real
b2mag
[[
gnu
::
unused
]],
UInt
idx1
,
UInt
idx2
,
UInt
*
next_i
)
{
if
(
b1mag
>
0.0
)
{
*
next_i
=
idx1
;
return
0
;
}
else
{
*
next_i
=
idx2
;
return
1
;
}
}
/* -------------------------------------------------------------------------- */
inline
Real
CADD_template
::
DistanceBetweenTwoNodes
(
Real
x1
,
Real
y1
,
Real
z1
,
Real
x2
,
Real
y2
,
Real
z2
)
{
Real
dist
=
sqrt
((
x1
-
x2
)
*
(
x1
-
x2
)
+
(
y1
-
y2
)
*
(
y1
-
y2
)
+
(
z1
-
z2
)
*
(
z1
-
z2
));
return
dist
;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif
/* __LIBMULTISCALE_TEMPLATE_CADD_HH__ */
Event Timeline
Log In to Comment