Page MenuHomec4science
Files F72840424

compute_barnett.cc
No OneTemporary
Actions

File Metadata

Created
Wed, Jul 17, 07:57

compute_barnett.cc
View Options

/**
* @file compute_dd_displacement.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Moseley Philip Arthur <philip.moseley@epfl.ch>
*
* @date Thu Nov 21 22:11:13 2013
*
* @brief This compute allows to compute the displacement field on a set of
* points due to dislocation segments
*
* @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 "compute_barnett.hh"
#include "domain_dd_interface.hh"
#include "domain_multiscale.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "lm_communicator.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
ComputeBarnett::ComputeBarnett(const LMID &id) : LMObject(id) {
this->removeInput("input1");
this->createInput("segments");
this->createInput("points");
this->createArrayOutput("displacement");
}
/* -------------------------------------------------------------------------- */
ComputeBarnett::~ComputeBarnett() {}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2, typename Vec3, typename Vec4>
Vector<3> computeClosurePoint(Vec1 &&A, Vec2 &&B, Vec3 &&b, Vec4 &&slip_plane) {
LM_ASSERT(slip_plane.dot(b) < 1e-10,
"invalid slip plane: " << slip_plane.dot(b));
// Marc Fivel & Christophe Depres,
// Philosophical Magazine, Volume 94, Issue 28, 2014
// The closure point is defined as the intersection point
// between arbitrary vector of u to the current slip plane.
// The slip plane can be defined by cross product of burgers vector
// and vector of AB
Vector<3> O = {0.0, 0.0, 0.0};
Vector<3> u = {1.0, 1.0, 1.0}; // arbitrary vector from the origin point O
// http: // mathworld.wolfram.com/Line-PlaneIntersection.html
Vector<3> slip_normal = (B - A).cross(b);
Vector<3> slip_direction;
// print slip_normal
// print np.linalg.norm(slip_normal)
if (slip_normal.norm() < 1e-15) {
slip_direction = b.cross(slip_plane);
B = A + slip_direction;
}
Vector<3> X1 = A + b;
Vector<3> X2 = A;
Vector<3> X3 = B;
Vector<3> X4 = O;
Vector<3> X5 = u;
Matrix<4> nom_mat = Matrix<4>::Zero();
nom_mat.block<1, 4>(0, 0).array() = 1.;
nom_mat.block<3, 1>(1, 0).array() = X1;
nom_mat.block<3, 1>(1, 1).array() = X2;
nom_mat.block<3, 1>(1, 2).array() = X3;
nom_mat.block<3, 1>(1, 3).array() = X4;
Real nom = nom_mat.determinant();
Matrix<4> &denom_mat = nom_mat;
denom_mat(0, 3) = 0.;
denom_mat.block<3, 1>(1, 3).array() = X5 - X4;
Real denom = denom_mat.determinant();
// print denom_mat
// print denom
Real t = -nom / denom;
auto C = X4 + (X5 - X4) * t;
return C;
}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2, typename Vec3, typename Vec4>
Real SolidAngleSegmentBarnett(Vec1 &vec_la, Vec2 &vec_lb, Vec3 &vec_lc,
Vec4 &N) {
DUMP("vec_la:" << vec_la, DBG_DETAIL);
DUMP("vec_lb:" << vec_lb, DBG_DETAIL);
DUMP("vec_lc:" << vec_lc, DBG_DETAIL);
auto correct_cos = [&](Real &&val) -> Real {
if (val > 1.)
return 1.;
if (val < -1.)
return -1;
return val;
};
auto val_a = std::acos(correct_cos(vec_lb.dot(vec_lc)));
auto val_b = std::acos(correct_cos(vec_la.dot(vec_lc)));
auto val_c = std::acos(correct_cos(vec_la.dot(vec_lb)));
auto val_s = .5 * (val_a + val_b + val_c);
DUMP("val_a:" << val_a, DBG_DETAIL);
DUMP("val_b:" << val_b, DBG_DETAIL);
DUMP("val_c:" << val_c, DBG_DETAIL);
DUMP("val_s:" << val_s, DBG_DETAIL);
DUMP("val_a:" << val_a << " " << std::tan(0.5 * val_a), DBG_DETAIL);
DUMP("val_s - val_a:" << val_s - val_a << " "
<< std::tan(0.5 * (val_s - val_a)),
DBG_DETAIL);
DUMP("val_s - val_b:" << val_s - val_b << " "
<< std::tan(0.5 * (val_s - val_b)),
DBG_DETAIL);
DUMP("val_s - val_c :" << val_s - val_c << " "
<< std::tan(0.5 * (val_s - val_c)),
DBG_DETAIL);
auto tan2e_4 =
(std::tan(0.5 * val_s) * std::tan(0.5 * (val_s - val_a)) *
std::tan(0.5 * (val_s - val_b)) * std::tan(0.5 * (val_s - val_c)));
if (std::abs(tan2e_4) < 1e-12)
tan2e_4 = 0.;
LM_ASSERT(tan2e_4 >= 0., "internal error: " << tan2e_4);
DUMP("tan2e_4:" << tan2e_4, DBG_DETAIL);
DUMP("std::sqrt(tan2e_4):" << std::sqrt(tan2e_4), DBG_DETAIL);
Real val_e = 4.0 * std::atan(std::sqrt(tan2e_4));
if (val_e > 2.0 * M_PI)
val_e = 0.;
if (val_e < 0.0)
val_e = 0.;
DUMP("val_e:" << val_e, DBG_DETAIL);
DUMP("la dot N:" << vec_la.dot(N), DBG_DETAIL);
Real res = -1. * std::copysign(1., vec_la.dot(N)) * val_e;
DUMP("omega_abc:" << res, DBG_DETAIL);
return res;
}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2, typename Vec3>
auto computeProjectedPoint(Vec1 &&A, Vec2 &&C, Vec3 &&burg) {
// Marc Fivel &Christophe Depres,
// Philosophical Magazine, Volume 94, Issue 28, 2014
auto X1 = A + burg * 1e+8;
auto X2 = A - burg * 1e+8;
auto X0 = C;
auto AC = X0 - X1;
auto AB = X2 - X1;
auto n = AB / AB.norm();
Real dotACn = AC.dot(n);
auto P = X1 + dotACn * n;
return P;
}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2, typename Vec3, typename Vec4,
typename Vec5>
Vector<3> segmentDisplacementsOfAtom(Vec1 &&A, Vec2 &&B, Vec3 &&C, Vec4 &&X,
Vec5 &&burg, Real poisson) {
// computes the displacement field at point X of
// a triangular loop made of 3 segments AB, BC, CA
// burg: the burgers vector
// pois: the poisson ratio
DUMP("A: " << A, DBG_DETAIL);
DUMP("B: " << B, DBG_DETAIL);
DUMP("C: " << C, DBG_DETAIL);
Vector<3> vec_ab = B - A;
Vector<3> vec_bc = C - B;
Vector<3> vec_ca = C - A;
DUMP("AB: " << vec_ab, DBG_DETAIL);
DUMP("BC: " << vec_bc, DBG_DETAIL);
DUMP("CA: " << vec_ca, DBG_DETAIL);
Real ab_norm = vec_ab.norm();
Real bc_norm = vec_bc.norm();
Real ca_norm = vec_ca.norm();
// Zero length segment; nothing to do.
if (ab_norm < 1e-20 or bc_norm < 1e-20 or ca_norm < 1e-20)
return Vector<3>::Zero();
// compute t_ab, t_bc, t_ca
auto t_ab = vec_ab / ab_norm;
auto t_bc = vec_bc / bc_norm;
auto t_ca = vec_ca / ca_norm;
// compute R_A, R_B, R_C
Vector<3> vec_ra = A - X;
auto vec_rb = B - X;
auto vec_rc = C - X;
DUMP("vec_ra:" << vec_ra, DBG_DETAIL);
DUMP("vec_rb:" << vec_rb, DBG_DETAIL);
DUMP("vec_rc:" << vec_rc, DBG_DETAIL);
// compute lambda_A, lambda_B, lambda_C
auto ra_norm = vec_ra.norm();
auto rb_norm = vec_rb.norm();
auto rc_norm = vec_rc.norm();
DUMP("ra_norm:" << ra_norm, DBG_DETAIL);
DUMP("rb_norm:" << rb_norm, DBG_DETAIL);
DUMP("rc_norm:" << rc_norm, DBG_DETAIL);
auto vec_la = vec_ra / ra_norm;
auto vec_lb = vec_rb / rb_norm;
auto vec_lc = vec_rc / rc_norm;
DUMP("vec_la:" << vec_la, DBG_DETAIL);
DUMP("vec_lb:" << vec_lb, DBG_DETAIL);
DUMP("vec_lc:" << vec_lc, DBG_DETAIL);
// Calculate the solid angle.
// There will be no contribution to solid angle if AB and BC are colinear.
auto N = vec_ab.cross(vec_bc);
auto n_norm = N.norm();
Real omega_abc;
if (n_norm < 1e-20)
omega_abc = 0.;
else {
N /= n_norm;
omega_abc = SolidAngleSegmentBarnett(vec_la, vec_lb, vec_lc, N);
}
DUMP("omega_abc:" << omega_abc, DBG_DETAIL);
// Calculate Barnett
// fAB = (b x tAB) ln[(Rb / Ra).((1 + lB.tAB) / (1 + lA.tAB))]
auto c_ab = std::log(rb_norm / ra_norm *
((1 + vec_lb.dot(t_ab)) / (1 + vec_la.dot(t_ab))));
auto c_bc = std::log(rc_norm / rb_norm *
((1 + vec_lc.dot(t_bc)) / (1 + vec_lb.dot(t_bc))));
auto c_ca = std::log(ra_norm / rc_norm *
((1 + vec_la.dot(t_ca)) / (1 + vec_lc.dot(t_ca))));
auto f_ab = c_ab * burg.cross(t_ab);
auto f_bc = c_bc * burg.cross(t_bc);
auto f_ca = c_ca * burg.cross(t_ca);
// Calculate Barnett gAB = [b.(lA x lB)](lA + lB) / (1 + lA.lB)
auto g_ab = burg.dot(vec_la.cross(vec_lb)) / (1 + vec_la.dot(vec_lb)) *
(vec_la + vec_lb);
auto g_bc = burg.dot(vec_lb.cross(vec_lc)) / (1 + vec_lb.dot(vec_lc)) *
(vec_lb + vec_lc);
auto g_ca = burg.dot(vec_lc.cross(vec_la)) / (1 + vec_lc.dot(vec_la)) *
(vec_lc + vec_la);
// Update the displacements.
auto coeff1 = omega_abc / (4.0 * M_PI);
DUMP(burg, DBG_DETAIL);
auto coeff3 = 1. / (8. * M_PI * (1. - poisson));
auto coeff2 = (1.0 - 2.0 * poisson) * coeff3;
DUMP("coeff1 = " << coeff1, DBG_DETAIL);
DUMP("coeff2 = " << coeff2, DBG_DETAIL);
DUMP("coeff3 = " << coeff3, DBG_DETAIL);
// U = -c2 * burg - c2 * (fAB + fBC + fCA) + c3 * (gAB + gBC + gCA)
auto U = (-coeff1 * burg - coeff2 * (f_ab + f_bc + f_ca) +
coeff3 * (g_ab + g_bc + g_ca));
DUMP("U" << U, DBG_DETAIL);
return U;
}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2, typename Vec3, typename Vec4,
typename Vec5>
Vector<3> computeAnalyticDisplacements(Vec1 &&A, Vec2 &&B, Vec3 &&X,
Vec4 &&burg, Vec5 &&slip_plane,
Real poisson) {
// Compute closure point for given AB segment
auto &&C = computeClosurePoint(A, B, burg, slip_plane);
DUMP("ClosurePoint=" << C << " " << A << " " << B, DBG_DETAIL);
// fivel's original method
// Compute point P1, projected point of C on line[A, A + b]
auto P1 = computeProjectedPoint(A, C, burg);
DUMP("P1=" << P1 << " " << A << " " << B, DBG_DETAIL);
// Compute point P2, projected point of C on line[B, B + b]
auto P2 = computeProjectedPoint(B, C, burg);
DUMP("P2=" << P2 << " " << A << " " << B, DBG_DETAIL);
//
// Compute displacement from the 2 triangles identified
auto U = segmentDisplacementsOfAtom(P2, A, B, X, burg, poisson);
U += segmentDisplacementsOfAtom(P1, A, P2, X, burg, poisson);
// barnett's original method
// U = segmentDisplacementsOfAtom(C, A, B, X, burg, pois)
if (isNaN(U) or isInf(U))
LM_FATAL("computed an invalid value a position: " << X << " " << U);
DUMP("For X=" << X << " U=" << U, DBG_DETAIL);
return U;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerDD, typename Vec1, typename Vec2>
Vector<3> computeDisloDisplacements(ContainerDD &edges, Vec1 &&X,
Vec2 &&slip_plane, Real poisson) {
auto &nodes = edges.getContainerNodes();
Vector<3> disp = Vector<3>::Zero();
for (auto &elem : edges.getContainerElems()) {
auto &&nds = elem.nodes();
Vector<3> A = nodes.get(nds[0]).position();
Vector<3> B = nodes.get(nds[1]).position();
Vector<3> burgers = elem.burgers();
disp += computeAnalyticDisplacements(A, B, X, burgers, slip_plane, poisson);
}
return disp;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerDD, typename ContainerPoints>
std::enable_if_t<ContainerDD::Dim != 3 or ContainerPoints::Dim != 3>
ComputeBarnett::build(ContainerDD &, ContainerPoints &) {
LM_ASSERT(ContainerDD::Dim == 3, "barnett compute requires 3d.");
LM_ASSERT(ContainerPoints::Dim == 3, "barnett compute requires 3d.");
}
/* -------------------------------------------------------------------------- */
template <typename ContainerDD, typename ContainerPoints>
std::enable_if_t<ContainerDD::Dim == 3 and ContainerPoints::Dim == 3>
ComputeBarnett::build(ContainerDD &cont_dd, ContainerPoints &cont_points) {
CommGroup group_dd = cont_dd.getCommGroup();
CommGroup group_points = cont_points.getCommGroup();
if (group_dd.size() > 1 or group_points.size() > 1)
LM_FATAL("barnett compute is not parallel yet");
// this->copyContainerInfo(cont_points);
this->clear();
this->getOutputAsArray().resize(0, 3);
for (auto &&p : cont_points) {
Vector<3> disp = computeDisloDisplacements(cont_dd, p.position0(),
this->slip_plane, poisson);
this->getOutputAsArray().push_back(disp);
}
DUMP(this->getOutputAsArray().size(), DBG_DETAIL);
// this->name_computed.clear();
// this->name_computed.push_back("barnett_dispX");
// this->name_computed.push_back("barnett_dispY");
// this->name_computed.push_back("barnett_dispZ");
// if (!allPositions)
// allPositions = &cont.allGatherAllData();
// // Calculate the displacements.
// std::vector<Real> allDisplacements(allPositions->size());
// DomainInterface *domDD = DomainMultiScale::getManager().getObject(dom);
// if (DomainDDInterface *_ptr = dynamic_cast<DomainDDInterface *>(domDD))
// _ptr->computeDisplacements(*allPositions, allDisplacements);
// else
// LM_FATAL("Domain passed is not of DD type");
// // FilterID id = compute_list[0];
// // ComputeInterface * my_compute
// // =
// //
// dynamic_cast<ComputeInterface*>(FilterManager::getManager().getObject(id));
// // Distribute results such that each processor has full
// // displacements only for all of its original points.
// UInt localSize = cont.size();
// CommGroup group = this->getCommGroup();
// Communicator &comm = Communicator::getCommunicator();
// MPI_Comm mpiComm = comm.getMpiGroup(group);
// UInt np = comm.getNBprocsOnGroup(group);
// std::vector<Real> displs(np * localSize);
// std::vector<int> &scnts(cont.rcnts), sdispls(np);
// std::vector<int> rcnts(np, localSize), rdispls(np);
// for (unsigned i = 0; i < np; ++i) {
// sdispls[i] = (i == 0) ? 0 : scnts[i - 1];
// rdispls[i] = i * localSize;
// }
// MPI_Alltoallv(&allDisplacements[0], &scnts[0], &sdispls[0], MPI_DOUBLE,
// &displs[0], &rcnts[0], &rdispls[0], MPI_DOUBLE, mpiComm);
// // Output displacements corresponding to input points.
// this->empty();
// this->setDim(3);
// for (unsigned i = 0; i < localSize; ++i) {
// Real sum(0.0);
// for (unsigned p = 0; p < np; ++p)
// sum += displs[i + rdispls[p]];
// this->add(sum);
// }
}
/* -------------------------------------------------------------------------- */
/* LMDESC BARNETT
Calculate DD displacements at fieldpoints defined by INPUT.
*/
/* LMEXAMPLE
COMPUTE ddDisp BARNETT INPUT points DOMAIN dd
*/
inline void ComputeBarnett::declareParams() {
ComputeInterface::declareParams();
/* LMKEYWORD POISSON
Specify Poisson's ratio of the material
*/
this->parseKeyword("POISSON", this->poisson, 0.3);
/* LMKEYWORD SLIP_PLANE
Specify the slip plane to consider
*/
this->parseKeyword("SLIP_PLANE", this->slip_plane);
}
/* -------------------------------------------------------------------------- */
void ComputeBarnett::compute_make_call() {
DUMP(this->getID() << ": Compute", DBG_INFO);
this->build<dispatch>(this->getInput("segments"), this->getInput("points"));
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__

Event Timeline

c4science · Help