template_cadd.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Jul 3, 01:38

template_cadd.cc
View Options

	/**
	* @file template_cadd.cc
	*
	* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
	* @author Till Junge <till.junge@epfl.ch>
	* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
	*
	* @date Thu Sep 18 14:51:07 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/>.
	*
	*/

	#define TIMER
	/* -------------------------------------------------------------------------- */
	#include "template_cadd.hh"
	#include "lib_bridging.hh"
	#include "lm_common.hh"
	#include <algorithm>
	#include <assert.h>
	#include <cmath>
	#include <fstream>
	#include <iomanip>
	#include <unistd.h>
	#include <utility>
	#include <vector>

	#include "mpi.h"
	#include <stdio.h>
	#include <stdlib.h>
	#include <string.h>

	#include "compute_extract.hh"
	#include "container_mesh.hh"
	#include "geometry_manager.hh"
	#include "reference_manager_interface.hh"

	#include "communicator.hh"
	#include "domain_dd_interface.hh"
	#include "ref_point_data.hh"

	/* -------------------------------------------------------------------------- */
	__BEGIN_LIBMULTISCALE__

	/* -------------------------------------------------------------------------- */

	CADD_template::CADD_template(const std::string &name)
	: LMObject(name), CouplingDDMD(name),
	pad_atoms(createOutput("pad-" + name)) {
	this->detect_freq = 50;
	this->num_close_loop_in_infinity = 0;
	this->close_in_infinity = false;
	this->multiple_loops = false;
	this->detect_dxa = false;
	this->detect_partial = false;
	this->detect_simple_bool = false;
	}

	/* -------------------------------------------------------------------------- */

	CADD_template::~CADD_template() { clearAll(); }

	/* -------------------------------------------------------------------------- */
	void CADD_template::mergeDuplicateNodes() {
	LM_TOIMPLEMENT;
	// for (auto &&nd : this->domA.getContainer().getContainerNodes()) {
	// LM_TOIMPLEMENT;
	// // UInt nbrs = nd.getNbrOfNeighs();
	// UInt nbrs = std::numeric_limits<UInt>::max();
	// UInt idx = nd.getIndex();

	// if (this->md_cube->contains(nd.position()) && nbrs == 2) {
	// UInt idx2 = nd.getIndex();
	// for (auto &&nd2 [[gnu::unused]] :
	// this->domA.getContainer().getContainerNodes()) {
	// LM_TOIMPLEMENT;
	// // auto dX = nd.position() - nd2.position();
	// // Real dist = dX.norm();
	// auto dist = std::numeric_limits<Real>::max();
	// if ((idx != idx2) && (dist < 1e-2)) {
	// this->domA.freenodearm(idx);
	// this->domA.freenode(idx);

	// this->domA.freenodearm(idx2);
	// this->domA.freenode(idx2);
	// }
	// }
	// }
	// }
	}

	/* -------------------------------------------------------------------------- */
	template <typename DomainDD, typename DomainMD>
	void CADD_template::init(DomainDD &domDD, DomainMD &domA) {

	// fabricating the geometries
	auto md_geom = GeometryManager::getManager().getGeometry(this->md_geom);
	this->md_cube = std::dynamic_pointer_cast<Cube>(md_geom);

	auto dd_geom = GeometryManager::getManager().getGeometry(this->dd_geom);
	this->dd_cube = std::dynamic_pointer_cast<Cube>(dd_geom);

	Vector<3u> xmin(this->md_cube->getXmin());
	Vector<3u> xmax(this->md_cube->getXmax());

	for (UInt i : {0, 1, 2}) {
	xmin[i] -= this->pad_thickness + this->distance_ddline_from_pad;
	xmax[i] += this->pad_thickness + this->distance_ddline_from_pad;
	}

	this->md_cube->setXmin(xmin);
	this->md_cube->setXmax(xmax);

	// constructing some flags
	this->detect_simple_bool = false;
	simple_direction = 3;

	for (auto &&i : {0, 1, 2}) {

	if (this->close_direction[i] != 0)
	this->close_in_infinity = true;

	if (this->detect_simple_direction[i] != 0) {
	this->detect_simple_bool = true;
	simple_direction = i;
	}
	}

	if (this->aniso) {
	Real c11 = this->elastic_coefs[0];
	Real c12 = this->elastic_coefs[1];
	Real c22 = this->elastic_coefs[4];
	Real c66 = this->elastic_coefs[12];
	Real lambda2 = sqrt(c11 / c22);

	this->c11bar = sqrt(c11 * c12);
	this->lambda = sqrt(lambda2);
	this->phi = 0.5 * acos((c12 * c12 + 2.0 * c12 * c66 - c11bar * c11bar) /
	(2.0 * c11bar * c66));
	if (this->volterra)
	LM_FATAL("cannot be used at the same time with volterra and aniso");
	}

	// sending infos on parameters to paradis
	Real maxSeg, burgMag;
	if (this->is_in_dd) {

	auto domDD_xmin = domDD.getDomainBoundingBox().getXmin();
	auto domDD_xmax = domDD.getDomainBoundingBox().getXmax();
	if ((this->dd_cube->getXmin() - domDD_xmin).norm() > 1e-7 \|\|
	(this->dd_cube->getXmax() - domDD_xmax).norm() > 1e-7)
	DUMP("[WARNING] Check DD box lengths defined Libmultiscale configuration "
	"and ParaDis control files",
	DBG_MESSAGE);

	LM_TOIMPLEMENT;
	// this->burgMag0 = domDD.getBurgersMagnitude() * 1e+10; // in angstrom unit

	int my_rank_in_DDprocs = this->comm_group_dd.getMyRank();

	if (my_rank_in_DDprocs == 0) {
	UInt mdNBprocs = this->comm_group_atomic.size();
	for (UInt i = 0; i < mdNBprocs; ++i) {
	this->comm_group_atomic.send(&maxSeg, 1, i, "maximum length of dd");
	this->comm_group_atomic.send(&burgMag, 1, i,
	"burgers vector magnitude");
	}
	}
	}

	// receive information on parameters to paradis
	if (this->is_in_atomic) {
	this->comm_group_dd.receive(&maxSeg, 1, 0, "maximum length of dd segments");
	this->comm_group_dd.receive(&burgMag, 1, 0, "burgers vector magnitude");
	// this->maxSeg0 = maxSeg;
	this->burgMag0 = burgMag;
	}
	}

	/* -------------------------------------------------------------------------- */
	template <typename DomainDD, typename DomainMD>
	void CADD_template::coupling(CouplingStage stage, DomainDD &domDD,
	DomainMD &domA) {

	//! Main coupling calls
	if (stage == COUPLING_STEP3) {
	STARTTIMER("applyDetection");
	this->applyDetection(domDD, domA);
	STOPTIMER("applyDetection");

	STARTTIMER("applyTemplateToPad");
	this->applyTemplateToPad(domDD, domA);
	STOPTIMER("applyTemplateToPad");
	return;
	}
	}

	/* ------------------------------------------------------------------------ */

	template <typename DomainDD, typename DomainMD>
	void CADD_template::applyDetection(DomainDD &domDD, DomainMD &domA) {

	if (current_step % this->detect_freq != 0)
	return;

	auto &all_group = Communicator::getCommunicator().getGroup("all");

	if (this->is_in_atomic) {

	STARTTIMER("detect_cores_by_compute");
	auto local_dislocations = this->detect_cores_by_compute();
	STOPTIMER("detect_cores_by_compute");

	ContainerGenericDDMesh<3> dislocations("detected_dislocations");
	comm_group_A.gather(local_dislocations, dislocations, 0,
	"Gather the detected dislocations");

	if (this->comm_group_atomic.getMyRank() == 0) {
	auto root_proc_DD = this->comm_group_dd[0].getAbsoluteRank();
	all_group.send(dislocations, root_proc_DD,
	"Send the detected segments to root of group_DD");
	}
	}

	if (this->is_in_dd) {

	STARTTIMER("receiveDetectedSegments");
	auto root_proc_MD = this->comm_group_atomic[0].getAbsoluteRank();
	ContainerGenericDDMesh<3> dislocations("detected_dislocations");

	all_group.receive(dislocations, root_proc_MD,
	"Receive the detected segments from root of group_MD");

	STOPTIMER("receiveDetectedSegments");

	STARTTIMER("moveSlavedDDNodes");
	/* loop over the dislocation nodes and treat the pinned ones
	in the coupling region, as they have to be moved to match
	the detected ones
	*/
	for (auto nd : domDD.getContainer().getContainerNodes()) {

	if (nd.getConstraint() != DomainParaDiS::_fixed)
	continue;
	auto &&dd_node = nd.position();
	if (not this->md_cube->contains(nd.position()))
	continue;

	Vector<3> new_position = Vector<3>::Zero();
	UInt counter = 0.0;
	for (auto &&detected_node : dislocations.getContainerNodes()) {
	auto &&detected_pos = detected_node.position();
	Real dist = (detected_pos - nd.position).norm();
	if (dist < this->position_tolerance) {
	new_position += detected_pos;
	counter += 1.0;
	}
	if (counter == 0.0)
	LM_FATAL(
	"cannot find the a detected node close by: " << nd.position());

	new_position /= counter;
	nd.position() = new_position;
	}
	nd.getContstraint() = DomainParaDiS::_free;
	}
	}
	STOPTIMER("moveSlavedDDNodes");

	STARTTIMER("DDremesh_in_CADD3d");
	domDD.remesh();
	STOPTIMER("DDremesh_in_CADD3d");

	// release or constrain again after coupling step
	for (auto &&nd : domA.getContainer().getContainerNodes()) {

	if (!detect_simple_bool) {
	if (this->md_cube->contains(nd.position())) {
	nd.constrain();
	} else {
	LM_TOIMPLEMENT;
	// UInt nbrs = nd.getNbrOfNeighs();
	// if (nbrs == 2)
	// nd.unconstrain();
	}
	}
	if (domA.is_it_bound()) {
	if (this->dd_cube->contains(nd.position())) {
	nd.constrain();
	}
	}
	}
	}

	/* -------------------------------------------------------------------------- */
	ContainerGenericDDMesh<3> CADD_template::detect_cores_by_compute() {

	LMID id = compute_list[0];
	auto &&my_compute = FilterManager::getManager().getObject(id);
	my_compute.compute();
	ContainerArray<Real> dislo_pos;

	Vector<3u> xmin(this->md_cube->getXmin());
	Vector<3u> xmax(this->md_cube->getXmax());
	Real pt = this->distance_ddline_from_pad;

	xmin -= pt;
	xmax += pt;

	Cube dd_cube_extended(3);
	dd_cube_extended.setXmin(xmin);
	dd_cube_extended.setXmax(xmax);

	auto segments = my_compute.evalOutput<ContainerGenericDDMesh<3>>();
	if (this->detect_simple_bool)
	return segments;

	FilterGeometry geom_filter("deteced_and_filtered");
	geom_filter.setParam("GEOMETRY", dd_cube_extended);
	geom_filter.compute(segments);
	return geom_filter.evalOutput<ContainerGenericDDMesh<3>>();
	}

	/* --------------------------------------------------------------------------

	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

	-------------------------------------------------------------------------- */
	template <typename DomainDD, typename DomainA>
	void CADD_template::applyTemplateToPad(DomainDD &domDD, DomainA &domA) {

	// set to zero force and velocity of the pad
	auto pad = GeometryManager::getManager().getGeometry(pad_geom);
	for (auto &&at : domDD.getContainer()) {
	if (pad->contains(at.position())) {
	at.velocity() = Vector<3>::Zero();
	at.force() = Vector<3>::Zero();
	}
	}

	// if not the right moment, do nothing
	if (current_step % detect_freq != 0)
	return;

	// retreive the dd_points (communication)
	std::vector<JaehyunDDPoint> dd_points = this->sendDDSegments2MD(domA);

	if (this->is_in_atomic) {

	STARTTIMER("applyDisplacement");

	UInt num_dislo = 0;
	std::vector<std::vector<Real>> dislocations;
	// Here, we analyze dd_points data to find out
	// where and how the cores exist

	if ((this->aniso) \|\| (this->volterra)) {

	for (auto &dd_pt : dd_points) {
	for (UInt j = 0; j < 2; ++j) {
	auto neigh = dd_pt.neighs[j];
	auto burg = neigh.burg * this->burgMag0;
	Real bmag =
	burg[0] * cos(this->angle_between_b_and_x / 180.0 * M_PI) +
	burg[2] * sin(this->angle_between_b_and_x / 180.0 * M_PI);
	bmag = roundf(bmag * 10.0) / 10.0;
	if (bmag / fabs(bmag) > 0) {
	Real idx = dd_pt.idx;
	auto &ddX = dd_pt.X;
	Real nbi = neigh.idx;
	if ((idx >= 0.0) && (nbi < 0.0)) {
	num_dislo += 1;
	std::vector<Real> coord_vector;
	coord_vector.push_back(ddX[0]);
	coord_vector.push_back(ddX[1]);
	coord_vector.push_back(ddX[2]);
	dislocations.push_back(coord_vector);
	}
	}
	}
	}
	DUMP(num_dislo << " straight dislocations will be imposed", DBG_MESSAGE);
	for (UInt i = 0; i < num_dislo; ++i) {
	DUMP("(" << dislocations[i][0] << ", " << dislocations[i][1] << ", "
	<< dislocations[i][2] << ")",
	DBG_MESSAGE);
	}
	}

	auto md_ranges = this->md_cube->getXmax() - this->md_cube->getXmin();
	Real mdxlen = md_ranges[0];
	Real mdzlen = md_ranges[2];
	UInt count_true_disp = 0;

	for (auto &&at : this->pad_atoms.get<DomainA::ContainerSubset>()) {
	auto &&disp = at.displacement();
	auto &&pos0 = at.position0();
	auto &&pos0_pbc = at.position0();
	auto &&pos0_pbc0 = at.position0();

	if (fabs(disp[2]) > mdzlen / 2.) {
	Real signZ = disp[2] / fabs(disp[2]);
	pos0_pbc[2] = pos0[2] + mdzlen * signZ;
	count_true_disp += 1;
	}

	if (fabs(disp[0]) > mdxlen / 2.) {
	Real signX = disp[0] / fabs(disp[0]);
	pos0_pbc[0] = pos0[0] + mdxlen * signX;
	}

	for (UInt i = 0; i < 3; ++i) {
	pos0_pbc0[i] = pos0_pbc[i];
	}

	if (this->aniso)
	functionifAniso(num_dislo, pos0_pbc, pos0_pbc0, dislocations, at);

	else if (this->volterra)
	functionifVolterra(num_dislo, pos0_pbc, pos0_pbc0, dislocations, at);

	else
	functionElse(dd_points, pos0_pbc, at);
	}

	STOPTIMER("applyDisplacement");
	// this->domC.updatePAD();
	} // if(is_in_atomic ...)
	}

	/* --------------------------------------------------------------------------
	*/
	void CADD_template::ComputeAnalyticDisplacements(
	Vector<3> &analytic_disp, Vector<3> &B [[gnu::unused]],
	Vector<3> &A [[gnu::unused]], Vector<3> &burg [[gnu::unused]],
	Vector<3> &pos0
	[[gnu::unused]]) { // Compute closure point for given AB segment
	Vector<3> C [[gnu::unused]] = Vector<3>::Zero();
	LM_TOIMPLEMENT;
	// this->domA.computeClosurePoint(B, A, burg, C);

	// Compute point P1, projected point of C on line consisted to A and Ab
	Vector<3> P1 [[gnu::unused]] = Vector<3>::Zero();
	LM_TOIMPLEMENT;
	// this->domA.computeProjectedPoint(B, C, burg, P1);

	// Compute point P2, projected point of C on line consisted to B and Bb
	Vector<3> P2 [[gnu::unused]] = Vector<3>::Zero();
	LM_TOIMPLEMENT;
	// this->domA.computeProjectedPoint(A, C, burg, P2);

	Vector<3> dispPointU1 = Vector<3>::Zero();
	LM_TOIMPLEMENT;
	// this->domA.computeBarnett(P2, B, A, burg, dispPointU1, pos0,
	// this->pois0);

	analytic_disp += dispPointU1;

	Vector<3> dispPointU2 = Vector<3>::Zero();
	LM_TOIMPLEMENT;
	// this->domA.computeBarnett(P1, B, P2, burg, dispPointU2, pos0,
	// this->pois0);
	analytic_disp += dispPointU2;
	}

	/* --------------------------------------------------------------------------
	*/
	void CADD_template::ComputeOffsetParallelCurve(
	std::vector<JaehyunDDPoint> &dd_points) {

	for (auto dd_pt : dd_points) {
	std::vector<Real> angles;
	Real count = 2.0;
	auto node_X = dd_pt.X;

	UInt slip_case = 3;
	Real avgangle = 0.;

	// computes some angle made by the arm in the slip plane
	for (auto &&neigh : dd_pt.neighs) {
	auto &next_X = neigh.X;
	// auto next_normal = neigh.normal;
	Real angle_between_neighbor = 0.0;
	Vector<3> segment = next_X - node_X;

	// should construct a rotation which puts
	// the normal as the direction Z
	LM_TOIMPLEMENT;
	angle_between_neighbor = atan2(segment[0], segment[1]);
	angles.push_back(angle_between_neighbor);
	avgangle += angle_between_neighbor;
	}
	avgangle /= dd_pt.neighs.size();

	Vector<3> dX1 = {0.0, 0.0, 0.0};
	Vector<3> dX2 = {0.0, 0.0, 0.0};
	computeOffset(node_X, avgangle, dX1, dX2, slip_case);

	if (count != 0.0) {
	if (isnan(dX1[0]) \|\| isnan(dX1[1]) \|\| isnan(dX1[2]) \|\| isnan(dX2[0]) \|\|
	isnan(dX2[1]) \|\| isnan(dX2[2]))
	LM_FATAL("wrong dX1 and dX2");
	}

	angles[0] = radian2degree(angles[0]);
	angles[1] = radian2degree(angles[1]);
	avgangle = (angles[0] + angles[1]) / count;
	findCorrectOffsets(avgangle, dX1, dX2, slip_case);

	dd_pt.tmp_X = node_X + dX1;
	dd_pt.tmp2_X = node_X + dX2;
	}
	return;
	}

	/* --------------------------------------------------------------------------
	*/
	void CADD_template::computeOffset(Vector<3> &current_node [[gnu::unused]],
	Real ang, Vector<3> &dX1, Vector<3> &dX2,
	UInt slip_case) {
	Real rmax = this->core_rcut;
	Real slope = tan(ang);

	Real val1 = sqrt(rmax * rmax / (1.0 + slope * slope));
	Real val2 = sqrt(rmax * rmax - rmax * rmax / (1.0 + slope * slope));
	Real val3 = -sqrt(rmax * rmax / (1.0 + slope * slope));
	Real val4 = -sqrt(rmax * rmax - rmax * rmax / (1.0 + slope * slope));

	if (slip_case == 0) {
	dX1[0] = 0.0;
	dX1[1] = val1;
	dX1[2] = val2;
	dX2[0] = 0.0;
	dX2[1] = val3;
	dX2[2] = val4;
	} else if (slip_case == 1) {
	dX1[0] = val1;
	dX1[1] = 0.0;
	dX1[2] = val2;
	dX2[0] = val3;
	dX2[1] = 0.0;
	dX2[2] = val4;
	} else if (slip_case == 2) {
	dX1[0] = val1;
	dX1[1] = val2;
	dX1[2] = 0.0;
	dX2[0] = val3;
	dX2[1] = val4;
	dX2[2] = 0.0;
	} else
	LM_FATAL("cannot define angle with neighbor");
	return;
	}

	/* --------------------------------------------------------------------------
	*/

	void CADD_template::findCorrectOffsets(Real ang, Vector<3> &dX1, Vector<3> &dX2,
	UInt slip_case) {
	Real dev = 1e+5;
	Real saved_i = 1.0;
	Real saved_k = 1.0;

	for (int i = -1; i < 3;) {
	for (int k = -1; k < 3;) {
	Real tmp = -1.0;
	if (slip_case == 0)
	tmp = fabs(radian2degree(atan2(dX2[1] * (Real)i, dX2[2] * (Real)k)) -
	ang);
	else if (slip_case == 1)
	tmp = fabs(radian2degree(atan2(dX2[2] * (Real)i, dX2[0] * (Real)k)) -
	ang);
	else if (slip_case == 2)
	tmp = fabs(radian2degree(atan2(dX2[0] * (Real)i, dX2[1] * (Real)k)) -
	ang);
	else
	LM_FATAL("unknown slip_case");

	if (tmp < dev) {
	dev = tmp;
	saved_i = (Real)i;
	saved_k = (Real)k;
	}
	k += 2;
	}
	i += 2;
	}

	if (slip_case == 0) {
	dX1[0] = dX1[0];
	dX1[1] = dX1[1] * saved_i;
	dX1[2] = dX1[2] * saved_k;

	dX2[0] = dX2[0];
	dX2[1] = dX2[1] * saved_i;
	dX2[2] = dX2[2] * saved_k;
	} else if (slip_case == 1) {
	dX1[0] = dX1[0] * saved_k;
	dX1[1] = dX1[1];
	dX1[2] = dX1[2] * saved_i;

	dX2[0] = dX2[0] * saved_k;
	dX2[1] = dX2[1];
	dX2[2] = dX2[2] * saved_i;
	} else if (slip_case == 2) {
	dX1[0] = dX1[0] * saved_i;
	dX1[1] = dX1[1] * saved_k;
	dX1[2] = dX1[2];

	dX2[0] = dX2[0] * saved_i;
	dX2[1] = dX2[1] * saved_k;
	dX2[2] = dX2[2];
	}

	return;
	}

	/* --------------------------------------------------------------------------
	*/

	bool CADD_template::in_side_ellipse(Real pos[], Real rmax, Real rmin,
	Real *factor) {
	Real eq =
	((pos[0] / rmax) * (pos[0] / rmax) + (pos[1] / rmin) * (pos[1] / rmin));
	Real val = 1.0 - eq;
	if (val >= 0.0) {
	*factor = sqrt(eq);
	return true;
	} else {
	*factor = sqrt(eq);
	return false;
	}
	}

	/* --------------------------------------------------------------------------
	*/
	bool CADD_template::is_it_infinite_connection(UInt index0 [[gnu::unused]],
	UInt *index0_infinite
	[[gnu::unused]]) {

	LM_TOIMPLEMENT;
	// int int_index0_infinite;
	// bool found = this->domA.getPairedIndex(index0, &int_index0_infinite);
	// if (found) {
	// *index0_infinite = (UInt)int_index0_infinite;
	// return true;
	// } else
	// return false;
	}

	/* --------------------------------------------------------------------------
	*/

	template <typename RefAtom>
	void CADD_template::shear_strain(RefAtom at, Vector<3> &disp) {
	Real eps_yx = this->shear_stresses[1] / this->mu0;
	Real eps_yz = this->shear_stresses[3] / this->mu0;

	Real eps_xy = this->shear_stresses[0] / this->mu0;
	Real eps_zy = this->shear_stresses[2] / this->mu0;

	disp[0] = eps_yx * at.position0()[1];
	disp[1] = eps_xy * at.position0()[0] + eps_zy * at.position0()[2];
	disp[2] = eps_yz * at.position0()[1];
	}

	/* --------------------------------------------------------------------------
	*/
	template <typename RefAtom>
	void CADD_template::tilt_strain(RefAtom at, Vector<3> &disp) {
	// Real LX = this->md_box_lengths[0];
	Real LY = this->md_box_lengths[1];
	Real LZ = this->md_box_lengths[2];

	// Real LY = fabs(this->md_yrange[1]) + fabs(this->md_yrange[0]);
	// Real LZ = fabs(this->md_zrange[1]) + fabs(this->md_zrange[0]);

	Real eps_xy = this->tilt_stresses[0] / this->mu0;
	Real eps_xz = this->tilt_stresses[1] / this->mu0;
	Real eps_yz = this->tilt_stresses[2] / this->mu0;

	disp[0] = -1.0 * eps_xy * (at.position0()[1] + LY / 2.0) -
	1.0 * eps_xz * (at.position0()[2] + LZ / 2.0);
	// disp[1] = -1.0eps_yz(at.position0(2) + LZ/2.0); // full dislocation
	disp[1] = -1.0 * eps_yz * (at.position0()[2] + LZ / 2.0); // half dislocation
	disp[2] = 0.0;
	}

	/* --------------------------------------------------------------------------
	*/

	template <typename DomainA>
	void CADD_template::close_dislocation_loop(
	std::vector<JaehyunDDPoint> &dd_points,
	std::vector<JaehyunDDPoint> &dd_points_global, DomainA &domA) {

	bool added_point_in_infinity(false);
	Real maxdist = domA.getMaxSeg();

	UInt count = 0;
	for (auto &dd_pt_global : dd_points_global) {

	Real idx_real = dd_pt_global.idx;
	Real num_neigh = dd_pt_global.num_neigh;

	auto crr_X = dd_pt_global.X;
	auto &dd_pt = dd_points[count];

	dd_pt.idx = idx_real;
	dd_pt.num_neigh = num_neigh;
	dd_pt.X = crr_X;

	for (UInt n = 0; n < num_neigh; ++n) {
	auto neigh = dd_pt_global.neighs[n];
	Real nbr_idx = neigh.idx;
	auto &nbr_X = neigh.X;

	auto burg = neigh.burg;
	auto normal = neigh.normal;

	auto diff = nbr_X - crr_X;
	Real dist = diff.norm();

	UInt nbr_idx0_infinite = 0;
	Real nbr_idx_infinite_real = 0.0;

	bool found(false);
	UInt idx0 = (UInt)dd_pt_global.idx0;

	if ((this->close_in_infinity) &&
	(is_it_infinite_connection(idx0, &nbr_idx0_infinite)) &&
	(num_neigh == 1)) {

	// connectivity over an infinity boundary for single neighbors nodes:
	// frank-read sources

	Vector<3> nbr_X_infinite;

	for (auto &dd_pt_global2 : dd_points_global) {
	Real tmp = dd_pt_global2.idx0;
	UInt nbnbr = dd_pt_global2.num_neigh;

	if (((tmp - nbr_idx0_infinite) < 1e-2) && (nbnbr == 1)) {
	nbr_idx_infinite_real = dd_pt_global2.idx;
	nbr_X_infinite = dd_pt_global2.X;
	found = true;
	}
	}

	if (!found)
	LM_FATAL("not found a node infinitely connected for idx0 " << idx0);

	DUMP("infinity points are considered between "
	<< idx0 << " (" << crr_X << ") and " << nbr_idx0_infinite
	<< " (" << nbr_X_infinite << ") in ("
	<< this->close_direction[0] << ", " << this->close_direction[1]
	<< ", " << this->close_direction[2] << ") direction",
	DBG_MESSAGE);

	dd_pt.num_neigh = 2;
	dd_pt.neighs[0].idx = nbr_idx;
	dd_pt.neighs[0].X = nbr_X;
	dd_pt.neighs[0].burg = burg;
	dd_pt.neighs[0].normal = normal;

	dd_pt.neighs[1].idx = -1.0 * idx_real - 10000.0;
	for (UInt i : {0, 1, 2}) {
	dd_pt.neighs[1].X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
	}
	dd_pt.neighs[1].burg = burg * -1.;
	dd_pt.neighs[1].normal = normal;

	auto &dd_pt_plus1 = *((&dd_pt) + 1);

	dd_pt_plus1.idx = -1.0 * idx_real - 10000.0;
	dd_pt_plus1.num_neigh = 2;

	for (UInt i : {0, 1, 2}) {
	dd_pt_plus1.X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
	}
	dd_pt_plus1.neighs[0].idx = idx_real;
	dd_pt_plus1.neighs[0].X = crr_X;
	dd_pt_plus1.neighs[0].burg = burg;
	dd_pt_plus1.neighs[0].normal = normal;
	dd_pt_plus1.neighs[1].idx = -1.0 * nbr_idx_infinite_real - 10000.0;
	for (UInt i : {0, 1, 2}) {
	dd_pt_plus1.neighs[1].X[i] =
	nbr_X_infinite[i] + close_direction[i] * 1E+3;
	}
	dd_pt_plus1.neighs[1].burg = burg * -1.;
	dd_pt_plus1.neighs[1].normal = normal;
	added_point_in_infinity = true;
	} else if ((close_in_infinity) && (num_neigh == 2) &&
	(dist > 5.0 * maxdist)) {
	// connectivity over an infinity boundary for Real neighbors nodes:
	// straight dislocation
	dd_pt.neighs[n].idx = -1.0 * dd_pt.idx - 10000.0;
	for (UInt i : {0, 1, 2}) {
	dd_pt.neighs[n].X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
	}
	dd_pt.neighs[n].burg = burg;
	dd_pt.neighs[n].normal = normal;

	auto &dd_pt_plus1 = *((&dd_pt) + 1);

	dd_pt_plus1.idx = -1.0 * dd_pt.idx - 10000.0;
	dd_pt_plus1.num_neigh = 2;
	for (auto i : {0, 1, 2}) {
	dd_pt_plus1.X[i] = crr_X[i] + this->close_direction[i] * 1E+3;
	}

	dd_pt_plus1.idx0 = dd_pt.idx;
	dd_pt_plus1.neighs[0].X = crr_X;
	dd_pt_plus1.neighs[0].burg = burg * -1.;
	dd_pt_plus1.neighs[0].normal = normal;

	dd_pt_plus1.neighs[1].idx = -1.0 * nbr_idx - 10000.0;
	for (auto i : {0, 1, 2}) {
	dd_pt_plus1.neighs[1].X[i] =
	nbr_X[i] + this->close_direction[i] * 1E+3;
	}
	dd_pt_plus1.neighs[1].burg = burg;
	dd_pt_plus1.neighs[1].normal = normal;

	DUMP("Implicit infinity connection between "
	<< idx0 << " (" << dd_pt.X << ") and " << dd_pt_plus1.X << " ("
	<< nbr_X << ").",
	DBG_MESSAGE);

	added_point_in_infinity = true;
	} else if (num_neigh == 1) {
	// && (this->detect_partial)) {
	// straight connection between two single neighbor node: partial
	// detection
	dd_pt.neighs[0].idx = nbr_idx;
	dd_pt.neighs[0].X = nbr_X;
	dd_pt.neighs[0].burg = burg;
	dd_pt.neighs[0].normal = normal;

	Real temp_val = 1e+4;
	// IteratorDDNodes nodes_loop3 =
	// this->domA.getContainerNodes().getIterator();
	// Real tmpx,tmpy,tmpz;

	for (auto &dd_pt_global2 : dd_points_global) {
	UInt num_neigh0 = dd_pt_global2.num_neigh;
	auto diff = dd_pt_global2.X - crr_X;
	Real dist0 = diff.norm();

	if ( // has to be terminated in MD domain
	(num_neigh0 == 1) &&
	// avoid self selection
	(0.1 < dist0) &&
	// find the node minimum distance
	(dist0 < temp_val) &&
	// avoid chosing the node in different slip plane
	(fabs(dd_pt_global2.X[1] - crr_X[1]) < 1.0)) {

	temp_val = dist;
	dd_pt.neighs[0].idx = dd_pt_global2.idx;
	dd_pt.neighs[0].X = dd_pt_global2.X;

	DUMP("Straight connection between the two fixed DD nodes "
	<< "(" << dd_pt.X << ") and "
	<< "(" << dd_pt_global2.X << ").",
	DBG_MESSAGE);
	}
	}

	dd_pt.num_neigh = 2.0;
	dd_pt.neighs[1].burg = -1.0 * dd_pt.neighs[0].burg;
	dd_pt.neighs[1].normal = dd_pt.neighs[0].normal;

	} else {
	// normal connection
	dd_pt.neighs[n].idx = nbr_idx;
	dd_pt.neighs[n].X = nbr_X;
	dd_pt.neighs[n].burg = burg;
	dd_pt.neighs[n].normal = normal;
	}
	}
	if (added_point_in_infinity) {
	count += 2;
	added_point_in_infinity = false;
	} else {
	count += 1;
	}
	}
	return;
	}

	/* ------------------------------------------------------------------------ */

	template <typename DomainA>
	UInt CADD_template::CollectingLocalDDpoints(
	std::vector<JaehyunDDPoint> &dd_points0_local, DomainA &domA) {

	bool problem_in_neigh = false;

	int count = 0;
	for (auto &node : domA.getContainer().getContainerNodes()) {
	Real idx;
	idx = (Real)node.getIndex();

	Real idx0;
	idx0 = (Real)node.getIndex0();

	LM_TOIMPLEMENT;
	// int num_neigh = node.getNbrOfNeighs();
	int num_neigh = std::numeric_limits<int>::max();

	Vector<3> crr_X = {roundf(node.position()[0] * 10.0) / 10.0,
	roundf(node.position()[1] * 10.0) / 10.0,
	roundf(node.position()[2] * 10.0) / 10.0};
	if (num_neigh >= 3) {
	DUMP("[collectingLocalDDpoints] problem in DD node idx"
	<< idx << " (" << crr_X << "): number of neighbors "
	<< num_neigh,
	DBG_MESSAGE);
	num_neigh = 2;
	problem_in_neigh = true;
	}
	dd_points0_local[count].idx = idx;
	dd_points0_local[count].num_neigh = (Real)num_neigh;
	dd_points0_local[count].X = crr_X;
	dd_points0_local[count].idx0 = (Real)idx0;

	for (int n = 0; n < num_neigh; ++n) {
	LM_TOIMPLEMENT;
	// Real nbr_idx = (Real)node.getNeighborNodeIndex(n);

	// Vector<Dim> nbr_X = {
	// roundf(node.getNeighborNodeCoord(n, 0) * 10.0) / 10.0,
	// roundf(node.getNeighborNodeCoord(n, 1) * 10.0) / 10.0,
	// roundf(node.getNeighborNodeCoord(n, 2) * 10.0) / 10.0};

	// Vector<Dim> burg = {node.burgersVectorX(n), node.burgersVectorY(n),
	// node.burgersVectorZ(n)};

	// Vector<Dim> normal = {node.normalVectorX(n), node.normalVectorY(n),
	// node.normalVectorZ(n)};

	// auto &neigh = dd_points0_local[count].neighs[n];

	// neigh.idx = nbr_idx;
	// neigh.burg = burg;
	// neigh.normal = normal;
	}
	}

	if (this->dd_points_out) {

	// Communicator &comm = Communicator::getCommunicator();
	UInt my_rank_in_DD = this->comm_group_dd.getMyRank();
	std::stringstream fname0;
	fname0 << this->directory_dd_points << "/local_dd_points0_step"
	<< current_step << ".txt.proc" << my_rank_in_DD;
	std::ofstream myfile0;
	for (auto &dd_pt0 : dd_points0_local) {
	if (!myfile0.is_open()) {
	myfile0.open(fname0.str().c_str());
	}
	myfile0 << dd_pt0 << "\t";
	}
	myfile0.close();
	}
	if (problem_in_neigh)
	return 1;
	else
	return 0;
	}

	/* --------------------------------------------------------------------------
	*/
	bool CADD_template::coarse_dislocation_loop(
	std::vector<JaehyunDDPoint> &dd_points,
	std::vector<JaehyunDDPoint> &dd_points0) {

	std::vector<UInt> indices;
	Vector<3> burg = {0.0, 0.0, 0.0};
	Vector<3> norm = {0.0, 0.0, 0.0};
	Real ci = 0.;
	UInt count = 0;

	UInt master_count = 0;
	for (UInt i = 0; i < dd_points0.size(); ++i) {
	if (!(std::find(indices.begin(), indices.end(), i) != indices.end())) {

	auto &dd_pt = dd_points0[i];
	Real ci0 = dd_pt.idx;
	auto cp0 = dd_pt.X;

	if (master_count >= dd_points0.size()) {
	DUMP("[coarse_dislocation_loop] exit 1", DBG_MESSAGE);
	return false;
	}

	auto &dd_pt_master = dd_points[master_count];
	auto &dd_pt_master_1 = dd_points[master_count - 1];
	auto &dd_pt_master_2 = dd_points[master_count - count];

	dd_pt_master.idx = ci0; // idx
	dd_pt_master.num_neigh = 2.0; // #nbr
	dd_pt_master.X = cp0; // x-coord

	master_count += 1;
	count += 1;

	ci = ci0;
	Vector<3> cp = cp0;

	indices.push_back(i);

	bool loop_start(true);
	bool not_finish(true);

	while (not_finish) {
	if (!loop_start) {
	if (is_it_closed(cp, cp0, ci, ci0)) {
	not_finish = false;

	if (master_count - 1 >= dd_points0.size()) {
	DUMP("[coarse_dislocation_loop] exit 2", DBG_MESSAGE);
	return false;
	}

	dd_pt_master_1.neighs[0].idx = ci0; // idx-nbr1
	dd_pt_master_1.neighs[0].X = cp0; // x-coord-nbr1
	dd_pt_master_1.neighs[0].burg = burg; // bx-nbr1
	dd_pt_master_1.neighs[0].normal = norm; // bx-nbr1

	dd_pt_master_2.neighs[1].idx = dd_pt_master_1.idx; // idx-nbr2
	dd_pt_master_2.neighs[1].X = dd_pt_master_1.X; // x-coord-nbr2
	dd_pt_master_2.neighs[1].burg = -1.0 * burg; // bx-nbr2
	dd_pt_master_2.neighs[1].normal = norm; // bx-nbr2
	count = 0;
	break;
	}
	}
	UInt count_inside = 0;
	bool straight(true);

	while (straight) {
	auto n1 = dd_pt.neighs[0].X;
	Real b1mag = dd_pt.neighs[0].burg[0] *
	cos(this->angle_between_b_and_x / 180.0 * M_PI) +
	dd_pt.neighs[0].burg[2] *
	sin(this->angle_between_b_and_x / 180.0 * M_PI);
	Real n1idx_real = dd_pt.neighs[1].idx;
	UInt n1idx = (UInt)n1idx_real;

	auto n2 = dd_pt.neighs[1].X;
	Real b2mag = dd_pt.neighs[1].burg[0] *
	cos(this->angle_between_b_and_x / 180.0 * M_PI) +
	dd_pt.neighs[1].burg[2] *
	sin(this->angle_between_b_and_x / 180.0 * M_PI);
	Real n2idx_real = dd_pt.neighs[2].idx;
	UInt n2idx = (UInt)n2idx_real;

	UInt saved_i = i;
	UInt tmp = find_one_direction(n1, n2, b1mag, b2mag, n1idx, n2idx, &i);
	if ((tmp == 0) && (master_count == 1)) {
	burg = 1.0 * dd_points0[saved_i].neighs[0].burg;
	} else if ((tmp == 1) && (master_count == 1)) {
	burg = -1.0 * dd_points0[saved_i].neighs[0].burg;
	}
	norm = dd_points0[saved_i].neighs[0].normal;

	// check whether this local connection is important or not
	Real area_in_yplane1 = AreaTriangle(cp, n1, n2, 1);

	Vector<3> last_registered_node = dd_pt_master_1.X;
	Vector<3> next_node = dd_pt.X;

	Real area_in_yplane2 =
	AreaTriangle(last_registered_node, cp, next_node, 1);

	if ((area_in_yplane1 > this->min_area) \|\|
	(area_in_yplane2 > this->min_area) \|\|
	(fabs(n2[1] - cp[1]) > 1.0) \|\|
	(fabs(n1[1] - cp[1]) > 1.0)) { // \|\| (dist > this->maxSeg0*10.0)){

	if (master_count - 1 >= dd_points0.size() - 1) {
	DUMP("[coarse_dislocation_loop] exit 3", DBG_MESSAGE);
	return false;
	}
	dd_pt_master.neighs[1].idx = dd_pt_master_1.idx; // idx-nbr2
	dd_pt_master.neighs[1].X = dd_pt_master_1.X; // x-coord-nbr2
	dd_pt_master.neighs[1].burg = -1.0 * burg; // bx-nbr2
	dd_pt_master.neighs[1].normal = norm; // nx-nbr2
	dd_pt_master_1.neighs[0].idx = ci; // idx-nbr1
	dd_pt_master_1.neighs[0].X = cp; // x-coord-nbr1
	dd_pt_master_1.neighs[0].burg = 1.0 * burg; // bx-nbr1
	dd_pt_master_1.neighs[0].normal = norm; // bx-nbr1
	if (!loop_start) {
	dd_pt_master.idx = ci; // idx
	dd_pt_master.num_neigh = 2.0; // #nbr
	dd_pt_master.X = cp; // x-coord
	master_count += 1;
	count += 1;
	}
	}
	loop_start = false;
	indices.push_back(i);

	ci = dd_pt.idx;
	cp = dd_pt.X;
	count_inside += 1;

	if (is_it_closed(cp, cp0, ci, ci0))
	break;
	if (count_inside >= dd_points0.size()) {
	DUMP("[coarse_dislocation_loop] exit 4", DBG_MESSAGE);
	return false;
	}
	}
	if (master_count >= dd_points0.size()) {
	DUMP("[coarse_dislocation_loop] exit 5", DBG_MESSAGE);
	return false;
	}
	}
	if (master_count >= dd_points0.size()) {
	DUMP("[coarse_dislocation_loop] exit 6", DBG_MESSAGE);
	return false;
	}
	}
	if (master_count >= dd_points0.size()) {
	DUMP("[coarse_dislocation_loop] exit 7", DBG_MESSAGE);
	return false;
	}
	}

	return true;
	}

	/* ------------------------------------------------------------------------ */
	template <typename DomainA> void CADD_template::SetPadAtoms(DomainA &domA) {
	UInt registered = 0;
	int my_rank_in_MDprocs = this->comm_group_atomic.getMyRank();
	LM_TOIMPLEMENT;
	// this->pad_atoms.clear();
	auto ref_manager = domA.getContainer().getRefManager();
	auto pad = GeometryManager::getManager().getGeometry(pad_geom);

	for (auto &at : domA.getContainer()) {
	if (pad->contains(at.position())) {
	LM_TOIMPLEMENT;
	// this->pad_atoms.push_back(at);
	registered += 1;
	}
	}
	this->comm_group_atomic.reduce(&registered, 1, "atoms: registered as PAD",
	OP_SUM);
	ref_manager->addSubSet("pad atoms", this->pad_atoms);
	if (my_rank_in_MDprocs == 0) {
	DUMP(registered << " atoms are registered as pad atoms", DBG_MESSAGE);
	}
	}

	/* --------------------------------------------------------------------------
	*/
	void CADD_template::find_indices_for_neighbors(
	std::vector<JaehyunDDPoint> &dd_points) {
	for (auto &dd_pt : dd_points) {
	if (dd_pt.num_neigh != 2)
	LM_FATAL("this is not a proper dislocation network");

	for (UInt j = 0; j < 2; ++j) {
	auto neigh = dd_pt.neighs[j];
	Real next_idx = neigh.idx;
	Vector<3> next_X = neigh.X;

	UInt k = 0;
	for (auto &dd_pt2 : dd_points) {
	Real node_idx = dd_pt2.idx;
	Vector<3> normal = dd_pt2.X;
	Real dist = (normal - next_X).norm();

	if ((dist < 1e-3) && (fabs(node_idx - next_idx) <= 1e-3)) {
	neigh.idx = (Real)k;
	// break;
	++k;
	}
	}
	}
	}
	}

	// -------------------------------------------------------------------------
	// */
	template <typename RefAtom>
	void CADD_template::applyDisplacement(RefAtom at, Vector<3> &dislo_disp) {

	auto disp_strain = Vector<3>::Zero();
	auto disp_tilt = Vector<3>::Zero();

	shear_strain(at, disp_strain);
	tilt_strain(at, disp_tilt);

	at.position() = at.position0() + dislo_disp - disp_strain + disp_tilt;
	}
	/* ------------------------------------------------------------------------ */
	void CADD_template::coordinate_for_2d_template(Real angle, Real local_coord[],
	Real global_coord[] [
	[gnu::unused]],
	Real template_coord[]) {
	template_coord[0] =
	local_coord[0] * local_coord[0] / fabs(local_coord[0] + 1e-8) *
	sin(angle * M_PI / 180.0) +
	local_coord[2] * local_coord[2] / fabs(local_coord[2] + 1e-8) *
	cos(angle * M_PI / 180.0);
	template_coord[1] =
	local_coord[1] * local_coord[1] / fabs(local_coord[1] + 1e-8);
	}

	/////// This part below should be removed later
	////////////////////////////////////
	/* ------------------------------------------------------------------------ */
	inline Real CADD_template::eval_q_coef(Real *X) {

	Real inside = X[0] * X[0] +
	2.0 * X[0] * X[1] * this->lambda * cos(this->phi) +
	X[1] * X[1] * this->lambda * this->lambda;
	return sqrt(inside);
	}

	/* ------------------------------------------------------------------------ */
	inline Real CADD_template::eval_t_coef(Real *X) {
	Real inside = X[0] * X[0] -
	2.0 * X[0] * X[1] * this->lambda * cos(this->phi) +
	X[1] * X[1] * this->lambda * this->lambda;
	return sqrt(inside);
	}

	/* ------------------------------------------------------------------------ */
	inline Real CADD_template::evalArctan(Real X, Real Y) {
	if ((X >= 0.0) && (Y > 0.0)) {
	return atan(Y / X);
	} else if ((X < 0.0) && (Y >= 0.0)) {
	return atan(-Y / -X) + M_PI;
	} else if ((X < 0.0) && (Y < 0.0)) {
	return atan(-Y / -X) - M_PI;
	} else if ((X >= 0.0) && (Y < 0.0)) {
	return atan(Y / X);
	} else {
	LM_FATAL("unexpected case");
	}
	}

	/* ------------------------------------------------------------------------ */
	inline Real CADD_template::computeBeDisplacementAniso(Real b_edge, Real *X) {

	Real q = eval_q_coef(X);
	Real t = eval_t_coef(X);
	Real c12 = this->elastic_coefs[1];
	Real c66 = this->elastic_coefs[12];

	Real first = 0.0;
	first += evalArctan((X[0] + X[1] * this->lambda * cos(this->phi)),
	(X[1] * this->lambda * sin(this->phi)));
	first += evalArctan((X[0] - X[1] * this->lambda * cos(this->phi)),
	(X[1] * this->lambda * sin(this->phi)));

	Real second = (this->c11bar * this->c11bar - c12 * c12) /
	(2.0 * this->c11bar * c66 * sin(2.0 * this->phi)) *
	floor(log(q / t) * 1e+5) / 1e+5;

	Real res = b_edge / (4.0 * M_PI) * (first + second);
	return res;
	}

	/* ------------------------------------------------------------------------ */
	inline Real CADD_template::computeNDisplacementAniso(Real b_edge, Real *X) {
	Real q = eval_q_coef(X);
	Real t = eval_t_coef(X);
	Real c12 = this->elastic_coefs[1];

	Real first =
	(this->c11bar - c12) * cos(this->phi) * floor(log(q * t) * 1e+5) / 1e+5;
	Real second =
	(this->c11bar + c12) * sin(this->phi) *
	atan2((X[1] * X[1] * this->lambda * this->lambda * sin(2.0 * this->phi)),
	(X[0] * X[0] -
	this->lambda * this->lambda * X[1] * X[1] * cos(2.0 * this->phi)));
	Real res = -this->lambda * b_edge /
	(4.0 * M_PI * this->c11bar * sin(2.0 * this->phi)) *
	(first - second);
	return res;
	}

	/* ------------------------------------------------------------------------ */

	inline Real CADD_template::computeLDisplacementAniso(Real b_screw, Real *X) {
	Real c44 = this->elastic_coefs[9];
	Real c45 = this->elastic_coefs[10];
	Real c55 = this->elastic_coefs[11];
	Real res = b_screw / 2.0 / M_PI *
	atan2(sqrt(c44 * c55 - c45 * c45) * X[1], c44 * X[0] - c45 * X[1]);
	return res;
	}

	/* --------------------------------------------------------------------------
	*/
	inline Real CADD_template::computeBeDisplacement(Real b_edge, Real *X) {
	Real x_2 = X[0] * X[0];
	Real y_2 = X[1] * X[1];
	Real nu = this->pois0;

	Real epsilon = 1e-15;
	Real res = b_edge / (2.0 * M_PI) *
	(atan2(X[1], X[0]) +
	X[0] * X[1] / (2.0 * (1 - nu) * (x_2 + y_2 + epsilon)));

	return res;
	}

	/* --------------------------------------------------------------------------
	*/
	inline Real CADD_template::computeNDisplacement(Real b_edge, Real *X) {
	Real x_2 = X[0] * X[0];
	Real y_2 = X[1] * X[1];
	Real epsilon = 1e-15;
	Real nu = this->pois0;
	// Real b_edge2 = b_edge*b_edge;
	// Real res = -b_edge/(2.0M_PI)((1-2.0nu)/(4.0(1-nu))*(log((x_2 +
	// y_2+epsilon)/b_edge2))
	// + (x_2 - y_2)/
	// (4.0(1-nu)(x_2 + y_2+epsilon))
	// );
	// return res;
	return -b_edge / (2.0 * M_PI) *
	((1 - 2.0 * nu) / (4.0 * (1 - nu)) * log(x_2 + y_2 + epsilon) +
	(x_2 - y_2) / (4.0 * (1 - nu) * (x_2 + y_2 + epsilon)));
	}

	/* --------------------------------------------------------------------------
	*/

	inline Real CADD_template::computeLDisplacement(Real b_screw, Real *X) {
	Real res = b_screw / (2.0 * M_PI) * (atan2(X[1], X[0]));
	return res;
	}

	/* --------------------------------------------------------------------------
	*/
	/* LMDESC CADDTEMPLATE
	CADDTEMPLATE is a simple CADD-like coupler which can deal with straing egde
	dislocations normal to the interface (i.e., 1;3400;0cno elastic problem)
	*/
	/* LMEXAMPLE
	COUPLING_CODE
	*/

	void CADD_template::declareParams() {
	/* LMKEYWORD PAD_GEOM
	Identifier of the pad geom
	*/
	this->parseKeyword("PAD_GEOM", this->pad_geom);
	/* LMKEYWORD DETECT_DXA
	Boolean to notify input from DXA results
	*/
	this->parseTag("DETECT_DXA", this->detect_dxa, false);
	/* LMKEYWORD MULTIPLE_LOOPS
	Boolean to notify multiple loops
	*/
	this->parseTag("MULTIPLE_LOOPS", this->multiple_loops, false);
	/* LMKEYWORD DETECT_PARTIAL
	Boolean to notify input for partially detection dislocation
	*/
	this->parseTag("DETECT_PARTIAL", this->detect_partial, false);
	/* LMKEYWORD PRINT_DD_POINTS
	Boolean to print internal data named dd\_points and its offset lines
	*/
	this->parseTag("PRINT_DD_POINTS", this->dd_points_out, false);
	/* LMKEYWORD COARSEN_DD
	Boolean to decide to coarse dislocation loop before apply Barnett fields
	*/
	this->parseTag("COARSEN_DD", this->coarsen_dd, false);
	/* LMKEYWORD MIN_AREA
	Specifies the minimum area to be coarsed out
	*/
	this->parseKeyword("MIN_AREA", this->min_area, 0.0);
	/* LMKEYWORD DIRECTORY_DD_POINTS
	fies the directory dd points are written
	*/
	this->parseKeyword("DIRECTORY_DD_POINTS", this->directory_dd_points, "/.");
	/* LMKEYWORD DD_BOUNDARY
	Boolean to set DD boundary
	*/
	this->parseTag("DD_BOUNDARY", this->dd_boundary, false);
	/* LMKEYWORD SLIP_PLANE
	Specifies the normal vector to the slip plane
	*/
	this->parseVectorKeyword("SLIP_PLANE", 3, this->normal_quant);
	/* LMKEYWORD ANGLE_BETWEEN_B_AND_X
	Specifies the angle between Burgers vector and X(11-2) axis
	*/
	this->parseKeyword("ANGLE_BETWEEN_B_AND_X", this->angle_between_b_and_x);
	/* LMKEYWORD MD_GEOM
	Identifier of the MD geometry
	*/
	this->parseKeyword("MD_GEOM", this->md_geom);
	/* LMKEYWORD DD_GEOM
	Identifier of the DD geometry
	*/
	this->parseKeyword("DD_GEOM", this->dd_geom);
	/* LMKEYWORD UPDATE_FREQUENCY
	Specifies the frequency at which detection computations are performed
	*/
	this->parseKeyword("UPDATE_FREQUENCY", this->detect_freq, 50);
	/* LMKEYWORD DETECT_COMPUTE
	Add computes to input to filter
	*/
	this->parseKeyword("DETECT_COMPUTE", compute_list);
	/* LMKEYWORD NUM_TEMPLATES
	Specifies the number of templates provided.
	*/
	this->parseKeyword("NUM_TEMPLATES", this->num_templates);
	/* LMKEYWORD TEMPLATE_FILE
	Specifies the template file.
	*/
	this->parseKeyword("TEMPLATE_FILE", this->templatefilename_list);
	/* LMKEYWORD CORE_RCUT
	Specifies the radius of the region where core template applied
	*/
	this->parseKeyword("CORE_RCUT", this->core_rcut);
	/* LMKEYWORD RATIO_OF_PURE_CORE_TEMPLATE
	Specifies the ratio of region where only pure core template is applied
	with respect to ellipse defined by RMAX and RMIN
	where
	*/
	this->parseKeyword("RATIO_OF_PURE_CORE_TEMPLATE", this->factor);
	/* LMKEYWORD SHEAR_STRESSES
	Specifies the shear stersses with respect to Y plane.
	*/
	this->parseVectorKeyword("SHEAR_STRESSES", 4, this->shear_stresses);
	/* LMKEYWORD TILT_STRESSES
	Specifies the shear stersses imposed by tilting.
	*/
	this->parseVectorKeyword("TILT_STRESSES", 3, this->tilt_stresses);
	/* LMKEYWORD SHEAR_MODULUS
	Specifies the shear modulus in Burgers vector direction.
	mega-Pascal unit
	*/
	this->parseKeyword("SHEAR_MODULUS", this->mu0);
	/* LMKEYWORD POISSON_RATIO
	Specifies the shear modulus in Burgers vector direction.
	mega-Pascal unit
	*/
	this->parseKeyword("POISSON_RATIO", this->pois0);

	/* LMKEYWORD DOMAIN_DD
	Specifies the dd domain id
	*/
	this->parseKeyword("DOMAIN_DD", this->domdd);

	/* LM_KEYWORD DISTANCE_DDLINE_FROM_PAD
	Specify the distance from pad (inner) boudaries to DD constrained nodes.
	*/
	this->parseKeyword("DISTANCE_DDLINE_FROM_PAD",
	this->distance_ddline_from_pad);
	/* LMKEYWORD PAD_THICKNESS
	Specify the thickness of pad zone
	*/
	this->parseKeyword("PAD_THICKNESS", this->pad_thickness);

	// In case of (0, 0, 1):
	// point1 point2
	// ----------------------
	// \| \| z
	// \| \| ^
	// \| \| \|
	// o-o-o-x x-o-o-o-o +--->x
	/* LMKEYWORD CLOSE_IN_INFINITY
	Specifies the direction in order to close
	dislocation loop
	*/
	this->parseVectorKeyword("CLOSE_IN_INFINITY", 3, this->close_direction,
	VEC_DEFAULTS(0, 0, 0));
	/* LMKEYWORD NUBER_INFINITY_CONNECTION
	Specifies the direction in order to close
	dislocation loop
	*/
	this->parseKeyword("NUBER_INFINITY_CONNECTION",
	this->num_close_loop_in_infinity, 0);
	/* LMKEYWORD MD_BOX_LENGTHS
	Specifies the (original) MD box size
	to specify correct tilt stress values
	*/
	this->parseVectorKeyword("MD_BOX_LENGTHS", 3, this->md_box_lengths);
	/* LMKEYWORD DETECT_SIMPLE_DIRECTION
	Specifies the (original) MD box size
	to specify correct tilt stress values
	*/
	this->parseVectorKeyword("DETECT_SIMPLE_DIRECTION", 3,
	this->detect_simple_direction,
	VEC_DEFAULTS(0, 0, 0));
	// /* LM**KEYWORD NB_SLIPS
	// Specifies the number of slip planes
	// */
	// this->parseKeyword("NB_SLIPS", this->nb_slips);
	// /* LM**KEYWORD LEN_INTER_SLIPS
	// Specifies length between two slip planes
	// */
	// this->parseKeyword("LEN_INTER_SLIPS", this->len_inter_slips);

	/* LMKEYWORD ANISO
	Boolean to set anisotropicity:
	[1] Elastic Strain Fields and Dislocation Mobility, Volume 31, J.Lothe
	Dislocations in anisotropic media
	[2] Plane-strain discrete dislocation plasticity incorporating
	anisotropic
	elasticity IJSS 48(2), January 2011
	*/
	this->parseTag("ANISO", this->aniso, false);
	/* LMKEYWORD VOLTERRA
	TODO
	*/
	this->parseTag("VOLTERRA", this->volterra, false);

	/* LMKEYWORD POSITION_TOLERANCE
	The position of the nodes found from the detection
	have to be imposed to the current dislocation nodes.
	This parameter specify the distance tolerance for which nodes
	in the previous configuration are mapped to the detected ones.
	*/
	this->parseKeyword("POSITION_TOLERANCE", this->position_tolerance);

	/* LMKEYWORD ELASTIC_COEFS
	Specify the thirteen elastic coefficients
	of the constitutive matrix of anisotropic media in the following order:
	C11, C12, C13, C21, C22, C23, C31, C32, C33, C44, C45, C55, C66
	*/
	this->parseVectorKeyword("ELASTIC_COEFS", 13, this->elastic_coefs,
	VEC_DEFAULTS(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0));

	LM_TOIMPLEMENT;
	// forward parsable of the core template needs to be done
	// this->addSubParsableObject(this->core_template);
	}

	/* ---------------------------------------------------------------------- */

	template <typename DomainA>
	std::vector<JaehyunDDPoint> CADD_template::sendDDSegments2MD(DomainA &domA) {

	UInt data_length = 33;
	std::vector<JaehyunDDPoint> dd_points;
	std::vector<JaehyunDDPoint> dd_points0;
	UInt num_dd_pts;
	UInt num_dd_pts0;

	if (this->is_in_dd) {
	LM_TOIMPLEMENT;
	// UInt my_rank_in_DD = this->comm_group_dd.groupRank(lm_my_proc_id);
	// UInt nbDDprocs = this->comm_group_dd.getNBprocsOnGroup();

	// declare data array for (unclosed) dd points in a local processor
	UInt num_dd_pts0_local = domA.getNBNodes();
	std::vector<JaehyunDDPoint> dd_points0_local(num_dd_pts0_local);
	// collecting dd points
	UInt collecting = 0;
	collecting = CollectingLocalDDpoints(dd_points0_local, domA);
	this->comm_group_dd.reduce(
	&collecting, 1,
	"boolean to detect any errors occurring in collecting local nodes",
	OP_SUM);

	// evaluate total number of points
	UInt num_dd_pts0_global = 0;
	LM_TOIMPLEMENT;
	// MPI_Reduce(&num_dd_pts0_local, &num_dd_pts0_global, 1, MPI_INT,
	// MPI_SUM,
	// comm.realRank(0, ),
	// comm.getMpiGroup(this->comm_group_dd));

	// declare data array for (unclosed) total dd points only at 0 processor
	std::vector<JaehyunDDPoint> dd_points0_global;

	LM_TOIMPLEMENT;
	// if (my_rank_in_DD == 0) {
	// (unclosed) total dd points
	// dd_points0_global.resize(num_dd_pts0_global);
	//}

	// declare data list for number of dd nodal points for all the DDD
	// processors
	std::vector<int> num_dd_pts0_local_list;

	LM_TOIMPLEMENT;
	// if (my_rank_in_DD == 0)
	// num_dd_pts0_local_list.resize(nbDDprocs);

	// gather the total local data length into the declared data list
	LM_TOIMPLEMENT;
	// int local_length = num_dd_pts0_local * local_data_length;
	// MPI_Gather(&local_length, 1, MPI_INT, &num_dd_pts0_local_list[0], 1,
	// MPI_INT, comm.realRank(0, this->comm_group_dd),
	// comm.getMpiGroup(this->comm_group_dd));

	// declear data array for starting number
	std::vector<int> lengths_dd_points_local;
	LM_TOIMPLEMENT;
	// UInt totlen = 0;
	// if (my_rank_in_DD == 0) {
	// lengths_dd_points_local.resize(nbDDprocs);
	// lengths_dd_points_local[0] = 0;
	// totlen += num_dd_pts0_local_list[0];
	// for (UInt i = 1; i < nbDDprocs; i++) {
	// lengths_dd_points_local[i] = lengths_dd_points_local[i - 1] +
	// totlen; totlen = num_dd_pts0_local_list[i];
	// }
	// }

	// gather local dd points to global points
	// inputs: dd_points0_local
	// local_length
	// num_dd_pts0_local_list
	// lengths_dd_points_local

	LM_TOIMPLEMENT;
	// MPI_Gatherv(&dd_points0_local[0], // local dd points
	// local_length, // local length
	// MPI_REAL,
	// &dd_points0_global[0], // global dd points
	// &num_dd_pts0_local_list[0], // number of length list
	// &lengths_dd_points_local[0], // starting number list
	// MPI_REAL, comm.realRank(0, this->comm_group_dd),
	// comm.getMpiGroup(this->comm_group_dd));

	LM_TOIMPLEMENT;
	int my_rank_in_DD = std::numeric_limits<int>::max();
	if (my_rank_in_DD == 0) {

	UInt num_pairs_in_infinity = domA.getNBpairsInfinite();
	UInt num_dd_pts0 = num_dd_pts0_global;
	if (this->close_in_infinity)
	num_dd_pts0 += (num_pairs_in_infinity + num_close_loop_in_infinity) * 2;

	// (closed) total dd nodal points
	dd_points0.resize(num_dd_pts0);

	close_dislocation_loop(dd_points0, dd_points0_global, domA);

	find_indices_for_neighbors(dd_points0);

	num_dd_pts = num_dd_pts0;
	// (closed) coarsened dd nodal points
	dd_points.resize(num_dd_pts);

	STARTTIMER("coarsening");
	if ((coarsen_dd) and (collecting == 0)) {
	bool coarsen = coarse_dislocation_loop(dd_points, dd_points0);
	if (!coarsen) {
	// in case of fail, copy all dd_points0 into dd_points
	DUMP("[WARNING] Fail to coarse dislocation loop. Total "
	<< num_dd_pts0 << " points",
	DBG_MESSAGE);
	for (UInt i = 0; i < num_dd_pts; ++i) {
	for (UInt j = 0; j < data_length; ++j) {
	dd_points[i * data_length + j] = dd_points0[i * data_length + j];
	}
	}
	} else {
	// in case of success,
	DUMP("Coarse_dislocation_loop end. Reduced from "
	<< num_dd_pts0 << " points to " << num_dd_pts << " points",
	DBG_MESSAGE);
	find_indices_for_neighbors(dd_points);
	}
	} else {
	DUMP("[WARNING] Coarse_dislocation_loop is skipped. coarsen_dd: "
	<< coarsen_dd << " collecting: " << collecting << " Total "
	<< num_dd_pts0 << " points",
	DBG_MESSAGE);
	for (UInt i = 0; i < num_dd_pts; ++i) {
	for (UInt j = 0; j < data_length; ++j) {
	dd_points[i * data_length + j] = dd_points0[i * data_length + j];
	}
	}
	}

	STOPTIMER("coarsening");
	ComputeOffsetParallelCurve(dd_points);
	// CorrectIntersetingPoints(num_dd_pts,dd_points);

	UInt mdNBprocs = this->comm_group_atomic.size();
	for (UInt i = 0; i < mdNBprocs; ++i) {
	// int real_rank = this->comm_group_atomic.realRank(i);
	int my_rank_in_group = this->comm_group_atomic.getMyRank();
	this->comm_group_atomic.send(&num_dd_pts, 1, my_rank_in_group,
	"number of dd points");
	this->comm_group_atomic.send(&num_dd_pts0, 1, my_rank_in_group,
	"number of dd points0");
	this->comm_group_atomic.send(&(dd_points[0]), num_dd_pts0 * data_length,
	my_rank_in_group, "dd_points (vector)");
	}
	}

	// BARRIER
	this->comm_group_dd.synchronize();
	}

	if (this->is_in_atomic) {

	STARTTIMER("SetPadAtoms");
	SetPadAtoms(domA);
	STOPTIMER("SetPadAtoms");

	STARTTIMER("ReceiveDDdata");
	this->comm_group_dd.receive(&num_dd_pts, 1, 0, "number of dd points");
	this->comm_group_dd.receive(&num_dd_pts0, 1, 0, "number of dd points0");

	dd_points.resize(num_dd_pts0 * data_length);

	this->comm_group_dd.receive(&(dd_points[0]), num_dd_pts0 * data_length, 0,
	"dd_points (vector)");
	STOPTIMER("ReceiveDDdata");

	std::vector<UInt> indices;
	// UInt count_pad_core = 0;
	// UInt count_pad_around = 0;
	// UInt count_far_field = 0;
	// UInt count_true_disp = 0;
	}

	return dd_points;
	}

	/* ---------------------------------------------------------------------- */

	template <typename RefAtom>
	void CADD_template::functionifAniso(
	UInt num_dislo, Vector<3> &pos0_pbc, Vector<3> &pos0_pbc0,
	std::vector<std::vector<Real>> &dislocations, RefAtom &at) {

	Vector<3> analytic_Disp = {0.0, 0.0, 0.0};
	for (UInt nd = 0; nd < num_dislo; ++nd) {
	for (UInt pd = 0; pd < 2; ++pd) {
	if (pd == 0) {
	// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
	pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
	pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] - 1.0);
	} else if (pd == 1) {
	// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
	pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
	pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] + 1.0);
	}

	Real rotation_matrix[] = {0.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0};
	Real pos0_pbc_rotate[] = {0.0, 0.0, 0.0};
	for (UInt i = 0; i < 3; ++i) {
	pos0_pbc_rotate[i] = 0.0;
	for (UInt j = 0; j < 3; ++j) {
	pos0_pbc_rotate[i] += rotation_matrix[3 * i + j] * pos0_pbc[j];
	}
	}

	Real b_screw =
	this->burgMag0 * cos(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;
	Real b_edge =
	this->burgMag0 * sin(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;

	Real analytic_disp[] = {0.0, 0.0, 0.0};
	Real dispV[] = {0.0, 0.0, 0.0};
	analytic_disp[0] = computeBeDisplacementAniso(b_edge, pos0_pbc_rotate);
	analytic_disp[1] = computeNDisplacementAniso(b_edge, pos0_pbc_rotate);
	analytic_disp[2] = computeLDisplacementAniso(b_screw, pos0_pbc_rotate);

	for (UInt i = 0; i < 3; ++i) {
	for (UInt j = 0; j < 3; ++j) {
	analytic_Disp[i] += rotation_matrix[3 * j + i] * analytic_disp[j];
	}
	}
	for (UInt i = 0; i < 3; ++i)
	analytic_Disp[i] += dispV[i];
	}
	}
	applyDisplacement(at, analytic_Disp);
	}

	/* ------------------------------------------------------------------------ */
	template <typename RefAtom>
	void CADD_template::functionifVolterra(
	UInt num_dislo, Vector<3> &pos0_pbc, Vector<3> &pos0_pbc0,
	std::vector<std::vector<Real>> &dislocations, RefAtom &at) {

	Vector<3> analytic_Disp = {0.0, 0.0, 0.0};

	for (UInt nd = 0; nd < num_dislo; ++nd) {

	for (UInt pd = 0; pd < 2; ++pd) {
	if (pd == 0) {
	// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
	pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
	pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] - 1.0);
	} else if (pd == 1) {
	// pos0_pbc[0] = pos0_pbc0[0] - dislocations[nd][0];
	pos0_pbc[1] = pos0_pbc0[1] - dislocations[nd][1];
	pos0_pbc[2] = pos0_pbc0[2] - (dislocations[nd][2] + 1.0);
	}

	Real rotation_matrix[] = {0.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0};
	Real pos0_pbc_rotate[] = {0.0, 0.0, 0.0};
	for (UInt i = 0; i < 3; ++i) {
	pos0_pbc_rotate[i] = 0.0;
	for (UInt j = 0; j < 3; ++j) {
	pos0_pbc_rotate[i] += rotation_matrix[3 * i + j] * pos0_pbc[j];
	}
	}

	Real b_screw =
	this->burgMag0 * cos(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;
	Real b_edge =
	this->burgMag0 * sin(this->angle_between_b_and_x / 180.0 * M_PI) / 2.;
	Real analytic_disp[] = {0.0, 0.0, 0.0};
	Real dispV[] = {0.0, 0.0, 0.0};
	analytic_disp[0] = computeBeDisplacement(b_edge, pos0_pbc_rotate);
	analytic_disp[1] = computeNDisplacement(b_edge, pos0_pbc_rotate);
	analytic_disp[2] = computeLDisplacement(b_screw, pos0_pbc_rotate);

	for (UInt i = 0; i < 3; ++i) {
	for (UInt j = 0; j < 3; ++j) {
	dispV[i] += rotation_matrix[3 * j + i] * analytic_disp[j];
	}
	}
	for (UInt i = 0; i < 3; ++i)
	analytic_Disp[i] += dispV[i];
	}
	}
	applyDisplacement(at, analytic_Disp);
	}
	/* ------------------------------------------------------------------------ */

	template <typename RefAtom>
	void CADD_template::functionElse(std::vector<JaehyunDDPoint> &dd_points,
	Vector<3> &pos0_pbc, RefAtom &at) {

	UInt count_pad_core = 0;
	UInt count_far_field = 0;
	UInt count_pad_around = 0;

	Vector<3> barnett_disp = {0.0, 0.0, 0.0};
	Vector<3> dislo_disp = {0.0, 0.0, 0.0};
	Vector<3> temp_disp = {0.0, 0.0, 0.0};

	for (auto &dd_pt : dd_points) {
	// bool selected = false;
	Vector<3> disp_for_segment = {0.0, 0.0, 0.0};
	auto &A = dd_pt.X;
	auto &neigh1 = dd_pt.neighs[0];
	auto &neigh2 = dd_pt.neighs[2];
	auto &pts1 = neigh1.X;
	auto &pts2 = neigh2.X;
	for (UInt j = 0; j < 2; ++j) {
	auto &neigh = dd_pt.neighs[j];
	Vector<3> burg = neigh.burg * this->burgMag0;
	Real bmag = burg[0] * cos(this->angle_between_b_and_x / 180.0 * M_PI) +
	burg[2] * sin(this->angle_between_b_and_x / 180.0 * M_PI);
	bmag = roundf(bmag * 10.0) / 10.0;

	if (bmag / fabs(bmag) > 0) {
	auto &B = neigh.X;
	auto &normal = neigh.normal;
	UInt slip_case = 3;
	if (fabs(normal[0]) - 1.0 == 0.0)
	slip_case = 0;
	else if (fabs(normal[1]) - 1.0 == 0.0)
	slip_case = 1;
	else if (fabs(normal[2]) - 1.0 == 0.0)
	slip_case = 2;
	else {
	DUMP("[warning] cannot define slip case", DBG_MESSAGE);
	slip_case = 1;
	}

	ComputeAnalyticDisplacements(disp_for_segment, B, A, burg, pos0_pbc);

	UInt idx_number = (UInt)neigh.idx;
	auto &pts3 = dd_points[idx_number].tmp_X;
	auto &pts4 = dd_points[idx_number].tmp2_X;
	if (PointInSquare(pts1, pts2, pts3, pts4, pos0_pbc, slip_case)) {
	Real distance = distance_between_line_and_point3D(pos0_pbc, A, B);
	if (distance <= this->core_rcut) {
	bool success(false);
	success = ComputeTemplateDisplacements(temp_disp, A, B, burg,
	normal, pos0_pbc);
	if (!success) {
	barnett_disp += disp_for_segment;
	} else {
	if (distance / this->core_rcut <= this->factor) {
	++count_pad_core;
	barnett_disp += temp_disp + disp_for_segment;
	} else {
	Real temp_ratio =
	(this->core_rcut - distance) /
	(this->core_rcut - this->core_rcut * this->factor);
	++count_pad_around;
	barnett_disp += temp_disp * temp_ratio + disp_for_segment;
	}
	}
	} else {
	barnett_disp += disp_for_segment;
	++count_far_field;
	}

	} else {
	barnett_disp += disp_for_segment;
	++count_far_field;
	} // if(PointInSquare(pts ... )
	} // if(bmag/abs(bmag)>0)
	} // for (UInt j=0; j < num_neigh; ++j){
	} // for (UInt i=0; i < num_dd_pts;
	dislo_disp = barnett_disp;
	applyDisplacement(at, dislo_disp);
	}
	/* ------------------------------------------------------------------------ */

	DECLARE_COUPLER_INIT_MAKE_CALL(CADD_template, domDD, domMD)
	/* ------------------------------------------------------------------------ */

	__END_LIBMULTISCALE__

template_cadd.ccNo OneTemporaryActions

File Metadata

template_cadd.ccView Options

Event Timeline

template_cadd.cc
No OneTemporary
Actions

template_cadd.cc
View Options