Diffusion LibMultiScale (features/coupler_component) Browse

coupler_compon…

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__ */

libmultiscale/src/coupling/template_cadd.hhc7c5deee9c21features/coupler_compon…

template_cadd.hhView Options

libmultiscale/src/coupling/template_cadd.hh
c7c5deee9c21features/coupler_compon…

template_cadd.hh
View Options