diff --git a/src/filter/compute_dislodetection.cc b/src/filter/compute_dislodetection.cc
index 225f9dc..7136b0d 100644
--- a/src/filter/compute_dislodetection.cc
+++ b/src/filter/compute_dislodetection.cc
@@ -1,992 +1,992 @@
/**
* @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_dislocdetn.hh"
+#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 = allocAndGetCastedOutput<ContainerGenericMesh<3>>(
"detected_segments:" + this->getID());
/* ------------------------------------------------------------------------ */
// 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);
std::vector<Real> pos_buffer;
this->gatherAllData(pos_buffer, positions, 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...
{
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.add(*nodes_it);
// Add "regular" nodes
for (UInt i = 0; i < Nsegs - Nsegs_res; ++i) {
std::advance(nodes_it, nodes_buffer_inc);
nodes.add(*nodes_it);
}
--nodes_buffer_inc;
// Add "residual" nodes
for (UInt i = 0; i < Nsegs_res; ++i) {
std::advance(nodes_it, nodes_buffer_inc);
nodes.add(*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.add(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.add(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.add(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(std::vector<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(static_cast<std::vector<Real> &>(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__