Page MenuHomec4science
Files F64662505

domain_paradis.cc
No OneTemporary
Actions

File Metadata

Created
Tue, May 28, 12:59

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(libmultiscale::Vector<3u> &A,
libmultiscale::Vector<3u> &B,
libmultiscale::Vector<3u> &C,
libmultiscale::Vector<3u> &burg,
libmultiscale::Vector<3u> &dispPointU,
libmultiscale::Vector<3u> &dispPointX,
libmultiscale::Real pois);
void ComputeClosurePoint(libmultiscale::Vector<3u> &A,
libmultiscale::Vector<3u> &B,
libmultiscale::Vector<3u> &burg,
libmultiscale::Vector<3u> &C);
void ComputeProjectedPoint(libmultiscale::Vector<3u> &A,
libmultiscale::Vector<3u> &C,
libmultiscale::Vector<3u> &burg,
libmultiscale::Vector<3u> &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(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::computeDisplacements(std::vector<Real> &positions,
// std::vector<Real> &displacements) {
// if (positions.empty()) {
// displacements.clear();
// return;
// }
// if (!ParaDiSHome::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);
// ParaDiSHome::paradis_ptr->dispPointCount = size / 3;
// ParaDiSHome::paradis_ptr->dispPointX = &positions[0];
// ParaDiSHome::paradis_ptr->dispPointU = &displacements[0];
// // Calculate the displacements.
// PointDisplacements(ParaDiSHome::paradis_ptr);
// }
/* -------------------------------------------------------------------------- */
// void DomainParaDiS::computeBarnett(Vector<3u> &A, Vector<3u> &B, Vector<3u>
// &C,
// Vector<3u> &burg, Vector<3u> &dispPointU,
// Vector<3u> &dispPointX, Real pois) {
// SegmentDisplacementsOfAtom(A, B, C, burg, dispPointU, dispPointX, pois);
// }
// /* --------------------------------------------------------------------------
// */
// void DomainParaDiS::computeClosurePoint(Vector<3u> &A, Vector<3u> &B,
// Vector<3u> &burg, Vector<3u> &C) {
// ComputeClosurePoint(A, B, burg, C);
// }
// /* --------------------------------------------------------------------------
// */
// void DomainParaDiS::computeProjectedPoint(Vector<3u> &A, Vector<3u> &C,
// Vector<3u> &burg, Vector<3u> &P) {
// ComputeProjectedPoint(A, C, burg, P);
// }
// /* --------------------------------------------------------------------------
// */
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 (strcmp(ParaDiSHome::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 = "
<< 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);
printf(" +++ Beginning %d load-balancing steps at %s", initialDLBCycles,
asctime(localtime(&tp)));
}
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);
printf(" +++ Completed load-balancing steps at %s",
asctime(localtime(&tp)));
}
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
// ContainerParaDiSDisplacement displs = this->getContainerDisplacements();
// DUMP("Number of ParaDiS displacement fieldpoints: "<<displs.size(),
// DBG_INFO);
// End debug tests.
if (ParaDiSHome::paradis_ptr->myDomain == 0)
std::cout << "END: Initialize ParaDiS on " << nbprocs << " processors."
<< std::endl;
}
/* -------------------------------------------------------------------------- */
void DomainParaDiS::performStep1() {
throw;
ParaDiSHome::paradis_ptr->param->timeLMSCycle = 0.;
// if (current_step % timestep_reset_frequency == 0) {
// ParaDiSHome::paradis_ptr->param->nextDT = this->timestep;
// }
// while (ParaDiSHome::paradis_ptr->param->timeLMSCycle <
// ParaDiSHome::paradis_ptr->param->timeMax) {
// STARTTIMER("ParadisStep");
// ParadisStep(ParaDiSHome::paradis_ptr);
// TimerClearAll(ParaDiSHome::paradis_ptr);
// STOPTIMER("ParadisStep");
// }
// ParaDiSHome::paradis_ptr->cycleLMS++;
// this->getOutput("segments")->incRelease();
// if (this->boundary) {
// auto c =
// GeometryManager::getManager().getGeometry(this->geom)->getBoundingBox();
// auto &&maxval = c.getXmax();
// auto &&minval = c.getXmin();
// // Geometry* dd =
// GeometryManager::getManager().getGeometry(this->geom);
// // Cube & dd_cube = dynamic_cast<Cube&>(*dd);
// // to mimic strong grain boundary...
// // Real xmin = ParaDiSHome::paradis_ptr->param->xBoundMin + 10.0;
// // Real xmax = ParaDiSHome::paradis_ptr->param->xBoundMax - 10.0;
// // Real ymin = ParaDiSHome::paradis_ptr->param->yBoundMin + 10.0;
// // Real ymax = ParaDiSHome::paradis_ptr->param->yBoundMax - 10.0;
// // Real zmin = ParaDiSHome::paradis_ptr->param->zBoundMin + 10.0;
// // Real zmax = ParaDiSHome::paradis_ptr->param->zBoundMax - 10.0;
// // Real xmin = ParaDiSHome::paradis_ptr->param->xBoundMin;
// // Real xmax = ParaDiSHome::paradis_ptr->param->xBoundMax;
// // Real ymin = ParaDiSHome::paradis_ptr->param->yBoundMin;
// // Real ymax = ParaDiSHome::paradis_ptr->param->yBoundMax;
// // Real zmin = ParaDiSHome::paradis_ptr->param->zBoundMin;
// // Real zmax = ParaDiSHome::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->container_dd) {
// 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::updateGradient() {
// // ParaDiSHome::paradis_ptr->param->timeLMSCycle = 0.;
// // while (ParaDiSHome::paradis_ptr->param->timeLMSCycle <
// // ParaDiSHome::paradis_ptr->param->timeMax){
// // ParadisStep(ParaDiSHome::paradis_ptr);
// // TimerClearAll(ParaDiSHome::paradis_ptr);
// // }
// // ParaDiSHome::paradis_ptr->cycleLMS++;
// // this->dd_container.incRelease();
}
// void DomainParaDiS::performStep3() {}
// void DomainParaDiS::performStep4() {}
// void DomainParaDiS::performStep5() {}
/* -------------------------------------------------------------------------- */
void DomainParaDiS::imageForce(const double *position, const double *normal,
double *force) {
/*computes the image forces on a interface given by it's normal vector and
position
needed for basic CADD
*/
unsigned int n_dim = 3;
double totStress[3][3];
GetFieldPointStress(ParaDiSHome::paradis_ptr, position[0], position[1],
position[2], totStress);
for (unsigned int i = 0; i < n_dim; ++i) {
force[i] = 0.;
for (unsigned int j = 0; j < n_dim; ++j) {
force[i] += totStress[i][j] * normal[j];
}
}
} // CHANGE THIS?
/* -------------------------------------------------------------------------- */
/* Configure segment connecting two nodes */
void setNodeArm(Node_t *nh_src, Node_t *nh_dst, const UInt idxArmSrc,
const UInt idxArmDst, const Real *burgersVec,
const Real *gpNormVec) {
/* Source node */
nh_src->nbrTag[idxArmSrc].domainID = nh_dst->myTag.domainID;
nh_src->nbrTag[idxArmSrc].index = nh_dst->myTag.index;
nh_src->burgX[idxArmSrc] = burgersVec[0];
nh_src->burgY[idxArmSrc] = burgersVec[1];
nh_src->burgZ[idxArmSrc] = burgersVec[2];
nh_src->nx[idxArmSrc] = gpNormVec[0];
nh_src->ny[idxArmSrc] = gpNormVec[1];
nh_src->nz[idxArmSrc] = gpNormVec[2];
/* Distributive node */
nh_dst->nbrTag[idxArmDst].domainID = nh_src->myTag.domainID;
nh_dst->nbrTag[idxArmDst].index = nh_src->myTag.index;
nh_dst->burgX[idxArmDst] = -burgersVec[0];
nh_dst->burgY[idxArmDst] = -burgersVec[1];
nh_dst->burgZ[idxArmDst] = -burgersVec[2];
nh_dst->nx[idxArmDst] = gpNormVec[0];
nh_dst->ny[idxArmDst] = gpNormVec[1];
nh_dst->nz[idxArmDst] = gpNormVec[2];
}
/* -------------------------------------------------------------------------- */
/* Orthogonal projection of a node onto the glide plane */
inline void projectNodeOnGP(Node_t *nh_src, // source node
Node_t *nh_dst, // inheriting node
const UInt idxArmSrc) // arm index of source node
// connecting to the
// inheriting node
{
Real alpha = nh_src->nx[idxArmSrc] * (nh_src->x - nh_dst->x) +
nh_src->ny[idxArmSrc] * (nh_src->y - nh_dst->y) +
nh_src->nz[idxArmSrc] * (nh_src->z - nh_dst->z);
nh_dst->x += alpha * nh_src->nx[idxArmSrc];
nh_dst->y += alpha * nh_src->ny[idxArmSrc];
nh_dst->z += alpha * nh_src->nz[idxArmSrc];
}
/* -------------------------------------------------------------------------- */
/* Add hybrid segments to an existing DD line crossing an interface */
// NOTE: test version --> the method assumes that the mesh struct carries only 1
// dislocation!!
// NOTE: it is further assumed that the node of hybrid DD line in the continuum
// connecting to a node in the atomistic domain remains the nearest node
// to the interface during a simulation!
// NOTE: there is no particular check performed if "continuum" DD nodes
// previously connected to a DD node in the atomistic domain retain the
// this connectivity!!!
// NOTE:
// void DomainParaDiS::addHybSegments(
// std::vector<RefPointData<3>> &positions, // DD nodes in an atomistic
// domain
// std::vector<RefGenericElem<3>>
// &conn, // connectivity between DD nodes in the atomistic domain
// Real threshold) // max. distance between two adjacent DD nodes possibly
// // generating a hybrid segment
// {
// UInt nb_points = positions.size();
// Real dist;
// Real pos_ato[3], pos_con[3];
// Real burgersVec[3], gpNormVec[3];
// // Containers/Lists
// std::vector<UInt> numNbrs;
// // ParaDis
// Node_t *nh, *nh_conInt1 = NULL, *nh_conInt2 = NULL;
// // ParaDis - Iterators
// std::vector<Node_t *>::iterator nh_it;
// /* ------------------------------------------------------------------------
// */
// // Delete all existing DD segments in the atomistic domain
// // NOTE: the FreeNode() function doesn't remove the segments associated
// with
// // the node to delete
// // --> more precisely, they are not deallocated
// // --> this is of no particular issue since one can call
// // ReallocNodeArms() in step 3)
// for (nh_it = created_nodes.begin(); nh_it != created_nodes.end(); ++nh_it)
// {
// FreeNode(ParaDiSHome::paradis_ptr, (*nh_it)->myTag.index);
// }
// created_nodes.clear();
// /* ------------------------------------------------------------------------
// */
// /* 1) Import new nodes in ParaDis */
// for (UInt i = 0; i < nb_points; ++i) {
// nh = GetNewNativeNode(ParaDiSHome::paradis_ptr);
// auto &&pos = positions[i].position();
// nh->x = pos[0];
// nh->y = pos[1];
// nh->z = pos[2];
// nh->native = 1;
// nh->constraint =
// PINNED_NODE; // nodes in an atomistic domain remain always fixed
// nh->oldvX = 0.0;
// nh->oldvY = 0.0;
// nh->oldvZ = 0.0;
// AssignNodeToCell(ParaDiSHome::paradis_ptr, nh);
// created_nodes.push_back(nh);
// }
// numNbrs.resize(created_nodes.size());
// /* ------------------------------------------------------------------------
// */
// /* 2) Find nearest neighbor(s) in the continuum */
// // 2.1) First node
// nh = created_nodes[0];
// pos_ato[0] = nh->x;
// pos_ato[1] = nh->y;
// pos_ato[2] = nh->z;
// for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
// ++nh_it) {
// // Position of "continuum" node
// pos_con[0] = (*nh_it)->x;
// pos_con[1] = (*nh_it)->y;
// pos_con[2] = (*nh_it)->z;
// // Compute the distance between the continuum and the
// dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
// (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
// (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
// // If distance is smaller than a prescribed threshold
// // --> mark the node as a neighbor of the first "atomistic" DD node
// if (dist < threshold) {
// nh_conInt1 = *nh_it;
// numNbrs[0] = 1;
// break;
// }
// }
// // 2.2) Last node
// nh = created_nodes[nb_points - 1];
// pos_ato[0] = nh->x;
// pos_ato[1] = nh->y;
// pos_ato[2] = nh->z;
// for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
// ++nh_it) {
// // Position of "continuum" node
// pos_con[0] = (*nh_it)->x;
// pos_con[1] = (*nh_it)->y;
// pos_con[2] = (*nh_it)->z;
// // Compute the distance between the continuum and the
// dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
// (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
// (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
// // If distance is smaller than a prescribed threshold
// // --> mark the node as a neighbor of the last "atomistic" DD node
// if (dist < threshold) {
// nh_conInt2 = *nh_it;
// numNbrs[nb_points - 1] = 1;
// break;
// }
// }
// /* ------------------------------------------------------------------------
// */
// /* 3) Import connectivity and allocate memory for segment info */
// UInt nb_elem = conn.size();
// for (UInt i = 0; i < nb_elem; ++i) {
// RefGenericElem<3> &el = conn[i];
// LM_ASSERT(el.conn.size() == 2,
// "we cannot add other things than segments to DD model");
// for (UInt n = 0; n < 2; ++n) {
// UInt nd = el.conn[n];
// ++numNbrs[nd];
// }
// }
// for (UInt i = 0; i < nb_points; ++i) {
// ReallocNodeArms(created_nodes[i], numNbrs[i]);
// // AllocNodeArms(newNode, numNbrs[i]);
// numNbrs[i] = 0;
// }
// /* ------------------------------------------------------------------------
// */
// /* 4) Connect atomistic and continuum DD line */
// // 4.1) First interface node in the continuum connects to the atomistic DD
// // line
// if ((nh_conInt1 != NULL) && (nh_conInt2 == NULL)) {
// // Get the Burgers vector of the dislocation
// burgersVec[0] = nh_conInt1->burgX[0];
// burgersVec[1] = nh_conInt1->burgY[0];
// burgersVec[2] = nh_conInt1->burgZ[0];
// // Get the normal vector of the glide plane
// gpNormVec[0] = nh_conInt1->nx[0];
// gpNormVec[1] = nh_conInt1->ny[0];
// gpNormVec[2] = nh_conInt1->nz[0];
// // Re-allocate memory for the hybrid segment if necessary (--> should
// only
// // be done once before the first iteration)
// if (nh_conInt1->numNbrs < 2) {
// ReallocNodeArms(nh_conInt1, 2);
// nh_conInt1->constraint = UNCONSTRAINED;
// }
// // Establish connectivity between nodes in the continuum and the
// atomistic
// // domain ...
// setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
// // // ... and set homogeneous Dirichlet conditions on the other node
// // (--> u=0 <=> pinned node)
// // created_nodes[nb_points-1]->constraint = PINNED_NODE;
// // Project atomistic DD node on the glide plane
// projectNodeOnGP(nh_conInt1, created_nodes[0], 1);
// ++numNbrs[0];
// }
// // 4.2) Last interface node in the continuum connects to the atomistic DD
// line
// else if ((nh_conInt1 == NULL) && (nh_conInt2 != NULL)) {
// // Get the Burgers vector of the dislocation
// burgersVec[0] = nh_conInt2->burgX[0];
// burgersVec[1] = nh_conInt2->burgY[0];
// burgersVec[2] = nh_conInt2->burgZ[0];
// // Get the normal vector of the glide plane
// gpNormVec[0] = nh_conInt2->nx[0];
// gpNormVec[1] = nh_conInt2->ny[0];
// gpNormVec[2] = nh_conInt2->nz[0];
// // Re-allocate memory for the hybrid segment if necessary (--> should
// only
// // be done once before the first iteration)
// if (nh_conInt2->numNbrs < 2) {
// ReallocNodeArms(nh_conInt2, 2);
// nh_conInt2->constraint = UNCONSTRAINED;
// }
// // Establish connectivity between nodes in the continuum and the
// atomistic
// // domain ...
// setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
// gpNormVec);
// // // ... and set homogeneous Dirichlet conditions on the other node
// // (--> u=0 <=> pinned node)
// // created_nodes[0]->constraint = PINNED_NODE;
// // Project atomistic DD node on the glide plane
// projectNodeOnGP(nh_conInt2, created_nodes[nb_points - 1], 1);
// ++numNbrs[nb_points - 1];
// }
// // 4.3) Both interfaces nodes in the continuum connects to the atomistic DD
// // line
// else if ((nh_conInt1 != NULL) && (nh_conInt2 != NULL)) {
// // Get the Burgers vector of the dislocation
// burgersVec[0] = nh_conInt1->burgX[0];
// burgersVec[1] = nh_conInt1->burgY[0];
// burgersVec[2] = nh_conInt1->burgZ[0];
// // Get the normal vector of the glide plane
// gpNormVec[0] = nh_conInt1->nx[0];
// gpNormVec[1] = nh_conInt1->ny[0];
// gpNormVec[2] = nh_conInt1->nz[0];
// // Re-allocate memory for the hybrid segment if necessary (--> should
// only
// // be done once before the first iteration)
// if (nh_conInt1->numNbrs < 2) {
// ReallocNodeArms(nh_conInt1, 2);
// nh_conInt1->constraint = UNCONSTRAINED;
// }
// if (nh_conInt2->numNbrs < 2) {
// ReallocNodeArms(nh_conInt2, 2);
// nh_conInt2->constraint = UNCONSTRAINED;
// }
// // Establish connectivity between nodes in the continuum and the
// atomistic
// // domain
// setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
// setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
// gpNormVec);
// // Project atomistic DD nodes on the glide plane
// projectNodeOnGP(nh_conInt1, created_nodes[0], 1);
// projectNodeOnGP(nh_conInt1, created_nodes[nb_points - 1], 1);
// ++numNbrs[0];
// ++numNbrs[nb_points - 1];
// } else {
// LM_FATAL(std::endl
// << "Hybrid dislocation failed to connect!!" << std::endl);
// }
// /* ------------------------------------------------------------------------
// */
// /* 5) Now connect the remainder of the hybrid dislocation line in the
// * atomistic domain */
// for (UInt i = 0; i < nb_elem; ++i) {
// // Get connectivity
// RefGenericElem<3> &el = conn[i];
// UInt idxNhSrc = el.conn[0];
// UInt idxNhDst = el.conn[1];
// UInt &idxArmSrc = numNbrs[idxNhSrc];
// UInt &idxArmDst = numNbrs[idxNhDst];
// // Update arm information for each node of segment i
// setNodeArm(created_nodes[idxNhSrc], created_nodes[idxNhDst], idxArmSrc,
// idxArmDst, burgersVec, gpNormVec);
// // Project atomistic DD node on the glide plane
// projectNodeOnGP(created_nodes[idxNhSrc], created_nodes[idxNhDst],
// idxArmSrc);
// ++idxArmSrc;
// ++idxArmDst;
// }
// }
/* --------------------------------------------------------------------------
*/
UInt DomainParaDiS::getNBNodes() {
return this->container_dd.getContainerNodes().size();
}
/* --------------------------------------------------------------------------
*/
Real DomainParaDiS::getShearModulus() {
return ParaDiSHome::paradis_ptr->param->shearModulus;
}
/* --------------------------------------------------------------------------
*/
VectorProxy<3u> DomainParaDiS::getXmin() {
return VectorProxy<3u>(
std::forward_as_tuple(ParaDiSHome::paradis_ptr->param->xBoundMin,
ParaDiSHome::paradis_ptr->param->yBoundMin,
ParaDiSHome::paradis_ptr->param->zBoundMin));
}
/* --------------------------------------------------------------------------
*/
VectorProxy<3u> DomainParaDiS::getXmax() {
return VectorProxy<3u>(
std::forward_as_tuple(ParaDiSHome::paradis_ptr->param->xBoundMax,
ParaDiSHome::paradis_ptr->param->yBoundMax,
ParaDiSHome::paradis_ptr->param->zBoundMax));
}
/* --------------------------------------------------------------------------
*/
void DomainParaDiS::freenodearm(UInt index) {
Node_t node = *ParaDiSHome::paradis_ptr->nodeKeys[index];
FreeNodeArms(&node);
}
/* --------------------------------------------------------------------------
*/
void DomainParaDiS::freenode(UInt index) {
Node_t node = *ParaDiSHome::paradis_ptr->nodeKeys[index];
FreeNodeArms(&node);
}
/* --------------------------------------------------------------------------
*/
bool DomainParaDiS::is_it_bound() {
// ContainerParaDiSNode::iterator dbgNodes =
// this->getContainerNodes().begin();
// UInt counter = 0;
// for(RefParaDiSNode node=dbgNodes.getFirst(); !dbgNodes.end();
// node=dbgNodes.getNext()){
// ++counter;
// }
return this->boundary;
} // CHANGE THIS?
/* --------------------------------------------------------------------------
*/
Cube DomainParaDiS::getDomainBoundingBox() const {
Cube c;
c.getXmin() = VectorView<3>(ParaDiSHome::paradis_ptr->param->minCoordinates);
c.getXmax() = VectorView<3>(ParaDiSHome::paradis_ptr->param->maxCoordinates);
return c;
}
/* --------------------------------------------------------------------------
*/
Cube DomainParaDiS::getSubDomainBoundingBox() const { LM_TOIMPLEMENT; }
/* --------------------------------------------------------------------------
*/
void DomainParaDiS::setDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }
/* --------------------------------------------------------------------------
*/
void DomainParaDiS::setSubDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }
/* ------------------------------------------------------------------------ */
/* LMDESC PARADIS
This section is used to configure the ParaDiS plugin
*/
/* LMEXAMPLE
Section MultiScale AtomsUnits\\
...\\
GEOMETRY 1 CUBE BBOX -10 10 -10 10 -10 10
MODEL PARADIS dislomodel \\
...\\
endSection\\ \\
Section PARADIS:dislomodel RealUnits \\
NUMDLBSTEPS 0 \\
CONTROLFILE frank_read_src.ctrl \\
DATAFILE frank_read_src.data \\
TIMESTEP 1e-10 \\
DOMAIN_GEOMETRY 1 \\
endSection
*/
/* LMHERITATE domain_dd */
void DomainParaDiS::declareParams() {
DomainDD<ContainerNodesParaDiS, ContainerElemsParaDiS, 3>::declareParams();
/* LMKEYWORD TIMESTEP
Set the maximum time step used for ParaDiS
*/
this->parseKeyword("TIMESTEP", timestep);
/* LMKEYWORD TIMESTEP_RESET_FREQ
Because of collisions, ParaDiS tends to decrease the timestep.
At the provided frequency, LibMultiScale will attempt to reset the
timestep to the desired value (specified by the keyword TIMESTEP)
*/
this->parseKeyword("TIMESTEP_RESET_FREQ", this->timestep_reset_frequency,
100u);
/* LMKEYWORD NUMDLBSTEPS
Number of dynamic load balancing steps
*/
this->parseKeyword("NUMDLBSTEPS", numDLBCycles);
/* LMKEYWORD CONTROLFILE
Path to the ParaDiS control file
*/
this->parseKeyword("CONTROLFILE", control_filename);
/* LMKEYWORD DATAFILE
Path to the ParaDis data file
*/
this->parseKeyword("DATAFILE", data_filename);
/* LMKEYWORD MOBILITYFILE
Mobility input acquired from MD
*/
this->parseKeyword("MOBILITYFILE", mobility_filename, "");
/* LMKEYWORD BOUNDARY
Constraint boundary to mimic grain boundary
*/
this->parseTag("BOUNDARY", this->boundary, false);
this->makeItOptional("DOMAIN_GEOMETRY");
}
/* --------------------------------------------------------------------------
*/
__END_LIBMULTISCALE__

Event Timeline

c4science · Help