Page MenuHomec4science
Files F72840426

compute_cadd_template.cc
No OneTemporary
Actions

File Metadata

Created
Wed, Jul 17, 07:57

compute_cadd_template.cc
View Options

#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.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.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 [[gnu::unused]],
ElementType &segment
[[gnu::unused]]) {
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);
LM_TOIMPLEMENT;
// for (UInt i = 0; i < Dim; ++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 [[gnu::unused]],
UInt *indices_closest_angles [[gnu::unused]],
Real *closest_angles [[gnu::unused]]) {
// 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);
this->build<dispatch>(this->getInput("points"), this->getInput("segments"));
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__

Event Timeline

c4science · Help