diff --git a/src/coupling/point_association.cc b/src/coupling/point_association.cc
index 300e9e1..21c1ef9 100644
--- a/src/coupling/point_association.cc
+++ b/src/coupling/point_association.cc
@@ -1,588 +1,588 @@
/**
* @file point_association.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Thu Jul 24 14:21:58 2014
*
* @brief Coupling tool to associate point-like DOFs
*
* @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/>.
*
*/
//#define TIMER
/* -------------------------------------------------------------------------- */
#include "point_association.hh"
#include "filter_geometry.hh"
#include "geometry_manager.hh"
#include "lib_bridging.hh"
#include "lib_continuum.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "reference_manager_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
PointAssociation<ContainerA, ContainerB, Dim>::PointAssociation(
const std::string &name, ContainerA &contA, ContainerB &contB)
: LMObject(name), DofAssociation<ContainerA, ContainerB, Dim>(name, contA,
contB),
pointsAlist("pointsAlist:" + name), pointsBlist("pointsBlist:" + name),
unmatchedA("unmatchedA:" + name), unmatchedB("unmatchedB:" + name),
assoc_found(0), tolerance(1e-5), path(""), weight_path(""),
weights("weights:" + this->getID()),
weightsA("weightsA:" + this->getID()),
weightsB("weightsB:" + this->getID()),
true_disp("ComputeTrueDisplacement:" + this->getID()) {
if (this->in_group_A && this->contA.hasRefManager())
this->contA.getRefManager()->addSubSet(this->getID() + ":pointsA",
*this->pointsAlist);
if (this->in_group_B && this->contB.hasRefManager())
this->contB.getRefManager()->addSubSet(this->getID() + ":pointsB",
*this->pointsBlist);
/* registering computes for outer world */
LM_TOIMPLEMENT;
// FilterManager::getManager().addObject(&this->pointsAlist, false);
// FilterManager::getManager().addObject(&this->pointsBlist, false);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
PointAssociation<ContainerA, ContainerB, Dim>::~PointAssociation() {
if (this->grid) {
delete this->grid;
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::init() {
if (this->contA.hasRefManager() && this->contB.hasRefManager())
LM_FATAL("Not sure what to do if both sides are able to migrate particles");
// Build point lists for both domains
if (this->in_group_A) {
this->buildPointList(this->contA, this->unmatchedA);
this->local_points = this->unmatchedA->size();
this->buildPositions(this->unmatchedA.getOutput());
}
if (this->in_group_B)
this->buildPointList(this->contB, this->unmatchedB);
DUMP(this->getID() << " : Generating Communication Scheme", DBG_INFO_STARTUP);
MPI_Barrier(MPI_COMM_WORLD);
DUMP(this->getID() << " : exchange coarse geometries", DBG_DETAIL);
this->exchangeGeometries(this->unmatchedA, this->unmatchedB);
MPI_Barrier(MPI_COMM_WORLD);
DUMP(this->getID() << " : exchange points positions", DBG_DETAIL);
this->exchangePositions();
MPI_Barrier(MPI_COMM_WORLD);
if (this->in_group_B) {
// set the size and assign zero (for help in debugging)
this->field_buffer.assign(Dim * (this->pointsAlist->size()), 0);
this->fillGrid();
this->associate();
}
MPI_Barrier(MPI_COMM_WORLD);
DUMP(this->getID() << " : ExchangeAssociation", DBG_DETAIL);
std::vector<UInt> managed;
this->exchangeAssociationInformation(managed);
if (!this->in_group_A || !this->in_group_B) {
if (this->in_group_A) {
LM_ASSERT(managed.size() == this->local_points * this->nb_zone_B,
"internal error " << managed.size() << " != "
<< this->local_points * this->nb_zone_B);
std::vector<std::vector<UInt>> unassociated;
this->createDuoVectorA("A", managed, unassociated);
for (UInt i = 0; i < this->nb_zone_B; ++i) {
LM_ASSERT(unassociated[i].size() == 0,
"internal error : "
<< " for proc " << i << " " << unassociated[i].size()
<< " points were already associated");
}
}
if (this->in_group_B)
this->createDuoVectorB("B");
}
if (this->in_group_A) {
this->filterPointListAForUnmatched();
LM_ASSERT(this->unmatchedA->size() == this->pointsAlist->size(),
"internal error " << this->unmatchedA->size()
<< " != " << this->pointsAlist->size());
LM_TOIMPLEMENT;
UInt total_points_from_file;
// UInt total_points_from_file = this->unmatchedA->getNbPointsInFile();
UInt total_points = this->local_points;
this->comm_group_A.allReduce(&total_points, 1, "verif association", OP_SUM);
if (total_points != total_points_from_file) {
if (total_points > total_points_from_file) {
LM_FATAL("This should never happen, not even with a faulty mesh. If "
<< "you get this message, you're looking at a serious bug.");
}
// in this case, need to figure out which points are missing
LM_FATAL("problem in A association " << total_points
<< " != " << total_points_from_file);
}
}
if (this->in_group_B) {
this->filterPointListBForUnmatched();
this->local_points = this->pointsBlist->size();
this->filterAssoc();
LM_ASSERT(this->unmatchedB->size() == this->pointsBlist->size(),
"internal error");
#ifndef LM_OPTIMIZED
LM_TOIMPLEMENT;
UInt total_points_from_file;
// UInt total_points_from_file = this->unmatchedB.getNbPointsInFile();
UInt total_points = this->local_points;
this->comm_group_B.allReduce(&total_points, 1, "verif association", OP_SUM);
LM_ASSERT(total_points == total_points_from_file,
"problem in B association " << total_points
<< " != " << total_points_from_file);
#endif
}
checkPositionCoherency();
if (this->in_group_A && this->contA.hasRefManager()) {
try {
this->contA.getRefManager()->attachObject(this->getDuoVector(),
*this->pointsAlist);
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
} catch (LibMultiScaleException &e) {
}
}
if (this->in_group_B && this->contB.hasRefManager()) {
try {
this->contB.getRefManager()->attachObject(this->getDuoVector(),
*this->pointsBlist);
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
} catch (LibMultiScaleException &e) {
}
}
// DUMP(this->getID() << ": bounding box of " <<
// this->pointsAlist.getID() << " is "
// << this->pointsAlist.getBoundingBox(),DBG_MESSAGE);
// DUMP(this->getID() << ": bounding box of " <<
// this->pointsBlist.getID() << " is "
// << this->pointsBlist.getBoundingBox(),DBG_MESSAGE);
// generate the weights
if (weight_path != "") {
this->computeWeights();
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::computeWeights() {
ComputeExtract c_positions("ComputeExtract:" + this->getID());
c_positions.setParam("FIELD", _position0);
c_positions.buildManual(pointsAlist);
this->weights.setParam("FILENAME", this->weight_path);
this->weights.init();
this->weights.buildManual(c_positions);
if (this->in_group_A && !this->in_group_B)
LM_FATAL("not a parallel routine yet");
if (this->pointsAlist->size() != this->pointsBlist->size())
LM_FATAL("inconsistency here");
std::vector<Real> denominator;
// weightsA.resize(this->pointsAlist.nbElem());
// weightsB.resize(this->pointsAlist.nbElem());
weightsA.assign(this->pointsAlist->size(), 0.);
weightsB.assign(this->pointsAlist->size(), 0.);
// compute nominator and denominator
UInt index = 0;
for (auto &&point : this->pointsAlist) {
Real mass = point.mass();
weightsA[index] = weights[index] * mass;
++index;
}
index = 0;
for (auto &&point : this->pointsBlist) {
Real mass = point.mass();
weightsB[index] = (1. - weights[index]) * mass;
++index;
}
for (UInt index = 0; index < weights.size(); ++index) {
Real denom = weightsA[index] + weightsB[index];
weightsA[index] /= denom;
weightsB[index] /= denom;
}
LM_TOIMPLEMENT;
// FilterManager::getManager().addObject(&this->weights, false);
// FilterManager::getManager().addObject(&this->weightsA, false);
// FilterManager::getManager().addObject(&this->weightsB, false);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
template <typename container, typename PointList>
void PointAssociation<ContainerA, ContainerB, Dim>::buildPointList(
container &cont, PointList &pointList) {
casted_component<FilterGeometry>::cast<typename container::ContainerSubset>
filter_geom("FilterGeom:" + this->getID());
filter_geom.setParam("GEOMETRY", this->geomID);
filter_geom.setParam("DONOTRESETBOUNDINGBOX", true);
filter_geom.setParam("DO_NOT_REGISTER", true);
filter_geom.setParam("BASED_ON_NODES", true);
filter_geom.setParam("DOF_TYPE", dt_local_master);
filter_geom.init();
filter_geom.buildManual(cont);
DUMP("constructed filtergeom with " << filter_geom->size() << " points ",
DBG_INFO_STARTUP);
pointList.setParam("PATH", this->path);
pointList.setParam("TOL", this->tolerance);
pointList.init();
pointList.buildManual(filter_geom);
DUMP("constructed filter_point with " << pointList->size() << " points ",
DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::fillGrid() {
// create my grid
if (this->grid) {
this->deleteGrid();
}
Geometry &geom = *GeometryManager::getManager().getGeometry(this->geomID);
Cube cube = geom.getBoundingBox();
this->grid = new SpatialGridPoint<UInt, Dim>(cube, this->grid_division);
UInt index = 0;
for (auto &&point : this->unmatchedB) {
auto position = point.position0();
// I'ts important to use addAtom here instead of addElement, because here
// we're sensitive to atoms on the edge of grid_boxes
this->grid->addAtom(index, position);
++index;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename V1, typename V2>
inline bool coincidence(V1 &pointA, V2 &pointB, const Real &tolerance) {
Real distance = 0;
for (UInt dir = 0; dir < Dim; ++dir) {
distance += fabs(pointA[dir] - pointB[dir]);
}
return distance < tolerance;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename V1, typename V2>
inline Real temp_coincidence(V1 &pointA, V2 &pointB) {
Real distance = 0;
for (UInt dir = 0; dir < Dim; ++dir) {
distance += fabs(pointA[dir] - pointB[dir]);
}
return distance;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::associate() {
this->assoc.assign(this->local_points, lm_uint_max);
LM_ASSERT(this->positions.size() == this->local_points * Dim,
"inconsistency " << this->positions.size()
<< " != " << this->local_points * Dim);
// Loop over all positions of A world
bool found;
Real min_dist;
UInt not_found_counter = 0;
Real *posA_ptr = &this->positions[0];
for (UInt indexA = 0; indexA < this->local_points;
++indexA, posA_ptr += Dim) {
// extract the points of pointsBlist which are candidates for association
VectorView<Dim> positionA(posA_ptr);
LM_TOIMPLEMENT;
// auto &&subsetPointsB = this->grid->findSet(positionA);
found = false;
min_dist = lm_real_max;
// loop the candidates
LM_TOIMPLEMENT;
// for (auto &&indexB : subsetPointsB) {
// auto &&positionB = unmatchedB.get(indexB).position0();
// Real dist = temp_coincidence<Dim>(positionA, positionB);
// if (dist < min_dist)
// min_dist = dist;
// if (coincidence<Dim>(positionA, positionB, this->tolerance)) {
// this->assoc[indexA] = indexB;
// ++this->assoc_found;
// found = true;
// if (indexA < 10)
// DUMP("associate point " << positionA << ", " << indexB
// << " with position " << indexA,
// DBG_DETAIL);
// break;
// }
// }
if (!found) {
++not_found_counter;
{
std::stringstream sstr;
sstr << not_found_counter << ": couldn't associate point (";
for (UInt i = 0; i < Dim - 1; ++i)
sstr << positionA[i] << ", ";
sstr << positionA[Dim - 1] << "), "
<< "the minimal distance was " << min_dist;
DUMP(sstr.str(), DBG_INFO);
}
}
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::clearAll() {
this->assoc_found = 0;
this->assoc.clear();
delete this->grid;
pointsAlist->clear();
pointsBlist->clear();
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB,
Dim>::filterPointListAForUnmatched() {
this->pointsAlist->clear();
for (UInt i = 0; i < this->assoc.size(); ++i) {
if (this->assoc[i] != UINT_MAX) {
DUMP("for container, atom " << i << " becomes "
<< this->pointsAlist->size(),
DBG_ALL);
- this->pointsAlist->add(this->unmatchedA->get(i));
+ this->pointsAlist->push_back(this->unmatchedA->get(i));
} else {
DUMP("atom filtered " << i << " " << this->unmatchedA->get(i),
DBG_DETAIL);
}
}
this->pointsAlist->copyContainerInfo(unmatchedA);
LM_TOIMPLEMENT;
// this->pointsAlist->copyReleaseInfo(unmatchedA);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB,
Dim>::filterPointListBForUnmatched() {
this->pointsBlist->clear();
UInt counter = 0;
for (UInt i = 0; i < this->assoc.size(); ++i) {
if (this->assoc[i] != UINT_MAX) {
DUMP("for container, atom " << i << " becomes "
<< this->pointsBlist->size(),
DBG_ALL);
- this->pointsBlist->add(unmatchedB->get(this->assoc[i]));
+ this->pointsBlist->push_back(unmatchedB->get(this->assoc[i]));
this->assoc[i] = counter;
++counter;
} else {
DUMP("atom filtered " << i, DBG_DETAIL);
}
}
this->pointsBlist->getBoundingBox() = unmatchedB->getBoundingBox();
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::checkPositionCoherency() {
// try {
// std::stringstream sstr;
// sstr << this->getID() << "-duovector-" << current_step << "-" <<
// current_stage;
// this->getDuoVector().print(sstr.str());
// }
// catch (LibMultiScaleException & e){
// }
// extract the position of local atoms
if (this->in_group_A)
this->buildPositions(this->pointsAlist.getOutput());
// transport the positions to B side
this->distributeVectorA2B("temp_positions", this->positions, Dim);
if (this->in_group_B) {
UInt point_idx = 0;
for (auto &&point : this->pointsBlist) {
auto &&local_pos = point.position0();
std::stringstream sstr;
for (UInt i = 0; i < Dim; ++i)
sstr << this->positions[Dim * point_idx + i] << " ";
sstr << " and ";
for (UInt i = 0; i < Dim; ++i)
sstr << local_pos[i] << " ";
std::string _tmp = sstr.str();
DUMP("for point " << point_idx << " compare " << _tmp, DBG_DETAIL);
for (UInt k = 0; k < Dim; ++k) {
if (fabs(this->positions[Dim * point_idx + k] - local_pos[k]) > 1e-5) {
LM_FATAL(this->getID()
<< ": inconsistency in the communication scheme "
<< this->positions[Dim * point_idx + k]
<< " != " << local_pos[k] << " " << point_idx);
}
}
++point_idx;
}
}
DUMP("position coherency OK", DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::updateForMigration() {
DUMP("update migration for " << this->getID(), DBG_INFO);
DUMP("update local_points " << this->getID(), DBG_INFO);
STARTTIMER("syncMigration resize");
// update local number of points
if (this->in_group_A)
this->local_points = this->pointsAlist->size();
// update local number of points
if (this->in_group_B)
this->local_points = this->pointsBlist->size();
this->positions.resize(Dim * this->local_points);
STOPTIMER("syncMigration resize");
DUMP("synchMigration " << this->getID(), DBG_INFO);
STARTTIMER("syncMigration duo");
if ((this->in_group_A && this->contA.hasRefManager()) ||
(this->in_group_B && !this->contB.hasRefManager()))
this->synchronizeMigration(this->comm_group_A, this->comm_group_B);
if ((this->in_group_B && this->contB.hasRefManager()) ||
(this->in_group_A && !this->contA.hasRefManager()))
this->synchronizeMigration(this->comm_group_B, this->comm_group_A);
STOPTIMER("syncMigration duo");
DUMP("update release " << this->getID(), DBG_INFO);
LM_TOIMPLEMENT; // to review
// pointsAlist.setRelease(this->contA.getRelease());
// pointsBlist.setRelease(this->contB.getRelease());
if (check_coherency) {
DUMP("check position coherency " << this->getID(), DBG_INFO);
STARTTIMER("check position coherency");
this->checkPositionCoherency();
STOPTIMER("check position coherency");
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC PointAssociation
This coupler computes the association between point degrees of freedom of
different domains
*/
/* LMHERITATE ComputeTrueDisplacement
*/
template <typename ContainerA, typename ContainerB, UInt Dim>
void PointAssociation<ContainerA, ContainerB, Dim>::declareParams() {
DofAssociation<ContainerA, ContainerB, Dim>::declareParams();
this->addSubParsableObject(true_disp);
/* LMKEYWORD TOL
max distance between two points to consider them coincident
*/
this->parseKeyword("TOL", tolerance);
/* LMKEYWORD PATH
file path to the filteratom file
*/
this->parseKeyword("PATH", path);
/* LMKEYWORD PBC
Specify if periodic boundary conditions should be considered for each
direction
*/
this->parseVectorKeyword("PBC", spatial_dimension, this->pbc);
/* LMKEYWORD CHECK_COHERENCY
Perform a systematic check of the communication scheme.
\textbf{Be careful, it is extremely computationally expensive}
*/
this->parseKeyword("CHECK_COHERENCY", check_coherency, false);
/* LMKEYWORD WEIGHT_PATH
file path to the python script to generate the weights
*/
this->parseKeyword("WEIGHT_PATH", weight_path);
}
/* -------------------------------------------------------------------------- */
DECLARE_BRIDGING_ATOMIC_CONTINUUM_TEMPLATE(PointAssociation)
DECLARE_BRIDGING_CONTINUUM_ATOMIC_TEMPLATE(PointAssociation)
__END_LIBMULTISCALE__
diff --git a/src/coupling/template_cadd.cc b/src/coupling/template_cadd.cc
index b0b29ec..3babde5 100644
--- a/src/coupling/template_cadd.cc
+++ b/src/coupling/template_cadd.cc
@@ -1,2075 +1,2075 @@
/**
* @file template_cadd.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Thu Sep 18 14:51:07 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/>.
*
*/
#define TIMER
/* -------------------------------------------------------------------------- */
#include "template_cadd.hh"
#include "lib_bridging.hh"
#include "lm_common.hh"
#include <algorithm>
#include <assert.h>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <unistd.h>
#include <utility>
#include <vector>
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "compute_extract.hh"
#include "container_mesh.hh"
#include "geometry_manager.hh"
#include "reference_manager_interface.hh"
#include "communicator.hh"
#include "domain_dd_interface.hh"
#include "ref_point_data.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::mergeDuplicateNodes() {
for (auto &&nd : this->domA.getContainer().getContainerNodes()) {
LM_TOIMPLEMENT;
// UInt nbrs = nd.getNbrOfNeighs();
UInt nbrs = std::numeric_limits<UInt>::max();
UInt idx = nd.getIndex();
if (this->md_cube->contains(nd.position()) && nbrs == 2) {
UInt idx2 = nd.getIndex();
for (auto &&nd2[[gnu::unused]] : this->domA.getContainer().getContainerNodes()) {
LM_TOIMPLEMENT;
//auto dX = nd.position() - nd2.position();
// Real dist = dX.norm();
auto dist = std::numeric_limits<Real>::max();
if ((idx != idx2) && (dist < 1e-2)) {
this->domA.freenodearm(idx);
this->domA.freenode(idx);
this->domA.freenodearm(idx2);
this->domA.freenode(idx2);
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::init() {
// auto &comm = Communicator::getCommunicator();
// fabricating the geometries
// IteratorDDNodes nodes = this->domA.getContainerNodes().getIterator();
auto md_geom = GeometryManager::getManager().getGeometry(this->md_geom);
this->md_cube = std::dynamic_pointer_cast<Cube>(md_geom);
auto dd_geom = GeometryManager::getManager().getGeometry(this->dd_geom);
this->dd_cube = std::dynamic_pointer_cast<Cube>(dd_geom);
Vector<3u> xmin(this->md_cube->getXmin());
Vector<3u> xmax(this->md_cube->getXmax());
for (UInt i : {0, 1, 2}) {
xmin[i] -= this->pad_thickness + this->distance_ddline_from_pad;
xmax[i] += this->pad_thickness + this->distance_ddline_from_pad;
}
this->md_cube->setXmin(xmin);
this->md_cube->setXmax(xmax);
// constructing some flags
this->detect_simple_bool = false;
simple_direction = 3;
for (auto &&i : {0, 1, 2}) {
if (this->close_direction[i] != 0)
this->close_in_infinity = true;
if (this->detect_simple_direction[i] != 0) {
this->detect_simple_bool = true;
simple_direction = i;
}
}
if (this->aniso) {
Real c11 = this->elastic_coefs[0];
Real c12 = this->elastic_coefs[1];
Real c22 = this->elastic_coefs[4];
Real c66 = this->elastic_coefs[12];
Real lambda2 = sqrt(c11 / c22);
this->c11bar = sqrt(c11 * c12);
this->lambda = sqrt(lambda2);
this->phi = 0.5 * acos((c12 * c12 + 2.0 * c12 * c66 - c11bar * c11bar) /
(2.0 * c11bar * c66));
if (this->volterra)
LM_FATAL("cannot be used at the same time with volterra and aniso");
}
// sending infos on parameters to paradis
Real maxSeg, burgMag;
if (this->is_in_dd) {
Vector<Dim> domA_xmin = this->domA.getXmin();
Vector<Dim> domA_xmax = this->domA.getXmax();
if ((this->dd_cube->getXmin() - domA_xmin).norm() > 1e-7 ||
(this->dd_cube->getXmax() - domA_xmax).norm() > 1e-7)
DUMP("[WARNING] Check DD box lengths defined Libmultiscale configuration "
"and ParaDis control files",
DBG_MESSAGE);
this->maxSeg0 = this->domA.getMaxSeg();
this->minSeg0 = this->domA.getMinSeg();
this->burgMag0 =
this->domA.getBurgersMagnitude() * 1e+10; // in angstrom unit
int my_rank_in_DDprocs = this->comm_group_dd.getMyRank();
if (my_rank_in_DDprocs == 0) {
UInt mdNBprocs = this->comm_group_atomic.size();
for (UInt i = 0; i < mdNBprocs; ++i) {
this->comm_group_atomic.send(&maxSeg, 1, i, "maximum length of dd");
this->comm_group_atomic.send(&burgMag, 1, i,
"burgers vector magnitude");
}
}
}
// receive information on parameters to paradis
if (this->is_in_atomic) {
this->comm_group_dd.receive(&maxSeg, 1, 0, "maximum length of dd segments");
this->comm_group_dd.receive(&burgMag, 1, 0, "burgers vector magnitude");
this->maxSeg0 = maxSeg;
this->burgMag0 = burgMag;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::coupling(CouplingStage stage) {
//! Main coupling calls
if (stage == COUPLING_STEP3) {
STARTTIMER("applyDetection");
this->applyDetection();
STOPTIMER("applyDetection");
STARTTIMER("applyTemplateToPad");
this->applyTemplateToPad();
STOPTIMER("applyTemplateToPad");
return;
}
}
/* --------------------------------------------------------------------------
applyDetection does following jobs in two kinds of processors (MD and DD)
MD: 1. compute MD dislocation coordinates (DXA or user-defined compute)
2. send to DD processor
DD: 1. (at every step) constrain DD nodes located in MD domain
3. receive the MD dislocation coordinates from MD processor
4. move DD nodes to the closest MD segements
5. register the remained MD cores as new DD nodes
6. check topology - remesh
-------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::applyDetection() {
UInt num_dislo = 0;
std::vector<Vector<Dim>> dislocations;
// auto &comm = Communicator::getCommunicator();
if (this->is_in_atomic) {
if (current_step % this->detect_freq == 0) {
UInt my_rank_in_MD = this->comm_group_atomic.getMyRank();
STARTTIMER("detect_cores_by_compute");
this->detect_cores_by_compute(num_dislo, dislocations);
STOPTIMER("detect_cores_by_compute");
if (my_rank_in_MD == 0) {
UInt nbDDprocs = this->comm_group_dd.size();
for (UInt i = 0; i < nbDDprocs; ++i) {
// int real_rank = this->comm_group_dd.realRank(i);
// int rank_in_group = this->comm_group_dd.groupRank(real_rank);
this->comm_group_dd.send(&num_dislo, 1, i,
"number of dislocation cores");
this->comm_group_dd.send(&dislocations[0], num_dislo * 3, i,
"coordinates of dislocation cores");
}
}
}
}
if (this->is_in_dd) {
if (current_step % this->detect_freq == 0) {
STARTTIMER("receiveMDcores");
this->comm_group_atomic.receive(&num_dislo, 1, 0,
"number of dislocation cores");
std::vector<Vector<Dim>> MDs(num_dislo);
this->comm_group_atomic.receive(&MDs[0], num_dislo * 3, 0,
"coordinates of dislocation cores");
if (MDs.size() == 0)
DUMP("No atoms received!", DBG_MESSAGE);
STOPTIMER("receiveMDcores");
Real burgMag = this->domA.getBurgersMagnitude();
burgMag = 1e+10 * burgMag; // in angstrom unit
Real avgSeg0 = (this->minSeg0 + this->maxSeg0) / 2.;
STARTTIMER("moveSlavedDDNodes");
for (auto nd : this->domA.getContainer().getContainerNodes()) {
LM_TOIMPLEMENT;
// UInt nbrs = nd.getNbrOfNeighs();
UInt nbrs = std::numeric_limits<UInt>::max();
Vector<Dim> dd_node = nd.position();
if (((nd.getConstraint() == 7) && (nbrs == 2) &&
this->md_cube->contains(nd.position())) ||
((nd.getConstraint() == 7) && (nbrs == 2) &&
(detect_simple_bool))) {
LM_TOIMPLEMENT;
// Vector<Dim> nbr1({nd.getNeighborNodeCoord(0, 0),
// nd.getNeighborNodeCoord(0, 1),
// nd.getNeighborNodeCoord(0, 2)});
// Vector<Dim> nbr2({nd.getNeighborNodeCoord(1, 0),
// nd.getNeighborNodeCoord(1, 1),
// nd.getNeighborNodeCoord(1, 2)});
// not simple detection ????
if (!detect_simple_bool) {
Vector<Dim> newX = Vector<Dim>::Zero();
Real newC = 0.0;
for (auto &&pos[[gnu::unused]] : MDs) {
LM_TOIMPLEMENT;
// Real dist = (pos - dd_node).norm();
// Real dist1 = (pos - nbr1).norm();
// Real dist2 = (pos - nbr2).norm();
// if (((dist < avgSeg0) && (dist < dist1) && (dist < dist2)) ||
// ((dist1 < avgSeg0) && (dist2 < avgSeg0))) {
// newX += pos;
// newC += 1.0;
// }
}
if (newC != 0.0) {
newX /= newC;
nd.position()[0] = newX[0];
nd.position()[2] = newX[2];
}
nd.unconstrain();
}
// simple detection ????
else {
Real value = -99999.9;
for (auto &&pos : MDs) {
Real dist = this->maxSeg0 * 3.0;
dist = fabs(dd_node[simple_direction] - pos[simple_direction]);
if (dist < avgSeg0 * 2.0)
value = pos[simple_direction];
}
if (value == -99999.9) {
DUMP("DD node (" << dd_node << ")", DBG_MESSAGE);
for (auto &&pos : MDs) {
DUMP("MD core "
<< " (" << pos << ")",
DBG_MESSAGE);
}
LM_FATAL("cannot find the value");
}
nd.position()[simple_direction] = value;
}
} else if ((nd.getConstraint() == 7) && (nbrs == 1) &&
(detect_partial)) {
bool selected(false);
for (auto &&pos : MDs) {
Real dist = sqrt((pos - dd_node).dot(pos - dd_node));
if (dist < avgSeg0) {
nd.position()[0] = pos[0];
nd.position()[2] = pos[2];
avgSeg0 = dist;
selected = true;
}
}
if (selected == false)
DUMP("cannot find the value in partial detection", DBG_MESSAGE);
}
}
if (!detect_simple_bool) {
// merging duplicate nodes
}
STOPTIMER("moveSlavedDDNodes");
STARTTIMER("DDremesh_in_CADD3d");
this->domA.remesh();
STOPTIMER("DDremesh_in_CADD3d");
// release or constrain again after coupling step
for (auto &&nd : this->domA.getContainer().getContainerNodes()) {
if (!detect_simple_bool) {
if (this->md_cube->contains(nd.position())) {
nd.constrain();
} else {
LM_TOIMPLEMENT;
// UInt nbrs = nd.getNbrOfNeighs();
// if (nbrs == 2)
// nd.unconstrain();
}
}
if (this->domA.is_it_bound()) {
if (this->dd_cube->contains(nd.position())) {
nd.constrain();
}
}
}
} // if(current_step % ....)
} // if(is_in_dd) {
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::detect_cores_by_compute(
UInt &num_dislo, std::vector<Vector<Dim>> &dislocations) {
LMID id = compute_list[0];
LM_TOIMPLEMENT;
// ComputeInterface *my_compute = dynamic_cast<ComputeInterface *>(
// FilterManager::getManager().getObject(id));
// if (!my_compute)
// LM_FATAL("invalid compute for detect cores " << compute_list[0]);
// my_compute->build();
// UInt n_nodes_xyz = my_compute->size();
// UInt count = 0;
// UInt num_dislo0 = 0;
std::vector<Vector<Dim>> dislo_pos;
Vector<3u> xmin(this->md_cube->getXmin());
Vector<3u> xmax(this->md_cube->getXmax());
Real pt = this->distance_ddline_from_pad;
for (UInt i : {0, 1, 2}) {
xmin[i] -= pt;
xmax[i] += pt;
}
Cube dd_cube_extended(Dim);
dd_cube_extended.setXmin(xmin);
dd_cube_extended.setXmax(xmax);
LM_TOIMPLEMENT;
// for (count = 0; count < n_nodes_xyz;) {
// Vector<Dim> dd_pos = {my_compute->get(count), my_compute->get(count + 1),
// my_compute->get(count + 2)};
// if (dd_cube_extended.contains(dd_pos) || (this->detect_simple_bool)) {
// dislo_pos.push_back(dd_pos);
// dislocations.push_back(dd_pos);
// num_dislo0 += 1;
// }
// count += 3;
// }
// num_dislo = num_dislo0;
// my_compute->empty();
return;
}
/* --------------------------------------------------------------------------
applyTemplateToPad does following jobs in two kinds of processors (MD and DD)
7. (at every step) set velocity of PAD atoms to zero
(initially it was done first)
DD: 1. close DD dislocation structures
2. compute number of closed DD dislocation loops
2. send to MD processor (with array form)
MD: 1. receive the informations about DD structures from DD processor
2. compute offset parallel curves in order to decide template regions
for each atom...
for each loop...
3. compute Barnett displacement fields
4. compute template displacments when the atom is in the parallel
curves.
5. update total displacements
6. apply total displacements
-------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::applyTemplateToPad() {
// set to zero force and velocity of the pad
auto pad = GeometryManager::getManager().getGeometry(pad_geom);
for (auto &&at : this->domC.getContainer()) {
if (pad->contains(at.position())) {
at.velocity() = Vector<Dim>::Zero();
at.force() = Vector<Dim>::Zero();
}
}
// if not the right moment, do nothing
if (current_step % detect_freq != 0)
return;
// retreive the dd_points (communication)
std::vector<JaehyunDDPoint> dd_points = this->sendDDSegments2MD();
if (this->is_in_atomic) {
STARTTIMER("applyDisplacement");
UInt num_dislo = 0;
std::vector<std::vector<Real>> dislocations;
// Here, we analyze dd_points data to find out
// where and how the cores exist
if ((this->aniso) || (this->volterra)) {
for (auto &dd_pt : dd_points) {
for (UInt j = 0; j < 2; ++j) {
auto neigh = dd_pt.neighs[j];
auto burg = neigh.burg * this->burgMag0;
Real bmag =
burg[0] * cos(this->angle_between_b_and_x / 180.0 * M_PI) +
burg[2] * sin(this->angle_between_b_and_x / 180.0 * M_PI);
bmag = roundf(bmag * 10.0) / 10.0;
if (bmag / fabs(bmag) > 0) {
Real idx = dd_pt.idx;
auto &ddX = dd_pt.X;
Real nbi = neigh.idx;
if ((idx >= 0.0) && (nbi < 0.0)) {
num_dislo += 1;
std::vector<Real> coord_vector;
coord_vector.push_back(ddX[0]);
coord_vector.push_back(ddX[1]);
coord_vector.push_back(ddX[2]);
dislocations.push_back(coord_vector);
}
}
}
}
DUMP(num_dislo << " straight dislocations will be imposed", DBG_MESSAGE);
for (UInt i = 0; i < num_dislo; ++i) {
DUMP("(" << dislocations[i][0] << ", " << dislocations[i][1] << ", "
<< dislocations[i][2] << ")",
DBG_MESSAGE);
}
}
auto md_ranges = this->md_cube->getXmax() - this->md_cube->getXmin();
Real mdxlen = md_ranges[0];
Real mdzlen = md_ranges[2];
UInt count_true_disp = 0;
// IteratorAtomsSubset it = this->pad_atoms.getIterator();
for (auto &&at : this->pad_atoms) {
Vector<Dim> disp = at.displacement();
Vector<Dim> pos0 = at.position0();
Vector<Dim> pos0_pbc = at.position0();
Vector<Dim> pos0_pbc0 = at.position0();
if (fabs(disp[2]) > mdzlen / 2.) {
Real signZ = disp[2] / fabs(disp[2]);
pos0_pbc[2] = pos0[2] + mdzlen * signZ;
count_true_disp += 1;
}
if (fabs(disp[0]) > mdxlen / 2.) {
Real signX = disp[0] / fabs(disp[0]);
pos0_pbc[0] = pos0[0] + mdxlen * signX;
}
for (UInt i = 0; i < 3; ++i) {
pos0_pbc0[i] = pos0_pbc[i];
}
if (this->aniso)
functionifAniso(num_dislo, pos0_pbc, pos0_pbc0, dislocations, at);
else if (this->volterra)
functionifVolterra(num_dislo, pos0_pbc, pos0_pbc0, dislocations, at);
else
functionElse(dd_points, pos0_pbc, at);
}
STOPTIMER("applyDisplacement");
// this->domC.updatePAD();
} // if(is_in_atomic ...)
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::ComputeAnalyticDisplacements(
Vector<Dim> &analytic_disp, Vector<Dim> &B, Vector<Dim> &A,
Vector<Dim> &burg, Vector<Dim> &pos0) {
// Compute closure point for given AB segment
Vector<Dim> C[[gnu::unused]] = Vector<Dim>::Zero();
LM_TOIMPLEMENT;
//this->domA.computeClosurePoint(B, A, burg, C);
// Compute point P1, projected point of C on line consisted to A and Ab
Vector<Dim> P1[[gnu::unused]] = Vector<Dim>::Zero();
LM_TOIMPLEMENT;
//this->domA.computeProjectedPoint(B, C, burg, P1);
// Compute point P2, projected point of C on line consisted to B and Bb
Vector<Dim> P2[[gnu::unused]] = Vector<Dim>::Zero();
LM_TOIMPLEMENT;
//this->domA.computeProjectedPoint(A, C, burg, P2);
Vector<Dim> dispPointU1 = Vector<Dim>::Zero();
LM_TOIMPLEMENT;
//this->domA.computeBarnett(P2, B, A, burg, dispPointU1, pos0, this->pois0);
analytic_disp += dispPointU1;
Vector<Dim> dispPointU2 = Vector<Dim>::Zero();
LM_TOIMPLEMENT;
//this->domA.computeBarnett(P1, B, P2, burg, dispPointU2, pos0, this->pois0);
analytic_disp += dispPointU2;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::ComputeOffsetParallelCurve(
std::vector<JaehyunDDPoint> &dd_points) {
for (auto dd_pt : dd_points) {
std::vector<Real> angles;
Real count = 2.0;
auto node_X = dd_pt.X;
UInt slip_case = 3;
Real avgangle = 0.;
// computes some angle made by the arm in the slip plane
for (auto &&neigh : dd_pt.neighs) {
auto &next_X = neigh.X;
//auto next_normal = neigh.normal;
Real angle_between_neighbor = 0.0;
Vector<Dim> segment = next_X - node_X;
// should construct a rotation which puts
// the normal as the direction Z
LM_TOIMPLEMENT;
angle_between_neighbor = atan2(segment[0], segment[1]);
angles.push_back(angle_between_neighbor);
avgangle += angle_between_neighbor;
}
avgangle /= dd_pt.neighs.size();
Vector<Dim> dX1 = {0.0, 0.0, 0.0};
Vector<Dim> dX2 = {0.0, 0.0, 0.0};
computeOffset(node_X, avgangle, dX1, dX2, slip_case);
if (count != 0.0) {
if (isnan(dX1[0]) || isnan(dX1[1]) || isnan(dX1[2]) || isnan(dX2[0]) ||
isnan(dX2[1]) || isnan(dX2[2]))
LM_FATAL("wrong dX1 and dX2");
}
angles[0] = radian2degree(angles[0]);
angles[1] = radian2degree(angles[1]);
avgangle = (angles[0] + angles[1]) / count;
findCorrectOffsets(avgangle, dX1, dX2, slip_case);
dd_pt.tmp_X = node_X + dX1;
dd_pt.tmp2_X = node_X + dX2;
}
return;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::computeOffset(
Vector<Dim> ¤t_node, Real ang, Vector<Dim> &dX1, Vector<Dim> &dX2,
UInt slip_case) {
Real rmax = this->core_rcut;
Real slope = tan(ang);
Real val1 = sqrt(rmax * rmax / (1.0 + slope * slope));
Real val2 = sqrt(rmax * rmax - rmax * rmax / (1.0 + slope * slope));
Real val3 = -sqrt(rmax * rmax / (1.0 + slope * slope));
Real val4 = -sqrt(rmax * rmax - rmax * rmax / (1.0 + slope * slope));
if (slip_case == 0) {
dX1[0] = 0.0;
dX1[1] = val1;
dX1[2] = val2;
dX2[0] = 0.0;
dX2[1] = val3;
dX2[2] = val4;
} else if (slip_case == 1) {
dX1[0] = val1;
dX1[1] = 0.0;
dX1[2] = val2;
dX2[0] = val3;
dX2[1] = 0.0;
dX2[2] = val4;
} else if (slip_case == 2) {
dX1[0] = val1;
dX1[1] = val2;
dX1[2] = 0.0;
dX2[0] = val3;
dX2[1] = val4;
dX2[2] = 0.0;
} else
LM_FATAL("cannot define angle with neighbor");
return;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::findCorrectOffsets(
Real ang, Vector<Dim> &dX1, Vector<Dim> &dX2, UInt slip_case) {
Real dev = 1e+5;
Real saved_i = 1.0;
Real saved_k = 1.0;
for (int i = -1; i < 3;) {
for (int k = -1; k < 3;) {
Real tmp = -1.0;
if (slip_case == 0)
tmp = fabs(radian2degree(atan2(dX2[1] * (Real)i, dX2[2] * (Real)k)) -
ang);
else if (slip_case == 1)
tmp = fabs(radian2degree(atan2(dX2[2] * (Real)i, dX2[0] * (Real)k)) -
ang);
else if (slip_case == 2)
tmp = fabs(radian2degree(atan2(dX2[0] * (Real)i, dX2[1] * (Real)k)) -
ang);
else
LM_FATAL("unknown slip_case");
if (tmp < dev) {
dev = tmp;
saved_i = (Real)i;
saved_k = (Real)k;
}
k += 2;
}
i += 2;
}
if (slip_case == 0) {
dX1[0] = dX1[0];
dX1[1] = dX1[1] * saved_i;
dX1[2] = dX1[2] * saved_k;
dX2[0] = dX2[0];
dX2[1] = dX2[1] * saved_i;
dX2[2] = dX2[2] * saved_k;
} else if (slip_case == 1) {
dX1[0] = dX1[0] * saved_k;
dX1[1] = dX1[1];
dX1[2] = dX1[2] * saved_i;
dX2[0] = dX2[0] * saved_k;
dX2[1] = dX2[1];
dX2[2] = dX2[2] * saved_i;
} else if (slip_case == 2) {
dX1[0] = dX1[0] * saved_i;
dX1[1] = dX1[1] * saved_k;
dX1[2] = dX1[2];
dX2[0] = dX2[0] * saved_i;
dX2[1] = dX2[1] * saved_k;
dX2[2] = dX2[2];
}
return;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
bool CADD_template<DomainDD, DomainMD, Dim>::in_side_ellipse(Real pos[],
Real rmax,
Real rmin,
Real *factor) {
Real eq =
((pos[0] / rmax) * (pos[0] / rmax) + (pos[1] / rmin) * (pos[1] / rmin));
Real val = 1.0 - eq;
if (val >= 0.0) {
*factor = sqrt(eq);
return true;
} else {
*factor = sqrt(eq);
return false;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
bool CADD_template<DomainDD, DomainMD, Dim>::is_it_infinite_connection(
UInt index0, UInt *index0_infinite) {
int int_index0_infinite;
bool found = this->domA.getPairedIndex(index0, &int_index0_infinite);
if (found) {
*index0_infinite = (UInt)int_index0_infinite;
return true;
} else
return false;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::shear_strain(RefAtom at,
Vector<Dim> &disp) {
Real eps_yx = this->shear_stresses[1] / this->mu0;
Real eps_yz = this->shear_stresses[3] / this->mu0;
Real eps_xy = this->shear_stresses[0] / this->mu0;
Real eps_zy = this->shear_stresses[2] / this->mu0;
disp[0] = eps_yx * at.position0()[1];
disp[1] = eps_xy * at.position0()[0] + eps_zy * at.position0()[2];
disp[2] = eps_yz * at.position0()[1];
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::tilt_strain(RefAtom at,
Vector<Dim> &disp) {
// Real LX = this->md_box_lengths[0];
Real LY = this->md_box_lengths[1];
Real LZ = this->md_box_lengths[2];
// Real LY = fabs(this->md_yrange[1]) + fabs(this->md_yrange[0]);
// Real LZ = fabs(this->md_zrange[1]) + fabs(this->md_zrange[0]);
Real eps_xy = this->tilt_stresses[0] / this->mu0;
Real eps_xz = this->tilt_stresses[1] / this->mu0;
Real eps_yz = this->tilt_stresses[2] / this->mu0;
disp[0] = -1.0 * eps_xy * (at.position0()[1] + LY / 2.0) -
1.0 * eps_xz * (at.position0()[2] + LZ / 2.0);
// disp[1] = -1.0*eps_yz*(at.position0(2) + LZ/2.0); // full dislocation
disp[1] = -1.0 * eps_yz * (at.position0()[2] + LZ / 2.0); // half dislocation
disp[2] = 0.0;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::close_dislocation_loop(
std::vector<JaehyunDDPoint> &dd_points,
std::vector<JaehyunDDPoint> &dd_points_global) {
bool added_point_in_infinity(false);
Real maxdist = this->domA.getMaxSeg();
UInt count = 0;
for (auto &dd_pt_global : dd_points_global) {
Real idx_real = dd_pt_global.idx;
Real num_neigh = dd_pt_global.num_neigh;
auto crr_X = dd_pt_global.X;
auto &dd_pt = dd_points[count];
dd_pt.idx = idx_real;
dd_pt.num_neigh = num_neigh;
dd_pt.X = crr_X;
for (UInt n = 0; n < num_neigh; ++n) {
auto neigh = dd_pt_global.neighs[n];
Real nbr_idx = neigh.idx;
auto &nbr_X = neigh.X;
auto burg = neigh.burg;
auto normal = neigh.normal;
auto diff = nbr_X - crr_X;
Real dist = diff.norm();
UInt nbr_idx0_infinite = 0;
Real nbr_idx_infinite_real = 0.0;
bool found(false);
UInt idx0 = (UInt)dd_pt_global.idx0;
if ((this->close_in_infinity) &&
(is_it_infinite_connection(idx0, &nbr_idx0_infinite)) &&
(num_neigh == 1)) {
// connectivity over an infinity boundary for single neighbors nodes:
// frank-read sources
Vector<Dim> nbr_X_infinite;
for (auto &dd_pt_global2 : dd_points_global) {
Real tmp = dd_pt_global2.idx0;
UInt nbnbr = dd_pt_global2.num_neigh;
if (((tmp - nbr_idx0_infinite) < 1e-2) && (nbnbr == 1)) {
nbr_idx_infinite_real = dd_pt_global2.idx;
nbr_X_infinite = dd_pt_global2.X;
found = true;
}
}
if (!found)
LM_FATAL("not found a node infinitely connected for idx0 " << idx0);
DUMP("infinity points are considered between "
<< idx0 << " (" << crr_X << ") and " << nbr_idx0_infinite
<< " (" << nbr_X_infinite << ") in ("
<< this->close_direction[0] << ", " << this->close_direction[1]
<< ", " << this->close_direction[2] << ") direction",
DBG_MESSAGE);
dd_pt.num_neigh = 2;
dd_pt.neighs[0].idx = nbr_idx;
dd_pt.neighs[0].X = nbr_X;
dd_pt.neighs[0].burg = burg;
dd_pt.neighs[0].normal = normal;
dd_pt.neighs[1].idx = -1.0 * idx_real - 10000.0;
for (UInt i : {0, 1, 2}) {
dd_pt.neighs[1].X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
}
dd_pt.neighs[1].burg = burg * -1.;
dd_pt.neighs[1].normal = normal;
auto &dd_pt_plus1 = *((&dd_pt) + 1);
dd_pt_plus1.idx = -1.0 * idx_real - 10000.0;
dd_pt_plus1.num_neigh = 2;
for (UInt i : {0, 1, 2}) {
dd_pt_plus1.X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
}
dd_pt_plus1.neighs[0].idx = idx_real;
dd_pt_plus1.neighs[0].X = crr_X;
dd_pt_plus1.neighs[0].burg = burg;
dd_pt_plus1.neighs[0].normal = normal;
dd_pt_plus1.neighs[1].idx = -1.0 * nbr_idx_infinite_real - 10000.0;
for (UInt i : {0, 1, 2}) {
dd_pt_plus1.neighs[1].X[i] =
nbr_X_infinite[i] + close_direction[i] * 1E+3;
}
dd_pt_plus1.neighs[1].burg = burg * -1.;
dd_pt_plus1.neighs[1].normal = normal;
added_point_in_infinity = true;
} else if ((close_in_infinity) && (num_neigh == 2) &&
(dist > 5.0 * maxdist)) {
// connectivity over an infinity boundary for Real neighbors nodes:
// straight dislocation
dd_pt.neighs[n].idx = -1.0 * dd_pt.idx - 10000.0;
for (UInt i : {0, 1, 2}) {
dd_pt.neighs[n].X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
}
dd_pt.neighs[n].burg = burg;
dd_pt.neighs[n].normal = normal;
auto &dd_pt_plus1 = *((&dd_pt) + 1);
dd_pt_plus1.idx = -1.0 * dd_pt.idx - 10000.0;
dd_pt_plus1.num_neigh = 2;
for (auto i : {0, 1, 2}) {
dd_pt_plus1.X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
}
dd_pt_plus1.idx0 = dd_pt.idx;
dd_pt_plus1.neighs[0].X = crr_X;
dd_pt_plus1.neighs[0].burg = burg * -1.;
dd_pt_plus1.neighs[0].normal = normal;
dd_pt_plus1.neighs[1].idx = -1.0 * nbr_idx - 10000.0;
for (auto i : {0, 1, 2}) {
dd_pt_plus1.neighs[1].X[i] =
nbr_X[i] + this->close_direction[i] * 1E+3;
}
dd_pt_plus1.neighs[1].burg = burg;
dd_pt_plus1.neighs[1].normal = normal;
DUMP("Implicit infinity connection between "
<< idx0 << " (" << dd_pt.X << ") and " << dd_pt_plus1.X << " ("
<< nbr_X << ").",
DBG_MESSAGE);
added_point_in_infinity = true;
} else if (num_neigh == 1) {
// && (this->detect_partial)) {
// straight connection between two single neighbor node: partial
// detection
dd_pt.neighs[0].idx = nbr_idx;
dd_pt.neighs[0].X = nbr_X;
dd_pt.neighs[0].burg = burg;
dd_pt.neighs[0].normal = normal;
Real temp_val = 1e+4;
// IteratorDDNodes nodes_loop3 =
// this->domA.getContainerNodes().getIterator();
// Real tmpx,tmpy,tmpz;
for (auto &dd_pt_global2 : dd_points_global) {
UInt num_neigh0 = dd_pt_global2.num_neigh;
auto diff = dd_pt_global2.X - crr_X;
Real dist0 = diff.norm();
if ( // has to be terminated in MD domain
(num_neigh0 == 1) &&
// avoid self selection
(0.1 < dist0) &&
// find the node minimum distance
(dist0 < temp_val) &&
// avoid chosing the node in different slip plane
(fabs(dd_pt_global2.X[1] - crr_X[1]) < 1.0)) {
temp_val = dist;
dd_pt.neighs[0].idx = dd_pt_global2.idx;
dd_pt.neighs[0].X = dd_pt_global2.X;
DUMP("Straight connection between the two fixed DD nodes "
<< "(" << dd_pt.X << ") and "
<< "(" << dd_pt_global2.X << ").",
DBG_MESSAGE);
}
}
dd_pt.num_neigh = 2.0;
dd_pt.neighs[1].burg = -1.0 * dd_pt.neighs[0].burg;
dd_pt.neighs[1].normal = dd_pt.neighs[0].normal;
} else {
// normal connection
dd_pt.neighs[n].idx = nbr_idx;
dd_pt.neighs[n].X = nbr_X;
dd_pt.neighs[n].burg = burg;
dd_pt.neighs[n].normal = normal;
}
}
if (added_point_in_infinity) {
count += 2;
added_point_in_infinity = false;
} else {
count += 1;
}
}
return;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
UInt CADD_template<DomainDD, DomainMD, Dim>::CollectingLocalDDpoints(
std::vector<JaehyunDDPoint> &dd_points0_local) {
bool problem_in_neigh = false;
int count = 0;
for (auto &node : this->domA.getContainer().getContainerNodes()) {
Real idx;
idx = (Real)node.getIndex();
Real idx0;
idx0 = (Real)node.getIndex0();
LM_TOIMPLEMENT;
//int num_neigh = node.getNbrOfNeighs();
int num_neigh = std::numeric_limits<int>::max();
Vector<Dim> crr_X = {roundf(node.position()[0] * 10.0) / 10.0,
roundf(node.position()[1] * 10.0) / 10.0,
roundf(node.position()[2] * 10.0) / 10.0};
if (num_neigh >= 3) {
DUMP("[collectingLocalDDpoints] problem in DD node idx"
<< idx << " (" << crr_X << "): number of neighbors "
<< num_neigh,
DBG_MESSAGE);
num_neigh = 2;
problem_in_neigh = true;
}
dd_points0_local[count].idx = idx;
dd_points0_local[count].num_neigh = (Real)num_neigh;
dd_points0_local[count].X = crr_X;
dd_points0_local[count].idx0 = (Real)idx0;
for (int n = 0; n < num_neigh; ++n) {
LM_TOIMPLEMENT;
//Real nbr_idx = (Real)node.getNeighborNodeIndex(n);
// Vector<Dim> nbr_X = {
// roundf(node.getNeighborNodeCoord(n, 0) * 10.0) / 10.0,
// roundf(node.getNeighborNodeCoord(n, 1) * 10.0) / 10.0,
// roundf(node.getNeighborNodeCoord(n, 2) * 10.0) / 10.0};
// Vector<Dim> burg = {node.burgersVectorX(n), node.burgersVectorY(n),
// node.burgersVectorZ(n)};
// Vector<Dim> normal = {node.normalVectorX(n), node.normalVectorY(n),
// node.normalVectorZ(n)};
// auto &neigh = dd_points0_local[count].neighs[n];
// neigh.idx = nbr_idx;
// neigh.burg = burg;
// neigh.normal = normal;
}
}
if (this->dd_points_out) {
// Communicator &comm = Communicator::getCommunicator();
UInt my_rank_in_DD = this->comm_group_dd.getMyRank();
std::stringstream fname0;
fname0 << this->directory_dd_points << "/local_dd_points0_step"
<< current_step << ".txt.proc" << my_rank_in_DD;
std::ofstream myfile0;
for (auto &dd_pt0 : dd_points0_local) {
if (!myfile0.is_open()) {
myfile0.open(fname0.str().c_str());
}
myfile0 << dd_pt0 << "\t";
}
myfile0.close();
}
if (problem_in_neigh)
return 1;
else
return 0;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
bool CADD_template<DomainDD, DomainMD, Dim>::coarse_dislocation_loop(
std::vector<JaehyunDDPoint> &dd_points,
std::vector<JaehyunDDPoint> &dd_points0) {
std::vector<UInt> indices;
Vector<Dim> burg = {0.0, 0.0, 0.0};
Vector<Dim> norm = {0.0, 0.0, 0.0};
Real ci = 0.;
UInt count = 0;
UInt master_count = 0;
for (UInt i = 0; i < dd_points0.size(); ++i) {
if (!(std::find(indices.begin(), indices.end(), i) != indices.end())) {
auto &dd_pt = dd_points0[i];
Real ci0 = dd_pt.idx;
auto cp0 = dd_pt.X;
if (master_count >= dd_points0.size()) {
DUMP("[coarse_dislocation_loop] exit 1", DBG_MESSAGE);
return false;
}
auto &dd_pt_master = dd_points[master_count];
auto &dd_pt_master_1 = dd_points[master_count - 1];
auto &dd_pt_master_2 = dd_points[master_count - count];
dd_pt_master.idx = ci0; // idx
dd_pt_master.num_neigh = 2.0; // #nbr
dd_pt_master.X = cp0; // x-coord
master_count += 1;
count += 1;
ci = ci0;
Vector<Dim> cp = cp0;
indices.push_back(i);
bool loop_start(true);
bool not_finish(true);
while (not_finish) {
if (!loop_start) {
if (is_it_closed(cp, cp0, ci, ci0)) {
not_finish = false;
if (master_count - 1 >= dd_points0.size()) {
DUMP("[coarse_dislocation_loop] exit 2", DBG_MESSAGE);
return false;
}
dd_pt_master_1.neighs[0].idx = ci0; // idx-nbr1
dd_pt_master_1.neighs[0].X = cp0; // x-coord-nbr1
dd_pt_master_1.neighs[0].burg = burg; // bx-nbr1
dd_pt_master_1.neighs[0].normal = norm; // bx-nbr1
dd_pt_master_2.neighs[1].idx = dd_pt_master_1.idx; // idx-nbr2
dd_pt_master_2.neighs[1].X = dd_pt_master_1.X; // x-coord-nbr2
dd_pt_master_2.neighs[1].burg = -1.0 * burg; // bx-nbr2
dd_pt_master_2.neighs[1].normal = norm; // bx-nbr2
count = 0;
break;
}
}
UInt count_inside = 0;
bool straight(true);
while (straight) {
auto n1 = dd_pt.neighs[0].X;
Real b1mag = dd_pt.neighs[0].burg[0] *
cos(this->angle_between_b_and_x / 180.0 * M_PI) +
dd_pt.neighs[0].burg[2] *
sin(this->angle_between_b_and_x / 180.0 * M_PI);
Real n1idx_real = dd_pt.neighs[1].idx;
UInt n1idx = (UInt)n1idx_real;
auto n2 = dd_pt.neighs[1].X;
Real b2mag = dd_pt.neighs[1].burg[0] *
cos(this->angle_between_b_and_x / 180.0 * M_PI) +
dd_pt.neighs[1].burg[2] *
sin(this->angle_between_b_and_x / 180.0 * M_PI);
Real n2idx_real = dd_pt.neighs[2].idx;
UInt n2idx = (UInt)n2idx_real;
UInt saved_i = i;
UInt tmp = find_one_direction(n1, n2, b1mag, b2mag, n1idx, n2idx, &i);
if ((tmp == 0) && (master_count == 1)) {
burg = 1.0 * dd_points0[saved_i].neighs[0].burg;
} else if ((tmp == 1) && (master_count == 1)) {
burg = -1.0 * dd_points0[saved_i].neighs[0].burg;
}
norm = dd_points0[saved_i].neighs[0].normal;
// check whether this local connection is important or not
Real area_in_yplane1 = AreaTriangle(cp, n1, n2, 1);
Vector<Dim> last_registered_node = dd_pt_master_1.X;
Vector<Dim> next_node = dd_pt.X;
Real area_in_yplane2 =
AreaTriangle(last_registered_node, cp, next_node, 1);
if ((area_in_yplane1 > this->min_area) ||
(area_in_yplane2 > this->min_area) ||
(fabs(n2[1] - cp[1]) > 1.0) ||
(fabs(n1[1] - cp[1]) > 1.0)) { // || (dist > this->maxSeg0*10.0)){
if (master_count - 1 >= dd_points0.size() - 1) {
DUMP("[coarse_dislocation_loop] exit 3", DBG_MESSAGE);
return false;
}
dd_pt_master.neighs[1].idx = dd_pt_master_1.idx; // idx-nbr2
dd_pt_master.neighs[1].X = dd_pt_master_1.X; // x-coord-nbr2
dd_pt_master.neighs[1].burg = -1.0 * burg; // bx-nbr2
dd_pt_master.neighs[1].normal = norm; // nx-nbr2
dd_pt_master_1.neighs[0].idx = ci; // idx-nbr1
dd_pt_master_1.neighs[0].X = cp; // x-coord-nbr1
dd_pt_master_1.neighs[0].burg = 1.0 * burg; // bx-nbr1
dd_pt_master_1.neighs[0].normal = norm; // bx-nbr1
if (!loop_start) {
dd_pt_master.idx = ci; // idx
dd_pt_master.num_neigh = 2.0; // #nbr
dd_pt_master.X = cp; // x-coord
master_count += 1;
count += 1;
}
}
loop_start = false;
indices.push_back(i);
ci = dd_pt.idx;
cp = dd_pt.X;
count_inside += 1;
if (is_it_closed(cp, cp0, ci, ci0))
break;
if (count_inside >= dd_points0.size()) {
DUMP("[coarse_dislocation_loop] exit 4", DBG_MESSAGE);
return false;
}
}
if (master_count >= dd_points0.size()) {
DUMP("[coarse_dislocation_loop] exit 5", DBG_MESSAGE);
return false;
}
}
if (master_count >= dd_points0.size()) {
DUMP("[coarse_dislocation_loop] exit 6", DBG_MESSAGE);
return false;
}
}
if (master_count >= dd_points0.size()) {
DUMP("[coarse_dislocation_loop] exit 7", DBG_MESSAGE);
return false;
}
}
return true;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::SetPadAtoms() {
UInt registered = 0;
int my_rank_in_MDprocs = this->comm_group_atomic.getMyRank();
this->pad_atoms.clear();
auto ref_manager = this->domC.getContainer().getRefManager();
auto pad = GeometryManager::getManager().getGeometry(pad_geom);
for (auto &at : this->domC.getContainer()) {
if (pad->contains(at.position())) {
- this->pad_atoms.add(at);
+ this->pad_atoms.push_back(at);
registered += 1;
}
}
this->comm_group_atomic.reduce(®istered, 1, "atoms: registered as PAD",
OP_SUM);
ref_manager->addSubSet("pad atoms", this->pad_atoms);
if (my_rank_in_MDprocs == 0) {
DUMP(registered << " atoms are registered as pad atoms", DBG_MESSAGE);
}
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::find_indices_for_neighbors(
std::vector<JaehyunDDPoint> &dd_points) {
for (auto &dd_pt : dd_points) {
if (dd_pt.num_neigh != 2)
LM_FATAL("this is not a proper dislocation network");
for (UInt j = 0; j < 2; ++j) {
auto neigh = dd_pt.neighs[j];
Real next_idx = neigh.idx;
Vector<Dim> next_X = neigh.X;
UInt k = 0;
for (auto &dd_pt2 : dd_points) {
Real node_idx = dd_pt2.idx;
Vector<Dim> normal = dd_pt2.X;
Real dist = (normal - next_X).norm();
if ((dist < 1e-3) && (fabs(node_idx - next_idx) <= 1e-3)) {
neigh.idx = (Real)k;
// break;
++k;
}
}
}
}
}
// /*
// --------------------------------------------------------------------------
// */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::applyDisplacement(
RefAtom at, Vector<Dim> &dislo_disp) {
Vector<Dim> disp_strain = Vector<Dim>::Zero();
Vector<Dim> disp_tilt = Vector<Dim>::Zero();
shear_strain(at, disp_strain);
tilt_strain(at, disp_tilt);
at.position() = at.position0() + dislo_disp - disp_strain + disp_tilt;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::coordinate_for_2d_template(
Real angle, Real local_coord[], Real global_coord[],
Real template_coord[]) {
template_coord[0] =
local_coord[0] * local_coord[0] / fabs(local_coord[0] + 1e-8) *
sin(angle * M_PI / 180.0) +
local_coord[2] * local_coord[2] / fabs(local_coord[2] + 1e-8) *
cos(angle * M_PI / 180.0);
template_coord[1] =
local_coord[1] * local_coord[1] / fabs(local_coord[1] + 1e-8);
}
/////// This part below should be removed later
////////////////////////////////////
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real CADD_template<DomainDD, DomainMD, Dim>::eval_q_coef(Real *X) {
Real inside = X[0] * X[0] +
2.0 * X[0] * X[1] * this->lambda * cos(this->phi) +
X[1] * X[1] * this->lambda * this->lambda;
return sqrt(inside);
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real CADD_template<DomainDD, DomainMD, Dim>::eval_t_coef(Real *X) {
Real inside = X[0] * X[0] -
2.0 * X[0] * X[1] * this->lambda * cos(this->phi) +
X[1] * X[1] * this->lambda * this->lambda;
return sqrt(inside);
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real CADD_template<DomainDD, DomainMD, Dim>::evalArctan(Real X, Real Y) {
if ((X >= 0.0) && (Y > 0.0)) {
return atan(Y / X);
} else if ((X < 0.0) && (Y >= 0.0)) {
return atan(-Y / -X) + M_PI;
} else if ((X < 0.0) && (Y < 0.0)) {
return atan(-Y / -X) - M_PI;
} else if ((X >= 0.0) && (Y < 0.0)) {
return atan(Y / X);
} else {
LM_FATAL("unexpected case");
}
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real
CADD_template<DomainDD, DomainMD, Dim>::computeBeDisplacementAniso(Real b_edge,
Real *X) {
Real q = eval_q_coef(X);
Real t = eval_t_coef(X);
Real c12 = this->elastic_coefs[1];
Real c66 = this->elastic_coefs[12];
Real first = 0.0;
first += evalArctan((X[0] + X[1] * this->lambda * cos(this->phi)),
(X[1] * this->lambda * sin(this->phi)));
first += evalArctan((X[0] - X[1] * this->lambda * cos(this->phi)),
(X[1] * this->lambda * sin(this->phi)));
Real second = (this->c11bar * this->c11bar - c12 * c12) /
(2.0 * this->c11bar * c66 * sin(2.0 * this->phi)) *
floor(log(q / t) * 1e+5) / 1e+5;
Real res = b_edge / (4.0 * M_PI) * (first + second);
return res;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real
CADD_template<DomainDD, DomainMD, Dim>::computeNDisplacementAniso(Real b_edge,
Real *X) {
Real q = eval_q_coef(X);
Real t = eval_t_coef(X);
Real c12 = this->elastic_coefs[1];
Real first =
(this->c11bar - c12) * cos(this->phi) * floor(log(q * t) * 1e+5) / 1e+5;
Real second =
(this->c11bar + c12) * sin(this->phi) *
atan2((X[1] * X[1] * this->lambda * this->lambda * sin(2.0 * this->phi)),
(X[0] * X[0] -
this->lambda * this->lambda * X[1] * X[1] * cos(2.0 * this->phi)));
Real res = -this->lambda * b_edge /
(4.0 * M_PI * this->c11bar * sin(2.0 * this->phi)) *
(first - second);
return res;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real
CADD_template<DomainDD, DomainMD, Dim>::computeLDisplacementAniso(Real b_screw,
Real *X) {
Real c44 = this->elastic_coefs[9];
Real c45 = this->elastic_coefs[10];
Real c55 = this->elastic_coefs[11];
Real res = b_screw / 2.0 / M_PI *
atan2(sqrt(c44 * c55 - c45 * c45) * X[1], c44 * X[0] - c45 * X[1]);
return res;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real
CADD_template<DomainDD, DomainMD, Dim>::computeBeDisplacement(Real b_edge,
Real *X) {
Real x_2 = X[0] * X[0];
Real y_2 = X[1] * X[1];
Real nu = this->pois0;
Real epsilon = 1e-15;
Real res = b_edge / (2.0 * M_PI) *
(atan2(X[1], X[0]) +
X[0] * X[1] / (2.0 * (1 - nu) * (x_2 + y_2 + epsilon)));
return res;
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real
CADD_template<DomainDD, DomainMD, Dim>::computeNDisplacement(Real b_edge,
Real *X) {
Real x_2 = X[0] * X[0];
Real y_2 = X[1] * X[1];
Real epsilon = 1e-15;
Real nu = this->pois0;
// Real b_edge2 = b_edge*b_edge;
// Real res = -b_edge/(2.0*M_PI)*((1-2.0*nu)/(4.0*(1-nu))*(log((x_2 +
// y_2+epsilon)/b_edge2))
// + (x_2 - y_2)/
// (4.0*(1-nu)*(x_2 + y_2+epsilon))
// );
// return res;
return -b_edge / (2.0 * M_PI) *
((1 - 2.0 * nu) / (4.0 * (1 - nu)) * log(x_2 + y_2 + epsilon) +
(x_2 - y_2) / (4.0 * (1 - nu) * (x_2 + y_2 + epsilon)));
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
inline Real
CADD_template<DomainDD, DomainMD, Dim>::computeLDisplacement(Real b_screw,
Real *X) {
Real res = b_screw / (2.0 * M_PI) * (atan2(X[1], X[0]));
return res;
}
/* --------------------------------------------------------------------------
*/
/* LMDESC CADDTEMPLATE
CADDTEMPLATE is a simple CADD-like coupler which can deal with straing egde
dislocations normal to the interface (i.e., 1;3400;0cno elastic problem)
*/
/* LMEXAMPLE
COUPLING_CODE
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::declareParams() {
/* LMKEYWORD PAD_GEOM
Identifier of the pad geom
*/
this->parseKeyword("PAD_GEOM", this->pad_geom);
/* LMKEYWORD DETECT_DXA
Boolean to notify input from DXA results
*/
this->parseTag("DETECT_DXA", this->detect_dxa, false);
/* LMKEYWORD MULTIPLE_LOOPS
Boolean to notify multiple loops
*/
this->parseTag("MULTIPLE_LOOPS", this->multiple_loops, false);
/* LMKEYWORD DETECT_PARTIAL
Boolean to notify input for partially detection dislocation
*/
this->parseTag("DETECT_PARTIAL", this->detect_partial, false);
/* LMKEYWORD PRINT_DD_POINTS
Boolean to print internal data named dd\_points and its offset lines
*/
this->parseTag("PRINT_DD_POINTS", this->dd_points_out, false);
/* LMKEYWORD COARSEN_DD
Boolean to decide to coarse dislocation loop before apply Barnett fields
*/
this->parseTag("COARSEN_DD", this->coarsen_dd, false);
/* LMKEYWORD MIN_AREA
Specifies the minimum area to be coarsed out
*/
this->parseKeyword("MIN_AREA", this->min_area, 0.0);
/* LMKEYWORD DIRECTORY_DD_POINTS
fies the directory dd points are written
*/
this->parseKeyword("DIRECTORY_DD_POINTS", this->directory_dd_points, "/.");
/* LMKEYWORD DD_BOUNDARY
Boolean to set DD boundary
*/
this->parseTag("DD_BOUNDARY", this->dd_boundary, false);
/* LMKEYWORD SLIP_PLANE
Specifies the normal vector to the slip plane
*/
this->parseVectorKeyword("SLIP_PLANE", 3, this->normal_quant);
/* LMKEYWORD ANGLE_BETWEEN_B_AND_X
Specifies the angle between Burgers vector and X(11-2) axis
*/
this->parseKeyword("ANGLE_BETWEEN_B_AND_X", this->angle_between_b_and_x);
/* LMKEYWORD MD_GEOM
Identifier of the MD geometry
*/
this->parseKeyword("MD_GEOM", this->md_geom);
/* LMKEYWORD DD_GEOM
Identifier of the DD geometry
*/
this->parseKeyword("DD_GEOM", this->dd_geom);
/* LMKEYWORD UPDATE_FREQUENCY
Specifies the frequency at which detection computations are performed
*/
this->parseKeyword("UPDATE_FREQUENCY", this->detect_freq, 50);
/* LMKEYWORD DETECT_COMPUTE
Add computes to input to filter
*/
this->parseKeyword("DETECT_COMPUTE", compute_list);
/* LMKEYWORD NUM_TEMPLATES
Specifies the number of templates provided.
*/
this->parseKeyword("NUM_TEMPLATES", this->num_templates);
/* LMKEYWORD TEMPLATE_FILE
Specifies the template file.
*/
this->parseKeyword("TEMPLATE_FILE", this->templatefilename_list);
/* LMKEYWORD CORE_RCUT
Specifies the radius of the region where core template applied
*/
this->parseKeyword("CORE_RCUT", this->core_rcut);
/* LMKEYWORD RATIO_OF_PURE_CORE_TEMPLATE
Specifies the ratio of region where only pure core template is applied
with respect to ellipse defined by RMAX and RMIN
where
*/
this->parseKeyword("RATIO_OF_PURE_CORE_TEMPLATE", this->factor);
/* LMKEYWORD SHEAR_STRESSES
Specifies the shear stersses with respect to Y plane.
*/
this->parseVectorKeyword("SHEAR_STRESSES", 4, this->shear_stresses);
/* LMKEYWORD TILT_STRESSES
Specifies the shear stersses imposed by tilting.
*/
this->parseVectorKeyword("TILT_STRESSES", 3, this->tilt_stresses);
/* LMKEYWORD SHEAR_MODULUS
Specifies the shear modulus in Burgers vector direction.
mega-Pascal unit
*/
this->parseKeyword("SHEAR_MODULUS", this->mu0);
/* LMKEYWORD POISSON_RATIO
Specifies the shear modulus in Burgers vector direction.
mega-Pascal unit
*/
this->parseKeyword("POISSON_RATIO", this->pois0);
/* LMKEYWORD DOMAIN_DD
Specifies the dd domain id
*/
this->parseKeyword("DOMAIN_DD", this->domdd);
/* LM_KEYWORD DISTANCE_DDLINE_FROM_PAD
Specify the distance from pad (inner) boudaries to DD constrained nodes.
*/
this->parseKeyword("DISTANCE_DDLINE_FROM_PAD",
this->distance_ddline_from_pad);
/* LMKEYWORD PAD_THICKNESS
Specify the thickness of pad zone
*/
this->parseKeyword("PAD_THICKNESS", this->pad_thickness);
// In case of (0, 0, 1):
// point1 point2
// *----------------------*
// | | z
// | | ^
// | | |
// o-o-o-x x-o-o-o-o +--->x
/* LMKEYWORD CLOSE_IN_INFINITY
Specifies the direction in order to close
dislocation loop
*/
this->parseVectorKeyword("CLOSE_IN_INFINITY", 3, this->close_direction,
VEC_DEFAULTS(0, 0, 0));
/* LMKEYWORD NUBER_INFINITY_CONNECTION
Specifies the direction in order to close
dislocation loop
*/
this->parseKeyword("NUBER_INFINITY_CONNECTION",
this->num_close_loop_in_infinity, 0);
/* LMKEYWORD MD_BOX_LENGTHS
Specifies the (original) MD box size
to specify correct tilt stress values
*/
this->parseVectorKeyword("MD_BOX_LENGTHS", 3, this->md_box_lengths);
/* LMKEYWORD DETECT_SIMPLE_DIRECTION
Specifies the (original) MD box size
to specify correct tilt stress values
*/
this->parseVectorKeyword("DETECT_SIMPLE_DIRECTION", 3,
this->detect_simple_direction,
VEC_DEFAULTS(0, 0, 0));
// /* LM**KEYWORD NB_SLIPS
// Specifies the number of slip planes
// */
// this->parseKeyword("NB_SLIPS", this->nb_slips);
// /* LM**KEYWORD LEN_INTER_SLIPS
// Specifies length between two slip planes
// */
// this->parseKeyword("LEN_INTER_SLIPS", this->len_inter_slips);
/* LMKEYWORD ANISO
Boolean to set anisotropicity:
[1] Elastic Strain Fields and Dislocation Mobility, Volume 31, J.Lothe
Dislocations in anisotropic media
[2] Plane-strain discrete dislocation plasticity incorporating
anisotropic
elasticity IJSS 48(2), January 2011
*/
this->parseTag("ANISO", this->aniso, false);
/* LMKEYWORD VOLTERRA
TODO
*/
this->parseTag("VOLTERRA", this->volterra, false);
/* LMKEYWORD ELASTIC_COEFS
Specify the thirteen elastic coefficients
of the constitutive matrix of anisotropic media in the following order:
C11, C12, C13, C21, C22, C23, C31, C32, C33, C44, C45, C55, C66
*/
this->parseVectorKeyword("ELASTIC_COEFS", 13, this->elastic_coefs,
VEC_DEFAULTS(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0));
LM_TOIMPLEMENT;
// forward parsable of the core template needs to be done
// this->addSubParsableObject(this->core_template);
}
/* ---------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
std::vector<JaehyunDDPoint>
CADD_template<DomainDD, DomainMD, Dim>::sendDDSegments2MD() {
// Communicator &comm = Communicator::getCommunicator();
// UInt local_data_length = 26;
UInt data_length = 33;
std::vector<JaehyunDDPoint> dd_points;
std::vector<JaehyunDDPoint> dd_points0;
UInt num_dd_pts;
UInt num_dd_pts0;
if (this->is_in_dd) {
LM_TOIMPLEMENT;
// UInt my_rank_in_DD = this->comm_group_dd.groupRank(lm_my_proc_id);
// UInt nbDDprocs = this->comm_group_dd.getNBprocsOnGroup();
// declare data array for (unclosed) dd points in a local processor
UInt num_dd_pts0_local = this->domA.getNBNodes();
std::vector<JaehyunDDPoint> dd_points0_local(num_dd_pts0_local);
// collecting dd points
UInt collecting = 0;
collecting = CollectingLocalDDpoints(dd_points0_local);
this->comm_group_dd.reduce(
&collecting, 1,
"boolean to detect any errors occurring in collecting local nodes",
OP_SUM);
// evaluate total number of points
UInt num_dd_pts0_global = 0;
LM_TOIMPLEMENT;
// MPI_Reduce(&num_dd_pts0_local, &num_dd_pts0_global, 1, MPI_INT, MPI_SUM,
// comm.realRank(0, ),
// comm.getMpiGroup(this->comm_group_dd));
// declare data array for (unclosed) total dd points only at 0 processor
std::vector<JaehyunDDPoint> dd_points0_global;
LM_TOIMPLEMENT;
// if (my_rank_in_DD == 0) {
// (unclosed) total dd points
// dd_points0_global.resize(num_dd_pts0_global);
//}
// declare data list for number of dd nodal points for all the DDD
// processors
std::vector<int> num_dd_pts0_local_list;
LM_TOIMPLEMENT;
// if (my_rank_in_DD == 0)
// num_dd_pts0_local_list.resize(nbDDprocs);
// gather the total local data length into the declared data list
LM_TOIMPLEMENT;
// int local_length = num_dd_pts0_local * local_data_length;
// MPI_Gather(&local_length, 1, MPI_INT, &num_dd_pts0_local_list[0], 1,
// MPI_INT, comm.realRank(0, this->comm_group_dd),
// comm.getMpiGroup(this->comm_group_dd));
// declear data array for starting number
std::vector<int> lengths_dd_points_local;
LM_TOIMPLEMENT;
// UInt totlen = 0;
// if (my_rank_in_DD == 0) {
// lengths_dd_points_local.resize(nbDDprocs);
// lengths_dd_points_local[0] = 0;
// totlen += num_dd_pts0_local_list[0];
// for (UInt i = 1; i < nbDDprocs; i++) {
// lengths_dd_points_local[i] = lengths_dd_points_local[i - 1] + totlen;
// totlen = num_dd_pts0_local_list[i];
// }
// }
// gather local dd points to global points
// inputs: dd_points0_local
// local_length
// num_dd_pts0_local_list
// lengths_dd_points_local
LM_TOIMPLEMENT;
// MPI_Gatherv(&dd_points0_local[0], // local dd points
// local_length, // local length
// MPI_REAL,
// &dd_points0_global[0], // global dd points
// &num_dd_pts0_local_list[0], // number of length list
// &lengths_dd_points_local[0], // starting number list
// MPI_REAL, comm.realRank(0, this->comm_group_dd),
// comm.getMpiGroup(this->comm_group_dd));
LM_TOIMPLEMENT;
int my_rank_in_DD = std::numeric_limits<int>::max();
if (my_rank_in_DD == 0) {
UInt num_pairs_in_infinity = this->domA.getNBpairsInfinite();
UInt num_dd_pts0 = num_dd_pts0_global;
if (this->close_in_infinity)
num_dd_pts0 += (num_pairs_in_infinity + num_close_loop_in_infinity) * 2;
// (closed) total dd nodal points
dd_points0.resize(num_dd_pts0);
close_dislocation_loop(dd_points0, dd_points0_global);
find_indices_for_neighbors(dd_points0);
num_dd_pts = num_dd_pts0;
// (closed) coarsened dd nodal points
dd_points.resize(num_dd_pts);
STARTTIMER("coarsening");
if ((coarsen_dd) and (collecting == 0)) {
bool coarsen = coarse_dislocation_loop(dd_points, dd_points0);
if (!coarsen) {
// in case of fail, copy all dd_points0 into dd_points
DUMP("[WARNING] Fail to coarse dislocation loop. Total "
<< num_dd_pts0 << " points",
DBG_MESSAGE);
for (UInt i = 0; i < num_dd_pts; ++i) {
for (UInt j = 0; j < data_length; ++j) {
dd_points[i * data_length + j] = dd_points0[i * data_length + j];
}
}
} else {
// in case of success,
DUMP("Coarse_dislocation_loop end. Reduced from "
<< num_dd_pts0 << " points to " << num_dd_pts << " points",
DBG_MESSAGE);
find_indices_for_neighbors(dd_points);
}
} else {
DUMP("[WARNING] Coarse_dislocation_loop is skipped. coarsen_dd: "
<< coarsen_dd << " collecting: " << collecting << " Total "
<< num_dd_pts0 << " points",
DBG_MESSAGE);
for (UInt i = 0; i < num_dd_pts; ++i) {
for (UInt j = 0; j < data_length; ++j) {
dd_points[i * data_length + j] = dd_points0[i * data_length + j];
}
}
}
STOPTIMER("coarsening");
ComputeOffsetParallelCurve(dd_points);
// CorrectIntersetingPoints(num_dd_pts,dd_points);
UInt mdNBprocs = this->comm_group_atomic.size();
for (UInt i = 0; i < mdNBprocs; ++i) {
// int real_rank = this->comm_group_atomic.realRank(i);
int my_rank_in_group = this->comm_group_atomic.getMyRank();
this->comm_group_atomic.send(&num_dd_pts, 1, my_rank_in_group,
"number of dd points");
this->comm_group_atomic.send(&num_dd_pts0, 1, my_rank_in_group,
"number of dd points0");
this->comm_group_atomic.send(&(dd_points[0]), num_dd_pts0 * data_length,
my_rank_in_group, "dd_points (vector)");
}
}
// BARRIER
this->comm_group_dd.synchronize();
}
if (this->is_in_atomic) {
STARTTIMER("SetPadAtoms");
SetPadAtoms();
STOPTIMER("SetPadAtoms");
STARTTIMER("ReceiveDDdata");
this->comm_group_dd.receive(&num_dd_pts, 1, 0, "number of dd points");
this->comm_group_dd.receive(&num_dd_pts0, 1, 0, "number of dd points0");
dd_points.resize(num_dd_pts0 * data_length);
this->comm_group_dd.receive(&(dd_points[0]), num_dd_pts0 * data_length, 0,
"dd_points (vector)");
STOPTIMER("ReceiveDDdata");
std::vector<UInt> indices;
// UInt count_pad_core = 0;
// UInt count_pad_around = 0;
// UInt count_far_field = 0;
// UInt count_true_disp = 0;
}
return dd_points;
}
/* ---------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::functionifAniso(
UInt num_dislo, Vector<Dim> &pos0_pbc, Vector<Dim> &pos0_pbc0,
std::vector<std::vector<Real>> &dislocations, RefAtom &at) {
Vector<Dim> analytic_Disp = {0.0, 0.0, 0.0};
for (UInt nd = 0; nd < num_dislo; ++nd) {
for (UInt pd = 0; pd < 2; ++pd) {
if (pd == 0) {
// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] - 1.0);
} else if (pd == 1) {
// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] + 1.0);
}
Real rotation_matrix[] = {0.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0};
Real pos0_pbc_rotate[] = {0.0, 0.0, 0.0};
for (UInt i = 0; i < 3; ++i) {
pos0_pbc_rotate[i] = 0.0;
for (UInt j = 0; j < 3; ++j) {
pos0_pbc_rotate[i] += rotation_matrix[3 * i + j] * pos0_pbc[j];
}
}
Real b_screw =
this->burgMag0 * cos(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;
Real b_edge =
this->burgMag0 * sin(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;
Real analytic_disp[] = {0.0, 0.0, 0.0};
Real dispV[] = {0.0, 0.0, 0.0};
analytic_disp[0] = computeBeDisplacementAniso(b_edge, pos0_pbc_rotate);
analytic_disp[1] = computeNDisplacementAniso(b_edge, pos0_pbc_rotate);
analytic_disp[2] = computeLDisplacementAniso(b_screw, pos0_pbc_rotate);
for (UInt i = 0; i < 3; ++i) {
for (UInt j = 0; j < 3; ++j) {
analytic_Disp[i] += rotation_matrix[3 * j + i] * analytic_disp[j];
}
}
for (UInt i = 0; i < 3; ++i)
analytic_Disp[i] += dispV[i];
}
}
applyDisplacement(at, analytic_Disp);
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::functionifVolterra(
UInt num_dislo, Vector<Dim> &pos0_pbc, Vector<Dim> &pos0_pbc0,
std::vector<std::vector<Real>> &dislocations, RefAtom &at) {
Vector<Dim> analytic_Disp = {0.0, 0.0, 0.0};
for (UInt nd = 0; nd < num_dislo; ++nd) {
for (UInt pd = 0; pd < 2; ++pd) {
if (pd == 0) {
// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] - 1.0);
} else if (pd == 1) {
// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] + 1.0);
}
Real rotation_matrix[] = {0.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0};
Real pos0_pbc_rotate[] = {0.0, 0.0, 0.0};
for (UInt i = 0; i < 3; ++i) {
pos0_pbc_rotate[i] = 0.0;
for (UInt j = 0; j < 3; ++j) {
pos0_pbc_rotate[i] += rotation_matrix[3 * i + j] * pos0_pbc[j];
}
}
Real b_screw =
this->burgMag0 * cos(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;
Real b_edge =
this->burgMag0 * sin(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;
Real analytic_disp[] = {0.0, 0.0, 0.0};
Real dispV[] = {0.0, 0.0, 0.0};
analytic_disp[0] = computeBeDisplacement(b_edge, pos0_pbc_rotate);
analytic_disp[1] = computeNDisplacement(b_edge, pos0_pbc_rotate);
analytic_disp[2] = computeLDisplacement(b_screw, pos0_pbc_rotate);
for (UInt i = 0; i < 3; ++i) {
for (UInt j = 0; j < 3; ++j) {
dispV[i] += rotation_matrix[3 * j + i] * analytic_disp[j];
}
}
for (UInt i = 0; i < 3; ++i)
analytic_Disp[i] += dispV[i];
}
}
applyDisplacement(at, analytic_Disp);
}
/* --------------------------------------------------------------------------
*/
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD, DomainMD, Dim>::functionElse(
std::vector<JaehyunDDPoint> &dd_points, Vector<Dim> &pos0_pbc,
RefAtom &at) {
UInt count_pad_core = 0;
UInt count_far_field = 0;
UInt count_pad_around = 0;
Vector<Dim> barnett_disp = {0.0, 0.0, 0.0};
Vector<Dim> dislo_disp = {0.0, 0.0, 0.0};
Vector<Dim> temp_disp = {0.0, 0.0, 0.0};
for (auto &dd_pt : dd_points) {
// bool selected = false;
Vector<Dim> disp_for_segment = {0.0, 0.0, 0.0};
auto &A = dd_pt.X;
auto &neigh1 = dd_pt.neighs[0];
auto &neigh2 = dd_pt.neighs[2];
auto &pts1 = neigh1.X;
auto &pts2 = neigh2.X;
for (UInt j = 0; j < 2; ++j) {
auto &neigh = dd_pt.neighs[j];
Vector<Dim> burg = neigh.burg * this->burgMag0;
Real bmag = burg[0] * cos(this->angle_between_b_and_x / 180.0 * M_PI) +
burg[2] * sin(this->angle_between_b_and_x / 180.0 * M_PI);
bmag = roundf(bmag * 10.0) / 10.0;
if (bmag / fabs(bmag) > 0) {
auto &B = neigh.X;
auto &normal = neigh.normal;
UInt slip_case = 3;
if (fabs(normal[0]) - 1.0 == 0.0)
slip_case = 0;
else if (fabs(normal[1]) - 1.0 == 0.0)
slip_case = 1;
else if (fabs(normal[2]) - 1.0 == 0.0)
slip_case = 2;
else {
DUMP("[warning] cannot define slip case", DBG_MESSAGE);
slip_case = 1;
}
ComputeAnalyticDisplacements(disp_for_segment, B, A, burg, pos0_pbc);
UInt idx_number = (UInt)neigh.idx;
auto &pts3 = dd_points[idx_number].tmp_X;
auto &pts4 = dd_points[idx_number].tmp2_X;
if (PointInSquare(pts1, pts2, pts3, pts4, pos0_pbc, slip_case)) {
Real distance = distance_between_line_and_point3D(pos0_pbc, A, B);
if (distance <= this->core_rcut) {
bool success(false);
success = ComputeTemplateDisplacements(temp_disp, A, B, burg,
normal, pos0_pbc);
if (!success) {
barnett_disp += disp_for_segment;
} else {
if (distance / this->core_rcut <= this->factor) {
++count_pad_core;
barnett_disp += temp_disp + disp_for_segment;
} else {
Real temp_ratio =
(this->core_rcut - distance) /
(this->core_rcut - this->core_rcut * this->factor);
++count_pad_around;
barnett_disp += temp_disp * temp_ratio + disp_for_segment;
}
}
} else {
barnett_disp += disp_for_segment;
++count_far_field;
}
} else {
barnett_disp += disp_for_segment;
++count_far_field;
} // if(PointInSquare(pts ... )
} // if(bmag/abs(bmag)>0)
} // for (UInt j=0; j < num_neigh; ++j){
} // for (UInt i=0; i < num_dd_pts;
dislo_disp = barnett_disp;
applyDisplacement(at, dislo_disp);
}
/* -------------------------------------------------------------------------- */
DECLARE_DD_ATOMIC_TEMPLATE(CADD_template);
__END_LIBMULTISCALE__
diff --git a/src/filter/compute_cadd_template.cc b/src/filter/compute_cadd_template.cc
index d4a7410..26999b2 100644
--- a/src/filter/compute_cadd_template.cc
+++ b/src/filter/compute_cadd_template.cc
@@ -1,429 +1,429 @@
#include "compute_cadd_template.hh"
#include "factory_multiscale.hh"
#include "lib_dd.hh"
#include "lib_dumper.hh"
#include "lib_filter.hh"
#include "lib_md.hh"
#include "lib_stimulation.hh"
#include "lm_common.hh"
#include <Eigen/Dense>
#include <algorithm>
#include <fstream>
#include <pybind11/eigen.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
void ComputeCADDTemplate::loadNPZ(std::string &filename) {
py::object load = py::module::import("numpy").attr("load");
auto npz = load(filename);
auto field_names = npz.attr("files").cast<std::list<std::string>>();
auto fields = npz.cast<py::dict>().cast<std::map<std::string, py::object>>();
// std::for_each(field_names.begin(), field_names.end(),
// [](auto &&v) { std::cout << v << std::endl; });
auto triangles =
fields["triangles"].cast<Eigen::Array<UInt, Eigen::Dynamic, 3>>();
auto points = fields["points"].cast<Eigen::Array<Real, Eigen::Dynamic, 2>>();
auto corrections =
fields["correction"].cast<Eigen::Array<Real, Eigen::Dynamic, 3>>();
auto &nodes = mesh.getContainerNodes();
RefPointData<2> p;
for (UInt i = 0; i < points.rows(); ++i) {
p = VectorView<2>(&points(i, 0));
- nodes.add(p);
+ nodes.push_back(p);
}
for (UInt i = 0; i < corrections.rows(); ++i) {
auto corr = VectorView<3>(&corrections(i, 0));
this->corrections.push_back(corr);
}
auto &elems = mesh.getContainerElems();
for (UInt i = 0; i < triangles.rows(); ++i) {
RefGenericElem<2> el;
for (UInt n = 0; n < 3; ++n)
el.conn.push_back(triangles(i, n));
- elems.add(el);
+ elems.push_back(el);
}
}
/* -------------------------------------------------------------------------- */
ComputeCADDTemplate::ComputeCADDTemplate(const LMID &name)
: LMObject(name), mesh(this->getID() + ":mesh") {
this->removeInput("input1");
this->createInput("points");
this->createInput("segments");
this->inversion = TOP_RIGHT;
this->search_dim = 1;
this->tol = -1;
this->r_cut = -1;
this->factor[0] = this->factor[1] = 1;
};
/* -------------------------------------------------------------------------- */
ComputeCADDTemplate::~ComputeCADDTemplate() {}
/* -------------------------------------------------------------------------- */
inline bool is_in_triangle(Vector<2> &a, Vector<2> &b, Vector<2> &c,
Vector<2> p) {
Vector<2> v0;
Vector<2> v1;
Vector<2> v2;
Real u, v;
v0 = c - a;
v1 = b - a;
v2 = p - a;
// Compute dot products
Real dot00 = v0.dot(v0);
Real dot01 = v0.dot(v1);
Real dot02 = v0.dot(v2);
Real dot11 = v1.dot(v1);
Real dot12 = v1.dot(v2);
// Compute barycentric coordinates
Real invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01);
u = (dot11 * dot02 - dot01 * dot12) * invDenom;
v = (dot00 * dot12 - dot01 * dot02) * invDenom;
// Check if point is in triangle
return (u >= -1e-8) && (v >= -1e-8) && (u + v < 1.0 + 1e-8);
}
/* -------------------------------------------------------------------------- */
// returns the number of the element an atom is in
// or -1 if it's outside the template
inline int ComputeCADDTemplate::findElement(Vector<Dim> &at_coords) {
auto &nodes = this->mesh.getContainerNodes();
int i = 0;
for (auto &&elem : this->mesh.getContainerElems()) {
auto node1 = nodes[elem.localIndexes()[0]].position0();
auto node2 = nodes[elem.localIndexes()[1]].position0();
auto node3 = nodes[elem.localIndexes()[2]].position0();
if (is_in_triangle(node1, node2, node3, at_coords.block<2, 1>(0, 0))) {
return i;
}
++i;
}
return UINT_MAX;
}
/* -------------------------------------------------------------------------- */
void ComputeCADDTemplate::init() { loadNPZ(filename); }
/* -------------------------------------------------------------------------- */
/* LMDESC CADD_TEMPLATE
TODO
*/
/* LMEXAMPLE COMPUTE template CADD_TEMPLATE INPUT md TODO
*/
void ComputeCADDTemplate::declareParams() {
ComputeInterface::declareParams();
/* LMKEYWORD FILENAME
Set the filename containing a core template
*/
this->parseKeyword("FILENAME", filename);
// /* LMKEYWORD ANGLE_BETWEEN_B_AND_X
// Specifies the angle between Burgers vector and X(11-2) axis
// */
// this->parseKeyword("ANGLE_BETWEEN_B_AND_X", angle_between_b_and_x);
}
/* -------------------------------------------------------------------------- */
void ComputeCADDTemplate::coreDispAtom(Vector<Dim> &pos_core,
Vector<Dim> &disp) {
int elem_num = this->findElement(pos_core);
if (elem_num != -1) {
Vector<3> coords = {1.0, pos_core[0], pos_core[1]};
auto &elems = this->mesh.getContainerElems();
auto &nodes = this->mesh.getContainerNodes();
auto &elem = elems[elem_num];
auto &node1 = nodes[elem.localIndexes()[0]].position0();
auto &node2 = nodes[elem.localIndexes()[1]].position0();
auto &node3 = nodes[elem.localIndexes()[2]].position0();
Eigen::Matrix<Real, 4, 3> Su;
Su << 1., 0., 0., node1, node2, node3;
Matrix<3> SxInv;
auto &&tmp = Su * SxInv * coords;
disp += tmp.block<3, 1>(1, 0);
}
}
/* -------------------------------------------------------------------------- */
template <typename Points, typename Segment>
std::enable_if_t<Points::Dim == 3u>
ComputeCADDTemplate::build(Points &points, Segment &segment) {
for (auto &&point : points) {
computeTemplateDisplacements(point, segment);
}
}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2>
Matrix<3u> computeTemplateRotation(Vec1 &&_dd_line_dir,
Vec2 &&_slip_normal_dir) {
Vector<3u> dd_line_dir(_dd_line_dir);
Vector<3u> slip_normal_dir(_slip_normal_dir);
// construct a rotation where:
// - the normal to slip plane becomes z
// - the line direction becomes y
// - normal to the line in the slip plane becomes x
// checks is the normal is normal to the line direction
LM_ASSERT(std::abs(dd_line_dir.dot(slip_normal_dir)) < 1e-10,
"normal and line not defined properly");
// construct the first basis vector
Vector<3u> first_basis_vector = dd_line_dir.cross(slip_normal_dir);
Matrix<3u> Rmat = Matrix<3>::Zero();
Rmat.block<3, 1>(0, 0) = first_basis_vector;
Rmat.block<3, 1>(0, 1) = dd_line_dir;
Rmat.block<3, 1>(0, 2) = slip_normal_dir;
return Rmat.transpose();
}
/* -------------------------------------------------------------------------- */
template <typename Point, typename ElementType>
std::enable_if_t<!std::is_same<ElementType, RefGenericDDElem<3>>::value>
ComputeCADDTemplate::computeTemplateDisplacements(Point &p,
ElementType &segment) {
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template <typename Point, typename ElementType>
std::enable_if_t<std::is_same<ElementType, RefGenericDDElem<3>>::value>
ComputeCADDTemplate::computeTemplateDisplacements(Point &p,
ElementType &segment) {
auto &&burgers = segment.burgers();
auto &&normal = segment.normal();
LM_ASSERT(std::abs(burgers.norm()) > 1e-8,
"No found burgers vector for neighbors");
auto &&pos0 = p.position0();
auto conn = segment.localIndexes();
auto node1 = segment.getNode(0).position();
auto node2 = segment.getNode(1).position();
// computes the direction of the line
auto dd_line = node2 - node1;
// computes a unitary vector in the direction of the line
Vector<Dim> dd_line_direction = dd_line;
dd_line_direction.normalize();
// computes the relatice position of the atom
Vector<Dim> pos_atom = pos0 - node1;
// Rmat is a rotation of only the coordinates in the slip plane coordinates
// This is to be used to convert to/from the reference where the templates
// were computed
Matrix<Dim> Rmat = computeTemplateRotation(dd_line_direction, normal);
// project the atom position to the segment
Vector<Dim> projected_atom_pos =
pos_atom.dot(dd_line_direction) * dd_line_direction;
// vector from segment to atom position
Vector<Dim> vec_atom_2_dd_core = pos_atom - projected_atom_pos;
Vector<Dim> template_coords = Rmat * vec_atom_2_dd_core;
Vector<Dim> disp = Vector<Dim>::Zero();
this->coreDispAtom(template_coords, disp);
for (UInt i = 0; i < Dim; ++i)
- this->add(disp[i]);
+ this->push_back(disp[i]);
}
/* -------------------------------------------------------------------------- */
//****************
// part concrening the interpolation between angles
//****************
// character angle: angle between burgers vector and dislocation line
// Real dangle = fabs(burgers.dot(dd_line_direction)) / burgers.norm();
// Real dangle_deg =
// std::acos(dangle) * 180.0 / M_PI; // character angle (degree)
// Real distance = fabs(90.0 - dangle_deg);
// Real cangle = distance + 90.0; // change in a range between 90 and 180
// Real closest_angles[2] = {0.0, 0.0};
// UInt indices_closest_angles[2] = {1, 1};
// compute_angle_ratio(cangle, indices_closest_angles, closest_angles);
// DO NOT DELETE: displace elsewhere !!!!!
// contained the displacement of the two templates
// involved for angle interpolation
// std::vector<Vector<Dim>> d = {0.0, 0.0};
// DO NOT DELETE: displace elsewhere !!!!!
// loop two neighboring angle templates for interpolation
// for (int j = 0; j < 2; ++j) {
// Real tmp_angle = fabs(this->angle_between_b_and_x - cangle);
// tdisp = Rmat.transpose() * disp_template0;
// tdisp[0] = disp_template0[2] * cos(tmp_angle / 180.0 * M_PI) +
// disp_template0[0] * sin(tmp_angle / 180.0 * M_PI);
// tdisp[1] = disp_template0[1];
// tdisp[2] = disp_template0[0] * cos(tmp_angle / 180.0 * M_PI) -
// disp_template0[2] * sin(tmp_angle / 180.0 * M_PI);
// DO NOT DELETE: displace elsewhere !!!!!
// in between angles interpolation
// tdisp[0] = (dx[1] - dx[0]) / (nei_angles[1] - nei_angles[0]) *
// (cangle - nei_angles[0]) +
// dx[0];
// tdisp[1] = (dy[1] - dy[0]) / (nei_angles[1] - nei_angles[0]) *
// (cangle - nei_angles[0]) +
// dy[0];
// tdisp[2] = (dz[1] - dz[0]) / (nei_angles[1] - nei_angles[0]) *
// (cangle - nei_angles[0]) +
// dz[0];
// /* --------------------------------------------------------------------------
// */
// template <typename DomainDD, typename DomainMD, UInt Dim>
// bool CADD_template<DomainDD, DomainMD, Dim>::PointInSquare(
// Vector<Dim> &pts1, Vector<Dim> &pts2, Vector<Dim> &pts3, Vector<Dim>
// &pts4,
// Vector<Dim> &pts, UInt slip_case) {
// Real area0 = 0.0;
// Real area1 = AreaTriangle(pts1, pts3, pts4, slip_case);
// Real area2 = AreaTriangle(pts4, pts2, pts1, slip_case);
// Real area3 = AreaTriangle(pts1, pts2, pts3, slip_case);
// Real area4 = AreaTriangle(pts2, pts3, pts4, slip_case);
// area0 = (fabs(area1) + fabs(area2) + fabs(area3) + fabs(area4)) / 2.;
// Vector<2> int_X;
// if (intersectTwoD(pts1, pts3, pts2, pts4, int_X, slip_case)) {
// Real temp[3] = {0.0, 0.0, 0.0};
// for (UInt i = 0; i < 3; ++i) {
// temp[i] = pts2[i];
// pts2[i] = pts1[i];
// pts1[i] = temp[i];
// }
// }
// Real tot_area = 0.0;
// Real tot_area1 = AreaTriangle(pts, pts1, pts2, slip_case);
// Real tot_area2 = AreaTriangle(pts, pts3, pts1, slip_case);
// Real tot_area3 = AreaTriangle(pts, pts4, pts3, slip_case);
// Real tot_area4 = AreaTriangle(pts, pts2, pts4, slip_case);
// tot_area =
// fabs(tot_area1) + fabs(tot_area2) + fabs(tot_area3) + fabs(tot_area4);
// if (area0 + 1.0 < tot_area)
// return false;
// else
// return true;
// }
// template <typename _Input> void ComputeCADDTemplate<_Input>::init() {
// if (this->is_in_atomic) {
// for (UInt i = 0; i < this->num_templates; ++i) {
// // TwoDTemplate<RefAtom, IteratorAtomsSubset> *single_template =
// // new TwoDTemplate<RefAtom, IteratorAtomsSubset>();
// // this->core_templates.push_back(single_template);
// // this->core_templates[i]->init(this->templatefilename_list[i]);
// // }
// LM_TOIMPLEMENT;
// }
// }
// this function retreives the two closest angles
// this is done to make the interpolation between angles easier
// TODO: review move elsewhere
inline void compute_angle_ratio(Real angle, UInt *indices_closest_angles,
Real *closest_angles) {
// TO REVIEW;
// to enforce the specific core template
// if(this->use_specific_template) {
// angle = this->specific_template_angle;
// } else {
// if (angle < 90.0) {
// angle = 180.0 - angle;
// } else if ((180.0 < angle) && (angle <= 270.0)) {
// angle = 360.0 - angle;
// } else if ((270.0 < angle) && (angle <= 360.0)) {
// angle = angle - 180.0;
// }
// // }
// for (UInt j = 0; j < this->num_templates - 1; ++j) {
// two_angles[0] = this->core_templates[j]->get_theta();
// two_angles[1] = this->core_templates[j + 1]->get_theta();
// two_indices[0] = j;
// two_indices[1] = j + 1;
// if ((angle - two_angles[0]) * (angle - two_angles[1]) <= 0.0) {
// break;
// }
// }
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void ComputeCADDTemplate::compute_make_call() {
DUMP(this->getID() << ": Compute", DBG_INFO);
call_compute(*this, [&](auto &... args) { this->build(args...); },
this->getInput("points"), this->getInput("segments"));
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
diff --git a/src/filter/compute_dislodetection.cc b/src/filter/compute_dislodetection.cc
index 7136b0d..dd5b468 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_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);
+ 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.add(*nodes_it);
+ 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.add(*nodes_it);
+ 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.add(el_buffer);
+ 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.add(pt_buffer);
+ 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.add(el_buffer);
+ 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(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__
diff --git a/src/filter/filter_point.cc b/src/filter/filter_point.cc
index 5ed9059..e76a9c3 100644
--- a/src/filter/filter_point.cc
+++ b/src/filter/filter_point.cc
@@ -1,302 +1,302 @@
/**
* @file filter_point.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Tue Jul 29 17:03:26 2014
*
* @brief This allows to filter a set of DOFs by matching a list of coordinates
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "filter_point.hh"
#include "communicator.hh"
#include "compute_compatibility.hh"
#include "compute_interface.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "ref_point_data.hh"
#include <fstream>
#include <list>
#include <sstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename Ref>
inline bool findPointInSet(Ref &point, std::list<UInt> &indices,
const std::vector<Real> &positions,
const Real &tolerance) {
//static const UInt Dim = Ref::Dim;
std::list<UInt>::iterator index = indices.begin();
std::list<UInt>::iterator end_index = indices.end();
Real distance = 0;
do {
LM_TOIMPLEMENT;
//distance = point.computeInitialDistance(&positions.at(Dim * (*index)));
++index;
} while ((distance > tolerance) && (index != end_index));
if (index == end_index) {
return false;
} else {
indices.erase(--index);
return true;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void setPositionsCopy(const std::vector<Real> &positions, Cube cube,
Real tolerance, std::vector<Real> &positions_copy) {
auto &&xmax = cube.getXmax();
auto &&xmin = cube.getXmin();
for (UInt k = 0; k < Dim; ++k) {
xmax[k] += 2. * tolerance;
xmin[k] -= 2. * tolerance;
}
cube.setDimensions(xmin.cast<PhysicalScalar<Length>>(),
xmax.cast<PhysicalScalar<Length>>());
UInt _n = positions.size() / Dim;
for (UInt i = 0; i < _n; ++i) {
Vector<Dim> tmp(&positions[i * Dim]);
bool test = static_cast<Geometry &>(cube).template contains<Dim>(tmp);
if (test) {
for (UInt k = 0; k < Dim; ++k)
positions_copy.push_back(positions[i * Dim + k]);
}
}
}
/* -------------------------------------------------------------------------- */
UInt FilterPoint::getNbSharedPoints() {
LM_TOIMPLEMENT;
// UInt counter = 0;
// auto it = this->begin(this->dt);
// auto end = this->end(this->dt);
// for (; it != end; ++it) {
// auto &&point = *it;
// if ((point.getDOFType() & dt_local_master) == false) {
// counter++;
// }
// }
// return counter;
}
/* -------------------------------------------------------------------------- */
template <typename Cont> void FilterPoint::build(Cont &cont) {
auto &output =
allocAndGetCastedOutput<typename Cont::ContainerSubset>(this->getID());
std::vector<Real> positions_copy;
static const UInt Dim = Cont::Dim;
Cube cube(Dim);
auto it = cont.begin(this->dt);
auto end = cont.end(this->dt);
for (; it != end; ++it) {
auto &&point = *it;
auto &&pos = point.position0();
cube.extendBoundingBox(pos);
}
setPositionsCopy<Dim>(this->positions, cube, this->tolerance, positions_copy);
DUMP(this->getID() << ": positions_copy " << positions_copy.size() << " "
<< positions.size(),
DBG_INFO);
output.clear();
std::list<UInt> indices;
for (UInt i = 0; i < Dim; ++i) {
// This is a dummy point I add for the findPointInSet algo, cleaned later
positions_copy.push_back(lm_real_max);
}
UInt nb_points = positions_copy.size() / Dim;
for (UInt i = 0; i < nb_points; ++i) {
indices.push_back(i);
}
DUMP("number of indices = " << indices.size(), DBG_INFO);
it = cont.begin(this->dt);
for (; it != end; ++it) {
auto &&point = *it;
if (findPointInSet(point, indices, positions_copy, this->tolerance)) {
- output.add(point);
+ output.push_back(point);
}
}
// remove dummy point
indices.pop_back();
for (UInt i = 0; i < Dim; ++i) {
positions_copy.pop_back();
}
nb_points = output.size();
if (this->dt != dt_local_master) {
UInt shared = this->getNbSharedPoints();
DUMP("locally found " << nb_points << " total of which " << shared
<< " are shared.",
DBG_DETAIL);
nb_points -= shared;
}
// Communicator &comm = Communicator::getCommunicator();
CommGroup group = cont.getCommGroup();
group.allReduce(&nb_points, 1, "DeviationOperator reduce res", OP_SUM);
if (nb_points == 0) {
LM_FATAL(this->getID() << ": Could not match any point out of "
<< positions.size() / Dim);
}
if (nb_points != positions.size() / Dim) {
std::stringstream sstr;
sstr << this->getID() << "-found_pos-" << lm_my_proc_id;
ComputeCompatibility found_positions(sstr.str());
found_positions.name_computed.resize(Dim);
for (UInt k = 0; k < Dim; ++k) {
std::stringstream sstr;
sstr << "pos" << k;
found_positions.name_computed[k] = sstr.str();
}
for (UInt i = 0; i < output.getArray().size(); ++i) {
auto &&point = output.getArray()[i];
auto &&pos = point.position0();
found_positions.push_back(pos[0]);
found_positions.push_back(pos[1]);
found_positions.push_back(pos[2]);
}
found_positions.setCommGroup(group);
LM_TOIMPLEMENT;
// DumperText error_positions(sstr.str(), found_positions);
// error_positions.setParam("SEPARATOR", std::string(","));
// error_positions.setParam("RANK", comm.groupRank(lm_my_proc_id, group));
// error_positions.manualDump();
LM_FATAL(this->getID() << ": Could match " << nb_points - 1 << " out of "
<< positions.size() / Dim << " positions"
<< " locally I found " << output.size()
<< " points");
}
dimension = Dim;
}
/* -------------------------------------------------------------------------- */
FilterPoint::FilterPoint(const LMID &name)
: LMObject(name), path(""), tolerance(1e-5), dimension(0) {}
/* -------------------------------------------------------------------------- */
FilterPoint::~FilterPoint() {}
/* -------------------------------------------------------------------------- */
void FilterPoint::init() {
// yeah, could be done more elegantly
std::ifstream infile(this->path.c_str());
if (!infile.good()) {
LM_FATAL("The atom file " << path.c_str()
<< " could NOT be opened, rdstate = " << std::hex
<< infile.rdstate() << std::endl
<< " badbit = " << infile.bad() << std::endl
<< " failbit = " << infile.fail() << std::endl
<< " eofbit = " << infile.eof() << std::endl
<< " goodbit = " << infile.good() << std::endl);
}
char tmp;
while (infile.good()) {
tmp = infile.peek();
switch (tmp) {
case '#':
infile.ignore(1024, '\n');
break;
default:
while (infile.good()) {
this->positions.push_back(0);
infile >> this->positions.back();
}
break;
}
}
if (this->positions.size() == 0)
LM_FATAL("no valid positions found in file " << this->path.c_str());
this->positions.pop_back();
DUMP("found " << this->positions.size() / 3 << " points", DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
UInt FilterPoint::getNbPointsInFile() {
return this->positions.size() / dimension;
}
/* -------------------------------------------------------------------------- */
/* LMDESC POINT
This filter allows to extract a set of points and/or
elements contained in a provided point
*/
/* LMHERITANCE filter */
/* LMEXAMPLE FILTER fpoint POINT INPUT md PATH ./my-point-list.txt TOL 1e-5 */
void FilterPoint::declareParams() {
FilterInterface::declareParams();
/* LMKEYWORD PATH
Specify the file containing the points to filter
*/
this->parseKeyword("PATH", path);
/* LMKEYWORD TOL
Specify the file containing the points to filter
*/
this->parseKeyword("TOL", tolerance, 1e-5);
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_MAKE_CALL(FilterPoint);
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__