Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/dd/paradis/domain_paradis.cc
1d1b2315bf45master

Last Change
Raw File

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 "domain_paradis.hh"
#include "field.hh"
#include "geometry_manager.hh"
#include "lm_communicator.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 RealUnits\\
...\\
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__
c4science · Help