Homec4science
Diffusion LibMultiScale (features/cadd) Browse

libmultiscale/src/coupling/template_cadd.hh
d8e9d1cc530cfeatures/cadd

Last Change
Raw File

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_cadd_template.hh"
#include "container_mesh.hh"
#include "coupling_dd_md.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_parsable.hh"
#include "lm_parser.hh"
#include "ref_point_data.hh"
/* -------------------------------------------------------------------------- */
#include <string>
/* -------------------------------------------------------------------------- */
__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>
ContainerGenericDDMesh<3> sendDDSegments2MD(DomainA &domA);
template <typename Vec, typename RefAtom>
void functionifAniso(UInt num_dislo, Vec &pos0_pbc, Vec &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 Vec, typename RefAtom>
void functionElse(ContainerGenericDDMesh<3> &dd_points, Vec &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();
/**
Performs the following tasks on processors MD and DD
MD: 1. detect the MD dislocation segments (DXA or user-defined compute)
2. send to DD processor(s)
DD: 1. (at every step) fix the DD nodes located within MD domain
3. receive the detected dislocation segments from MD processors
4. move DD nodes to the closest detected segments
5. register the remaining detected nodes as new DD nodes
6. check topology - remesh
*/
//! Impose the position of detected segments to coupled DDD segments
template <typename DomainDD, typename DomainMD>
void applyDetection(DomainDD &domDD,
DomainMD &domA); // first coupling step
/**
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
*/
//! Impose the template onto the boundary atoms in the PAD
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(ContainerGenericDDMesh<3> &dd_points,
ContainerGenericDDMesh<3> &dd_points_global,
DomainA &domA);
bool coarse_dislocation_loop(ContainerGenericDDMesh<3> &dd_points,
ContainerGenericDDMesh<3> &dd_points0);
void find_indices_for_neighbors(ContainerGenericDDMesh<3> &dd_points);
void number_of_dislocation_loops(UInt num_dd_pts, UInt *num_loops,
ContainerGenericDDMesh<3> &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(ContainerGenericDDMesh<3> &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(ContainerGenericDDMesh<3> &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
ComputeCADDTemplate core_template;
// 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<ComputeInterface *> 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;
Quantity<Length> position_tolerance;
// 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__ */
c4science · Help