domain_paradis.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, May 27, 08:56

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 "field.hh"
	#include "geometry_manager.hh"
	/* -------------------------------------------------------------------------- */

	// ParaDiS includes
	extern "C" {
	// void TimerClearAll(Home_t *home);
	#include "Decomp.h"
	#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 GetFieldPointStress(Home_t *home, real8 x, real8 y, real8 z,
	real8 totStress[3][3]);

	void RecycleGhostNodes(Home_t *home);
	void CommSendGhosts(Home_t *home);
	void CommSendVelocity(Home_t *home);
	void TrapezoidTimestepIntegrator(Home_t *home);
	void ApplyDeltaStress(Home_t *home, real8 deltaStress[3][3]);
	void CommSendRemesh(Home_t *home);
	#ifdef FPES_ON
	#include <fpcontrol.h>
	#endif
	}

	#undef PARALLEL

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

	__BEGIN_LIBMULTISCALE__

	DomainParaDiS::DomainParaDiS(LMID ID, CommGroup &group)
	: LMObject(ID), DomainDD<ContainerNodesParaDiS, ContainerElemsParaDiS, 3>(
	ID, group) {}
	/* -------------------------------------------------------------------------- */

	DomainParaDiS::~DomainParaDiS() {}

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

	void DomainParaDiS::remesh() {
	RemeshAfterCoupling(ParaDiSHome::paradis_ptr);

	// //ParadisStep.c from line 336 to 418
	// InitTopologyExemptions(ParaDiSHome::paradis_ptr);
	// ClearOpList(ParaDiSHome::paradis_ptr);
	// SortNodesForCollision(ParaDiSHome::paradis_ptr);
	// SplitMultiNodes(ParaDiSHome::paradis_ptr);

	// CrossSlip(ParaDiSHome::paradis_ptr);
	// HandleCollisions(ParaDiSHome::paradis_ptr);
	// ProximityCollisions(ParaDiSHome::paradis_ptr);

	// #ifdef PARALLEL
	// #ifdef PARADIS_IN_LIBMULTISCALE
	// CommSendRemesh(ParaDiSHome::paradis_ptr);
	// #endif
	// #endif

	// FixRemesh(ParaDiSHome::paradis_ptr);

	// Remesh(ParaDiSHome::paradis_ptr);
	// Migrate(ParaDiSHome::paradis_ptr);

	// #ifdef PARALLEL
	// #ifdef PARADIS_IN_LIBMULTISCALE
	// RecycleGhostNodes(ParaDiSHome::paradis_ptr);
	// #endif
	// #endif

	// SortNativeNodes(ParaDiSHome::paradis_ptr);

	// #ifdef PARALLEL
	// #ifdef PARADIS_IN_LIBMULTISCALE
	// CommSendGhosts(ParaDiSHome::paradis_ptr);
	// #endif
	// #endif
	}

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

	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

	ParaDiSHome::paradis_ptr = InitHome();
	TimerInit(ParaDiSHome::paradis_ptr);
	TimerStart(ParaDiSHome::paradis_ptr, TOTAL_TIME);
	TimerStart(ParaDiSHome::paradis_ptr, INITIALIZE);

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

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

	if (nbprocs > 1)
	MPI_Barrier(ParaDiSHome::paradis_ptr->MPI_COMM_PARADIS);

	ParaDiSHome::paradis_ptr->param->timeMax = this->timestep;
	ParaDiSHome::paradis_ptr->cycle = ParaDiSHome::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 = ParaDiSHome::paradis_ptr->param->numDLBCycles;
	TimerStart(ParaDiSHome::paradis_ptr, INITIALIZE);
	if ((ParaDiSHome::paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
	time(&tp);
	DUMP(" +++ Beginning " << initialDLBCycles << " load-balancing steps at "
	<< asctime(localtime(&tp)),
	DBG_INFO_STARTUP);
	}
	while (ParaDiSHome::paradis_ptr->param->numDLBCycles > 0) {
	ParadisStep(ParaDiSHome::paradis_ptr);
	ParaDiSHome::paradis_ptr->cycle++;
	ParaDiSHome::paradis_ptr->param->numDLBCycles--;
	}
	if ((ParaDiSHome::paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
	time(&tp);
	DUMP(" +++ Completed load-balancing steps at " << asctime(localtime(&tp)),
	DBG_INFO_STARTUP);
	}
	TimerStop(ParaDiSHome::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(ParaDiSHome::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.
	*/
	ParaDiSHome::paradis_ptr->cycle = ParaDiSHome::paradis_ptr->param->cycleStart;
	ParaDiSHome::paradis_ptr->cycleLMS = 1;

	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

	if (ParaDiSHome::paradis_ptr->myDomain == 0)
	DUMP("END: Initialize ParaDiS on " << nbprocs << " processors.",
	DBG_INFO_STARTUP);

	if (ParaDiSHome::paradis_ptr->param->maxDT != this->timestep)
	DUMP("Careful the timestep set by paradis is not the same "
	"as the one described by libmultiscale: "
	<< this->timestep
	<< " != " << ParaDiSHome::paradis_ptr->param->maxDT,
	DBG_WARNING);

	ParaDiSHome::paradis_ptr->param->maxDT = this->timestep;
	ParaDiSHome::paradis_ptr->param->deltaTT = this->timestep;
	}

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

	void DomainParaDiS::updatePosition() {

	// This is the time during the internal time steping loop
	auto &timeLMSCycle = ParaDiSHome::paradis_ptr->param->timeLMSCycle;
	timeLMSCycle = 0.;
	// This is the maximal time jump
	auto &timeMax = ParaDiSHome::paradis_ptr->param->timeMax;

	/*
	Because of collisions, ParaDiS tends to decrease the timestep.
	At the provided frequency, LibMultiScale will attempt to reset the
	timestep to the desired value (specified by the keyword TIMESTEP)
	*/
	if (current_step % timestep_reset_frequency == 0) {
	ParaDiSHome::paradis_ptr->param->nextDT = this->timestep;
	}

	// Perform leap steps until timeMax is reached
	while (timeLMSCycle < timeMax) {
	STARTTIMER("ParadisStep");
	ParadisStep(ParaDiSHome::paradis_ptr);
	TimerClearAll(ParaDiSHome::paradis_ptr);
	STOPTIMER("ParadisStep");
	}

	// the global cycle counter is managed by LM: increment it
	ParaDiSHome::paradis_ptr->cycleLMS++;
	this->changeRelease();
	}
	/* -------------------------------------------------------------------------- */

	void DomainParaDiS::updateGradient() {
	auto home = ParaDiSHome::paradis_ptr;
	// Calculate the net charge tensor for each cell (includes global comm)
	CellCharge(home);
	// compute the forces on all nodes
	NodeForce(home, FULL);
	}
	/* -------------------------------------------------------------------------- */

	void DomainParaDiS::enforceCompatibility() {

	this->updateVelocity();

	// copy back the positions to the nodal representation
	auto &nodes = this->getContainer().getContainerNodes();

	for (auto &&[pos, nd] : zip(positions.rowwise(), nodes)) {
	nd.position() = pos;
	}

	auto home = ParaDiSHome::paradis_ptr;

	/* In some simulations, it is necessary to recalculate and distribute
	* the glide plane infromation for segments after the nodes have been
	* repositioned. Do so now if needed.
	*/
	ResetGlidePlanes(home);

	/* Increment the per-burgers vector density gain/loss with
	* changes for this cycle. This must be done immediately
	* after timestep integration!
	*
	* Note: This is currently only applicable to BCC simulations.
	*/
	GetDensityDelta(home);

	// Calculate the new plastic strain.
	DeltaPlasticStrain(home);

	/* The call to GenerateOutput will update the time and cycle counters,
	* determine if any output needs to be generated at this stage, and
	* call the appropriate I/O functions if necessary.
	*/
	GenerateOutput(home, STAGE_CYCLE);

	/* Before doing topological changes, set flags indicating any
	* nodes exempt from topological changes. These flags are used
	* in both splitting multi-arm nodes and collisions, so this
	* function should be invoked before either of those items are done.
	*/
	InitTopologyExemptions(home);

	/* Now do all the topological changes from segment interactions
	* (collisions, multinode splitting)... Clear the list of local
	* operations that will be sent to the remote domains for processsing,
	* then split any multi-arm nodes that need splitting, cross slip
	* nodes (as needed/allowed), handle all local collisions, then
	* send remote nodes the list of ops needed to keep their data in sync.
	*/
	ClearOpList(home);
	SortNodesForCollision(home);

	SplitMultiNodes(home);

	/* Call a generic cross-slip dispatch function that will call
	* (if necessary) the cross-slip function appropriate to the
	* type of material in use.
	*/
	CrossSlip(home);

	/* Search for dislocation segments in close proximity to each other
	* and if necessary handle any collision between them.
	*/
	HandleCollisions(home);

	CommSendRemesh(home);
	FixRemesh(home);

	/* Under certain circumstances, parallel topological changes can
	* create double links between nodes; links which can not be detected
	* until after FixRemesh() is called... so, a quick check has to be
	* done to clean up these potential double-links here, or they will
	* cause problems later on. Should only have to check nodes local
	* to this domain.
	*/
	Node_t *node;
	for (int i = 0; i < home->newNodeKeyPtr; i++) {
	if ((node = home->nodeKeys[i]) == (Node_t *)NULL)
	continue;
	(void)RemoveDoubleLinks(home, node, 0);
	node->flags &= ~NODE_CHK_DBL_LINK;
	}

	// Invoke mesh coarsen/refine
	Remesh(home);

	Real deltaStress[3][3];

	// Define load curve and calculate change in applied stress this cycle
	LoadCurve(home, deltaStress);

	/* This is only needed when we do force calcs for only
	* a subset of the nodes at the beginning of the timestep. It will
	* adjust the nodal forces based on the current delta stress and
	* recalculate the nodal velocities so we have more accurate values
	* when we enter the timestep integrator at the beginning of the next
	* cycle.
	*/
	ApplyDeltaStress(home, deltaStress);

	/* If necessary, use the current load data to generate a new
	* domain decomposition to rebalance the workload among the
	* processors.
	*/
	Rebalance(home, DLB_USE_WALLCLK_TIME);

	/* Send any nodes that have moved beyond the domain's
	* boundaries to the domain the node now belongs to.
	*/
	Migrate(home);

	// Recycle all the ghost nodes: move them back to the free Queue
	RecycleGhostNodes(home);

	// Sort the native nodes into their proper subcells.
	SortNativeNodes(home);

	// Communicate ghost cells to/from neighbor domains
	CommSendGhosts(home);

	#ifdef NAN_CHECK
	/* For debug only: Abort if any of the nodes have position or
	* velocity values that are NaNs or infinites. Be sure to do this
	* before we write the restart files so we don't get bad data
	* in the restart.
	*/
	CheckForNANS(home);
	#endif

	/* If memory debugging is enabled, run a consistency check on all
	* allocated memory blocks looking for corruption in any of the
	* block headers or trailers.
	*/
	#ifdef DEBUG_MEM
	ParadisMemCheck();
	#endif

	CheckMemUsage(home, "ParadisStep-complete");

	/* Zero out the count of force calculations done this cycle
	* so the load-balancing in the next step is based on accurate
	* values.
	*/
	home->cycleForceCalcCount = 0;

	// For Debugging and testing only...
	#if DEBUG_CHECK_FOR_ZERO_SEG
	CheckForEmptySimulation(home);
	#endif
	}
	/* ------------------------------------------------------------------------- */

	void DomainParaDiS::updateAcceleration() {
	auto &nodes = this->getContainer().getContainerNodes();
	if (nodes.size() == accelerations.rows())
	return;
	accelerations.resize(nodes.size(), 3);
	accelerations.setZero();
	}

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

	void DomainParaDiS::updateVelocity() {
	// it is wrongly named, as it updates velocity
	bool doAll = true;
	auto home = ParaDiSHome::paradis_ptr;
	auto mobIterError = CalcNodeVelocities(home, 0, doAll);
	if (mobIterError != 0) {
	DUMP("velocity could not be computed: " << mobIterError, DBG_WARNING);
	}
	CommSendVelocity(home);
	}

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

	ArrayView DomainParaDiS::primal() {

	auto &nodes = this->getContainer().getContainerNodes();

	positions.resize(nodes.size(), 3);
	for (auto &&[pos, nd] : zip(positions.rowwise(), nodes)) {
	pos = Vector<3>(nd.position());
	}
	return ArrayView(positions.data(), positions.rows(), positions.cols());
	}

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

	ArrayView DomainParaDiS::primalTimeDerivative() {
	this->updateVelocity();
	auto &nodes = this->getContainer().getContainerNodes();
	velocities.resize(nodes.size(), 3);

	for (auto &&[vel, nd] : zip(velocities.rowwise(), nodes)) {
	vel = Vector<3>(nd.velocity());
	}
	return ArrayView(velocities.data(), velocities.rows(), velocities.cols());
	}
	/* ------------------------------------------------------------------------ */
	ArrayView DomainParaDiS::acceleration() {
	this->updateAcceleration();
	return ArrayView(accelerations.data(), accelerations.rows(),
	accelerations.cols());
	}
	/* ------------------------------------------------------------------------ */

	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(ParaDiSHome::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(ParaDiSHome::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(ParaDiSHome::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(ParaDiSHome::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 ParaDiSHome::paradis_ptr->param->shearModulus;
	}
	/* ------------------------------------------------------------------------ */
	VectorProxy<3u> DomainParaDiS::getXmin() {
	return VectorProxy<3u>(
	std::forward_as_tuple(ParaDiSHome::paradis_ptr->param->xBoundMin,
	ParaDiSHome::paradis_ptr->param->yBoundMin,
	ParaDiSHome::paradis_ptr->param->zBoundMin));
	}
	/* ------------------------------------------------------------------------ */
	VectorProxy<3u> DomainParaDiS::getXmax() {
	return VectorProxy<3u>(
	std::forward_as_tuple(ParaDiSHome::paradis_ptr->param->xBoundMax,
	ParaDiSHome::paradis_ptr->param->yBoundMax,
	ParaDiSHome::paradis_ptr->param->zBoundMax));
	}
	/* ------------------------------------------------------------------------ */
	void DomainParaDiS::freenodearm(UInt index) {
	Node_t node = *ParaDiSHome::paradis_ptr->nodeKeys[index];
	FreeNodeArms(&node);
	}
	/* ------------------------------------------------------------------------ */
	void DomainParaDiS::freenode(UInt index) {
	Node_t node = *ParaDiSHome::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() const {

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

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

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

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

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

	/* 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 TIMESTEP_RESET_FREQ
	Because of collisions, ParaDiS tends to decrease the timestep.
	At the provided frequency, LibMultiScale will attempt to reset the
	timestep to the desired value (specified by the keyword TIMESTEP)
	*/
	this->parseKeyword("TIMESTEP_RESET_FREQ", this->timestep_reset_frequency,
	100u);
	/* 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

	Deprecated option: use a stimulator instead.
	*/
	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