Diffusion LibMultiScale (features/new_parralelism) Browse

new_parralelism

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__

libmultiscale/src/dd/paradis/domain_paradis.cc6e0d714926eefeatures/new_parralelism

domain_paradis.ccView Options

libmultiscale/src/dd/paradis/domain_paradis.cc
6e0d714926eefeatures/new_parralelism

domain_paradis.cc
View Options