Page MenuHomec4science
Files F90464621

template_cadd.hh
No OneTemporary
Actions

File Metadata

Created
Fri, Nov 1, 22:24

template_cadd.hh
View Options

/**
* @file template_cadd.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Wed Jan 15 17:00:43 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/>.
*
*/
#ifndef __LIBMULTISCALE_TEMPLATE_CADD_HH__
#define __LIBMULTISCALE_TEMPLATE_CADD_HH__
#include "bridging.hh"
#include "compute_compatibility.hh"
#include "coupling_dd_md.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_parsable.hh"
#include "lm_parser.hh"
#include <string>
#include "container_mesh.hh"
#include "ref_point_data.hh"
//#include "twod_template.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
struct JaehyunDDNeighbor {
Real idx; // 0
Vector<3> X; // 1,2,3
Vector<3> burg; // 4,5,6
Vector<3> normal; // 7,8,9
};
struct JaehyunDDPoint {
Real idx; // 0
Real num_neigh; // 1
Vector<3> X; // 2,3,4
Real idx0; // 5
std::vector<JaehyunDDNeighbor> neighs; // 6
Vector<3> tmp_X; // 27,28,29
Vector<3> tmp2_X; // 30,31,32
};
inline std::ostream &operator<<(std::ostream &os, JaehyunDDNeighbor &ref) {
os << "\t\tnbr_idx:" << ref.idx << std::endl;
os << "\t\tnbr_x:" << ref.X << std::endl;
os << "\t\tbx:" << ref.burg << std::endl;
os << "\t\tnx:" << ref.normal << std::endl;
return os;
}
inline std::ostream &operator<<(std::ostream &os, JaehyunDDPoint &ref) {
os << "JaehyunDD point:" << std::endl;
os << "\tidx:" << ref.idx << std::endl;
os << "\tnum_neigh:" << ref.num_neigh << std::endl;
os << "\tcrr_x:" << ref.X << std::endl;
os << "\tidx0:" << ref.idx0 << std::endl;
for (UInt i = 0; i < ref.num_neigh; ++i) {
os << ref.neighs[i];
}
return os;
}
/* -------------------------------------------------------------------------- */
class CADD_template : public CouplingDDMD {
public:
DECLARE_COUPLER(CADD_template, domDD, domMD);
template <typename DomainDD, typename DomainA>
void coupling(CouplingStage stage, DomainDD &domDD, DomainA &domA);
template <typename DomainDD, typename DomainA>
void init(DomainDD &domDD, DomainA &domA);
protected:
template <typename DomainA>
std::vector<JaehyunDDPoint> sendDDSegments2MD(DomainA &domA);
template <typename RefAtom>
void
functionifAniso(UInt num_dislo, Vector<3> &pos0_pbc, Vector<3> &pos0_pbc0,
std::vector<std::vector<Real>> &dislocations, RefAtom &at);
template <typename RefAtom>
void
functionifVolterra(UInt num_dislo, Vector<3> &pos0_pbc, Vector<3> &pos0_pbc0,
std::vector<std::vector<Real>> &dislocations, RefAtom &at);
template <typename RefAtom>
void functionElse(std::vector<JaehyunDDPoint> &dd_points, Vector<3> &pos0_pbc,
RefAtom &at);
inline void setGeometries(int pad_id, int atom_id, int geom_id,
int boundary_id);
inline virtual void clearAll();
inline Real compute_dd_theta(Real x3[], Real x4[]);
template <typename Vec>
inline bool intersectTwoD(Vec &p1, Vec &p2, Vec &p3, Vec &p4, Vector<2> &pts,
UInt slip_case);
template <typename Vec>
inline bool intersectThreeD(Vec &p1, Vec &p2, Vec &p3, Vec &p4, Vec &pts,
UInt slip_case);
inline bool intersectTwoD_infinite(Real p1[], Real p2[], Real p3[], Real p4[],
Real *ptsx, Real *ptsz);
template <typename Vec>
inline bool intersect(Vec &x1, Vec &x2, Vec &x3, Vec &x4, UInt slip_case);
inline Real ccw(Real Ax, Real Ay, Real Bx, Real By, Real Cx, Real Cy);
inline Real determinant(Real a, Real b, Real c, Real d, Real e, Real f,
Real g, Real h, Real i);
inline Real radian2degree(Real rad);
inline Real angle_on_plane(Real pos[], UInt axis);
// inline void RotationMatrixTemplateToRealLattices(Real theta, Real theta0,
// std::vector<Real>
// &resultR);
// inline void MatrixMultiplication(std::vector<Real> &mt, Real vec[],
// Real result[]);
// inline void CrossProduct(Real vec1[], Real vec2[], Real result[]);
inline Real distance_between_line_and_point(Real x0, Real z0, Real pts1[],
Real pts2[]);
inline void intersect_of_point_on_line_3d(Real p0[], Real p1[], Real p2[],
Real result[]);
template <typename Vec>
inline Real distance_between_line_and_point3D(Vec &pos0, Vec &pts1,
Vec &pts2);
template <typename Vec>
inline Real AreaTriangle(Vec &p1, Vec &p2, Vec &p3, UInt slip_case);
template <typename Vec>
inline bool is_it_closed(Vec cp, Vec &cp0, Real &ci, Real &ci0);
inline UInt next_index(UInt next_idx, UInt num_dd_pts0,
std::vector<UInt> indices);
template <typename Vec>
inline UInt find_one_direction(Vec &nb1, Vec &nb2, Real b1mag, Real b2mag,
UInt idx1, UInt idx2, UInt *next_i);
inline Real DistanceBetweenTwoNodes(Real x1, Real y1, Real z1, Real x2,
Real y2, Real z2);
void mergeDuplicateNodes();
//! main coupling methods
template <typename DomainDD, typename DomainMD>
void applyDetection(DomainDD &domDD,
DomainMD &domA); // first coupling step
template <typename DomainDD, typename DomainMD>
void applyTemplateToPad(DomainDD &domDD,
DomainMD &domA); // second coupling step
//! lower level auxiliary coupling methods
void detect_core();
// void detect_core_gui();
ContainerGenericDDMesh<3> detect_cores_by_compute();
template <typename RefAtom>
void applyTemplateToPoint(RefAtom &at, const Real pos_core[3]);
template <typename DomainA>
void close_dislocation_loop(std::vector<JaehyunDDPoint> &dd_points,
std::vector<JaehyunDDPoint> &dd_points_global,
DomainA &domA);
bool coarse_dislocation_loop(std::vector<JaehyunDDPoint> &dd_points,
std::vector<JaehyunDDPoint> &dd_points0);
void find_indices_for_neighbors(std::vector<JaehyunDDPoint> &dd_points);
void number_of_dislocation_loops(UInt num_dd_pts, UInt *num_loops,
std::vector<JaehyunDDPoint> &dd_points);
// void set_loop_number(Real node[], UInt num_dd_pts, Real **dd_points, UInt
// num_loops, UInt *i_index);
UInt startNewDDPoint(UInt num_dd_pts, std::vector<UInt> indices);
UInt findIndex_nbrnode(Real node[], UInt num_dd_pts, Real **dd_points);
// void position_of_atom(Real pos0[], int *position, Real *template_disp,
// Real *barnett_disp, UInt num_dd_pts, Real **dd_points, Real maxSeg,
// Real *f);
void center_of_circumcircle(Real pts1[], Real pts2[], Real pts3[], Real *cx,
Real *cy, Real *cz);
void coordinate_for_2d_template(Real distance, Real distance_from_dd_core[],
Real global_coord[], Real template_coords[]);
void compute_angle_ratio(Real angle, UInt *two_indices, Real *two_angles);
void computeNextPositionBasedOnCompute(std::vector<Real> computed_x,
std::vector<Real> computed_z, Real x,
Real z, Real mindist, Real maxdist,
Real *next_position);
void newNodesRegistration(std::vector<Real> NextDDXs,
std::vector<Real> NextDDYs,
std::vector<Real> NextDDZs);
template <typename RefDDNode>
void moveDDconstNode2MDcore(std::vector<Real> &computed_x,
std::vector<Real> &computed_y,
std::vector<Real> &computed_z, RefDDNode node);
template <typename RefDDNode, typename IteratorDDNodes>
void letAsFreeCombinedNode(RefDDNode node, IteratorDDNodes nodes);
void findPointInDDLine(std::vector<Real> ddx, std::vector<Real> ddy,
std::vector<Real> ddz, Real pts0[], Real pts1[],
Real pts2[], Real pts3[], std::vector<Real> &NextDDXs,
std::vector<Real> &NextDDYs,
std::vector<Real> &NextDDZs);
int evaluate_edge(Real x, Real y, Real z, Real pts1[], Real pts2[],
Real pts3[], Real pts4[], Real cent[]);
bool in_side_ellipse(Real pos[], Real rmax, Real rmin, Real *factor);
bool is_it_infinite_connection(UInt index0, UInt *index2);
template <typename RefAtom>
void applyDisplacement(RefAtom at, Vector<3> &total_disp);
template <typename DomainA> void SetPadAtoms(DomainA &domA);
void SetP0Atoms();
void barnett_disp_atom_all_segments(Real pos0[], Real barnett_disp[],
Real **dd_points, UInt num_dd_pts,
Real burgMag, Real pois,
UInt loop_number);
void compute_closure_point(Real **dd_points, UInt num_dd_pts,
UInt loop_number, Real A[], Real B[], Real b[],
Real C[]);
template <typename IteratorDDNodes>
void compute_pad_dd(UInt *Npad_dd, UInt *Npad_dd_column,
IteratorDDNodes nodes, Real md_xrange[], Real md_yrange[],
Real md_zrange[]);
void constructElems(std::vector<Real> nextDDXs, std::vector<Real> nextDDYs,
std::vector<Real> nextDDZs,
std::vector<RefGenericElem<3>> &els,
std::vector<Real> &normal_ptr,
std::vector<Real> &burgers_ptr);
template <typename Vec>
bool PointInSquare(Vec &pts1, Vec &pts2, Vec &pts3, Vec &pts4, Vec &pts,
UInt slip_case) {
LM_TOIMPLEMENT;
}
bool is_it_pad_atom(Real pos0[], Real pos[], Geometry *pad);
template <typename Vec>
bool ComputeTemplateDisplacements(Vec &template_disp, Vec &current_node,
Vec &next_node, Vec &burgers, Vec &normal,
Vec &pos0) {
LM_TOIMPLEMENT;
}
void ComputeAnalyticDisplacements(Vector<3> &barnett_disp, Vector<3> &B,
Vector<3> &A, Vector<3> &burg,
Vector<3> &pos0);
void ComputeBoundaryPointsInCoreRegion(Real pts1[], Real pts2[],
Real currnet_node[], Real center[]);
UInt FindNeighborIndex(UInt next_i, Real **dd_points, UInt num_dd_pts);
template <typename DomainA>
UInt CollectingLocalDDpoints(std::vector<JaehyunDDPoint> &dd_points0_local,
DomainA &domA);
Real *CollectingDDpoints();
template <typename RefAtom>
inline void shear_strain(RefAtom at, Vector<3> &disp);
template <typename RefAtom>
inline void tilt_strain(RefAtom at, Vector<3> &disp);
inline void
ComputeOffsetParallelCurve(std::vector<JaehyunDDPoint> &dd_points);
inline void computeOffset(Vector<3> &current_node, Real slope,
Vector<3> &pts1, Vector<3> &pts2, UInt slip_case);
inline void FindOffsetPoints(Real node_idx, Real *parallel_curve,
UInt num_dd_pts, Real *pts1, Real *pts2);
inline bool near_md_boundary(Real *A, Real *B, Real multiple);
inline void findCorrectOffsets(Real ang, Vector<3> &pts1, Vector<3> &pts2,
UInt slip_case);
inline void
CorrectIntersetingPoints(UInt num_pts,
Real **dd_points); //, Real *parallel_curve);
// inline void ComputeOffsetParallelCurve(UInt num_dd_pts, Real **dd_points,
// Real *parallel_curve);
// inline void computeOffset(Real *current_node, Real slope, Real *pts1, Real
// *pts2);
// inline void FindOffsetPoints(Real node_idx, Real *parallel_curve, UInt
// num_dd_pts, Real *normal);
inline void RelativePositionToPlane(Real *node, Real *normal, Real *pos,
Real *position);
// inline void findCorrectOffsets(Real ang, Real *pts1, Real *pts2);
// inline void CorrectIntersetingPoints(UInt num_pts, Real **dd_points, Real
// *parallel_curve);
inline void find_another_point_on_same_plane(Real *dd_node1, Real *burg,
Real *normal1, Real *normal2,
UInt num_dd_pts,
Real **dd_points, Real *dX1,
Real *dX2);
// inline Real BeDisplacement(Real b_edge, Real nu, Real X[]);
// inline Real NDisplacement(Real b_edge, Real nu, Real X[]);
// inline Real LDisplacement(Real b_screw, Real nu, Real X[]);
// inline void volterra_strain(Real pos[], Real disp[], Real pos0[]);
inline Real eval_q_coef(Real *X);
inline Real eval_t_coef(Real *X);
inline Real evalArctan(Real X, Real Y);
inline Real computeBeDisplacementAniso(Real b_edge, Real *X);
inline Real computeNDisplacementAniso(Real b_edge, Real *X);
inline Real computeLDisplacementAniso(Real b_screw, Real *X);
inline Real computeBeDisplacement(Real b_edge, Real *X);
inline Real computeNDisplacement(Real b_edge, Real *X);
inline Real computeLDisplacement(Real b_screw, Real *X);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
//! detected core positions
std::vector<std::vector<Real>> md_cores;
int num_cores;
int detect_freq;
Real cluster_radius;
//! The pad zone and DD nodes in atomistic domain
OutputContainer pad_atoms;
//! geometry ids for the pad zone and atomistic domain
LMID pad_geom;
// GeomID top_moving_geom;
// GeomID bot_moving_geom;
// GeomID slip_moving_geom;
// FilterID pad_filter;
int normal_quant[3];
Real center[3];
Real angle_between_b_and_x;
int close_direction[3];
int num_close_loop_in_infinity;
bool close_in_infinity;
bool multiple_loops;
bool current_position_pad;
bool coarsen_dd;
bool use_specific_template;
Real specific_template_angle;
Real min_area;
int detect_simple_direction[3];
bool detect_simple_bool;
UInt simple_direction;
LMID md_geom;
LMID dd_geom;
Vector<3> md_bbox_min, md_bbox_max;
std::shared_ptr<Cube> md_cube;
std::shared_ptr<Cube> dd_cube;
// Real top_moving_xrange[2], top_moving_yrange[2], top_moving_zrange[2];
// Real bot_moving_xrange[2], bot_moving_yrange[2], bot_moving_zrange[2];
// Real slip_moving_xrange[2],slip_moving_yrange[2],slip_moving_zrange[2];
Real md_box_lengths[3];
Real md_top_z, md_bot_z;
LMID detect_geom;
bool detect_dxa;
bool detect_partial;
bool dd_points_out;
bool dd_boundary;
UInt num_templates;
std::vector<std::string> templatefilename_list;
std::string directory_dd_points;
//! template related data
// std::vector<TwoDTemplate<RefAtom, IteratorAtomsSubset> *> core_templates;
// TwoDTemplate<RefAtom, IteratorAtomsSubset> core_template;
//! list of computes as added input
std::vector<LMID> compute_list;
// std::vector<FilterID> truedisp;
// std::vector<FilterID> barnett_compute;
//! computed values filtered out
std::vector<ComputeCompatibility *> out_computes;
Real epsilon[2];
Real core_rcut; // core template correction radius
Real rmin; // minor radius of the region where core template applied
Real factor;
Real mu0;
Real shear_stresses[4]; // # sigma_xy, sigma_yx, sigma_zy, sigma_yz
Real tilt_stresses[3]; // # xy, xz, yz
// UInt nb_slips;
// Real len_inter_slips;
// UInt node_index_list[8];
// UInt num_pairs_in_infinity;
LMID domdd;
// int volterra_replica;
// Real volterra_angle;
Real distance_ddline_from_pad;
Real pad_thickness;
// Real release_Xhi;
// Real release_Xlo;
// Real release_Zhi;
// Real release_Zlo;
Real pois0;
Real maxSeg0;
Real minSeg0;
Real burgMag0;
// Real maxDDY0;
// Real minDDY0;
bool aniso;
bool volterra;
Real elastic_coefs[13];
Real lambda;
Real phi;
Real c11bar;
};
/* -------------------------------------------------------------------------- */
template <>
inline UInt Parser::parse(dislo_orientation &orientation,
std::stringstream &line, UInt) {
int val;
UInt nb = parse(val, line);
orientation = static_cast<dislo_orientation>(val);
return nb;
}
/* -------------------------------------------------------------------------- */
inline void CADD_template::setGeometries(int pad_id, int atom_id, int geom_id,
[[gnu::unused]] int boundary_id) {
this->pad_geom = pad_id;
this->md_geom = atom_id;
this->detect_geom = geom_id;
}
/* -------------------------------------------------------------------------- */
inline void CADD_template::clearAll() {
// this->md_cores.clear();
// this->pad_atoms.empty();
// this->core_templates.clear();
this->out_computes.clear();
}
/* -------------------------------------------------------------------------- */
inline Real CADD_template::compute_dd_theta(Real x3[], Real x4[]) {
Real norm1, bx1, bz1, bx2, bz2, norm2, inner, angle_with_ONEZERO, angdist;
bx1 = 0.0;
bz1 = 0.0;
bx2 = 1.0;
bz2 = 0.0;
norm1 = 1.0;
norm2 = sqrt(pow((x4[0] - x3[0]), 2) + pow((x4[1] - x3[1]), 2));
inner = (bx2 - bx1) * (x4[0] - x3[0]) + (bz2 - bz1) * (x4[1] - x3[1]);
angle_with_ONEZERO = acos(inner / norm1 / norm2) * 180.0 / M_PI;
angdist = fabs(90.0 - abs(angle_with_ONEZERO));
return (180.0 - angdist);
}
/* -------------------------------------------------------------------------- */
template <typename Vec>
inline bool CADD_template::intersectTwoD(Vec &p1, Vec &p2, Vec &p3, Vec &p4,
Vector<2> &pts, UInt slip_case) {
if (intersect(p1, p2, p3, p4, slip_case)) {
UInt id0 = 0;
UInt id1 = 1;
if (slip_case == 0) {
id0 = 1;
id1 = 2;
} else if (slip_case == 1) {
id0 = 0;
id1 = 2;
} else if (slip_case == 2) {
id0 = 0;
id1 = 1;
}
Real det = (p1[id0] - p2[id0]) * (p3[id1] - p4[id1]) -
(p1[id1] - p2[id1]) * (p3[id0] - p4[id0]);
Real nominator1 =
(p1[id0] * p2[id1] - p1[id1] * p2[id0]) * (p3[id0] - p4[id0]) -
(p1[id0] - p2[id0]) * (p3[id0] * p4[id1] - p3[id1] * p4[id0]);
Real nominator2 =
(p1[id0] * p2[id1] - p1[id1] * p2[id0]) * (p3[id1] - p4[id1]) -
(p1[id1] - p2[id1]) * (p3[id0] * p4[id1] - p3[id1] * p4[id0]);
Real px = nominator1 / det;
Real py = nominator2 / det;
Real ratioid0;
if (fabs(p2[id0] - p1[id0]) < 1e-5)
ratioid0 = 0.0;
else
ratioid0 = roundf((px - p1[id0]) * 10.0) / 10.0 / (p2[id0] - p1[id0]);
Real ratioid1;
if (fabs(p2[id1] - p1[id1]) < 1e-5)
ratioid1 = 0.0;
else
ratioid1 = roundf((py - p1[id1]) * 10.0) / 10.0 / (p2[id1] - p1[id1]);
if ((ratioid0 <= 1.0) && (0.0 <= ratioid0) && (ratioid1 <= 1.0) &&
(0.0 <= ratioid1)) {
pts[0] = px;
pts[1] = py;
return true;
} else {
pts[0] = -1.0;
pts[1] = -1.0;
return false;
}
}
pts[0] = -1.0;
pts[1] = -1.0;
return false;
}
/* -------------------------------------------------------------------------- */
template <typename Vec>
inline bool CADD_template::intersectThreeD(Vec &p1, Vec &p2, Vec &p3, Vec &p4,
Vec &pts, UInt slip_case) {
// different slip plane
if ((fabs(p1[1] - p3[1]) > 1.0) || (fabs(p2[1] - p3[1]) > 1.0) ||
(fabs(p1[1] - p4[1]) > 1.0) || (fabs(p2[1] - p4[1]) > 1.0)) {
return false;
} else {
Vector<2> pts2D = {pts[0], pts[1]};
if (intersectTwoD(p1, p2, p3, p4, pts2D, slip_case)) {
pts[1] = p1[1];
return true;
} else {
return false;
}
}
return false;
}
/* -------------------------------------------------------------------------- */
inline bool CADD_template::intersectTwoD_infinite(Real p1[], Real p2[],
Real p3[], Real p4[],
Real *ptsx, Real *ptsz) {
Real x1 = p1[0];
Real y1 = p1[2];
Real x2 = p2[0];
Real y2 = p2[2];
Real x3 = p3[0];
Real y3 = p3[2];
Real x4 = p4[0];
Real y4 = p4[2];
Real denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
if (fabs(denom) <= 1e-3) {
*ptsx = -1.0;
*ptsz = -1.0;
return false;
} else {
*ptsx =
((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) /
denom;
*ptsz =
((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) /
denom;
return true;
}
}
/* -------------------------------------------------------------------------- */
template <typename Vec>
inline bool CADD_template::intersect(Vec &x1, Vec &x2, Vec &x3, Vec &x4,
UInt slip_case) {
if (slip_case == 0)
return ccw(x1[1], x1[2], x3[1], x3[2], x4[1], x4[2]) !=
ccw(x2[1], x2[2], x3[1], x3[2], x4[1], x4[2]) &&
ccw(x1[1], x1[2], x2[1], x2[2], x3[1], x3[2]) !=
ccw(x1[1], x1[2], x2[1], x2[2], x4[1], x4[2]);
else if (slip_case == 1)
return ccw(x1[0], x1[2], x3[0], x3[2], x4[0], x4[2]) !=
ccw(x2[0], x2[2], x3[0], x3[2], x4[0], x4[2]) &&
ccw(x1[0], x1[2], x2[0], x2[2], x3[0], x3[2]) !=
ccw(x1[0], x1[2], x2[0], x2[2], x4[0], x4[2]);
else if (slip_case == 2)
return ccw(x1[0], x1[1], x3[0], x3[1], x4[0], x4[1]) !=
ccw(x2[0], x2[1], x3[0], x3[1], x4[0], x4[1]) &&
ccw(x1[0], x1[1], x2[0], x2[1], x3[0], x3[1]) !=
ccw(x1[0], x1[1], x2[0], x2[1], x4[0], x4[1]);
return false;
}
/* -------------------------------------------------------------------------- */
inline Real CADD_template::ccw(Real Ax, Real Ay, Real Bx, Real By, Real Cx,
Real Cy) {
return (Cy - Ay) * (Bx - Ax) > (By - Ay) * (Cx - Ax);
}
/* -------------------------------------------------------------------------- */
inline Real CADD_template::determinant(Real a, Real b, Real c, Real d, Real e,
Real f, Real g, Real h, Real i) {
return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
}
/* -------------------------------------------------------------------------- */
inline Real CADD_template::radian2degree(Real rad) {
// switch radian to degree having a range (0,360)
Real ang = rad * 180.0 / M_PI;
if (ang <= 0.0)
return 360.0 + ang;
else
return ang;
}
/* -------------------------------------------------------------------------- */
inline Real CADD_template::angle_on_plane(Real pos[], UInt axis) {
if (axis == 0) {
return atan(pos[1] / pos[2]);
} else if (axis == 1) {
return atan(pos[2] / pos[0]);
} else if (axis == 2) {
return atan(pos[0] / pos[1]);
} else {
LM_FATAL("wrong axis");
}
}
// inline void CADD_template::RotationMatrixTemplateToRealLattices(Real theta,
// Real theta0,
// std::vector<Real>
// &resultR) {
// // |R[0] R[1] R[2]|
// // |R[3] R[4] R[5]|
// // |R[6] R[7] R[8]|
// resultR[0] = cos((theta - theta0) * M_PI / 180.0);
// resultR[2] = -sin((theta - theta0) * M_PI / 180.0);
// resultR[6] = sin((theta - theta0) * M_PI / 180.0);
// resultR[8] = cos((theta - theta0) * M_PI / 180.0);
// resultR[1] = this->normal_quant[0];
// resultR[4] = this->normal_quant[1];
// resultR[7] = this->normal_quant[2];
// resultR[3] = 0.0;
// resultR[5] = 0.0;
// }
// inline void CADD_template::MatrixMultiplication(std::vector<Real> &mt, Real
// vec[],
// Real result[]) {
// result[0] = mt[0] * vec[0] + mt[1] * vec[1] + mt[2] * vec[2];
// result[1] = mt[3] * vec[0] + mt[4] * vec[1] + mt[5] * vec[2];
// result[2] = mt[6] * vec[0] + mt[7] * vec[1] + mt[8] * vec[2];
// }
// inline void CADD_template::CrossProduct(Real vec1[], Real vec2[], Real
// result[]) {
// result[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
// result[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
// result[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
// }
/* -------------------------------------------------------------------------- */
inline Real CADD_template::distance_between_line_and_point(Real x0, Real z0,
Real pts1[],
Real pts2[]) {
Real x2 = pts2[0];
Real x1 = pts1[0];
Real z2 = pts2[2];
Real z1 = pts1[2];
Real dist = fabs((x1 - x2) * (z1 - z0) - (x1 - x0) * (z2 - z1)) /
sqrt((x2 - x1) * (x2 - x1) + (z2 - z1) * (z2 - z1));
return dist;
}
/* -------------------------------------------------------------------------- */
inline void CADD_template::intersect_of_point_on_line_3d(Real p0[], Real p1[],
Real p2[],
Real result[]) {
Real nom = (p1[0] - p0[0]) * (p2[0] - p1[0]) +
(p1[1] - p0[1]) * (p2[1] - p1[1]) +
(p1[2] - p0[2]) * (p2[2] - p1[2]);
Real denom = (p2[0] - p1[0]) * (p2[0] - p1[0]) +
(p2[1] - p1[1]) * (p2[1] - p1[1]) +
(p2[2] - p1[2]) * (p2[2] - p1[2]);
Real t = -1.0 * nom / (denom + 1e-8);
result[0] = p1[0] + (p2[0] - p1[0]) * t;
result[1] = p1[1] + (p2[1] - p1[1]) * t;
result[2] = p1[2] + (p2[2] - p1[2]) * t;
}
/* -------------------------------------------------------------------------- */
template <typename Vec>
inline Real CADD_template::distance_between_line_and_point3D(Vec &pos0,
Vec &pts1,
Vec &pts2) {
Real x0 = pos0[0];
Real y0 = pos0[1];
Real z0 = pos0[2];
Real u1 = x0 - pts1[0];
Real u2 = y0 - pts1[1];
Real u3 = z0 - pts1[2];
Real v1 = x0 - pts2[0];
Real v2 = y0 - pts2[1];
Real v3 = z0 - pts2[2];
Real cross[3] = {u2 * v3 - u3 * v2, u3 * v1 - u1 * v3, u1 * v2 - u2 * v1};
Real norm_cross =
sqrt(cross[0] * cross[0] + cross[1] * cross[1] + cross[2] * cross[2]);
Real r1 = pts2[0] - pts1[0];
Real r2 = pts2[1] - pts1[1];
Real r3 = pts2[2] - pts1[2];
Real denom = sqrt(r1 * r1 + r2 * r2 + r3 * r3);
Real dist = norm_cross / denom;
return dist;
}
/* -------------------------------------------------------------------------- */
template <typename Vec>
inline Real CADD_template::AreaTriangle(Vec &p1, Vec &p2, Vec &p3,
UInt slip_case) {
Real nom = 0.0;
if (slip_case == 0)
nom = p1[1] * (p2[2] - p3[2]) + p2[1] * (p3[2] - p1[2]) +
p3[1] * (p1[2] - p2[2]);
else if (slip_case == 1)
nom = p1[0] * (p2[2] - p3[2]) + p2[0] * (p3[2] - p1[2]) +
p3[0] * (p1[2] - p2[2]);
else if (slip_case == 2)
nom = p1[0] * (p2[1] - p3[1]) + p2[0] * (p3[1] - p1[1]) +
p3[0] * (p1[1] - p2[1]);
Real area = fabs(nom / 2.0) + 1e-8;
return area;
}
/* -------------------------------------------------------------------------- */
template <typename Vec>
inline bool CADD_template::is_it_closed(Vec cp, Vec &cp0, Real &ci, Real &ci0) {
Real dist = (cp - cp0).norm();
if ((dist < 1e-8) && (fabs(ci0 - ci) < 1e-5)) {
return true;
} else {
return false;
}
}
/* -------------------------------------------------------------------------- */
inline UInt CADD_template::next_index(UInt next_idx, UInt num_dd_pts0,
std::vector<UInt> indices) {
for (UInt i = 0; i < num_dd_pts0; ++i) {
if (!(std::find(indices.begin(), indices.end(), i) != indices.end())) {
return i;
}
}
return next_idx;
}
/* -------------------------------------------------------------------------- */
template <typename Vec>
inline UInt CADD_template::find_one_direction(
Vec &nb1 [[gnu::unused]], Vec &nb2 [[gnu::unused]], Real b1mag,
Real b2mag [[gnu::unused]], UInt idx1, UInt idx2, UInt *next_i) {
if (b1mag > 0.0) {
*next_i = idx1;
return 0;
} else {
*next_i = idx2;
return 1;
}
}
/* -------------------------------------------------------------------------- */
inline Real CADD_template::DistanceBetweenTwoNodes(Real x1, Real y1, Real z1,
Real x2, Real y2, Real z2) {
Real dist = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) +
(z1 - z2) * (z1 - z2));
return dist;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_TEMPLATE_CADD_HH__ */

Event Timeline

c4science · Help