Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/dd/paradis/domain_paradis.cc
4a06aa073a0blammps_2020

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 "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__
c4science · Help