Diffusion LibMultiScale (features/components) Browse

components

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(real8 A, real8 B, real8 C, real8 burg,
			real8 dispPointU, real8 dispPointX,
			real8 pois);
			void ComputeClosurePoint(real8 A, real8 B, real8 burg, real8 C);
			void ComputeProjectedPoint(real8 A, real8 C, real8 burg, real8 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(Real A, Real B, Real C, Real burg,
			Real dispPointU, Real dispPointX,
			Real pois) {
			SegmentDisplacementsOfAtom(A, B, C, burg, dispPointU, dispPointX, pois);
			}

			/* -------------------------------------------------------------------------- */
			void DomainPARADIS::computeClosurePoint(Real A, Real B, Real burg, Real C) {
			ComputeClosurePoint(A, B, burg, C);
			}
			/* -------------------------------------------------------------------------- */
			void DomainPARADIS::computeProjectedPoint(Real A, Real C, Real *burg,
			Real *P) {
			ComputeProjectedPoint(A, C, burg, P);
			}
			/* -------------------------------------------------------------------------- */

			void DomainPARADIS::init() {
			Communicator &comm = Communicator::getCommunicator();
			if (!comm.amIinGroup(this->getGroupID()))
			return;
			const UInt nbprocs = comm.getNBprocsOnGroup(this->getGroupID());

			// 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;

			const char *argv[8] = {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(), mobility_filename.c_str()};

			unsigned int argc = 8;

			/*
			* 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 (comm.groupRank(lm_my_proc_id, this->getGroupID()) == 0)
			std::cout << "BEGIN: Initialize ParaDiS on " << nbprocs << " processors."
			<< std::endl;
			Initialize(this->paradis_ptr, argc, const_cast<char **>(argv),
			comm.getMpiGroup(this->getGroupID()));

			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->getContainerNodes().setGrain(this->paradis_ptr);
			// this->getContainerDisplacements().setGrain(this->paradis_ptr);

			#ifndef LM_OPTIMIZED
			// Debug tests.
			unsigned int counter = 0;
			DUMP("Printing all initial ParaDiS nodes...", DBG_INFO);
			for (auto &&node : this->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;
			this->dd_container.setRelease(INITIAL_MODEL_RELEASE);
			}

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

			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->dd_container.incRelease();

			if (this->boundary) {
			Cube &c =
			GeometryManager::getManager().getGeometry(this->geom)->getBoundingBox();
			const Quantity<Length, 3> &maxval = c.getXmax();
			const Quantity<Length, 3> &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->getContainerNodes()) {
			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->getContainerNodes().size(); }
			/* -------------------------------------------------------------------------- */
			Real DomainPARADIS::getShearModulus() {
			return this->paradis_ptr->param->shearModulus;
			}
			/* -------------------------------------------------------------------------- */
			Real DomainPARADIS::getXmin(UInt dim) {
			if (dim == 0)
			return this->paradis_ptr->param->xBoundMin;
			else if (dim == 1)
			return this->paradis_ptr->param->yBoundMin;
			else if (dim == 2)
			return this->paradis_ptr->param->zBoundMin;
			else
			LM_FATAL("Unknown dimension");
			}
			/* -------------------------------------------------------------------------- */
			Real DomainPARADIS::getXmax(UInt dim) {
			if (dim == 0)
			return this->paradis_ptr->param->xBoundMax;
			else if (dim == 1)
			return this->paradis_ptr->param->yBoundMax;
			else if (dim == 2)
			return this->paradis_ptr->param->zBoundMax;
			else
			LM_FATAL("Unknown dimension");
			}
			/* -------------------------------------------------------------------------- */
			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?
			/* -------------------------------------------------------------------------- */

			/* 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<ContainerPARADISNode, DummySource, DummyObstacle,
			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);
			}

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

			__END_LIBMULTISCALE__

libmultiscale/src/dd/paradis/domain_paradis.cc4a06aa073a0bfeatures/components

domain_paradis.ccView Options

libmultiscale/src/dd/paradis/domain_paradis.cc
4a06aa073a0bfeatures/components

domain_paradis.cc
View Options