Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76001213
compute_dislodetection.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, Aug 5, 16:00
Size
34 KB
Mime Type
text/x-c++
Expires
Wed, Aug 7, 16:00 (2 d)
Engine
blob
Format
Raw Data
Handle
19646576
Attached To
rLIBMULTISCALE LibMultiScale
compute_dislodetection.cc
View Options
/**
* @file compute_dislocdetn.cc
*
* @author Max Hodapp <max.hodapp@epfl.ch>
*
* @date Thu Oct 29 2015
*
* @brief Dislocation detection for single crystals
*
* @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/>.
*
*/
/* TODO:
* - check 'Delaunay_Triangulation_3' constructor in CGAL
* --> try to avoid redundant memcopies to convert Reals into Points!!
* - add 2-D compatibility
* - add compatibility for 1) several dislocation lines, 2) dislocation loops
* - detection of partial dislocations
* - detection of dislocations in different crystals (e.g. bcc)
* - the method for MPI communication should be generalized s.t. they can be
* used by all classes which take possibly distributed data as input
*
*/
/* -------------------------------------------------------------------------- */
#include "compute_dislodetection.hh"
#include "lib_dumper.hh"
#include "lib_filter.hh"
#include "lib_md.hh"
#include "lib_stimulation.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
#define __TOL__ 1e-16
// numerical tolerance
/* -------------------------------------------------------------------------- */
// Forward declarations
// Methods are straight-forward modifications from the original CGAL
// implementation (see 'Fixed_alpha_shape_3.h')
template
<
typename
OutputIterator
>
OutputIterator
get_all_alpha_shape_edges
(
Ta
&
tet_alpha
,
OutputIterator
it
);
template
<
typename
OutputIterator
>
OutputIterator
get_all_alpha_shape_faces
(
Ta
&
tet_alpha
,
OutputIterator
it
);
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* IMPLEMENTATION - CONSTRUCTOR */
/* -------------------------------------------------------------------------- */
ComputeDisloDetection
::
ComputeDisloDetection
(
const
std
::
string
&
name
)
:
LMObject
(
name
),
meshTet
(
name
+
"-meshTet"
)
{
// Add 'meshTet' to filter manager (--> the object is now visible to the
// parser)
LM_TOIMPLEMENT
;
// FilterManager::getManager().addObject(&this->meshTet, false);
}
/* -------------------------------------------------------------------------- */
ComputeDisloDetection
::~
ComputeDisloDetection
()
{}
/* -------------------------------------------------------------------------- */
/* IMPLEMENTATION - PUBLIC MEMBER METHODS */
/* -------------------------------------------------------------------------- */
/* Main build */
template
<
typename
Cont
>
void
ComputeDisloDetection
::
build
(
Cont
&
cont
)
{
if
(
Cont
::
Dim
!=
3
)
{
LM_FATAL
(
"Dislocation detection only available in 3D!"
<<
std
::
endl
);
}
auto
&
output
=
allocOutput
<
ContainerGenericMesh
<
3
>>
(
"detected_segments"
,
"detected_segments"
);
/* ------------------------------------------------------------------------ */
// NOTE: Dislocation detection does not run in parallel
// ==> gather atomic positions from all other processes
/* ------------------------------------------------------------------------ */
ComputeExtract
positions
(
"ComputeExtract:"
+
this
->
getID
());
positions
.
setParam
(
"FIELD"
,
_position
);
positions
.
build
(
cont
);
ContainerArray
<
Real
>
pos_buffer
;
this
->
gatherAllData
(
pos_buffer
,
positions
.
evalArrayOutput
(),
0
);
// root process is 0
/* ------------------------------------------------------------------------ */
// Run dislocation detection only on the root process of the atomistic group
/* ------------------------------------------------------------------------ */
CommGroup
group
=
this
->
getCommGroup
();
if
(
group
.
size
()
>
1
)
{
if
(
group
.
getMyRank
()
==
0
)
{
DUMP
(
"Dislocation detection is currently not parallelized!"
,
DBG_MESSAGE
);
}
else
{
return
;
}
}
DUMP
(
"Starting dislocation detection ..."
,
DBG_MESSAGE
);
Real
norm2
[
3
];
// Containers/Lists
std
::
vector
<
Pt
>
pt_buffer
;
std
::
list
<
Vertex_handle
>
v_l
;
std
::
list
<
Cell_handle
>
c_l
;
/* ------------------------------------------------------------------------ */
// Delete previous mesh struct
output
.
getContainerNodes
().
clear
();
output
.
getContainerElems
().
clear
();
this
->
meshTet
.
getContainerNodes
().
clear
();
this
->
meshTet
.
getContainerElems
().
clear
();
/* ------------------------------------------------------------------------ */
/* Mesh atomistic domain */
// T tet;
//
// {
// std::vector< Real >::iterator pos_it;
// for ( pos_it = pos_buffer.begin(); pos_it != pos_buffer.end();
// std::advance(pos_it,3) )
// {
// tet.insert( Pt( pos_it[0], pos_it[1], pos_it[2] ) );
// }
// }
// Convert full position matrix to 'Point_3'
// --> CGAL seems to accept only nodes of this type
// NOTE: intuitively one would like to to do sth like
// T tet( pos_buffer.begin(), pos_buffer.end() )
// to avoid redundant copy operations and memory capacity
// --> however, the second method seems to be faster...
LM_TOIMPLEMENT
;
// {
// std::vector<Real>::iterator pos_it;
// for (pos_it = pos_buffer.begin(); pos_it != pos_buffer.end();
// std::advance(pos_it, 3)) {
// pt_buffer.push_back(Pt(pos_it[0], pos_it[1], pos_it[2]));
// }
// pos_buffer.clear();
// }
// Tetrahedralize the atomistic domain
T
tet
(
pt_buffer
.
begin
(),
pt_buffer
.
end
());
pt_buffer
.
clear
();
// Ceck validity of 'tet'
assert
(
tet
.
is_valid
());
// Generate the alpha shape from 'tet'
// NOTE: 'tet' will be destroyed
Ta
tet_alpha
(
tet
,
this
->
alpha
);
// Get iterators to the simplices of the alpha shape
tet_alpha
.
get_alpha_shape_vertices
(
std
::
back_inserter
(
v_l
),
Ta
::
REGULAR
);
tet_alpha
.
get_alpha_shape_cells
(
std
::
back_inserter
(
c_l
),
Ta
::
INTERIOR
);
DUMP
(
"-- Tetrahedralization successful! --"
,
DBG_MESSAGE
);
DUMP
(
"Vertices ... "
<<
v_l
.
size
(),
DBG_MESSAGE
);
DUMP
(
"Cells ...... "
<<
c_l
.
size
(),
DBG_MESSAGE
);
/* ------------------------------------------------------------------------ */
// Detect distorted tetrahedra
// Normalize the orientation matrix
norm2
[
0
]
=
sqrt
(
this
->
orient
[
0
]
*
this
->
orient
[
0
]
+
this
->
orient
[
1
]
*
this
->
orient
[
1
]
+
this
->
orient
[
2
]
*
this
->
orient
[
2
]);
norm2
[
1
]
=
sqrt
(
this
->
orient
[
3
]
*
this
->
orient
[
3
]
+
this
->
orient
[
4
]
*
this
->
orient
[
4
]
+
this
->
orient
[
5
]
*
this
->
orient
[
5
]);
norm2
[
2
]
=
sqrt
(
this
->
orient
[
6
]
*
this
->
orient
[
6
]
+
this
->
orient
[
7
]
*
this
->
orient
[
7
]
+
this
->
orient
[
8
]
*
this
->
orient
[
8
]);
this
->
orient
[
0
]
/=
norm2
[
0
];
this
->
orient
[
3
]
/=
norm2
[
1
];
this
->
orient
[
6
]
/=
norm2
[
2
];
this
->
orient
[
1
]
/=
norm2
[
0
];
this
->
orient
[
4
]
/=
norm2
[
1
];
this
->
orient
[
7
]
/=
norm2
[
2
];
this
->
orient
[
2
]
/=
norm2
[
0
];
this
->
orient
[
5
]
/=
norm2
[
1
];
this
->
orient
[
8
]
/=
norm2
[
2
];
// Assign reference bond vectors
this
->
setRefBondVec
();
// Run the detection
this
->
detect
(
output
,
tet_alpha
);
/* ------------------------------------------------------------------------ */
// Output tetMesh --> skip in test version since mesh cannot yet be visualized
// inside LM
this
->
outputMesh
(
v_l
,
c_l
);
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* IMPLEMENTATION - PRIVATE MEMBER METHODS */
/* -------------------------------------------------------------------------- */
/* Define the set of ideal reference vectors */
void
ComputeDisloDetection
::
setRefBondVec
()
{
if
(
this
->
lattice_t
==
std
::
string
(
"fcc"
))
{
if
(
this
->
disloc_t
==
std
::
string
(
"full"
))
{
this
->
Nref
=
18
;
this
->
bVecRef
=
new
Vec
[
18
];
// First nearest neighbors
this
->
bVecRef
[
0
]
=
Vec
(
0.5
,
0.5
,
0.
);
this
->
bVecRef
[
1
]
=
Vec
(
-
0.5
,
-
0.5
,
0.
);
this
->
bVecRef
[
2
]
=
Vec
(
0.5
,
0.
,
0.5
);
this
->
bVecRef
[
3
]
=
Vec
(
-
0.5
,
0.
,
-
0.5
);
this
->
bVecRef
[
4
]
=
Vec
(
0.
,
0.5
,
0.5
);
this
->
bVecRef
[
5
]
=
Vec
(
0.
,
-
0.5
,
-
0.5
);
this
->
bVecRef
[
6
]
=
Vec
(
-
0.5
,
0.5
,
0.
);
this
->
bVecRef
[
7
]
=
Vec
(
0.5
,
-
0.5
,
0.
);
this
->
bVecRef
[
8
]
=
Vec
(
-
0.5
,
0.
,
0.5
);
this
->
bVecRef
[
9
]
=
Vec
(
0.5
,
0.
,
-
0.5
);
this
->
bVecRef
[
10
]
=
Vec
(
0.
,
0.5
,
-
0.5
);
this
->
bVecRef
[
11
]
=
Vec
(
0.
,
-
0.5
,
0.5
);
// Second nearest neighbors
this
->
bVecRef
[
12
]
=
Vec
(
1.
,
0.
,
0.
);
this
->
bVecRef
[
13
]
=
Vec
(
-
1.
,
0.
,
0.
);
this
->
bVecRef
[
14
]
=
Vec
(
0.
,
1.
,
0.
);
this
->
bVecRef
[
15
]
=
Vec
(
0.
,
-
1.
,
0.
);
this
->
bVecRef
[
16
]
=
Vec
(
0.
,
0.
,
1.
);
this
->
bVecRef
[
17
]
=
Vec
(
0.
,
0.
,
-
1.
);
// Rotate the reference lattice to the crystals orientation
for
(
UInt
i
=
0
;
i
<
18
;
++
i
)
this
->
bVecRef
[
i
]
=
this
->
a0
*
Vec
(
this
->
orient
[
0
]
*
this
->
bVecRef
[
i
].
x
()
+
this
->
orient
[
1
]
*
this
->
bVecRef
[
i
].
y
()
+
this
->
orient
[
2
]
*
this
->
bVecRef
[
i
].
z
(),
this
->
orient
[
3
]
*
this
->
bVecRef
[
i
].
x
()
+
this
->
orient
[
4
]
*
this
->
bVecRef
[
i
].
y
()
+
this
->
orient
[
5
]
*
this
->
bVecRef
[
i
].
z
(),
this
->
orient
[
6
]
*
this
->
bVecRef
[
i
].
x
()
+
this
->
orient
[
7
]
*
this
->
bVecRef
[
i
].
y
()
+
this
->
orient
[
8
]
*
this
->
bVecRef
[
i
].
z
());
}
else
{
LM_TOIMPLEMENT
;
// detection of partial dislocations not yet supported
}
}
else
{
LM_TOIMPLEMENT
;
// other lattice types not yet supported
}
}
/* -------------------------------------------------------------------------- */
/* Compute the Burgers vector of a face */
// SHORT:
//
// Let a the faces of a tetrahedron be numbered 0,1,2,3. Then the Burgers vector
// for each face is defined as
//
// burgers_v_0 = v_5_ref - v_6_ref - v_4_ref,
// burgers_v_1 = v_2_ref + v_6_ref - v_3_ref,
// burgers_v_2 = v_3_ref - v_5_ref - v_1_ref,
// burgers_v_3 = v_1_ref + v_4_ref - v_2_ref,
//
// where the v_i_ref are the reference vectors corresponding to each of the 6
// edge vectors of a tetrahedron in the current configuration.
inline
Vec
ComputeDisloDetection
::
computeBurgersVec
(
const
Cell_handle
&
ch
,
// cell handle
const
UInt
idxFace
)
// face index
{
switch
(
idxFace
)
{
case
0
:
return
(
this
->
bVecRef
[
ch
->
info
()[
4
]]
-
this
->
bVecRef
[
ch
->
info
()[
5
]]
-
this
->
bVecRef
[
ch
->
info
()[
3
]]);
case
1
:
return
(
this
->
bVecRef
[
ch
->
info
()[
1
]]
+
this
->
bVecRef
[
ch
->
info
()[
5
]]
-
this
->
bVecRef
[
ch
->
info
()[
2
]]);
case
2
:
return
(
this
->
bVecRef
[
ch
->
info
()[
2
]]
-
this
->
bVecRef
[
ch
->
info
()[
4
]]
-
this
->
bVecRef
[
ch
->
info
()[
0
]]);
case
3
:
return
(
this
->
bVecRef
[
ch
->
info
()[
0
]]
+
this
->
bVecRef
[
ch
->
info
()[
3
]]
-
this
->
bVecRef
[
ch
->
info
()[
1
]]);
};
LM_FATAL
(
"Face index "
<<
idxFace
<<
" not valid!"
<<
std
::
endl
<<
"Range is (0, 1, 2, 3)."
<<
std
::
endl
);
}
/* -------------------------------------------------------------------------- */
/* Get remaining distorted faces of tetrahedron which has a distorted face */
// NOTE: dislocation junctions are currently not warranted!
UInt
ComputeDisloDetection
::
getDistortedFaces
(
Face
&
fh
,
// reference to 2nd distorted face
const
Cell_handle
&
ch
,
// cell handle
const
UInt
idxFace
)
// index of the detected face
{
UInt
ctr
=
0
;
// nbr of distorted faces
switch
(
idxFace
)
{
case
0
:
if
(
this
->
computeBurgersVec
(
ch
,
1
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
1
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
2
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
2
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
3
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
3
);
++
ctr
;
}
break
;
case
1
:
if
(
this
->
computeBurgersVec
(
ch
,
0
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
0
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
2
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
2
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
3
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
3
);
++
ctr
;
}
break
;
case
2
:
if
(
this
->
computeBurgersVec
(
ch
,
0
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
0
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
1
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
1
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
3
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
3
);
++
ctr
;
}
break
;
case
3
:
if
(
this
->
computeBurgersVec
(
ch
,
0
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
0
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
1
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
1
);
++
ctr
;
}
if
(
this
->
computeBurgersVec
(
ch
,
2
).
squared_length
()
>
__TOL__
)
{
fh
=
std
::
make_pair
(
ch
,
2
);
++
ctr
;
}
break
;
};
if
(
ctr
>
0
)
return
ctr
;
else
LM_FATAL
(
"Dislocation cannot end within a crystal"
<<
std
::
endl
);
}
/* -------------------------------------------------------------------------- */
/* Main method during a dislocation detection */
template
<
typename
Mesh
>
void
ComputeDisloDetection
::
detect
(
Mesh
&
mesh
,
Ta
&
tet_alpha
)
{
UInt
idx
,
size_buffer
,
Nlines
=
0
,
Nloops
=
0
;
std
::
pair
<
UInt
,
UInt
>
idxVh
;
Real
norm2
,
norm2_buffer
;
// LM
RefPointData
<
3
>
node
;
// LM - Containers/Lists
std
::
vector
<
RefPointData
<
3
>>
nodes_buffer
;
// CGAL
Pt
pt_centroid
;
Vec
eVec
;
Vertex_handle
vh1
,
vh2
;
Face
fh
;
Cell_handle
ch
;
// CGAL - Containers/Lists
std
::
list
<
Edge
>
eh_l
;
std
::
list
<
Face
>
fh_l
;
std
::
list
<
Face
>
fhR_l
;
std
::
list
<
Face
>
fhI_l
;
// CGAL - Iterators
std
::
list
<
Edge
>::
iterator
eh_it
;
std
::
list
<
Face
>::
iterator
fh_it
;
Ta
::
Cell_circulator
ch_circ
;
/* ------------------------------------------------------------------------ */
// Get list of all edges & faces from the alpha shape
get_all_alpha_shape_edges
(
tet_alpha
,
std
::
back_inserter
(
eh_l
));
get_all_alpha_shape_faces
(
tet_alpha
,
std
::
back_inserter
(
fh_l
));
// Get list of all faces which are on the boundary of the alpha shape
tet_alpha
.
get_alpha_shape_facets
(
std
::
back_inserter
(
fhR_l
),
Ta
::
REGULAR
);
// Get list of all faces which are in the interior of the alpha shape
tet_alpha
.
get_alpha_shape_facets
(
std
::
back_inserter
(
fhI_l
),
Ta
::
INTERIOR
);
// Iterate over all edges of the alpha shape
// An edge handle is a triple < Cell_handle, int int >, where the
// Cell_handle
// refers to a cell which contains the edge and the both integers are the
// indices
// of its vertices
for
(
eh_it
=
eh_l
.
begin
();
eh_it
!=
eh_l
.
end
();
++
eh_it
)
{
/* ********************************************************************** */
// 1) Compute the reference vector the current edge
// First obtain the vertex handles from the edge list ...
vh1
=
(
*
eh_it
).
first
->
vertex
((
*
eh_it
).
second
);
vh2
=
(
*
eh_it
).
first
->
vertex
((
*
eh_it
).
third
);
// ... then compute the vector p2-p1
eVec
=
vh2
->
point
()
-
vh1
->
point
();
// Compute the squared l^2-norm for all ev_ref(i)-ev and assigend the
// element
// of V_ref which minimzes it
// assert( Nref != 0 );
norm2_buffer
=
Vec
(
this
->
bVecRef
[
0
]
-
eVec
).
squared_length
();
idx
=
0
;
for
(
UInt
i
=
1
;
i
<
this
->
Nref
;
++
i
)
{
norm2
=
Vec
(
this
->
bVecRef
[
i
]
-
eVec
).
squared_length
();
if
(
norm2
<
norm2_buffer
)
{
norm2_buffer
=
norm2
;
idx
=
i
;
}
}
/* ********************************************************************** */
// 2) Store the corresponding reference vector in all cells incident to
// '*eit'
ch_circ
=
tet_alpha
.
incident_cells
((
*
eh_it
),
(
*
eh_it
).
first
);
ch
=
ch_circ
;
do
{
// Get the indices of each vertex
idxVh
.
first
=
ch_circ
->
index
(
vh1
);
idxVh
.
second
=
ch_circ
->
index
(
vh2
);
switch
(
idxVh
.
first
)
{
case
0
:
switch
(
idxVh
.
second
)
{
case
1
:
ch_circ
->
info
()[
0
]
=
idx
;
break
;
case
2
:
ch_circ
->
info
()[
1
]
=
idx
;
break
;
case
3
:
ch_circ
->
info
()[
2
]
=
idx
;
break
;
};
break
;
case
1
:
switch
(
idxVh
.
second
)
{
case
0
:
ch_circ
->
info
()[
0
]
=
idx
^
1
;
break
;
case
2
:
ch_circ
->
info
()[
3
]
=
idx
;
break
;
case
3
:
ch_circ
->
info
()[
4
]
=
idx
;
break
;
};
break
;
case
2
:
switch
(
idxVh
.
second
)
{
case
0
:
ch_circ
->
info
()[
1
]
=
idx
^
1
;
break
;
case
1
:
ch_circ
->
info
()[
3
]
=
idx
^
1
;
break
;
case
3
:
ch_circ
->
info
()[
5
]
=
idx
;
break
;
};
break
;
case
3
:
switch
(
idxVh
.
second
)
{
case
0
:
ch_circ
->
info
()[
2
]
=
idx
^
1
;
break
;
case
1
:
ch_circ
->
info
()[
4
]
=
idx
^
1
;
break
;
case
2
:
ch_circ
->
info
()[
5
]
=
idx
^
1
;
break
;
};
break
;
}
ch_circ
++
;
}
while
(
ch_circ
!=
ch
);
/* ********************************************************************** */
}
/* ------------------------------------------------------------------------ */
// 3) Dislocation line tracking
UInt
ctr
=
0
;
// counts the number of distorted faces per tetrahedron
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
while
(
fhR_l
.
size
()
>
0
)
{
fh_it
=
fhR_l
.
begin
();
// set iterator to the first element of the list
// Get the cell and the index of the opposed vertex of the current face
ch
=
(
*
fh_it
).
first
;
idx
=
(
*
fh_it
).
second
;
/* ********************************************************************** */
// If the Burgers vector is non-zero ==> mark both cells containing the
// current face
// --> the neighboring cell can simply be accessed via the index of the
// opposed vertex
if
(
this
->
computeBurgersVec
(
ch
,
idx
).
squared_length
()
>
__TOL__
)
{
// Check if current cell is an element of the alpha shape
// ==> if 'ch' is not interior --> set 'ch' to its opposed neighbor and
// update the index of the current face
if
(
tet_alpha
.
classify
(
ch
)
!=
Ta
::
INTERIOR
)
{
ch
=
ch
->
neighbor
(
idx
);
idx
=
tet_alpha
.
mirror_facet
(
*
fh_it
).
second
;
}
// Add the center of mass of the current cell to the list of nodes
pt_centroid
=
CGAL
::
centroid
(
Tet
(
ch
->
vertex
(
0
)
->
point
(),
ch
->
vertex
(
1
)
->
point
(),
ch
->
vertex
(
2
)
->
point
(),
ch
->
vertex
(
3
)
->
point
()));
// Push the center of mass to the set of nodes
node
.
position
()[
0
]
=
pt_centroid
.
x
();
node
.
position
()[
1
]
=
pt_centroid
.
y
();
node
.
position
()[
2
]
=
pt_centroid
.
z
();
nodes_buffer
.
push_back
(
node
);
// Face handle can now be savely removed from the list of boundary faces
fhR_l
.
remove
(
*
fh_it
);
// Compute the Burgers vectors of the 3 remaining faces
ctr
=
getDistortedFaces
(
fh
,
ch
,
idx
);
// Abort program if more than 2 faces with non-zero Burgers vector are
// found
if
(
ctr
>
1
)
LM_FATAL
(
std
::
endl
<<
"Number of faces with non-zero Burgers vector: "
<<
(
ctr
+
1
)
<<
std
::
endl
<<
"Dislocation junctions not supported yet!"
<<
std
::
endl
);
ctr
=
0
;
/* ******************************************************************** */
// Track the dislocation line through the atomistic domain as long as the
// subsequent face is in its interior
while
(
tet_alpha
.
classify
(
fh
)
!=
Ta
::
REGULAR
)
{
// Get the opposed neighbor and the index of the face thereof
ch
=
ch
->
neighbor
(
fh
.
second
);
idx
=
tet_alpha
.
mirror_facet
(
fh
).
second
;
// Add the center of mass of the current cell to the list of nodes
pt_centroid
=
CGAL
::
centroid
(
Tet
(
ch
->
vertex
(
0
)
->
point
(),
ch
->
vertex
(
1
)
->
point
(),
ch
->
vertex
(
2
)
->
point
(),
ch
->
vertex
(
3
)
->
point
()));
// Push the center of mass to the set of nodes
node
.
position
()[
0
]
=
pt_centroid
.
x
();
node
.
position
()[
1
]
=
pt_centroid
.
y
();
node
.
position
()[
2
]
=
pt_centroid
.
z
();
nodes_buffer
.
push_back
(
node
);
// Remove the face from the list of interior faces
// --> since there is no iterator pointing to 'fh', perform check if
// size is reduced ...
size_buffer
=
fhI_l
.
size
();
fhI_l
.
remove
(
fh
);
if
(
fhI_l
.
size
()
==
size_buffer
)
// ... else remove the face handle
// defined w.r.t. the opposed cell
fhI_l
.
remove
(
tet_alpha
.
mirror_facet
(
fh
));
// Compute the Burgers vectors of the 3 remaining faces
ctr
=
getDistortedFaces
(
fh
,
ch
,
idx
);
// Abort program if more than 2 faces with non-zero Burgers vector are
// found
if
(
ctr
>
1
)
LM_FATAL
(
std
::
endl
<<
"Number of faces with non-zero Burgers vector: "
<<
(
ctr
+
1
)
<<
std
::
endl
<<
"Dislocation junctions not supported yet!"
<<
std
::
endl
);
ctr
=
0
;
};
/* ******************************************************************** */
// Remove the last face from the list of boundary faces
size_buffer
=
fhR_l
.
size
();
fhR_l
.
remove
(
fh
);
if
(
fhR_l
.
size
()
==
size_buffer
)
fhR_l
.
remove
(
tet_alpha
.
mirror_facet
(
fh
));
/* ******************************************************************** */
// Add the DD nodes to the output container
if
(
nodes_buffer
.
size
()
>
1
)
outputDD
(
mesh
,
nodes_buffer
);
/* ******************************************************************** */
++
Nlines
;
// increase the number of detected dislocation lines by 1
}
/* ********************************************************************** */
else
{
fhR_l
.
remove
(
*
fh_it
);
// remove the current cell from the list
}
};
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
DUMP
(
"Detected dislocation lines: "
<<
Nlines
,
DBG_MESSAGE
);
// Test version --> throw error if more than one dislocation line is detected
if
(
Nlines
>
1
)
LM_FATAL
(
std
::
endl
<<
"Output of several dislocation lines is currently prohibited!"
<<
std
::
endl
);
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
// Track dislocation loops
while
(
fhI_l
.
size
()
>
0
)
{
fh_it
=
fhI_l
.
begin
();
// set iterator to the first element of the list
// Get the cell and the index of the opposed vertex of the current face
ch
=
(
*
fh_it
).
first
;
idx
=
(
*
fh_it
).
second
;
/* ********************************************************************** */
// If the Burgers vector is non-zero ==> mark both cells containing the
// current face
// --> the neighboring cell can simply be accessed via the index of the
// opposed vertex
if
(
this
->
computeBurgersVec
(
ch
,
idx
).
squared_length
()
>
__TOL__
)
{
LM_FATAL
(
std
::
endl
<<
"Dislocation loop detected."
<<
std
::
endl
<<
"Dislocation loops not supported yet!"
<<
std
::
endl
);
}
/* ********************************************************************** */
else
{
fhI_l
.
remove
(
*
fh_it
);
// remove the current cell from the list
}
};
/* ------------------------------------------------------------------------ */
DUMP
(
"Detected dislocation loops: "
<<
Nloops
,
DBG_MESSAGE
);
/* ------------------------------------------------------------------------ */
}
/* -------------------------------------------------------------------------- */
/* Outputs the DD mesh to the output container */
// NOTE: the point set is divided into sets of equal size. If 'Nnodes' is not
// divisible by 'nodesPerSeg' resdiual segments ared added with
// 'nodesPerSeg-1' nodes per segment
template
<
typename
Mesh
>
void
ComputeDisloDetection
::
outputDD
(
Mesh
&
mesh
,
std
::
vector
<
RefPointData
<
3
>>
&
nodes_buffer
)
{
UInt
Nnodes
,
nodes_buffer_inc
;
UInt
Nsegs
,
Nsegs_res
;
Vector
<
2
,
UInt
>
conn
;
// Containers/Lists
typename
std
::
vector
<
RefPointData
<
3
>>::
iterator
nodes_it
;
// LM
RefGenericElem
<
3
>
el_buffer
;
// LM - Containers/Lists
auto
&
nodes
=
mesh
.
getContainerNodes
();
auto
&
elems
=
mesh
.
getContainerElems
();
/* ------------------------------------------------------------------------ */
// Set initial connectivity
conn
[
0
]
=
nodes
.
size
()
-
1
;
conn
[
1
]
=
nodes
.
size
();
/* ------------------------------------------------------------------------ */
/* 1) Get the number of regular and residual segments */
Nnodes
=
nodes_buffer
.
size
()
-
1
;
// Set increment
if
(
this
->
nodesPerSeg
>
nodes_buffer
.
size
())
nodes_buffer_inc
=
Nnodes
;
else
nodes_buffer_inc
=
this
->
nodesPerSeg
-
1
;
Nsegs
=
UInt
(
Nnodes
/
nodes_buffer_inc
);
if
(
Nnodes
-
Nsegs
*
nodes_buffer_inc
>
0
)
{
++
Nsegs
;
Nsegs_res
=
nodes_buffer_inc
*
Nsegs
-
Nnodes
;
// If the number of residual segments is larger than the number of regular
// segments ...
if
(
Nsegs_res
>
Nsegs
)
{
nodes_buffer_inc
=
UInt
((
Nsegs
+
Nnodes
)
/
Nsegs
);
// ... choose an optimal increment
// (w.r.t. to the (largest possible)
// nbr of nodes per segment) s.t. the
// Nsegs == const.
// Update the number of resdiual segments
Nsegs_res
=
nodes_buffer_inc
*
Nsegs
-
Nnodes
;
}
}
else
Nsegs_res
=
0
;
/* ------------------------------------------------------------------------ */
/* 2) Add the nodes to the output container */
nodes_it
=
nodes_buffer
.
begin
();
nodes
.
push_back
(
*
nodes_it
);
// Add "regular" nodes
for
(
UInt
i
=
0
;
i
<
Nsegs
-
Nsegs_res
;
++
i
)
{
std
::
advance
(
nodes_it
,
nodes_buffer_inc
);
nodes
.
push_back
(
*
nodes_it
);
}
--
nodes_buffer_inc
;
// Add "residual" nodes
for
(
UInt
i
=
0
;
i
<
Nsegs_res
;
++
i
)
{
std
::
advance
(
nodes_it
,
nodes_buffer_inc
);
nodes
.
push_back
(
*
nodes_it
);
}
/* ------------------------------------------------------------------------ */
/* 3) Add the connectivity */
for
(
UInt
i
=
0
;
i
<
Nsegs
;
++
i
)
{
// Update connectivity
++
conn
[
0
];
++
conn
[
1
];
// Add connectivity to output container
el_buffer
.
setConnectivity
(
conn
);
elems
.
push_back
(
el_buffer
);
}
/* ------------------------------------------------------------------------ */
// Clear array of buffer nodes
nodes_buffer
.
clear
();
/* ------------------------------------------------------------------------ */
}
/* -------------------------------------------------------------------------- */
/* Outputs the tetrahedralization to an additional container */
void
ComputeDisloDetection
::
outputMesh
(
std
::
list
<
Vertex_handle
>
&
v_l
,
std
::
list
<
Cell_handle
>
&
c_l
)
{
UInt
ctr
=
0
;
Vector
<
4
,
UInt
>
conn
;
RefPointData
<
3
>
pt_buffer
;
RefGenericElem
<
3
>
el_buffer
;
// Containers/Lists
auto
&
nodes
=
this
->
meshTet
.
getContainerNodes
();
auto
&
elems
=
this
->
meshTet
.
getContainerElems
();
// Iterators
std
::
list
<
Vertex_handle
>::
iterator
v_it
;
std
::
list
<
Cell_handle
>::
iterator
c_it
;
/* ------------------------------------------------------------------------ */
// Add nodes
nodes
.
clear
();
for
(
v_it
=
v_l
.
begin
();
v_it
!=
v_l
.
end
();
++
v_it
)
{
// Get positions
pt_buffer
.
position
()[
0
]
=
(
*
v_it
)
->
point
().
x
();
pt_buffer
.
position
()[
1
]
=
(
*
v_it
)
->
point
().
y
();
pt_buffer
.
position
()[
2
]
=
(
*
v_it
)
->
point
().
z
();
nodes
.
push_back
(
pt_buffer
);
// Add corresponding index to info()
(
*
v_it
)
->
info
()
=
ctr
;
++
ctr
;
}
/* ------------------------------------------------------------------------ */
// Add elements
elems
.
clear
();
for
(
c_it
=
c_l
.
begin
();
c_it
!=
c_l
.
end
();
++
c_it
)
{
// Get connectivity via vertex info
conn
[
0
]
=
(
*
c_it
)
->
vertex
(
0
)
->
info
();
conn
[
1
]
=
(
*
c_it
)
->
vertex
(
1
)
->
info
();
conn
[
2
]
=
(
*
c_it
)
->
vertex
(
2
)
->
info
();
conn
[
3
]
=
(
*
c_it
)
->
vertex
(
3
)
->
info
();
el_buffer
.
setConnectivity
(
conn
);
elems
.
push_back
(
el_buffer
);
}
}
/* -------------------------------------------------------------------------- */
/* Collect container arrays of Reals by the root process */
// NOTE: this is a quick fix to gather all data on the root process since the
// dislocation detection is not parallelized yet!
// NOTE: the function was is mostly taken from 'compute.cc' and should be
// generalized since communications methods are confined to the compute
// interface...
void
ComputeDisloDetection
::
gatherAllData
(
ContainerArray
<
Real
>
&
data_recv
,
ContainerArray
<
Real
>
&
data_send
,
const
UInt
root_rank
)
{
CommGroup
group
=
this
->
getCommGroup
();
UInt
my_rank
=
group
.
getMyRank
();
UInt
nb_data
=
data_send
.
size
();
std
::
vector
<
UInt
>
nb_data_per_proc
;
// #ifndef LM_OPTIMIZED
// UInt total_data = this->getTotalNbData(root_rank);
// comm.synchronize(group);
// #endif //LM_OPTIMIZED
// gather_flag = true;
// gather_root_proc = root_rank;
if
(
root_rank
==
my_rank
)
{
// prepare the array to resizes the sizes
UInt
nb_procs
=
group
.
size
();
DUMP
(
"receive "
<<
nb_procs
<<
" procs"
,
DBG_INFO
);
DUMP
(
"my rank is "
<<
my_rank
,
DBG_INFO
);
nb_data_per_proc
.
resize
(
nb_procs
);
for
(
UInt
p
=
0
;
p
<
nb_procs
;
++
p
)
{
if
(
p
==
my_rank
)
nb_data_per_proc
[
p
]
=
nb_data
;
else
{
DUMP
(
"receive from proc "
<<
p
,
DBG_INFO
);
group
.
receive
(
&
nb_data_per_proc
[
p
],
1
,
p
,
"gatherData: receive number"
);
}
}
// compute total size of the gathered data
UInt
tot_size
=
0
;
for
(
UInt
p
=
0
;
p
<
nb_procs
;
++
p
)
{
DUMP
(
"nb_data_per_proc["
<<
p
<<
"]="
<<
nb_data_per_proc
[
p
],
DBG_INFO
);
tot_size
+=
nb_data_per_proc
[
p
];
}
// #ifndef LM_OPTIMIZED
// LM_ASSERT(total_data == tot_size, "mismatched of global sizes");
// #endif //LM_OPTIMIZED
// resize the receiving buffer
data_recv
.
resize
(
tot_size
);
// receive the data from the other processors
UInt
offset
=
0
;
for
(
UInt
p
=
0
;
p
<
nb_procs
;
++
p
)
{
// if p is my_rank/root_rank copy the data
if
(
p
==
my_rank
)
{
for
(
UInt
i
=
0
;
i
<
nb_data
;
++
i
)
{
data_recv
[
offset
+
i
]
=
(
data_send
)[
i
];
}
}
// else receive from distant proc
else
if
(
nb_data_per_proc
[
p
])
{
group
.
receive
(
&
data_recv
[
offset
],
nb_data_per_proc
[
p
],
p
,
"gatherData: receive data"
);
}
// increment the offset
offset
+=
nb_data_per_proc
[
p
];
}
}
else
{
DUMP
(
"send my nb_data = "
<<
nb_data
,
DBG_INFO
);
// send the amount of data I have
group
.
send
(
&
nb_data
,
1
,
root_rank
,
"gatherData: send number"
);
// if there is data to be sent : do so
if
(
nb_data
!=
0
)
group
.
send
(
data_send
,
root_rank
,
"gatherData: send data"
);
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC DISLODETECTION
Dislocation detection for single crystals based on Stukowski (2014)
*/
/* LMEXAMPLE
COMPUTE dxa DISLOCDETN INPUT Pa ALPHA 7.5
*/
inline
void
ComputeDisloDetection
::
declareParams
()
{
/* LMKEYWORD LATTICE_T
Lattice type (supported: \textit{fcc} - default: \textit{fcc})
*/
this
->
parseKeyword
(
"LATTICE_T"
,
this
->
lattice_t
,
"fcc"
);
/* LMKEYWORD A0
Lattice constant
*/
this
->
parseKeyword
(
"A0"
,
this
->
a0
);
/* LMKEYWORD ALPHA
Maximum squared circumradius for tetrahedrons (defines the alpha shape)
*/
this
->
parseKeyword
(
"ALPHA"
,
this
->
alpha
);
/* LMKEYWORD ORIENT
Orientation of the crystal (default: x [1 2 1], y [-1 0 1], z [1 -1 1])
*/
this
->
parseVectorKeyword
(
"ORIENT"
,
9
,
this
->
orient
,
VEC_DEFAULTS
(
1.
,
2.
,
1.
,
-
1.
,
0.
,
1.
,
1.
,
-
1.
,
1.
));
/* LMKEYWORD DISLOC_T
Dislocation type (supported: \textit{full} - default: \textit{full})
*/
this
->
parseKeyword
(
"DISLOC_T"
,
this
->
disloc_t
,
"full"
);
/* LMKEYWORD NODES_PER_SEG
Specify a number of nodes per segment to avoid artificial small segments
(default is 1).
{bf NOTE:} fix should be replaced with a curve-fitting procedure in future
versions
*/
this
->
parseKeyword
(
"NODES_PER_SEG"
,
this
->
nodesPerSeg
,
1u
);
}
/* -------------------------------------------------------------------------- */
/* BACKWARD IMPLEMENTATION */
/* -------------------------------------------------------------------------- */
template
<
typename
OutputIterator
>
OutputIterator
get_all_alpha_shape_edges
(
Ta
&
tet_alpha
,
OutputIterator
it
)
{
T
::
Finite_edges_iterator
eh_it
=
tet_alpha
.
finite_edges_begin
();
// use tet_alpha since tet will be destroyed!!
for
(;
eh_it
!=
tet_alpha
.
finite_edges_end
();
++
eh_it
)
{
if
((
tet_alpha
.
classify
(
*
eh_it
)
==
Ta
::
INTERIOR
)
||
(
tet_alpha
.
classify
(
*
eh_it
)
==
Ta
::
REGULAR
))
*
it
++
=
*
eh_it
;
}
return
it
;
}
template
<
typename
OutputIterator
>
OutputIterator
get_all_alpha_shape_faces
(
Ta
&
tet_alpha
,
OutputIterator
it
)
{
T
::
Finite_facets_iterator
fh_it
=
tet_alpha
.
finite_facets_begin
();
// use tet_alpha since tet will be destroyed!!
for
(;
fh_it
!=
tet_alpha
.
finite_facets_end
();
++
fh_it
)
{
if
((
tet_alpha
.
classify
(
*
fh_it
)
==
Ta
::
INTERIOR
)
||
(
tet_alpha
.
classify
(
*
fh_it
)
==
Ta
::
REGULAR
))
*
it
++
=
*
fh_it
;
}
return
it
;
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_MAKE_CALL
(
ComputeDisloDetection
)
/* -------------------------------------------------------------------------- */
#undef __TOL__
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
Event Timeline
Log In to Comment