domain_paradis.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, May 27, 08:32

domain_paradis.cc
View Options

	/**
	* @file domain_paradis.cc
	*
	* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
	* @author Till Junge <till.junge@epfl.ch>
	* @author Moseley Philip Arthur <philip.moseley@epfl.ch>
	* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
	* @author Max Hodapp <max.hodapp@epfl.ch>
	*
	* @date Fri Jul 11 15:47:44 2014
	*
	* @brief ParaDiS domain
	*
	* @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/>.
	*
	*/

	/* TODO:
	* - the function addHybSegments() is currently a rudimentary implementation
	* for simulations containing only one dislocations
	* ==> method has to be generalized and additional checks (e.g. if hybrid
	* segments are actually self-consistent, etc...) have to be performed to
	* improve robustness
	*
	*/
	/* -------------------------------------------------------------------------- */
	#define TIMER
	/* -------------------------------------------------------------------------- */
	#include "lm_common.hh"
	#define PARALLEL
	#include "communicator.hh"
	#include "domain_paradis.hh"
	#include "geometry_manager.hh"
	#include "stdio.h"
	/* -------------------------------------------------------------------------- */

	// ParaDiS includes
	extern "C" {
	// void TimerClearAll(Home_t *home);
	#include "Home.h"
	#include "Param.h"
	#include "Typedefs.h"
	/* -------------------------------------------------------------------------- */
	Home_t *InitHome();
	void ParadisStep(Home_t *home);
	void RemeshAfterCoupling(Home_t *home);
	void PointDisplacements(Home_t *home);
	void SegmentDisplacementsOfAtom(libmultiscale::Vector<3u> &A,
	libmultiscale::Vector<3u> &B,
	libmultiscale::Vector<3u> &C,
	libmultiscale::Vector<3u> &burg,
	libmultiscale::Vector<3u> &dispPointU,
	libmultiscale::Vector<3u> &dispPointX,
	libmultiscale::Real pois);
	void ComputeClosurePoint(libmultiscale::Vector<3u> &A,
	libmultiscale::Vector<3u> &B,
	libmultiscale::Vector<3u> &burg,
	libmultiscale::Vector<3u> &C);
	void ComputeProjectedPoint(libmultiscale::Vector<3u> &A,
	libmultiscale::Vector<3u> &C,
	libmultiscale::Vector<3u> &burg,
	libmultiscale::Vector<3u> &P);

	void GetFieldPointStress(Home_t *home, real8 x, real8 y, real8 z,
	real8 totStress[3][3]);
	#ifdef FPES_ON
	#include <fpcontrol.h>
	#endif
	}

	#undef PARALLEL

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

	__BEGIN_LIBMULTISCALE__

	void DomainParaDiS::remesh() {
	RemeshAfterCoupling(this->paradis_ptr);

	// //ParadisStep.c from line 336 to 418
	// InitTopologyExemptions(this->paradis_ptr);
	// ClearOpList(this->paradis_ptr);
	// SortNodesForCollision(this->paradis_ptr);
	// SplitMultiNodes(this->paradis_ptr);

	// CrossSlip(this->paradis_ptr);
	// HandleCollisions(this->paradis_ptr);
	// ProximityCollisions(this->paradis_ptr);

	// #ifdef PARALLEL
	// #ifdef PARADIS_IN_LIBMULTISCALE
	// CommSendRemesh(this->paradis_ptr);
	// #endif
	// #endif

	// FixRemesh(this->paradis_ptr);

	// Remesh(this->paradis_ptr);
	// Migrate(this->paradis_ptr);

	// #ifdef PARALLEL
	// #ifdef PARADIS_IN_LIBMULTISCALE
	// RecycleGhostNodes(this->paradis_ptr);
	// #endif
	// #endif

	// SortNativeNodes(this->paradis_ptr);

	// #ifdef PARALLEL
	// #ifdef PARADIS_IN_LIBMULTISCALE
	// CommSendGhosts(this->paradis_ptr);
	// #endif
	// #endif
	}

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

	void DomainParaDiS::computeDisplacements(std::vector<Real> &positions,
	std::vector<Real> &displacements) {
	if (positions.empty()) {
	displacements.clear();
	return;
	}
	if (!this->paradis_ptr) {
	return;
	}
	// Pass data to Paradis.
	unsigned int size = positions.size();
	LM_ASSERT(size % 3 == 0, "paradis computeDisplacements requires 3d.");
	if (displacements.size() != size)
	displacements.resize(size);
	this->paradis_ptr->dispPointCount = size / 3;
	this->paradis_ptr->dispPointX = &positions[0];
	this->paradis_ptr->dispPointU = &displacements[0];
	// Calculate the displacements.
	PointDisplacements(this->paradis_ptr);
	}
	/* -------------------------------------------------------------------------- */

	void DomainParaDiS::computeBarnett(Vector<3u> &A, Vector<3u> &B, Vector<3u> &C,
	Vector<3u> &burg, Vector<3u> &dispPointU,
	Vector<3u> &dispPointX, Real pois) {
	SegmentDisplacementsOfAtom(A, B, C, burg, dispPointU, dispPointX, pois);
	}

	/* -------------------------------------------------------------------------- */
	void DomainParaDiS::computeClosurePoint(Vector<3u> &A, Vector<3u> &B,
	Vector<3u> &burg, Vector<3u> &C) {
	ComputeClosurePoint(A, B, burg, C);
	}
	/* -------------------------------------------------------------------------- */
	void DomainParaDiS::computeProjectedPoint(Vector<3u> &A, Vector<3u> &C,
	Vector<3u> &burg, Vector<3u> &P) {
	ComputeProjectedPoint(A, C, burg, P);
	}
	/* -------------------------------------------------------------------------- */

	void DomainParaDiS::init() {

	auto &group = this->getCommGroup();
	if (!group.amIinGroup())
	return;
	const UInt nbprocs = group.size();

	// Mostly taken from ParaDiS source files Main.c and ParadisInit.c
	std::string dlb_flag = "-r";
	std::string datafile_flag = "-d";
	std::string ctrlfile_flag = "-c";
	std::string fake_argv0 = "IGNORED";
	std::string dlb_num;
	std::stringstream(dlb_num) << numDLBCycles;

	std::vector<const char *> argv({fake_argv0.c_str(), datafile_flag.c_str(),
	data_filename.c_str(), dlb_flag.c_str(),
	dlb_num.c_str(), ctrlfile_flag.c_str(),
	control_filename.c_str()});

	if (mobility_filename != "")
	argv.push_back(mobility_filename.c_str());

	unsigned int argc = argv.size();

	/*
	* On some systems, the getrusage() call made by Meminfo() to get
	* the memory resident set size does not work properly. In those
	* cases, the function will try to return the current heap size
	* instead. This initial call allows meminfo() to get a copy of
	* the original heap pointer so subsequent calls can calculate the
	* heap size by taking the difference of the original and current
	* heap pointers.
	*/
	Meminfo(&(this->memSize));

	/*
	* On linux systems (e.g. MCR) if built to have floating point exceptions
	* turned on, invoke macro to do so
	*/
	#ifdef FPES_ON
	unmask_std_fpes();
	#endif

	this->paradis_ptr = InitHome();
	TimerInit(this->paradis_ptr);
	TimerStart(this->paradis_ptr, TOTAL_TIME);
	TimerStart(this->paradis_ptr, INITIALIZE);

	if (group.getMyRank() == 0)
	std::cout << "BEGIN: Initialize ParaDiS on " << nbprocs << " processors."
	<< std::endl;
	Initialize(this->paradis_ptr, argc, const_cast<char **>(&argv[0]),
	group.getMPICommGroup());

	if (strcmp(this->paradis_ptr->param->timestepIntegrator,
	"trapezoidtimestep") != 0) {
	LM_FATAL("The timestep integrator for ParaDiS needs to be "
	<< "\"trapezoidtimestep\" in order for ParaDis to work with "
	<< "LibMultiscale. You specified \"timestepIntegrator = "
	<< this->paradis_ptr->param->timestepIntegrator
	<< "\" in the control file \"" << control_filename << "\".");
	}
	TimerStop(this->paradis_ptr, INITIALIZE);

	if (nbprocs > 1)
	MPI_Barrier(this->paradis_ptr->MPI_COMM_PARADIS);

	this->paradis_ptr->param->timeMax = this->timestep;
	this->paradis_ptr->cycle = this->paradis_ptr->param->cycleStart;

	/*
	* Perform the needed number (if any) of load-balance-only
	* steps before doing the main processing loop. These steps
	* perform only the minimal amount of stuff needed to
	* estimate per-process load, move boundaries and migrate
	* nodes among processsors to get a good initial balance.
	*/
	time_t tp;
	int initialDLBCycles = this->paradis_ptr->param->numDLBCycles;
	TimerStart(this->paradis_ptr, INITIALIZE);
	if ((this->paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
	time(&tp);
	printf(" +++ Beginning %d load-balancing steps at %s", initialDLBCycles,
	asctime(localtime(&tp)));
	}
	while (this->paradis_ptr->param->numDLBCycles > 0) {
	ParadisStep(this->paradis_ptr);
	this->paradis_ptr->cycle++;
	this->paradis_ptr->param->numDLBCycles--;
	}
	if ((this->paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
	time(&tp);
	printf(" +++ Completed load-balancing steps at %s",
	asctime(localtime(&tp)));
	}
	TimerStop(this->paradis_ptr, INITIALIZE);

	/*
	* Any time spent doing the initial DLB-only steps should
	* just be attributed to initialization time, so be sure to
	* reset the other timers before going into the main
	* computational loop
	*/
	TimerInitDLBReset(this->paradis_ptr);
	/*
	* The cycle number has been incremented during the initial
	* load-balance steps, and previous steps so reset it to the
	* proper startingvalue before entering the main processing loop.
	* In paradis, cycle uses 1-based indexing.
	*/
	this->paradis_ptr->cycle = this->paradis_ptr->param->cycleStart;
	this->paradis_ptr->cycleLMS = 1;

	this->container_dd.getContainerNodes().setGrain(this->paradis_ptr);
	this->container_dd.getContainerElems().setGrain(this->paradis_ptr);
	this->container_dd.getContainerElems().setNodes(
	this->container_dd.getContainerNodes());

	#ifndef LM_OPTIMIZED
	// Debug tests.
	unsigned int counter = 0;
	DUMP("Printing all initial ParaDiS nodes...", DBG_INFO);
	for (auto &&node : this->container_dd.getContainerNodes()) {
	DUMP("Processor " << lm_my_proc_id << ", Node " << ++counter << "("
	<< node.position() << ")",
	DBG_INFO);
	}
	#endif

	// ContainerParaDiSDisplacement displs = this->getContainerDisplacements();
	// DUMP("Number of ParaDiS displacement fieldpoints: "<<displs.size(),
	// DBG_INFO);
	// End debug tests.

	if (this->paradis_ptr->myDomain == 0)
	std::cout << "END: Initialize ParaDiS on " << nbprocs << " processors."
	<< std::endl;
	}

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

	void DomainParaDiS::performStep1() {
	this->paradis_ptr->param->timeLMSCycle = 0.;
	while (this->paradis_ptr->param->timeLMSCycle <
	this->paradis_ptr->param->timeMax) {
	STARTTIMER("ParadisStep");
	ParadisStep(this->paradis_ptr);
	TimerClearAll(this->paradis_ptr);
	STOPTIMER("ParadisStep");
	}
	this->paradis_ptr->cycleLMS++;
	this->getOutput()->incRelease();

	if (this->boundary) {
	auto c =
	GeometryManager::getManager().getGeometry(this->geom)->getBoundingBox();
	auto &&maxval = c.getXmax();
	auto &&minval = c.getXmin();

	// Geometry* dd = GeometryManager::getManager().getGeometry(this->geom);
	// Cube & dd_cube = dynamic_cast<Cube&>(*dd);

	// to mimic strong grain boundary...
	// Real xmin = this->paradis_ptr->param->xBoundMin + 10.0;
	// Real xmax = this->paradis_ptr->param->xBoundMax - 10.0;
	// Real ymin = this->paradis_ptr->param->yBoundMin + 10.0;
	// Real ymax = this->paradis_ptr->param->yBoundMax - 10.0;
	// Real zmin = this->paradis_ptr->param->zBoundMin + 10.0;
	// Real zmax = this->paradis_ptr->param->zBoundMax - 10.0;

	// Real xmin = this->paradis_ptr->param->xBoundMin;
	// Real xmax = this->paradis_ptr->param->xBoundMax;
	// Real ymin = this->paradis_ptr->param->yBoundMin;
	// Real ymax = this->paradis_ptr->param->yBoundMax;
	// Real zmin = this->paradis_ptr->param->zBoundMin;
	// Real zmax = this->paradis_ptr->param->zBoundMax;
	Real xmin = minval[0];
	Real xmax = maxval[0];
	Real ymin = minval[1];
	Real ymax = maxval[1];
	Real zmin = minval[2];
	Real zmax = maxval[2];
	for (auto &&node : this->container_dd) {
	if ((xmin > node.position()[0]) \|\| (node.position()[0] > xmax) \|\|
	(ymin > node.position()[1]) \|\| (node.position()[1] > ymax) \|\|
	(zmin > node.position()[2]) \|\| (node.position()[2] > zmax)) {
	node.slaving();
	node.constrain();
	}
	}
	}
	}
	/* -------------------------------------------------------------------------- */

	void DomainParaDiS::performStep2() {
	// this->paradis_ptr->param->timeLMSCycle = 0.;
	// while (this->paradis_ptr->param->timeLMSCycle <
	// this->paradis_ptr->param->timeMax){
	// ParadisStep(this->paradis_ptr);
	// TimerClearAll(this->paradis_ptr);
	// }
	// this->paradis_ptr->cycleLMS++;
	// this->dd_container.incRelease();
	}

	void DomainParaDiS::performStep3() {}
	void DomainParaDiS::performStep4() {}
	void DomainParaDiS::performStep5() {}

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

	void DomainParaDiS::addSegments(std::vector<RefPointData<3>> &positions,
	std::vector<RefGenericElem<3>> &conn,
	std::vector<Real> &burgers,
	std::vector<Real> &normals) {

	UInt nb_points = positions.size();

	std::vector<Node_t *> created_nodes;
	std::vector<UInt> numNbrs;

	for (UInt i = 0; i < nb_points; ++i) {
	Node_t *newNode = GetNewNativeNode(this->paradis_ptr);

	auto &&pos = positions[i].position();
	newNode->x = pos[0];
	newNode->y = pos[1];
	newNode->z = pos[2];

	newNode->native = 1;
	newNode->constraint = 7;

	newNode->oldvX = 0.0;
	newNode->oldvY = 0.0;
	newNode->oldvZ = 0.0;

	AssignNodeToCell(this->paradis_ptr, newNode);

	created_nodes.push_back(newNode);
	}

	numNbrs.resize(created_nodes.size());

	UInt nb_elem = conn.size();
	for (UInt i = 0; i < nb_elem; ++i) {
	RefGenericElem<3> &el = conn[i];
	LM_ASSERT(el.conn.size() == 2,
	"we cannot add other things than segments to DD model");
	for (UInt n = 0; n < 2; ++n) {
	UInt nd = el.conn[n];
	++numNbrs[nd];
	}
	}

	for (UInt i = 0; i < nb_points; ++i) {
	Node_t *newNode = created_nodes[i];
	AllocNodeArms(newNode, numNbrs[i]);
	numNbrs[i] = 0;
	}

	for (UInt i = 0; i < nb_elem; ++i) {
	RefGenericElem<3> &el = conn[i];
	UInt nd_src = el.conn[0];
	UInt nd_dst = el.conn[1];

	Node_t *srcNode = created_nodes[nd_src];
	Node_t *dstNode = created_nodes[nd_dst];

	UInt &current_arm_src = numNbrs[nd_src];
	UInt &current_arm_dst = numNbrs[nd_dst];

	srcNode->nbrTag[current_arm_src].domainID = dstNode->myTag.domainID;
	srcNode->nbrTag[current_arm_src].index = dstNode->myTag.index;

	srcNode->burgX[current_arm_src] = burgers[3 * i + 0];
	srcNode->burgY[current_arm_src] = burgers[3 * i + 1];
	srcNode->burgZ[current_arm_src] = burgers[3 * i + 2];

	srcNode->nx[current_arm_src] = normals[3 * i + 0];
	srcNode->ny[current_arm_src] = normals[3 * i + 1];
	srcNode->nz[current_arm_src] = normals[3 * i + 2];

	dstNode->nbrTag[current_arm_dst].domainID = srcNode->myTag.domainID;
	dstNode->nbrTag[current_arm_dst].index = srcNode->myTag.index;

	dstNode->burgX[current_arm_dst] = -1. * burgers[3 * i + 0];
	dstNode->burgY[current_arm_dst] = -1. * burgers[3 * i + 1];
	dstNode->burgZ[current_arm_dst] = -1. * burgers[3 * i + 2];

	dstNode->nx[current_arm_dst] = normals[3 * i + 0];
	dstNode->ny[current_arm_dst] = normals[3 * i + 1];
	dstNode->nz[current_arm_dst] = normals[3 * i + 2];

	++current_arm_src;
	++current_arm_dst;
	}

	for (UInt i = 0; i < nb_points; ++i) {
	if (numNbrs[i] == 1) {
	created_nodes[i]->constraint = PINNED_NODE;
	// Add all pinned nodes to the list of (possible) hybrid nodes
	nh_l_conInt.push_back(created_nodes[i]);
	}
	}
	}

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

	void DomainParaDiS::imageForce(const double position, const double normal,
	double *force) {
	/*computes the image forces on a interface given by it's normal vector and
	position
	needed for basic CADD
	*/
	unsigned int n_dim = 3;
	double totStress[3][3];
	GetFieldPointStress(this->paradis_ptr, position[0], position[1], position[2],
	totStress);
	for (unsigned int i = 0; i < n_dim; ++i) {
	force[i] = 0.;
	for (unsigned int j = 0; j < n_dim; ++j) {
	force[i] += totStress[i][j] * normal[j];
	}
	}

	} // CHANGE THIS?
	/* -------------------------------------------------------------------------- */

	/* Configure segment connecting two nodes */

	void setNodeArm(Node_t nh_src, Node_t nh_dst, const UInt idxArmSrc,
	const UInt idxArmDst, const Real *burgersVec,
	const Real *gpNormVec) {
	/* Source node */

	nh_src->nbrTag[idxArmSrc].domainID = nh_dst->myTag.domainID;
	nh_src->nbrTag[idxArmSrc].index = nh_dst->myTag.index;

	nh_src->burgX[idxArmSrc] = burgersVec[0];
	nh_src->burgY[idxArmSrc] = burgersVec[1];
	nh_src->burgZ[idxArmSrc] = burgersVec[2];

	nh_src->nx[idxArmSrc] = gpNormVec[0];
	nh_src->ny[idxArmSrc] = gpNormVec[1];
	nh_src->nz[idxArmSrc] = gpNormVec[2];

	/* Distributive node */

	nh_dst->nbrTag[idxArmDst].domainID = nh_src->myTag.domainID;
	nh_dst->nbrTag[idxArmDst].index = nh_src->myTag.index;

	nh_dst->burgX[idxArmDst] = -burgersVec[0];
	nh_dst->burgY[idxArmDst] = -burgersVec[1];
	nh_dst->burgZ[idxArmDst] = -burgersVec[2];

	nh_dst->nx[idxArmDst] = gpNormVec[0];
	nh_dst->ny[idxArmDst] = gpNormVec[1];
	nh_dst->nz[idxArmDst] = gpNormVec[2];
	}
	/* -------------------------------------------------------------------------- */

	/* Orthogonal projection of a node onto the glide plane */

	inline void projectNodeOnGP(Node_t *nh_src, // source node
	Node_t *nh_dst, // inheriting node
	const UInt idxArmSrc) // arm index of source node
	// connecting to the
	// inheriting node
	{
	Real alpha = nh_src->nx[idxArmSrc] * (nh_src->x - nh_dst->x) +
	nh_src->ny[idxArmSrc] * (nh_src->y - nh_dst->y) +
	nh_src->nz[idxArmSrc] * (nh_src->z - nh_dst->z);

	nh_dst->x += alpha * nh_src->nx[idxArmSrc];
	nh_dst->y += alpha * nh_src->ny[idxArmSrc];
	nh_dst->z += alpha * nh_src->nz[idxArmSrc];
	}
	/* -------------------------------------------------------------------------- */

	/* Add hybrid segments to an existing DD line crossing an interface */

	// NOTE: test version --> the method assumes that the mesh struct carries only 1
	// dislocation!!
	// NOTE: it is further assumed that the node of hybrid DD line in the continuum
	// connecting to a node in the atomistic domain remains the nearest node
	// to the interface during a simulation!
	// NOTE: there is no particular check performed if "continuum" DD nodes
	// previously connected to a DD node in the atomistic domain retain the
	// this connectivity!!!
	// NOTE:

	void DomainParaDiS::addHybSegments(
	std::vector<RefPointData<3>> &positions, // DD nodes in an atomistic domain
	std::vector<RefGenericElem<3>>
	&conn, // connectivity between DD nodes in the atomistic domain
	Real threshold) // max. distance between two adjacent DD nodes possibly
	// generating a hybrid segment
	{
	UInt nb_points = positions.size();
	Real dist;
	Real pos_ato[3], pos_con[3];
	Real burgersVec[3], gpNormVec[3];
	// Containers/Lists
	std::vector<UInt> numNbrs;
	// ParaDis
	Node_t nh, nh_conInt1 = NULL, *nh_conInt2 = NULL;
	// ParaDis - Iterators
	std::vector<Node_t *>::iterator nh_it;
	/* ------------------------------------------------------------------------ */
	// Delete all existing DD segments in the atomistic domain
	// NOTE: the FreeNode() function doesn't remove the segments associated with
	// the node to delete
	// --> more precisely, they are not deallocated
	// --> this is of no particular issue since one can call
	// ReallocNodeArms() in step 3)

	for (nh_it = created_nodes.begin(); nh_it != created_nodes.end(); ++nh_it) {
	FreeNode(paradis_ptr, (*nh_it)->myTag.index);
	}

	created_nodes.clear();
	/* ------------------------------------------------------------------------ */

	/* 1) Import new nodes in ParaDis */

	for (UInt i = 0; i < nb_points; ++i) {
	nh = GetNewNativeNode(this->paradis_ptr);

	auto &&pos = positions[i].position();
	nh->x = pos[0];
	nh->y = pos[1];
	nh->z = pos[2];

	nh->native = 1;
	nh->constraint =
	PINNED_NODE; // nodes in an atomistic domain remain always fixed

	nh->oldvX = 0.0;
	nh->oldvY = 0.0;
	nh->oldvZ = 0.0;

	AssignNodeToCell(this->paradis_ptr, nh);

	created_nodes.push_back(nh);
	}

	numNbrs.resize(created_nodes.size());
	/* ------------------------------------------------------------------------ */

	/* 2) Find nearest neighbor(s) in the continuum */

	// 2.1) First node
	nh = created_nodes[0];

	pos_ato[0] = nh->x;
	pos_ato[1] = nh->y;
	pos_ato[2] = nh->z;

	for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
	++nh_it) {
	// Position of "continuum" node
	pos_con[0] = (*nh_it)->x;
	pos_con[1] = (*nh_it)->y;
	pos_con[2] = (*nh_it)->z;
	// Compute the distance between the continuum and the
	dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
	(pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
	(pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
	// If distance is smaller than a prescribed threshold
	// --> mark the node as a neighbor of the first "atomistic" DD node
	if (dist < threshold) {
	nh_conInt1 = *nh_it;
	numNbrs[0] = 1;
	break;
	}
	}

	// 2.2) Last node
	nh = created_nodes[nb_points - 1];

	pos_ato[0] = nh->x;
	pos_ato[1] = nh->y;
	pos_ato[2] = nh->z;

	for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
	++nh_it) {
	// Position of "continuum" node
	pos_con[0] = (*nh_it)->x;
	pos_con[1] = (*nh_it)->y;
	pos_con[2] = (*nh_it)->z;
	// Compute the distance between the continuum and the
	dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
	(pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
	(pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
	// If distance is smaller than a prescribed threshold
	// --> mark the node as a neighbor of the last "atomistic" DD node
	if (dist < threshold) {
	nh_conInt2 = *nh_it;
	numNbrs[nb_points - 1] = 1;
	break;
	}
	}
	/* ------------------------------------------------------------------------ */

	/* 3) Import connectivity and allocate memory for segment info */

	UInt nb_elem = conn.size();
	for (UInt i = 0; i < nb_elem; ++i) {
	RefGenericElem<3> &el = conn[i];
	LM_ASSERT(el.conn.size() == 2,
	"we cannot add other things than segments to DD model");
	for (UInt n = 0; n < 2; ++n) {
	UInt nd = el.conn[n];
	++numNbrs[nd];
	}
	}

	for (UInt i = 0; i < nb_points; ++i) {
	ReallocNodeArms(created_nodes[i], numNbrs[i]);
	// AllocNodeArms(newNode, numNbrs[i]);
	numNbrs[i] = 0;
	}
	/* ------------------------------------------------------------------------ */

	/* 4) Connect atomistic and continuum DD line */

	// 4.1) First interface node in the continuum connects to the atomistic DD
	// line
	if ((nh_conInt1 != NULL) && (nh_conInt2 == NULL)) {
	// Get the Burgers vector of the dislocation
	burgersVec[0] = nh_conInt1->burgX[0];
	burgersVec[1] = nh_conInt1->burgY[0];
	burgersVec[2] = nh_conInt1->burgZ[0];
	// Get the normal vector of the glide plane
	gpNormVec[0] = nh_conInt1->nx[0];
	gpNormVec[1] = nh_conInt1->ny[0];
	gpNormVec[2] = nh_conInt1->nz[0];
	// Re-allocate memory for the hybrid segment if necessary (--> should only
	// be done once before the first iteration)
	if (nh_conInt1->numNbrs < 2) {
	ReallocNodeArms(nh_conInt1, 2);
	nh_conInt1->constraint = UNCONSTRAINED;
	}
	// Establish connectivity between nodes in the continuum and the atomistic
	// domain ...
	setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
	// // ... and set homogeneous Dirichlet conditions on the other node
	// (--> u=0 <=> pinned node)
	// created_nodes[nb_points-1]->constraint = PINNED_NODE;
	// Project atomistic DD node on the glide plane
	projectNodeOnGP(nh_conInt1, created_nodes[0], 1);

	++numNbrs[0];
	}
	// 4.2) Last interface node in the continuum connects to the atomistic DD line
	else if ((nh_conInt1 == NULL) && (nh_conInt2 != NULL)) {
	// Get the Burgers vector of the dislocation
	burgersVec[0] = nh_conInt2->burgX[0];
	burgersVec[1] = nh_conInt2->burgY[0];
	burgersVec[2] = nh_conInt2->burgZ[0];
	// Get the normal vector of the glide plane
	gpNormVec[0] = nh_conInt2->nx[0];
	gpNormVec[1] = nh_conInt2->ny[0];
	gpNormVec[2] = nh_conInt2->nz[0];
	// Re-allocate memory for the hybrid segment if necessary (--> should only
	// be done once before the first iteration)
	if (nh_conInt2->numNbrs < 2) {
	ReallocNodeArms(nh_conInt2, 2);
	nh_conInt2->constraint = UNCONSTRAINED;
	}
	// Establish connectivity between nodes in the continuum and the atomistic
	// domain ...
	setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
	gpNormVec);
	// // ... and set homogeneous Dirichlet conditions on the other node
	// (--> u=0 <=> pinned node)
	// created_nodes[0]->constraint = PINNED_NODE;
	// Project atomistic DD node on the glide plane
	projectNodeOnGP(nh_conInt2, created_nodes[nb_points - 1], 1);

	++numNbrs[nb_points - 1];
	}
	// 4.3) Both interfaces nodes in the continuum connects to the atomistic DD
	// line
	else if ((nh_conInt1 != NULL) && (nh_conInt2 != NULL)) {
	// Get the Burgers vector of the dislocation
	burgersVec[0] = nh_conInt1->burgX[0];
	burgersVec[1] = nh_conInt1->burgY[0];
	burgersVec[2] = nh_conInt1->burgZ[0];
	// Get the normal vector of the glide plane
	gpNormVec[0] = nh_conInt1->nx[0];
	gpNormVec[1] = nh_conInt1->ny[0];
	gpNormVec[2] = nh_conInt1->nz[0];
	// Re-allocate memory for the hybrid segment if necessary (--> should only
	// be done once before the first iteration)
	if (nh_conInt1->numNbrs < 2) {
	ReallocNodeArms(nh_conInt1, 2);
	nh_conInt1->constraint = UNCONSTRAINED;
	}
	if (nh_conInt2->numNbrs < 2) {
	ReallocNodeArms(nh_conInt2, 2);
	nh_conInt2->constraint = UNCONSTRAINED;
	}
	// Establish connectivity between nodes in the continuum and the atomistic
	// domain
	setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
	setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
	gpNormVec);
	// Project atomistic DD nodes on the glide plane
	projectNodeOnGP(nh_conInt1, created_nodes[0], 1);
	projectNodeOnGP(nh_conInt1, created_nodes[nb_points - 1], 1);

	++numNbrs[0];
	++numNbrs[nb_points - 1];
	} else {
	LM_FATAL(std::endl
	<< "Hybrid dislocation failed to connect!!" << std::endl);
	}
	/* ------------------------------------------------------------------------ */

	/* 5) Now connect the remainder of the hybrid dislocation line in the
	* atomistic domain */

	for (UInt i = 0; i < nb_elem; ++i) {
	// Get connectivity
	RefGenericElem<3> &el = conn[i];
	UInt idxNhSrc = el.conn[0];
	UInt idxNhDst = el.conn[1];
	UInt &idxArmSrc = numNbrs[idxNhSrc];
	UInt &idxArmDst = numNbrs[idxNhDst];
	// Update arm information for each node of segment i
	setNodeArm(created_nodes[idxNhSrc], created_nodes[idxNhDst], idxArmSrc,
	idxArmDst, burgersVec, gpNormVec);
	// Project atomistic DD node on the glide plane
	projectNodeOnGP(created_nodes[idxNhSrc], created_nodes[idxNhDst],
	idxArmSrc);

	++idxArmSrc;
	++idxArmDst;
	}
	}

	/* -------------------------------------------------------------------------- */
	UInt DomainParaDiS::getNBNodes() {
	return this->container_dd.getContainerNodes().size();
	}
	/* -------------------------------------------------------------------------- */
	Real DomainParaDiS::getShearModulus() {
	return this->paradis_ptr->param->shearModulus;
	}
	/* -------------------------------------------------------------------------- */
	VectorProxy<3u> DomainParaDiS::getXmin() {
	return VectorProxy<3u>(this->paradis_ptr->param->xBoundMin,
	this->paradis_ptr->param->yBoundMin,
	this->paradis_ptr->param->zBoundMin);
	}
	/* -------------------------------------------------------------------------- */
	VectorProxy<3u> DomainParaDiS::getXmax() {
	return VectorProxy<3u>(this->paradis_ptr->param->xBoundMax,
	this->paradis_ptr->param->yBoundMax,
	this->paradis_ptr->param->zBoundMax);
	}
	/* -------------------------------------------------------------------------- */
	void DomainParaDiS::freenodearm(UInt index) {
	Node_t node = *paradis_ptr->nodeKeys[index];
	FreeNodeArms(&node);
	}
	/* -------------------------------------------------------------------------- */
	void DomainParaDiS::freenode(UInt index) {
	Node_t node = *paradis_ptr->nodeKeys[index];
	FreeNodeArms(&node);
	}

	/* -------------------------------------------------------------------------- */
	bool DomainParaDiS::is_it_bound() {
	// ContainerParaDiSNode::iterator dbgNodes =
	// this->getContainerNodes().begin();
	// UInt counter = 0;
	// for(RefParaDiSNode node=dbgNodes.getFirst(); !dbgNodes.end();
	// node=dbgNodes.getNext()){
	// ++counter;
	// }
	return this->boundary;
	} // CHANGE THIS?

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

	Cube DomainParaDiS::getDomainBoundingBox() {

	Cube c;
	c.getXmin() = VectorView<3>(this->paradis_ptr->param->minCoordinates);
	c.getXmax() = VectorView<3>(this->paradis_ptr->param->maxCoordinates);
	return c;
	}
	/* -------------------------------------------------------------------------- */

	Cube DomainParaDiS::getSubDomainBoundingBox() { LM_TOIMPLEMENT; }
	/* -------------------------------------------------------------------------- */

	void DomainParaDiS::setDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }

	/* -------------------------------------------------------------------------- */
	void DomainParaDiS::setSubDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }

	// void getDomainDimensions(Real xmin, Real xmax) {
	// };

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

	/* LMDESC PARADIS
	This section is used to configure the ParaDiS plugin
	*/

	/* LMEXAMPLE

	Section MultiScale AtomsUnits\\
	...\\
	GEOMETRY 1 CUBE BBOX -10 10 -10 10 -10 10
	MODEL PARADIS dislomodel \\
	...\\
	endSection\\ \\

	Section PARADIS:dislomodel RealUnits \\
	NUMDLBSTEPS 0 \\
	CONTROLFILE frank_read_src.ctrl \\
	DATAFILE frank_read_src.data \\
	TIMESTEP 1e-10 \\
	DOMAIN_GEOMETRY 1 \\
	endSection
	*/

	/* LMHERITATE domain_dd */

	void DomainParaDiS::declareParams() {

	DomainDD<ContainerNodesParaDiS, ContainerElemsParaDiS, 3>::declareParams();

	/* LMKEYWORD TIMESTEP
	Set the maximum time step used for ParaDiS
	*/
	this->parseKeyword("TIMESTEP", timestep);
	/* LMKEYWORD NUMDLBSTEPS
	Number of dynamic load balancing steps
	*/
	this->parseKeyword("NUMDLBSTEPS", numDLBCycles);
	/* LMKEYWORD CONTROLFILE
	Path to the ParaDiS control file
	*/
	this->parseKeyword("CONTROLFILE", control_filename);
	/* LMKEYWORD DATAFILE
	Path to the ParaDis data file
	*/
	this->parseKeyword("DATAFILE", data_filename);
	/* LMKEYWORD MOBILITYFILE
	Mobility input acquired from MD
	*/
	this->parseKeyword("MOBILITYFILE", mobility_filename, "");
	/* LMKEYWORD BOUNDARY
	Constraint boundary to mimic grain boundary
	*/
	this->parseTag("BOUNDARY", this->boundary, false);

	this->makeItOptional("DOMAIN_GEOMETRY");
	}

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

	__END_LIBMULTISCALE__

domain_paradis.ccNo OneTemporaryActions

File Metadata

domain_paradis.ccView Options

Event Timeline

domain_paradis.cc
No OneTemporary
Actions

domain_paradis.cc
View Options