diff --git a/extra_packages/cadd b/extra_packages/cadd
index 720be54..a4885a3 160000
--- a/extra_packages/cadd
+++ b/extra_packages/cadd
@@ -1 +1 @@
-Subproject commit 720be54eb817a56ace376c2497fcff8db066a012
+Subproject commit a4885a39424cb59031298d7731401711ef1da022
diff --git a/packages/dumpers.cmake b/packages/dumpers.cmake
index 64b0975..f9ae649 100644
--- a/packages/dumpers.cmake
+++ b/packages/dumpers.cmake
@@ -1,98 +1,98 @@
#===============================================================================
# @file dumpers.cmake
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date Thu Sep 18 16:13:26 2014
#
# @brief Dumper package
#
# @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/>.
#
#===============================================================================
package_declare(dumpers NOT_OPTIONAL DESCRIPTION "dumpers package for LibMultiScale")
set(LIBMULTISCALE_DUMPER_HDRS_LIST_DUMPERS
dumper/dumper.hh
dumper/dumper_1d.hh
dumper/dumper_ekin.hh
dumper/dumper_epot.hh
dumper/dumper_grid.hh
dumper/dumper_lammps.hh
dumper/dumper_paraview.hh
dumper/dumper_restart.hh
dumper/dumper_text.hh
)
set(LIBMULTISCALE_DUMPERS_FILES
${LIBMULTISCALE_DUMPER_HDRS_LIST_DUMPERS}
dumper/dumper_interface.hh
dumper/dumper.cc
- dumper/dumper_dd_paraview.cc
+ dumper/dumper_paraview_dd.cc
dumper/dumper_paraview_atom.cc
dumper/dumper_paraview_mesh.cc
dumper/dumper_paraview.cc
dumper/dumper_1d.cc
dumper/dumper_ekin.cc
dumper/dumper_epot.cc
dumper/dumper_grid.cc
dumper/dumper_lammps.cc
dumper/dumper_restart.cc
dumper/dumper_text.cc
dumper/paraview_helper.hh
dumper/base64_reader.hh
dumper/base64_writer.hh
)
set(LIBMULTISCALE_DUMPER_REF_INPUT_LIST_DUMPERS
"((DumperParaview,PARAVIEW))"
"((Dumper1D,DUMPER1D))"
"((DumperEKin,EKIN))"
"((DumperEPot,EPOT))"
"((DumperGrid,GRID))"
"((DumperLammps,LAMMPS))"
"((DumperRestart,RESTART))"
)
set(LIBMULTISCALE_DUMPER_REAL_INPUT_LIST_DUMPERS
"((DumperText,TEXT))"
)
# ((DumperXYZ,XYZ))
# ((DumperVGroup,VGROUPE))
# ((DumperReduce,REDUCE))
# ((DumperReduce,REDUCE))
# ((DumperUser,USER))
# ((DumperGlobalStat,GLOBALSTAT))
# ((DumperContact,CONTACT))
# ((DumperNeighbor,NEIGHBOR))
# ((DumperCountPoints,COUNTPOINTS))
#LIBMULTISCALE_USE_FIG -> DumperFig,DUMPERFIG
#LIBMULTISCALE_USE_EPSN -> DumperEpsn,EPSN
package_declare_sources(dumpers ${LIBMULTISCALE_DUMPERS_FILES})
set_module_headers_from_package(dumper dumpers ${LIBMULTISCALE_DUMPER_HDRS_LIST_DUMPERS})
declare_boost_list_from_package(dumper_real_input dumpers ${LIBMULTISCALE_DUMPER_REAL_INPUT_LIST_DUMPERS})
-declare_boost_list_from_package(dumper_ref_input dumpers ${LIBMULTISCALE_DUMPER_REF_INPUT_LIST_DUMPERS})
\ No newline at end of file
+declare_boost_list_from_package(dumper_ref_input dumpers ${LIBMULTISCALE_DUMPER_REF_INPUT_LIST_DUMPERS})
diff --git a/src/common/lm_types.hh b/src/common/lm_types.hh
index dafbbe5..40040fc 100644
--- a/src/common/lm_types.hh
+++ b/src/common/lm_types.hh
@@ -1,343 +1,422 @@
/**
* @file lm_types.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Moseley Philip Arthur <philip.moseley@epfl.ch>
*
* @date Mon Sep 08 23:40:22 2014
*
* @brief This where the global types and enums of LibMultiScale are defined
*
* @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/>.
*
*/
#ifndef __LIBMULTISCALE_LM_TYPES_HH__
#define __LIBMULTISCALE_LM_TYPES_HH__
/* -------------------------------------------------------------------------- */
#include "lm_macros.hh"
#include <Eigen/Dense>
#include <ostream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
// error management
/* -------------------------------------------------------------------------- */
#define LM_EXIT_FAILURE 1
#define LM_EXIT_SUCCESS 0
/* -------------------------------------------------------------------------- */
// TYPEDEFS
/* -------------------------------------------------------------------------- */
typedef unsigned int UInt;
typedef double Real;
typedef Real REAL;
/* -------------------------------------------------------------------------- */
// ENUMS
/* -------------------------------------------------------------------------- */
enum DOFType {
dt_local_master = 1,
dt_local_slave = 2,
dt_local = 4,
dt_ghost = 8,
dt_all = 255
};
/* -------------------------------------------------------------------------- */
enum GeomType {
BALL = 1,
CUBE = 2,
INTER = 3,
SUB = 4,
ELLIPSOID = 5,
UNION = 6,
CYLINDER = 7,
CHAUDRON = 8,
CUBE_SURFACE = 9
};
/* -------------------------------------------------------------------------- */
enum COutputType { OUTPUT_REF = 1, OUTPUT_REAL = 2, OUTPUT_NULL = 32767 };
/* -------------------------------------------------------------------------- */
enum CouplingStage {
COUPLING_STEP1 = 1,
COUPLING_STEP2 = 2,
COUPLING_STEP3 = 4,
COUPLING_STEP4 = 8,
COUPLING_STEP5 = 16,
COUPLING_STEP6 = 32,
COUPLING_STEP7 = 64
};
/* -------------------------------------------------------------------------- */
enum Operator {
OP_NONE = 0,
OP_AVERAGE = 1,
OP_DENSITY = 2,
OP_DEVIATION = 4,
OP_MAX = 8,
OP_MIN = 16,
OP_SUM = 64
};
/* -------------------------------------------------------------------------- */
enum ActionType {
at_stimulator = 0,
at_dumper = 1,
at_filter = 2,
at_compute = 4,
at_uninitialised = 32767
};
/* -------------------------------------------------------------------------- */
enum LatticeType { _not_defined, mono_atom_1d, hex_2d, hex_2d_radial, fcc_3d };
template <LatticeType type> struct LatticeDim { static const UInt value; };
template <> struct LatticeDim<mono_atom_1d> { static const UInt value = 1; };
template <> struct LatticeDim<hex_2d> { static const UInt value = 2; };
template <> struct LatticeDim<hex_2d_radial> { static const UInt value = 2; };
template <> struct LatticeDim<fcc_3d> { static const UInt value = 3; };
/* -------------------------------------------------------------------------- */
enum IntegrationSchemeStage {
NONE_STEP = 0,
PRE_DUMP = 1,
PRE_STEP1 = 2,
PRE_STEP2 = 4,
PRE_STEP3 = 8,
PRE_STEP4 = 16,
PRE_STEP5 = 32,
PRE_STEP6 = 64,
PRE_STEP7 = 128,
PRE_FATAL = 256, // This stage only exists if a FATAL is triggered
ALL_STEP = 65536
};
/* -------------------------------------------------------------------------- */
class IntegrationSchemeMask {
public:
IntegrationSchemeMask(UInt mask) { stage = mask; }
IntegrationSchemeMask(const IntegrationSchemeStage &s) { stage = s; }
IntegrationSchemeMask() { stage = NONE_STEP; };
virtual ~IntegrationSchemeMask(){};
bool operator&(const IntegrationSchemeStage &s) const { return (stage & s); };
IntegrationSchemeMask &operator=(const IntegrationSchemeStage &s) {
stage = s;
return (*this);
};
IntegrationSchemeMask operator|=(const IntegrationSchemeStage &s) {
stage |= s;
return (*this);
};
IntegrationSchemeMask operator|=(const IntegrationSchemeMask &m) {
stage |= m.stage;
return (*this);
};
bool isNone() { return (stage == NONE_STEP); }
virtual void printself(std::ostream &stream, int indent = 0) const {
stream << "NONE_STEP";
if ((*this) & PRE_DUMP)
stream << "|PREDUMP";
if ((*this) & PRE_STEP1)
stream << "|PRE_STEP1";
if ((*this) & PRE_STEP2)
stream << "|PRE_STEP2";
if ((*this) & PRE_STEP3)
stream << "|PRE_STEP3";
if ((*this) & PRE_STEP4)
stream << "|PRE_STEP4";
if ((*this) & PRE_STEP5)
stream << "|PRE_STEP5";
if ((*this) & PRE_STEP6)
stream << "|PRE_STEP6";
if ((*this) & PRE_STEP7)
stream << "|PRE_STEP7";
if ((*this) & PRE_FATAL)
stream << "|PRE_FATAL";
if ((*this) & ALL_STEP)
stream << "|ALL_STEP";
};
private:
UInt stage;
};
inline std::ostream &operator<<(std::ostream &stream,
const IntegrationSchemeMask &_this) {
_this.printself(stream);
return stream;
}
inline std::ostream &operator<<(std::ostream &stream,
const IntegrationSchemeStage &_this) {
switch (_this) {
case NONE_STEP:
stream << "NONESTEP";
break;
case PRE_DUMP:
stream << "PREDUMP";
break;
case PRE_STEP1:
stream << "PRE_STEP1";
break;
case PRE_STEP2:
stream << "PRE_STEP2";
break;
case PRE_STEP3:
stream << "PRE_STEP3";
break;
case PRE_STEP4:
stream << "PRE_STEP4";
break;
case PRE_STEP5:
stream << "PRE_STEP5";
break;
case PRE_STEP6:
stream << "PRE_STEP6";
break;
case PRE_STEP7:
stream << "PRE_STEP7";
break;
case PRE_FATAL:
stream << "PRE_FATAL";
break;
case ALL_STEP:
stream << "ALL_STEP";
break;
}
return stream;
}
/* -------------------------------------------------------------------------- */
class CommGroup {
public:
static const int id_all;
static const int id_none;
static const int id_invalid;
CommGroup() { group_id = id_invalid; }
CommGroup(int id) { this->setID(id); }
bool operator==(CommGroup &grp) { return (grp.group_id == group_id); }
bool operator!=(CommGroup &grp) { return (grp.group_id != group_id); }
void setID(int id) { group_id = id; };
int getID() { return group_id; };
private:
int group_id;
};
/* -------------------------------------------------------------------------- */
enum PhysicalQuantity {
// physical quantities
Length,
Mass,
Energy,
Time,
MassDensity,
Force,
Pressure,
Temperature
};
/* -------------------------------------------------------------------------- */
// enum DDModelType {
// _DD2D,
// _PARADIS,
// };
/* -------------------------------------------------------------------------- */
/* describes the orientation of the edge dislo by indication which
* quadrant is shifted towards x=0*/
enum dislo_orientation {
TOP_RIGHT = 0,
TOP_LEFT = 1,
BOTTOM_RIGHT = 2,
BOTTOM_LEFT = 3
};
/* -------------------------------------------------------------------------- */
-template <UInt Dim, typename T=Real> using Vector = Eigen::Matrix<T, Dim, 1>;
-template <UInt Dim> using Array1D = Eigen::Array<Real, Dim, 1>;
-template <UInt Dim> using ArrayView1D = Eigen::Map<Eigen::Array<Real, Dim, 1>>;
-template <UInt Dim> using Matrix = Eigen::Matrix<Real, Dim, Dim>;
-template <UInt Dim> using VectorView = Eigen::Map<Eigen::Matrix<Real, Dim, 1>>;
+// template <UInt Dim, typename T=Real> using Vector = Eigen::Matrix<T, Dim,
+// 1u>;
+template <UInt Dim, typename T = Real>
+struct Vector : public Eigen::Matrix<T, Dim, 1u> {
+ using Eigen::Matrix<T, Dim, 1u>::Matrix;
+};
+// template <UInt Dim> using Array1D = Eigen::Array<Real, Dim, 1u>;
+template <UInt Dim> struct Array1D : public Eigen::Array<Real, Dim, 1u> {
+ using Eigen::Array<Real, Dim, 1u>::Array;
+ using Eigen::Array<Real, Dim, 1u>::operator=;
+};
+// template <UInt Dim> using ArrayView1D = Eigen::Map<Eigen::Array<Real, Dim,
+// 1u>>;
+template <UInt Dim>
+struct ArrayView1D : public Eigen::Map<Eigen::Array<Real, Dim, 1u>> {
+ using Eigen::Map<Eigen::Array<Real, Dim, 1u>>::Map;
+};
+// template <UInt Dim> using Matrix = Eigen::Matrix<Real, Dim, Dim>;
+template <UInt Dim> struct Matrix : public Eigen::Matrix<Real, Dim, Dim> {
+ using Eigen::Matrix<Real, Dim, Dim>::Matrix;
+};
+// template <UInt Dim> using VectorView = Eigen::Map<Eigen::Matrix<Real, Dim,
+// 1u>>;
+template <UInt Dim>
+struct VectorView : public Eigen::Map<Eigen::Matrix<Real, Dim, 1u>> {
+ using Eigen::Map<Eigen::Matrix<Real, Dim, 1u>>::Map;
+ using Eigen::Map<Eigen::Matrix<Real, Dim, 1u>>::Map::operator=;
+};
+// template <UInt Dim> using MatrixView = Eigen::Map<Eigen::Matrix<Real, Dim,
+// Dim>>;
template <UInt Dim>
-using MatrixView = Eigen::Map<Eigen::Matrix<Real, Dim, Dim>>;
+struct MatrixView : public Eigen::Map<Eigen::Matrix<Real, Dim, Dim>> {
+ using Eigen::Map<Eigen::Matrix<Real, Dim, Dim>>::Map;
+};
template <typename T> bool isNaN(T t) {
Eigen::ArrayXd a(t);
return a.isNaN().any();
}
template <typename T> bool isInf(T t) {
Eigen::ArrayXd a(t);
return a.isInf().any();
}
enum FieldType {
_position0 = 0,
_position = 1,
_displacement = 2,
_velocity = 3,
_force = 4,
_stress = 5,
_temperature = 6,
_grouprank = 7,
_strain = 8,
_epot = 9,
_applied_force = 10,
_mass = 11,
_tag = 12,
ft_uninitialised = 32767
};
template <UInt Dim, FieldType ft> struct FieldTraits {
using TensorType = Vector<Dim>;
using TensorViewType = VectorView<Dim>;
};
/* -------------------------------------------------------------------------- */
+// VectorProxy set of classes
+// helpful to wrap n variables (of the same type)
+// and manipulate them/convert them as an eigen vector
+
+template <typename T, int dim> struct vector_proxy_getter {
+ static double &item_getter(T &v, unsigned int i) {
+ if (dim == i)
+ return std::get<dim>(v);
+ return vector_proxy_getter<T, dim - 1>::item_getter(v, i);
+ }
+};
+
+template <typename T> struct vector_proxy_getter<T, -1> {
+ static double &item_getter(T &v, unsigned int i) {
+ return *new double(std::nan("1"));
+ }
+};
+
+template <typename T> struct vector_proxy {
+ vector_proxy(T tup) : tup(tup){};
+
+ auto &operator[](unsigned int i) {
+ return vector_proxy_getter<T, std::tuple_size<T>::value - 1>::item_getter(
+ this->tup, i);
+ }
+
+ T tup;
+};
+
+template <unsigned int Dim, typename... Args> struct real_tuple {
+ using type = typename real_tuple<Dim - 1, double &, Args...>::type;
+};
+
+template <typename... Args> struct real_tuple<0, Args...> {
+ using type = std::tuple<Args...>;
+};
+
+template <unsigned int Dim>
+struct VectorProxy : public vector_proxy<typename real_tuple<Dim>::type> {
+
+ using tuple_type = typename real_tuple<Dim>::type;
+
+ template <typename... Args>
+ VectorProxy(Args &... args) : vector_proxy<tuple_type>(tuple_type(args...)){};
+
+ operator Eigen::Matrix<double, Dim, 1u>() {
+ Eigen::Matrix<double, Dim, 1u> res;
+ for (unsigned int i = 0; i < Dim; ++i)
+ res[i] = (*this)[i];
+ return res;
+ }
+};
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
#endif /* __LIBMULTISCALE_LM_TYPES_HH__ */
diff --git a/src/dd/paradis/domain_paradis.cc b/src/dd/paradis/domain_paradis.cc
index 3f128a3..3fad5d0 100644
--- a/src/dd/paradis/domain_paradis.cc
+++ b/src/dd/paradis/domain_paradis.cc
@@ -1,912 +1,902 @@
/**
* @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();
+ 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 ¤t_arm_src = numNbrs[nd_src];
UInt ¤t_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();
+ 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");
+VectorProxy<3> DomainPARADIS::getXmin() {
+ return VectorProxy<3>(this->paradis_ptr->param->xBoundMin,
+ this->paradis_ptr->param->yBoundMin,
+ this->paradis_ptr->param->zBoundMin);
}
/* -------------------------------------------------------------------------- */
-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");
+VectorProxy<3> DomainPARADIS::getXmax() {
+ return VectorProxy<3>(this->paradis_ptr->param->xBoundMax,
+ this->paradis_ptr->param->yBoundMax,
+ this->paradis_ptr->param->zBoundMax);
}
/* -------------------------------------------------------------------------- */
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__
diff --git a/src/dd/paradis/domain_paradis.hh b/src/dd/paradis/domain_paradis.hh
index 575654a..83c44e0 100644
--- a/src/dd/paradis/domain_paradis.hh
+++ b/src/dd/paradis/domain_paradis.hh
@@ -1,218 +1,218 @@
/**
* @file domain_paradis.hh
*
* @author Till Junge <till.junge@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/>.
*
*/
#ifndef __LIBMULTISCALE_DOMAIN_PARADIS_HH__
#define __LIBMULTISCALE_DOMAIN_PARADIS_HH__
/* -------------------------------------------------------------------------- */
#include "container_array.hh"
#include "container_paradis.hh"
#include "domain_dd.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/**
* Class DummySource
* DD dummy container (defect fictive container ?)
*/
typedef ContainerPARADISNode DummySource;
typedef ContainerPARADISNode DummyObstacle;
/* -------------------------------------------------------------------------- */
/**
* Class DomainPARADIS
* domain to implement the plugin interface to PARADIS dislocation dynamics
*/
class DomainPARADIS
: public DomainDD<ContainerPARADISNode, DummySource, DummyObstacle, 3> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DomainPARADIS(DomainID ID, CommGroup GID)
: DomainDD<ContainerPARADISNode, DummySource, DummyObstacle, 3>(ID, GID) {
paradis_ptr = nullptr;
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
void performStep1();
void performStep2();
void performStep3();
void performStep4();
void performStep5();
void imageForce(const double *position, const double *normal, double *force);
void getSubDomainDimensions(Real *xmin, Real *xmax) { LM_TOIMPLEMENT; };
void getDomainDimensions(Real *xmin, Real *xmax) {
for (UInt dir = 0; dir < this->getDim(); ++dir) {
xmin[dir] = this->paradis_ptr->param->minCoordinates[dir];
xmax[dir] = this->paradis_ptr->param->maxCoordinates[dir];
}
};
//! compute displacement field for a set of positions
void computeDisplacements(std::vector<Real> &positions,
std::vector<Real> &displacements);
//! compute displacement field for a position with a given segment
void computeBarnett(Real *A, Real *B, Real *C, Real *burg, Real *dispPointU,
Real *dispPointX, Real pois);
void computeClosurePoint(Real *A, Real *B, Real *burg, Real *C);
void computeProjectedPoint(Real *A, Real *C, Real *burg, Real *P);
//! remesh
void remesh();
//! intialisation stage
void init();
//! parse the keywords
void declareParams();
//! add segments to the dd model
virtual void addSegments(std::vector<RefPointData<3>> &pos,
std::vector<RefGenericElem<3>> &conn,
std::vector<Real> &burgers,
std::vector<Real> &glide_normals);
UInt getNBNodes();
bool is_it_bound();
- Real getXmin(UInt dim);
- Real getXmax(UInt dim);
+ VectorProxy<3> getXmin();
+ VectorProxy<3> getXmax();
Real getShearModulus();
void freenodearm(UInt index);
void freenode(UInt index);
//! Add hybrid segments to the DD model
virtual void 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
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
//! set current step to a given value
void setCurrentStep(UInt ts) { LM_FATAL("to be implemented"); } // GUILLAUME?
//! get current step
UInt getCurrentStep() {
LM_FATAL("to be implemented");
return 0;
}; // CHANGE THIS
//! set timestep to a given value
void setTimeStep(Real ts){LM_FATAL("to be implemented")}; // CHANGE THIS
//! GetTimeStep
Real getTimeStep() { return this->timestep; }; // CHANGE THIS
//! return the stress component (what is i and j ?)
double &get_stress_component(int i, int j) {
int indices[3][3] = {{0, 5, 4}, {5, 1, 3}, {4, 3, 2}};
return this->paradis_ptr->param->appliedStress[indices[i][j]];
}
UInt getDim() { return 3; };
Real getPoisson() { return this->paradis_ptr->param->pois; };
Real getMaxSeg() { return this->paradis_ptr->param->maxSeg; };
Real getMinSeg() { return this->paradis_ptr->param->minSeg; };
Real getBurgersMagnitude() { return this->paradis_ptr->param->burgMag; };
UInt getNBpairsInfinite() { return this->paradis_ptr->param->numPairs; };
bool getPairedIndex(int index0, int *index0_infinite) {
UInt num_pairs = this->paradis_ptr->param->numPairs;
for (int i = 0; i < (int)num_pairs; ++i) {
Pair_t pair_ind;
pair_ind = *this->paradis_ptr->pairKeys[i];
if (pair_ind.index1 == index0) {
*index0_infinite = pair_ind.index2;
return true;
} else if (pair_ind.index2 == index0) {
*index0_infinite = pair_ind.index1;
return true;
}
}
return false;
};
/* ------------------------------------------------------------------------ */
/* Functions for AMELCG */
/* ------------------------------------------------------------------------ */
double getExternalWork() { return 0; }
void setExternalWork(double) {}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
bool boundary;
int memSize;
Home_t *paradis_ptr;
Quantity<Time> timestep;
UInt numDLBCycles;
std::string control_filename;
std::string data_filename;
// ParaDis - Containers/Lists
// list of ParaDis node handles of
// "continuum" DD nodes next to the
// artificial interface
std::vector<Node_t *> nh_l_conInt;
// list of ParaDis node handles confined
// to an atomistic domain
std::vector<Node_t *> created_nodes;
// NOTE: vector changes every iteration depending on a dislocation detection
// algorithm
std::string mobility_filename;
};
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_DOMAIN_PARADIS_HH__ */
diff --git a/src/dd/paradis/ref_node_paradis.hh b/src/dd/paradis/ref_node_paradis.hh
index 0560ab8..e1e1638 100644
--- a/src/dd/paradis/ref_node_paradis.hh
+++ b/src/dd/paradis/ref_node_paradis.hh
@@ -1,397 +1,397 @@
/**
* @file ref_node_paradis.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Max Hodapp <max.hodapp@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Mon Oct 28 19:23:14 2013
*
* @brief Nodal reference of ParaDiS
*
* @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/>.
*
*/
#ifndef __LIBMULTISCALE_REF_NODE_PARADIS_HH__
#define __LIBMULTISCALE_REF_NODE_PARADIS_HH__
/* -------------------------------------------------------------------------- */
// ParaDis includes
#define PARALLEL
#include <mpi.h>
extern "C" {
#include "Home.h"
#include "Typedefs.h"
}
#undef PARALLEL
/* -------------------------------------------------------------------------- */
#include "ref_node_dd.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
class DomainPARADIS;
class ComparatorRefPARADISNode;
/* -------------------------------------------------------------------------- */
/**
* Class RefPARADISNode
*
*/
class RefPARADISNode : public RefDDNode<3> {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
//! generic container of the original model container
typedef DomainPARADIS Domain;
//! node comparator class to be usable in maps
typedef ComparatorRefPARADISNode RefComparator;
static const UInt Dim = 3;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
RefPARADISNode() : fakestress(0) {}
RefPARADISNode(int index) : fakestress(0) { index_node = index; }
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
bool operator==(RefPARADISNode &ref) {
return (ref.index_node == index_node);
};
VectorView<3> position0() { return position(); }
VectorView<3> position() {
LM_TOIMPLEMENT;
// switch (i) {
// case 0:
// return paradis_ptr->nodeKeys[this->index_node]->x;
// break;
// case 1:
// return paradis_ptr->nodeKeys[this->index_node]->y;
// break;
// case 2:
// return paradis_ptr->nodeKeys[this->index_node]->z;
// break;
// }
- // double *p = NULL;
+ // Real *p = NULL;
// return *p;
}
- double &prev_position(unsigned int i) {
+ Real &prev_position(unsigned int i) {
switch (i) {
case 0:
return paradis_ptr->nodeKeys[this->index_node]->oldx;
break;
case 1:
return paradis_ptr->nodeKeys[this->index_node]->oldy;
break;
case 2:
return paradis_ptr->nodeKeys[this->index_node]->oldz;
break;
}
- double *p = NULL;
+ Real *p = NULL;
return *p;
}
VectorView<3> velocity() {
LM_TOIMPLEMENT;
// switch (i) {
// case 0:
// return paradis_ptr->nodeKeys[this->index_node]->vX;
// break;
// case 1:
// return paradis_ptr->nodeKeys[this->index_node]->vY;
// break;
// case 2:
// return paradis_ptr->nodeKeys[this->index_node]->vZ;
// break;
// }
- // double *p = NULL;
+ // Real *p = NULL;
// return *p;
}
- double &burgersVectorX(unsigned int n) {
+ Real &burgersVectorX(unsigned int n) {
return paradis_ptr->nodeKeys[this->index_node]->burgX[n];
}
- double &burgersVectorY(unsigned int n) {
+ Real &burgersVectorY(unsigned int n) {
return paradis_ptr->nodeKeys[this->index_node]->burgY[n];
}
- double &burgersVectorZ(unsigned int n) {
+ Real &burgersVectorZ(unsigned int n) {
return paradis_ptr->nodeKeys[this->index_node]->burgZ[n];
}
- double &normalVectorX(unsigned int n) {
+ Real &normalVectorX(unsigned int n) {
return paradis_ptr->nodeKeys[this->index_node]->nx[n];
}
- double &normalVectorY(unsigned int n) {
+ Real &normalVectorY(unsigned int n) {
return paradis_ptr->nodeKeys[this->index_node]->ny[n];
}
- double &normalVectorZ(unsigned int n) {
+ Real &normalVectorZ(unsigned int n) {
return paradis_ptr->nodeKeys[this->index_node]->nz[n];
}
int &numNbrsIndexNode(unsigned int idx) {
return paradis_ptr->nodeKeys[idx]->numNbrs;
}
- double &getNeighborNodeCoord(unsigned int n, unsigned int i) {
+ Real &getNeighborNodeCoord(unsigned int n, unsigned int i) {
UInt my_domain = paradis_ptr->myDomain;
Node_t current = *paradis_ptr->nodeKeys[this->index_node];
Tag_t *tag = current.nbrTag;
UInt tag_domain = tag[n].domainID;
int idx_for_neighbor = tag[n].index;
if (my_domain == tag_domain) {
switch (i) {
case 0:
return paradis_ptr->nodeKeys[idx_for_neighbor]->x;
break;
case 1:
return paradis_ptr->nodeKeys[idx_for_neighbor]->y;
break;
case 2:
return paradis_ptr->nodeKeys[idx_for_neighbor]->z;
break;
}
} else {
switch (i) {
case 0:
return paradis_ptr->remoteDomainKeys[tag_domain]
->nodeKeys[idx_for_neighbor]
->x;
break;
case 1:
return paradis_ptr->remoteDomainKeys[tag_domain]
->nodeKeys[idx_for_neighbor]
->y;
break;
case 2:
return paradis_ptr->remoteDomainKeys[tag_domain]
->nodeKeys[idx_for_neighbor]
->z;
break;
}
}
- double *p = NULL;
+ Real *p = NULL;
return *p;
}
int &getNeighborNodeIndex0(unsigned int n) {
Node_t current = *paradis_ptr->nodeKeys[this->index_node];
Tag_t *tag = current.nbrTag;
return tag[n].index0;
}
int &getNeighborNodeIndex(unsigned int n) {
Node_t current = *paradis_ptr->nodeKeys[this->index_node];
Tag_t *tag = current.nbrTag;
return tag[n].index;
}
int &getNeighborNodeSlave(unsigned int n) {
UInt my_domain = paradis_ptr->myDomain;
Node_t current = *paradis_ptr->nodeKeys[this->index_node];
Tag_t *tag = current.nbrTag;
UInt tag_domain = tag[n].domainID;
int idx_for_neighbor = tag[n].index;
if (my_domain == tag_domain)
return paradis_ptr->nodeKeys[idx_for_neighbor]->slave;
else
return paradis_ptr->remoteDomainKeys[tag_domain]
->nodeKeys[idx_for_neighbor]
->slave;
}
int &getNbrNeighborsOfNeighbor(unsigned int n) {
UInt my_domain = paradis_ptr->myDomain;
Node_t current = *paradis_ptr->nodeKeys[this->index_node];
Tag_t *tag = current.nbrTag;
UInt tag_domain = tag[n].domainID;
int idx_for_neighbor = tag[n].index;
if (my_domain == tag_domain)
return paradis_ptr->nodeKeys[idx_for_neighbor]->numNbrs;
else
return paradis_ptr->remoteDomainKeys[tag_domain]
->nodeKeys[idx_for_neighbor]
->numNbrs;
}
UInt &tag() { return this->index_node; }
int &getIndex0() {
return paradis_ptr->nodeKeys[this->index_node]->myTag.index0;
}
int &getIndex() {
return paradis_ptr->nodeKeys[this->index_node]->myTag.index;
}
VectorView<3> force() {
LM_TOIMPLEMENT;
// switch (i) {
// case 0:
// return paradis_ptr->nodeKeys[this->index_node]->fX;
// break;
// case 1:
// return paradis_ptr->nodeKeys[this->index_node]->fY;
// break;
// case 2:
// return paradis_ptr->nodeKeys[this->index_node]->fZ;
// break;
// }
- // double *p = NULL;
+ // Real *p = NULL;
// return *p;
}
MatrixView<3> stress() {
LM_TOIMPLEMENT;
// PROBLEM!!! ParaDiS only stores projected
// stresses \sigma \cdot b. ?????WHAT TO DO?
// return fakestress; // nodal_stress->at(index_node,i);
}
void getIndexes(int &n) { n = index_node; }
void setIndex(UInt index) { index_node = index; }
void setHome(Home_t *paradis_ptr) { this->paradis_ptr = paradis_ptr; };
// Get the number of neighbors
UInt getNbrOfNeighs() {
return UInt(this->paradis_ptr->nodeKeys[this->index_node]->numNbrs);
}
int getSlave() {
return this->paradis_ptr->nodeKeys[this->index_node]->slave;
}
+ //! constraint tag:
int getConstraint() {
return this->paradis_ptr->nodeKeys[this->index_node]->constraint;
}
void slaving() { paradis_ptr->nodeKeys[this->index_node]->slave = 1; }
void unslaving() { paradis_ptr->nodeKeys[this->index_node]->slave = 0; }
void constrain() { paradis_ptr->nodeKeys[this->index_node]->constraint = 7; }
void unconstrain() {
paradis_ptr->nodeKeys[this->index_node]->constraint = 0;
}
void setOldVel(unsigned int i, Real oldv) {
switch (i) {
case 0:
paradis_ptr->nodeKeys[this->index_node]->oldvX = oldv;
return;
case 1:
paradis_ptr->nodeKeys[this->index_node]->oldvY = oldv;
return;
case 2:
paradis_ptr->nodeKeys[this->index_node]->oldvZ = oldv;
return;
}
return;
}
void release(UInt nbr_index) {
paradis_ptr->nodeKeys[this->index_node]->constraint = 0;
UInt index_for_neighbor =
paradis_ptr->nodeKeys[this->index_node]->nbrTag[0].index;
paradis_ptr->nodeKeys[this->index_node]->oldvX =
0.5 * (paradis_ptr->nodeKeys[index_for_neighbor]->oldvX +
paradis_ptr->nodeKeys[nbr_index]->oldvX);
paradis_ptr->nodeKeys[this->index_node]->oldvY =
0.5 * (paradis_ptr->nodeKeys[index_for_neighbor]->oldvY +
paradis_ptr->nodeKeys[nbr_index]->oldvY);
paradis_ptr->nodeKeys[this->index_node]->oldvZ =
0.5 * (paradis_ptr->nodeKeys[index_for_neighbor]->oldvZ +
paradis_ptr->nodeKeys[nbr_index]->oldvZ);
InsertArm(paradis_ptr, paradis_ptr->nodeKeys[this->index_node],
¶dis_ptr->nodeKeys[nbr_index]->myTag,
-1.0 * (paradis_ptr->nodeKeys[this->index_node]->burgX[0]),
-1.0 * (paradis_ptr->nodeKeys[this->index_node]->burgY[0]),
-1.0 * (paradis_ptr->nodeKeys[this->index_node]->burgZ[0]),
paradis_ptr->nodeKeys[this->index_node]->nx[0],
paradis_ptr->nodeKeys[this->index_node]->ny[0],
paradis_ptr->nodeKeys[this->index_node]->nz[0], 0, 1);
paradis_ptr->nodeKeys[this->index_node]->numNbrs = 2;
}
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
friend class ContainerPARADISNode;
private:
Home_t *paradis_ptr;
unsigned int index_node;
unsigned int index0_node;
- double
- fakestress; // This is a temprary fix to make the dumpers work: FIXMETILL
+ Real fakestress; // This is a temprary fix to make the dumpers work: FIXMETILL
};
/* -------------------------------------------------------------------------- */
class ComparatorRefPARADISNode {
public:
bool operator()(RefPARADISNode &r1, RefPARADISNode &r2) const {
// return r1.index_atome < r2.index_atome;
LM_TOIMPLEMENT;
}
bool operator()(const RefPARADISNode &r1, const RefPARADISNode &r2) const {
LM_TOIMPLEMENT;
}
};
/* -------------------------------------------------------------------------- */
inline std::ostream &operator<<(std::ostream &os, RefPARADISNode &ref) {
int n;
ref.getIndexes(n);
os << "PARADISNode reference : " << n;
return os;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_REF_NODE_PARADIS_HH__ */
diff --git a/src/dd/ref_node_dd.hh b/src/dd/ref_node_dd.hh
index 980a74a..3718012 100644
--- a/src/dd/ref_node_dd.hh
+++ b/src/dd/ref_node_dd.hh
@@ -1,141 +1,142 @@
/**
* @file ref_node_dd.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Max Hodapp <max.hodapp@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Fri Nov 15 14:49:04 2013
*
* @brief Generic reference to DD nodes
*
* @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/>.
*
*/
#ifndef __LIBMULTISCALE_REF_NODE_DD_HH__
#define __LIBMULTISCALE_REF_NODE_DD_HH__
/* ----------------------------------------------------------------- */
#include "ref_point.hh"
/* ----------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* ----------------------------------------------------------------- */
template <UInt Dim> class RefDDNode;
/* ----------------------------------------------------------------- */
template <UInt Dim, FieldType ftype> class AccessorDDNodalDof {
public:
using TensorType = typename FieldTraits<Dim, ftype>::TensorViewType;
AccessorDDNodalDof(RefDDNode<Dim> &nd) : field(nullptr) { LM_TOIMPLEMENT; };
template <typename T> inline TensorType &operator=(const T &val) {
field = val;
return field;
}
template <typename T> inline T operator=(const T val[]) {
field = Eigen::Map<TensorType>(val);
return field;
}
template <typename T> inline T operator+=(const T &val) {
field += val;
return field;
}
inline Real &operator[](UInt index) { return field[index]; }
private:
TensorType field;
};
/* -------------------------------------------------------------------------- */
/**
* Class Refnode
*
*/
template <UInt Dim> class RefDDNode : public RefPoint<Dim> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
virtual ~RefDDNode(){};
RefDDNode(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
virtual VectorView<Dim> position() = 0;
- virtual Real &prev_position(UInt i) = 0;
+ // virtual Real &prev_position(UInt i) = 0;
virtual VectorView<Dim> position0() = 0;
- virtual Real &burgersVectorX(UInt n) = 0;
- virtual Real &burgersVectorY(UInt n) = 0;
- virtual Real &burgersVectorZ(UInt n) = 0;
+ virtual Real & burgersVectorX(UInt n) = 0;
+ virtual Real & burgersVectorY(UInt n) = 0;
+ virtual Real & burgersVectorZ(UInt n) = 0;
virtual VectorView<Dim> force() = 0;
virtual MatrixView<Dim> stress() = 0;
virtual Real getEPot() { LM_TOIMPLEMENT; };
- virtual Real &getNeighborNodeCoord(UInt n, UInt i) = 0;
- virtual int &getNeighborNodeIndex(UInt n) = 0;
- virtual int &getNeighborNodeSlave(UInt n) = 0;
- virtual int &getNbrNeighborsOfNeighbor(UInt n) = 0;
-
- virtual int &numNbrsIndexNode(UInt idx) = 0;
+ // virtual Real &getNeighborNodeCoord(UInt n, UInt i) = 0;
+ // virtual int &getNeighborNodeIndex(UInt n) = 0;
+ // virtual int &getNeighborNodeSlave(UInt n) = 0;
+ // virtual int &getNbrNeighborsOfNeighbor(UInt n) = 0;
+ // virtual int &numNbrsIndexNode(UInt idx) = 0;
+
Real &mass() {
LM_FATAL("this function does not exist for Dislocation dynamics models");
};
VectorView<Dim> velocity() {
LM_FATAL("this function does not exist for Dislocation dynamics models");
};
VectorView<Dim> displacement() {
LM_FATAL("this function does not exist for Dislocation dynamics models");
};
UInt getIndex() {
LM_FATAL("this function does not exist for Dislocation dynamics models");
};
+
template <FieldType ftype> AccessorDDNodalDof<Dim, ftype> field(){
AccessorDDNodalDof<Dim, ftype> acc(*this);
return acc;
}
// Get nodex index
virtual void getIndexes(int &n) = 0;
// Get nbr of neighbors
virtual UInt getNbrOfNeighs() = 0;
// Get index of the n-th neighbor
// virtual void getNeighIdx (int &idxNh_neigh, const int n ) = 0;
public:
static constexpr UInt dim = Dim;
};
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_REF_NODE_DD_HH__ */
diff --git a/src/dumper/dumper_dd_paraview.cc b/src/dumper/dumper_paraview_dd.cc
similarity index 100%
rename from src/dumper/dumper_dd_paraview.cc
rename to src/dumper/dumper_paraview_dd.cc
diff --git a/src/geometry/cube.cc b/src/geometry/cube.cc
index 6b65d08..fc48141 100644
--- a/src/geometry/cube.cc
+++ b/src/geometry/cube.cc
@@ -1,192 +1,192 @@
/**
* @file cube.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Mon Nov 25 17:17:19 2013
*
* @brief Cubic geometry
*
* @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/>.
*
*/
#include "cube.hh"
#include "geometry_manager.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
const Quantity<Length, 3> &Cube::getXmin() const { return this->xmin; }
/* -------------------------------------------------------------------------- */
const Quantity<Length, 3> &Cube::getXmax() const { return this->xmax; }
/* -------------------------------------------------------------------------- */
Cube &Geometry::getBoundingBox() {
Cube *bbox = NULL;
GeomID bboxID = id + ":bbox";
try {
bbox =
dynamic_cast<Cube *>(GeometryManager::getManager().getGeometry(bboxID));
GeometryManager::getManager().addGeometry(bbox);
} catch (LibMultiScaleException e) {
bbox = new Cube(dim, bboxID);
}
return *bbox;
}
/* -------------------------------------------------------------------------- */
Cube::Cube(UInt dim, GeomID id) : Geometry(dim, CUBE, id) {
for (UInt i = 0; i < 3; ++i) {
- ranges[2 * i] = 0;
- ranges[2 * i + 1] = 0;
+ xmin[i] = 0;
+ xmax[i] = 0;
hole[i] = 0;
}
empty = true;
}
/* -------------------------------------------------------------------------- */
void Cube::setXmin(UInt i, Real r) {
- this->ranges[2 * i + 0] = r;
+ this->xmin[i] = r;
this->init();
}
/* -------------------------------------------------------------------------- */
void Cube::setXmax(UInt i, Real r) {
- this->ranges[2 * i + 1] = r;
+ this->xmax[i] = r;
this->init();
}
/* -------------------------------------------------------------------------- */
const Real &Cube::Hole(UInt i) const { return hole[i]; }
/* -------------------------------------------------------------------------- */
void Cube::setHole(Real *X) {
for (UInt i = 0; i < 3; ++i) {
hole[i] = X[i];
}
}
/* -------------------------------------------------------------------------- */
void Cube::resetDimensions() {
for (UInt i = 0; i < 3; ++i) {
setXmax(i, 0.0);
setXmin(i, 0.0);
}
init();
}
/* -------------------------------------------------------------------------- */
Real *Cube::getSize() { return size; }
/* -------------------------------------------------------------------------- */
Real Cube::getSize(UInt i) { return size[i]; }
/* -------------------------------------------------------------------------- */
Real *Cube::getSizeHalf() { return size_half; }
/* -------------------------------------------------------------------------- */
Real *Cube::getInvSizeHalf() { return inv_size_half; }
/* -------------------------------------------------------------------------- */
Cube &Cube::getBoundingBox() { return *this; }
/* -------------------------------------------------------------------------- */
bool Cube::doIntersect(Geometry &geom) {
if (geom.getType() == Geometry::CUBE) {
Cube &c1 = *this;
Cube &c2 = dynamic_cast<Cube &>(geom);
DUMP("Intersect operation : " << c1.getXmin() << " " << c1.getXmax(),
DBG_DETAIL);
UInt test = true;
for (UInt i = 0; i < dim; ++i) {
if ((c1.getXmin()[i] <= c2.getXmin()[i]) &&
(c2.getXmin()[i] <= c1.getXmax()[i]))
test &= true;
else if ((c1.getXmin()[i] <= c2.getXmax()[i]) &&
(c2.getXmax()[i] <= c1.getXmax()[i]))
test &= true;
else if ((c2.getXmin()[i] <= c1.getXmin()[i]) &&
(c1.getXmin()[i] <= c2.getXmax()[i]))
test &= true;
else if ((c2.getXmin()[i] <= c1.getXmax()[i]) &&
(c1.getXmax()[i] <= c2.getXmax()[i]))
test &= true;
}
return test;
}
LM_FATAL("No implementation of intersection of cube with other geometries");
return false;
}
/* -------------------------------------------------------------------------- */
void Cube::init() {
for (UInt i = 0; i < 3; ++i) {
this->center[i] = getXmin()[i] + (getXmax()[i] - getXmin()[i]) / 2;
this->size[i] = getXmax()[i] - getXmin()[i];
this->size_half[i] = this->size[i] / 2.;
this->inv_size_half[i] = 1. / this->size_half[i];
xmin[i] = this->ranges[2 * i + 0];
xmax[i] = this->ranges[2 * i + 1];
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC CUBE
This geometry is the generalized description of a bounding box.
There is also a mechanism to generate holes
(as described by the HOLE keyword) inside of a plain bounding box.
*/
/* LMEXAMPLE
GEOMETRY segment CUBE BBOX -10 10 \\
GEOMETRY rectangle CUBE BBOX -10 10 -20 20 \\
GEOMETRY cube CUBE BBOX -10 10 -10 10 -10 10 \\
*/
void Cube::declareParams() {
/* LMKEYWORD BBOX
Give the bounding box coordinates. For example:\\
GEOMETRY full CUBE BBOX -L0/2 L0/2 -L1/2 L1/2 \\
\includegraphics[scale=.2]{images/geom-cube} \\
*/
this->parseVectorKeyword("BBOX", dim*2, ranges);
/* LMKEYWORD HOLE
In the bounding box description this provides the hole values
Below are provided a few self explanatory examples.
GEOMETRY 1d-hole CUBE BBOX -L0/2 L0/2 -L1/2 L1/2 HOLE L0/10 0 \\
\includegraphics[scale=.2]{images/geom-cube-hole1d} \\
GEOMETRY 2d-hole CUBE BBOX -L0/2 L0/2 -L1/2 L1/2 HOLE L0/10 L1/5 \\
\includegraphics[scale=.2]{images/geom-cube-hole2d} \\
*/
this->parseVectorKeyword("HOLE", dim, hole, VEC_DEFAULTS(0., 0., 0.));
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
diff --git a/src/geometry/cube.hh b/src/geometry/cube.hh
index bf791cb..97f6ff0 100644
--- a/src/geometry/cube.hh
+++ b/src/geometry/cube.hh
@@ -1,280 +1,292 @@
/**
* @file cube.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Nov 25 12:23:57 2013
*
* @brief Cubic geometry
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
#ifndef __LIBMULTISCALE_CUBE_HH__
#define __LIBMULTISCALE_CUBE_HH__
/* -------------------------------------------------------------------------- */
#include "geometry.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/**
* Class Cube
*
*/
class Cube : public virtual Geometry {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Cube(UInt dim = 3, GeomID id = "default geometry");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
virtual void declareParams();
virtual bool contains(Real my_x, Real my_y = 0, Real my_z = 0);
virtual void init();
/// function to print the contain of the class
virtual void printself(std::ostream &stream, UInt indent = 0) const {
Geometry::printself(stream, indent);
std::string tabu;
for (UInt i = 0; i < indent; ++i)
tabu += "\t";
stream << tabu << "CUBE geom : " << this->getXmin() << " "
<< this->getXmax() << " " << this->Hole(0) << " " << this->Hole(1)
<< " " << this->Hole(2) << std::endl;
};
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
const Quantity<Length, 3> &getXmax() const;
const Quantity<Length, 3> &getXmin() const;
+
+ template <UInt Dim>
+ void setXmax(Vector<Dim> xmax){
+ this->xmax = xmax;
+ this->init();
+ };
+
+ template <UInt Dim>
+ void setXmin(Vector<Dim> xmin) {
+ this->xmin = xmin;
+ this->init();
+ };
+
void setXmax(UInt i, Real r);
void setXmin(UInt i, Real r);
const Real &Hole(UInt i) const;
void setHole(Real *X);
void setDimensions(Real minx, Real maxx, Real miny, Real maxy, Real minz,
Real maxz);
void setDimensions(const Real *xmin, const Real *xmax);
void setDimensions(const Quantity<Length, 3> &xmin,
const Quantity<Length, 3> &xmax);
void resetDimensions();
void extendBoundingBox(Real x, Real y, Real z);
void extendBoundingBox(const Real *X);
void extendBoundingBox(Vector<1> &X);
void extendBoundingBox(Vector<2> &X);
void extendBoundingBox(Vector<3> &X);
void extendBoundingBox(VectorView<1> X);
void extendBoundingBox(VectorView<2> X);
void extendBoundingBox(VectorView<3> X);
Real *getSize();
Real getSize(UInt i);
Real *getSizeHalf();
Real *getInvSizeHalf();
Cube &getBoundingBox();
bool doIntersect(Geometry &geom);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
friend class GeomTools;
- Quantity<Length, 3 * 2> ranges;
+ Quantity<Length, 3*2> ranges;
+ Quantity<Length, 3> xmin;
+ Quantity<Length, 3> xmax;
Quantity<Length, 3> hole;
bool empty;
Real size[3];
Real size_half[3];
Real inv_size_half[3];
-
- Quantity<Length, 3> xmin;
- Quantity<Length, 3> xmax;
};
/* -------------------------------------------------------------------------- */
inline bool Cube::contains(Real my_x, Real my_y, Real my_z) {
rotate(my_x, my_y, my_z);
Real X1[3] = {my_x, my_y, my_z};
DUMP(X1[0] << " " << getXmax() << " " << getXmin(), DBG_ALL);
bool result = true;
for (UInt i = 0; i < dim; ++i) {
result &= (X1[i] <= getXmax()[i] && X1[i] >= getXmin()[i]);
}
if (result) {
bool result_hole[3] = {true, true, true};
for (UInt i = 0; i < dim; ++i) {
if (hole[i] != 0)
result_hole[i] =
!(X1[i] <= center[i] + hole[i] && X1[i] >= center[i] - hole[i]);
result &= result_hole[i];
}
}
for (UInt i = 0; i < dim; ++i) {
DUMP("result " << result << " " << X1[i] << " xmin " << getXmin()[i]
<< " xmax " << getXmax()[i],
DBG_ALL);
}
return result;
}
/* -------------------------------------------------------------------------- */
inline void Cube::extendBoundingBox(Vector<1> &x) {
Real X[3] = {x[0], 0., 0.};
extendBoundingBox(X);
}
inline void Cube::extendBoundingBox(Vector<2> &x) {
Real X[3] = {x[0], x[1], 0.};
extendBoundingBox(X);
}
inline void Cube::extendBoundingBox(Vector<3> &x) {
Real X[3] = {x[0], x[1], x[2]};
extendBoundingBox(X);
}
inline void Cube::extendBoundingBox(VectorView<1> x) {
Real X[3] = {x[0], 0., 0.};
extendBoundingBox(X);
}
inline void Cube::extendBoundingBox(VectorView<2> x) {
Real X[3] = {x[0], x[1], 0.};
extendBoundingBox(X);
}
inline void Cube::extendBoundingBox(VectorView<3> x) {
Real X[3] = {x[0], x[1], x[2]};
extendBoundingBox(X);
}
/* -------------------------------------------------------------------------- */
inline void Cube::extendBoundingBox(Real x, Real y, Real z) {
Real X[3] = {0., 0., 0.};
X[0] = x;
X[1] = y;
X[2] = z;
extendBoundingBox(X);
}
/* -------------------------------------------------------------------------- */
inline void Cube::extendBoundingBox(const Real *X) {
if (empty) {
setDimensions(X, X);
empty = false;
return;
}
Real new_xmin[3] = {0., 0., 0.};
Real new_xmax[3] = {0., 0., 0.};
for (UInt i = 0; i < dim; ++i) {
new_xmin[i] = std::min(getXmin()[i], X[i]);
new_xmax[i] = std::max(getXmax()[i], X[i]);
}
setDimensions(new_xmin, new_xmax);
}
/* -------------------------------------------------------------------------- */
inline void Cube::setDimensions(Real minx, Real maxx, Real miny, Real maxy,
Real minz, Real maxz) {
Real xmax[3];
Real xmin[3];
xmax[0] = maxx;
xmin[0] = minx;
xmax[1] = maxy;
xmin[1] = miny;
xmax[2] = maxz;
xmin[2] = minz;
setDimensions(xmin, xmax);
}
/* -------------------------------------------------------------------------- */
inline void Cube::setDimensions(const Real *xmin, const Real *xmax) {
for (UInt i = 0; i < this->dim; ++i) {
this->setXmax(i, xmax[i]);
this->setXmin(i, xmin[i]);
}
init();
}
/* -------------------------------------------------------------------------- */
inline void Cube::setDimensions(const Quantity<Length, 3> &xmin,
const Quantity<Length, 3> &xmax) {
for (UInt i = 0; i < this->dim; ++i) {
this->setXmax(i, xmax[i]);
this->setXmin(i, xmin[i]);
}
init();
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_CUBE_HH__ */
diff --git a/src/geometry/geometry.hh b/src/geometry/geometry.hh
index f4c42b1..37d17ef 100644
--- a/src/geometry/geometry.hh
+++ b/src/geometry/geometry.hh
@@ -1,252 +1,258 @@
/**
* @file geometry.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Jan 10 20:47:45 2014
*
* @brief Common mother of all geometries
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
#ifndef __LIBMULTISCALE_GEOMETRY_HH__
#define __LIBMULTISCALE_GEOMETRY_HH__
/* -------------------------------------------------------------------------- */
#include "lm_parsable.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class Cube;
/* -------------------------------------------------------------------------- */
/**
* Class Geometry
*
*/
class Geometry : public Parsable {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
enum GeomType {
BALL = 1,
CUBE = 2,
INTER = 3,
SUB = 4,
ELLIPSOID = 5,
UNION = 6,
CYLINDER = 7,
CHAUDRON = 8,
CUBE_SURFACE = 9
};
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Geometry(UInt Dim, GeomType Type, GeomID id);
virtual ~Geometry();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
//! return type of geometry
GeomType getType() const;
//! return dimention of geometry
UInt getDim() const;
//! return dimention of geometry
void setDim(UInt d);
//! set the center of the geometry
void setCenter(Real X, Real Y = 0, Real Z = 0);
//! set the rotation matrix from euler parameters
void setRotation(Real phi, Real theta, Real psi);
//! get the i-est coordinate of the center
Real getCenter(UInt i) const;
//! get the id of the geometry
GeomID getID() { return id; };
//! set the id of the geometry
void setID(GeomID &id) { this->id = id; };
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// function to print the contain of the class
virtual void printself(std::ostream &stream, UInt indent = 0) const;
//! init function
virtual void init() = 0;
//! function that return true if geometrie contaUInt poUInt x,y,z
virtual bool contains(Real x, Real y = 0, Real z = 0) = 0;
//! generic function that call the contains from a pointer
+ template <UInt Dim> inline bool contains(VectorProxy<Dim> &X) {
+ return this->contains<Dim>(Vector<Dim>(X));
+ }
template <UInt Dim> inline bool contains(Vector<Dim> &X) {
return this->contains<Dim>(X.data());
}
- template <UInt Dim> inline bool contains(VectorView<Dim> &X) {
+ template <UInt Dim> inline bool contains(VectorView<Dim> X) {
+ return this->contains<Dim>(X.data());
+ }
+ template <int Dim> inline bool contains(VectorView<Dim> X) {
return this->contains<Dim>(X.data());
}
template <UInt Dim> inline bool contains(Real *X) {
if (Dim == 1) {
return this->contains(X[0]);
} else if (Dim == 2) {
return this->contains(X[0], X[1]);
} else if (Dim == 3) {
return this->contains(X[0], X[1], X[2]);
} else {
LM_FATAL("Unknown dimension '" << Dim << "'");
return false;
}
}
//! apply rotation on a given poUInt to change coordinates system
void rotate(Real &x, Real &y, Real &z);
//! return the distance to the center
Real distToCenter(Real X, Real Y = 0, Real Z = 0);
//! return the distance to the center
Real distToCenter(Real *X);
//! return the unidimentional(X) vector to the center (see source code)
Real vecToCenter(UInt index, Real X);
//! return the normed vector coordinate to the center
Real normToCenter(UInt index, Real X, Real Y = 0.0, Real Z = 0.0);
//! return a bounding box of the current geometry
virtual Cube &getBoundingBox();
//! intersect with another geometry
virtual bool doIntersect(Geometry &geom) = 0;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
friend class GeomTools;
//! type of geom
GeomType type;
//! id of geom
GeomID id;
//! dimention
UInt dim;
//! center coordinates
Quantity<Length, 3> center;
//! center of rotation
Quantity<Length, 3> rotCenter;
//! rotation matrix by euler theorem with phi,teta,psi notation
/** see http://mathworld.wolfram.com/EulerAngles.html for more details */
Real rotMatrix[9];
};
/* -------------------------------------------------------------------------- */
inline void Geometry::rotate(Real &x, Real &y, Real &z) {
Real x2 = rotMatrix[0] * x + rotMatrix[3] * y + rotMatrix[6] * z;
Real y2 = rotMatrix[1] * x + rotMatrix[4] * y + rotMatrix[7] * z;
Real z2 = rotMatrix[2] * x + rotMatrix[5] * y + rotMatrix[8] * z;
x = x2;
y = y2;
z = z2;
}
/* -------------------------------------------------------------------------- */
inline Real Geometry::distToCenter(Real X, Real Y, Real Z) {
Real x[3] = {X, Y, Z};
return distToCenter(x);
}
/* -------------------------------------------------------------------------- */
inline Real Geometry::distToCenter(Real *X) {
Real res = 0.0;
for (UInt i = 0; i < 3; ++i) {
Real tmp = center[i] - X[i];
res += tmp * tmp;
}
return sqrt(res);
}
/* -------------------------------------------------------------------------- */
inline Real Geometry::vecToCenter(UInt index, Real X) {
Real res = X - center[index];
// LOG("vec au centre = " << res << " sur axe " << index);
return res;
}
/* -------------------------------------------------------------------------- */
inline Real Geometry::normToCenter(UInt index, Real X, Real Y, Real Z) {
Real res[3];
res[0] = X - center[0];
res[1] = Y - center[1];
res[2] = Z - center[2];
Real norm = sqrt(res[0] * res[0] + res[1] * res[1] + res[2] * res[2]);
if (!norm)
return 0.0;
DUMP("norm au centre = " << res[index] / norm << " sur axe " << index,
DBG_ALL);
return res[index] / norm;
}
/* -------------------------------------------------------------------------- */
/// standard output stream operator
inline std::ostream &operator<<(std::ostream &stream, Geometry &_this) {
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_GEOMETRY_HH__ */
diff --git a/src/md/lammps/domain_lammps.cc b/src/md/lammps/domain_lammps.cc
index 7946699..9274c95 100644
--- a/src/md/lammps/domain_lammps.cc
+++ b/src/md/lammps/domain_lammps.cc
@@ -1,1011 +1,1005 @@
/**
* @file domain_lammps.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Thu Jul 31 22:41:23 2014
*
* @brief This is the model wrapping LAMMPS
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
/* -------------------------------------------------------------------------- */
//#define CHECK_STABILITY
/* -------------------------------------------------------------------------- */
#include "domain_lammps.hh"
#include "filter_manager.hh"
#include "import_lammps.hh"
#include "lammps_common.hh"
#include "lm_common.hh"
#include "reference_manager_lammps.hh"
#include "trace_atom.hh"
#include <iomanip>
#include <iostream>
#include <iterator>
/* -------------------------------------------------------------------------- */
// LAMMPS includes
#include "angle.h"
#include "bond.h"
#include "change_box.h"
#include "dihedral.h"
#include "displace_box.h"
#include "fix_minimize.h"
#include "improper.h"
#include "integrate.h"
#include "kspace.h"
#include "min_cg.h"
#include "pair.h"
#include "variable.h"
#include "verlet.h"
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
#include "compute_lammps.hh"
#include "geometry_manager.hh"
#include "lm_parser.hh"
/* -------------------------------------------------------------------------- */
LAMMPS_NS::LAMMPS *lammps_main_object = NULL;
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <UInt Dim> void DomainLammps<Dim>::init() { this->initModel(); }
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::loadLammpsConfigFile() {
// forward LM variables to lammps
std::map<std::string, double> &vars = Parser::getAlgebraicVariables();
std::map<std::string, double>::iterator it = vars.begin();
std::map<std::string, double>::iterator end = vars.end();
while (it != end) {
std::string varname = it->first;
std::stringstream varval;
varval << std::setprecision(16) << std::scientific << it->second;
LAMMPS_NS::Variable *lmpvar =
static_cast<LAMMPS_NS::LAMMPS *>(this)->input->variable;
char *cvarname = const_cast<char *>(varname.c_str());
char cvarval[30];
varval.str().copy(cvarval, 29);
cvarval[varval.str().length()] = 0;
if (lmpvar->find(cvarname) >= 0)
LM_FATAL("cannot export LibMultiscale variable "
<< varname << " to LAMMPS since variable already exists");
DUMP(cvarname << " = " << cvarval, DBG_INFO);
lmpvar->add_equal(cvarname, cvarval);
++it;
}
std::string file_to_read = lammpsfile;
if (create_header_flag &&
Communicator::getCommunicator().groupRank(lm_my_proc_id,
this->getGroupID()) == 0) {
std::string unit_name;
if (code_unit_system == metal_unit_system)
unit_name = "metal";
else if (code_unit_system == atomic_unit_system)
unit_name = "real";
else if (code_unit_system == real_unit_system)
unit_name = "si";
else
LM_FATAL("unknown units");
file_to_read = "./generated_lammps.config";
std::ofstream file_lammps(file_to_read.c_str());
file_lammps << "units " << unit_name << std::endl;
file_lammps << "dimension " << Dim << std::endl;
file_lammps << "atom_style atomic" << std::endl;
file_lammps << "boundary ";
for (UInt i = 0; i < 3; ++i) {
file_lammps << boundaries[i] << " ";
}
file_lammps << std::endl;
file_lammps << "lattice " << lattice << " " << std::setprecision(15)
<< static_cast<Real>(lattice_size);
file_lammps << " orient x ";
for (UInt i = 0; i < 3; ++i)
file_lammps << lattice_orient_x[i] << " ";
if (Dim > 1) {
file_lammps << "orient y ";
for (UInt i = 0; i < 3; ++i)
file_lammps << lattice_orient_y[i] << " ";
}
if (Dim == 3) {
file_lammps << "orient z ";
for (UInt i = 0; i < 3; ++i)
file_lammps << lattice_orient_z[i] << " ";
}
file_lammps << "spacing ";
for (UInt i = 0; i < 3; ++i)
file_lammps << lattice_spacing[i] << " ";
file_lammps << "origin ";
for (UInt i = 0; i < Dim; ++i)
file_lammps << lattice_origin[i] << " ";
for (UInt i = Dim; i < 3; ++i)
file_lammps << 0.0 << " ";
file_lammps << std::endl;
file_lammps << "region box block ";
for (UInt i = 0; i < 6; ++i)
file_lammps << replica[i] << " ";
file_lammps << std::endl;
file_lammps << "create_box 1 box" << std::endl;
file_lammps << "create_atoms 1 box" << std::endl;
file_lammps << std::endl << std::endl;
if (triclinic == 1 && floorf(tilt[1] * 100 + 0.5) / 100 != 0.0) {
// double value;
std::string direction;
Real tilting = floorf(tilt[1] * 100 + 0.5) / 100;
if (tilt[0] == 0.0) {
// value = tilt[1]/lattice_size/lattice_spacing[0];
direction = "xy";
} else if (tilt[0] == 1.0) {
// value = tilt[1]/lattice_size/lattice_spacing[0];
direction = "xz";
} else if (tilt[0] == 2.0) {
// value = tilt[1]/lattice_size/lattice_spacing[0];
direction = "yz";
} else
LM_FATAL("tilt must be happend in xy(0), xz(1) or yz(2) axis!!!");
file_lammps << "change_box triclinic" << std::endl;
file_lammps << "displace_box all " << direction << " final " << tilting
<< " remap none units box" << std::endl;
// file_lammps << "displace_box all "<<direction<<" final "<<value<<"
// remap none"<< std::endl;
}
std::ifstream in(lammpsfile.c_str());
std::istreambuf_iterator<char> src(in.rdbuf());
std::istreambuf_iterator<char> end;
std::ostream_iterator<char> dest(file_lammps);
std::copy(src, end, dest);
}
DUMP("opening lammps config file " << file_to_read, DBG_INFO_STARTUP);
infile = fopen(file_to_read.c_str(), "rb");
if (infile == NULL)
LM_FATAL("cannot open " << file_to_read << " as a lammps config file !!!");
input->file();
DUMP("lammps file loaded", DBG_INFO);
// register lammps computes to libmultiscale system
UInt ncompute = this->modify->ncompute;
for (UInt i = 0; i < ncompute; ++i) {
ComputeLammps<LAMMPS_NS::Compute> *_ptr =
new ComputeLammps<LAMMPS_NS::Compute>(*(this->modify->compute[i]),
*this);
_ptr->setRelease(INITIAL_MODEL_RELEASE);
_ptr->setCommGroup(this->getGroupID());
FilterManager::getManager().addObject(_ptr);
}
// register lammps computes to libmultiscale system
UInt nfix = this->modify->nfix;
for (UInt i = 0; i < nfix; ++i) {
ComputeLammps<LAMMPS_NS::Fix> *_ptr =
new ComputeLammps<LAMMPS_NS::Fix>(*(this->modify->fix[i]), *this);
_ptr->setRelease(INITIAL_MODEL_RELEASE);
_ptr->setCommGroup(this->getGroupID());
FilterManager::getManager().addObject(_ptr);
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::performStep1() {
if (this->newGeomBox != invalidGeom &&
current_step >= this->when_change_box) {
this->changeBox();
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::changeBox() {
Cube &c = GeometryManager::getManager()
.getGeometry(this->newGeomBox)
->getBoundingBox();
this->setDomainDimensions(c.getXmin().getValues(), c.getXmax().getValues());
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::initModel() {
lammps_main_object = this;
ReferenceManagerLammps<Dim> *refMgr =
new ReferenceManagerLammps<Dim>(this->atoms);
this->atoms.setRefManager(refMgr);
refMgr->setMPIComm(
Communicator::getCommunicator().getMpiGroup(this->getGroupID()));
ImportLammpsInterface::createImportLammps(*refMgr);
if (this->getGeom().getDim() == 1)
LM_FATAL("LAMMPS cannot work in 1D");
ImportLammpsInterface::setGeom(&this->getGeom());
std::map<UInt, GeomID>::iterator it = geom_by_type.begin();
std::map<UInt, GeomID>::iterator end = geom_by_type.end();
while (it != end) {
Geometry *g = GeometryManager::getManager().getObject(it->second);
UInt itype = it->first;
ImportLammpsInterface::setGeom(g, itype);
++it;
}
loadLammpsConfigFile();
this->atoms.init();
update->whichflag = 0;
update->firststep = 0;
update->laststep = nb_step;
update->ntimestep = 0;
update->beginstep = 0;
update->endstep = update->laststep;
if (flag_units) {
force->ftm2v = code_unit_system.ft_m2v;
// force->nktv2p = 1;
force->mvv2e = code_unit_system.mvv2e;
}
if (this->domain->dimension != Dim)
LM_FATAL("mismatch dimension: check lammps config file");
if (std::string(this->atom->atom_style) != "atomic")
LM_FATAL("can only work with atomic style");
this->update_pad = false;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::initDOFs() {
{
UInt nall;
if (force->newton)
nall = atom->nlocal + atom->nghost;
else
nall = atom->nlocal;
Real **x = atom->x;
Real **x0 = atom->x0;
for (UInt i = 0; i < nall; i++) {
x0[i][0] = x[i][0];
x0[i][1] = x[i][1];
x0[i][2] = x[i][2];
}
}
if (this->shouldRestart()) {
this->restart_start = true;
this->readXMLFile(this->restartfile);
if (this->newGeomBox != invalidGeom &&
current_step >= this->when_change_box) {
this->changeBox();
this->performMigration(); // perform migration after change box to keep
// periodicity
this->when_change_box = current_step + nb_step; // to avoid repetitions
} else if (neighbor->decide()) {
this->performMigration();
this->resetBox();
}
this->performStep2();
LAMMPS_NS::LAMMPS::output->setup(false);
}
// conversion of masses from g/mol to Kg
if (flag_units) {
UnitsConverter uc;
uc.setOutUnits(code_unit_system);
if (strcmp(update->unit_style, "real") == 0)
uc.setInUnits(atomic_unit_system);
else if (strcmp(update->unit_style, "metal") == 0)
uc.setInUnits(metal_unit_system);
else
LM_FATAL("not known lammps unit system:"
<< " things need to be done to integrate");
uc.computeConversions();
for (int i = 1; i <= atom->ntypes; i++)
atom->mass[i] *= uc.getConversion<Mass>();
}
#ifdef TRACE_ATOM
for (int i = 0; i < atom->nlocal; ++i) {
Real X[3] = {atom->x0[i][0], atom->x0[i][1], atom->x0[i][2]};
SET_INTERNAL_TRACE_INDEX(X, i);
IF_TRACED(X, "traced atom has been detected at position "
<< internal_index);
}
VIEW_ATOM(RefLammps<Dim>);
#endif
this->atoms.setRelease(INITIAL_MODEL_RELEASE);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getEcin() { return 0.0; }
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getEpot() {
return update->minimize->pe_compute->compute_scalar();
}
/* -------------------------------------------------------------------------- */
extern "C" {
static UInt fargc = 1;
static char prog[20] = "lammps";
static char *___tmp = &prog[0];
static char **fargv = &___tmp;
//#include "../../common/openmp.hh"
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> DomainLammps<Dim>::~DomainLammps() {
if (this->atoms.hasRefManager())
delete (&this->atoms.getRefManager());
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
DomainLammps<Dim>::DomainLammps(DomainID ID, CommGroup GID)
: DomainAtomic<ContainerLammps<Dim>, Dim>(ID, GID),
LAMMPS(fargc, fargv, Communicator::getCommunicator().getMpiGroup(GID)),
flag_units(true), triclinic(0) {
// external_work = 0;
newGeomBox = invalidGeom;
when_change_box = 0;
lattice = "";
lattice_size = 0.;
for (UInt i = 0; i < 6; ++i)
replica[i] = 0;
create_header_flag = false;
restart_start = false;
for (UInt i = 0; i < 3; ++i) {
lattice_orient_x[i] = 0;
lattice_orient_y[i] = 0;
lattice_orient_z[i] = 0;
lattice_origin[i] = 0.;
}
lattice_orient_x[0] = 1;
lattice_orient_y[1] = 1;
lattice_orient_z[2] = 1;
lattice_spacing[0] = 1.0;
lattice_spacing[1] = 1.0;
lattice_spacing[2] = 1.0;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainLammps<Dim>::getSubDomainDimensions(Real *xmin, Real *xmax) {
for (UInt i = 0; i < 3; ++i) {
xmin[i] = domain->sublo[i];
xmax[i] = domain->subhi[i];
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainLammps<Dim>::getDomainDimensions(Real *xmin, Real *xmax) {
for (UInt i = 0; i < 3; ++i) {
xmin[i] = domain->boxlo[i];
xmax[i] = domain->boxhi[i];
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainLammps<Dim>::setDomainDimensions(const Real *xmin,
const Real *xmax) {
for (UInt i = 0; i < 3; ++i) {
domain->boxlo[i] = xmin[i];
domain->boxhi[i] = xmax[i];
}
if (neighbor->decide())
this->resetBox();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::resetBox() {
atom->setup();
if (domain->triclinic)
domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style)
neighbor->setup_bins();
comm->exchange();
if (atom->sortfreq > 0)
atom->sort();
comm->borders();
if (domain->triclinic)
domain->lamda2x(atom->nlocal + atom->nghost);
neighbor->build();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::performMigration() {
if (Communicator::getCommunicator().amIinGroup(this->getGroupID())) {
echangeAtomes();
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::echangeAtomes() {
VIEW_ATOM(RefLammps<Dim>);
UInt nflag = neighbor->decide();
// if (nflag == 0 && this->restart_start==false && this->update_pad == false)
// {
if (nflag == 0 && this->restart_start == false) {
timer->stamp();
comm->forward_comm();
timer->stamp(TIME_COMM);
} else {
if (modify->n_pre_exchange)
modify->pre_exchange();
if (domain->triclinic)
domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style)
neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
comm->borders();
if (domain->triclinic)
domain->lamda2x(atom->nlocal + atom->nghost);
timer->stamp(TIME_COMM);
if (modify->n_pre_neighbor)
modify->pre_neighbor();
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
this->restart_start = false;
this->update_pad = false;
}
if (this->when_change_box == current_time) {
// if one change box size, this part has to be fulfilled.
if (modify->n_pre_exchange)
modify->pre_exchange();
if (domain->triclinic)
domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style)
neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
comm->borders();
if (domain->triclinic)
domain->lamda2x(atom->nlocal + atom->nghost);
timer->stamp(TIME_COMM);
if (modify->n_pre_neighbor)
modify->pre_neighbor();
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
this->update_pad = false;
this->restart_start = false;
}
VIEW_ATOM(RefLammps<Dim>);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getFdotDir() {
// Real *f = atom->f[0];
Real fdotdirme = 0.0;
UInt n = update->minimize->nvec;
for (UInt i = 0; i < n; i++) {
// fdotdirme += f[i]*update->minimize->h[i];
}
Real fdotdirall;
MPI_Allreduce(&fdotdirme, &fdotdirall, 1, MPI_DOUBLE, MPI_SUM,
MPI_COMM_WORLD);
return fdotdirall;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getFMax() {
Real fmax = 0;
Real fmaxall;
Real *f = atom->f[0];
UInt n = update->minimize->nvec;
for (UInt i = 0; i < n; i++)
fmax = std::max(fabs(f[i]), fmax);
MPI_Allreduce(&fmax, &fmaxall, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
return fmaxall;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getDirMax() {
Real fmax = 0;
Real fmaxall;
// UInt n = update->minimize->nvec;
// for (UInt i = 0; i < n; i++) fmax = max(fabs(update->minimize->h[i]),fmax);
MPI_Allreduce(&fmax, &fmaxall, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
return fmaxall;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getFNorm2() {
Real fnorm2all = 0.0;
Real *f = atom->f[0];
Real fnorm2me = 0.0;
UInt n = update->minimize->nvec;
for (UInt i = 0; i < n; i++)
fnorm2me += f[i] * f[i];
MPI_Allreduce(&fnorm2me, &fnorm2all, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
return fnorm2all;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getFdotOldF() {
Real fdotoldf;
// Real *f = atom->f[0];
Real fdotoldfme = 0.0;
// UInt n = update->minimize->nvec;
// for (UInt i = 0; i < n; i++) fdotoldfme += f[i]*update->minimize->g[i];
MPI_Allreduce(&fdotoldfme, &fdotoldf, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
return fdotoldf;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::updateDirection(Real beta) {
for (int i = 0; i < update->minimize->nvec; i++) {
// update->minimize->h[i] = update->minimize->g[i] +
// beta*update->minimize->h[i];
// max_h = MAX(max_h,fabs(min_ptr->h[i]));
}
auto *geom_ptr = &this->getGeomConstrained();
if (geom_ptr) {
for (int j = 0; j < update->minimize->nvec / 3; ++j) {
Real x = atom->x0[j][0];
Real y = atom->x0[j][1];
Real z = atom->x0[j][2];
if (this->getGeomConstrained().contains(x, y, z)) {
// update->minimize->h[3*j] = 0;
// update->minimize->h[3*j+1] = 0;
// update->minimize->h[3*j+2] = 0;
}
// double norm =
// min_ptr->h[3*j]*min_ptr->h[3*j]
// + min_ptr->h[3*j+1]*min_ptr->h[3*j+1]
// + min_ptr->h[3*j+2]*min_ptr->h[3*j+2];
// norm = sqrt(norm);
// double factor = exp (-max_h/100/norm );
// if (norm < max_h/100) {
// min_ptr->h[3*j] *= factor;
// min_ptr->h[3*j+1] *= factor;
// min_ptr->h[3*j+2] *= factor;
// }
// if (norm < max_h/10) {
// min_ptr->h[3*j] = 0;
// min_ptr->h[3*j+1] = 0;
// min_ptr->h[3*j+2] = 0;
// }
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::saveForceVector() {
// Real *f = atom->f[0];
for (int i = 0; i < update->minimize->nvec; i++) {
// update->minimize->g[i] = f[i];
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainLammps<Dim>::displaceTowardsDirection(Real delta) {
// Real * x = atom->x[0];
// for (int i = 0; i < update->minimize->nvec; i++) {
// // std::cerr << "x[" << i << "/" << update->minimize->ndof << "](" <<
// x[i] << ") += " << delta << "*" << update->minimize->h[i] << std::endl;
// x[i] += delta*update->minimize->h[i];
// }
// double local_external_work = 0;
// double * v = lammps_main_object->atom->v[0];
// for (int i = 0 ; i < min_ptr->ndof ; ++i){
// local_external_work += delta*min_ptr->h[i]*v[i];
// }
// double global_external_work = 0;
// MPI_Allreduce(&local_external_work,&global_external_work,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
// external_work += global_external_work;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getExternalWork() {
// return external_work;
return 0;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::setExternalWork(Real work) {
// external_work = work;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> bool DomainLammps<Dim>::isPeriodic(UInt dir) {
switch (dir) {
case 0:
return this->domain->xperiodic;
break;
case 1:
return this->domain->yperiodic;
break;
case 2:
return this->domain->zperiodic;
break;
default:
LM_FATAL("called on non-dimension");
return false;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::force_clear(UInt vflag) {
UInt i;
// clear global force array
// nall includes ghosts only if either newton flag is set
UInt nall;
if (force->newton)
nall = atom->nlocal + atom->nghost;
else
nall = atom->nlocal;
Real **f = atom->f;
for (i = 0; i < nall; i++) {
f[i][0] = 0.0;
f[i][1] = 0.0;
f[i][2] = 0.0;
}
// if (integ.torqueflag) {
// Real **torque = atom->torque;
// for (i = 0; i < nall; i++) {
// torque[i][0] = 0.0;
// torque[i][1] = 0.0;
// torque[i][2] = 0.0;
// }
// }
// if (granflag) {
// Real **phia = atom->phia;
// for (i = 0; i < nall; i++) {
// phia[i][0] = 0.0;
// phia[i][1] = 0.0;
// phia[i][2] = 0.0;
// }
// }
// clear pair force array if using it this timestep to compute virial
// if (vflag == 2 && pairflag) {
// if (atom->nmax > maxpair) {
// maxpair = atom->nmax;
// memory->destroy_2d_Real_array(f_pair);
// f_pair = memory->create_2d_Real_array(maxpair,3,"verlet:f_pair");
// }
// for (i = 0; i < nall; i++) {
// f_pair[i][0] = 0.0;
// f_pair[i][1] = 0.0;
// f_pair[i][2] = 0.0;
// }
// }
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::setCurrentStep(UInt ts) {
current_step = ts;
// content of the Update::reset_timestep function
update->eflag_global = update->vflag_global = ts;
for (int i = 0; i < modify->ncompute; i++) {
modify->compute[i]->invoked_scalar = -1;
modify->compute[i]->invoked_vector = -1;
modify->compute[i]->invoked_array = -1;
modify->compute[i]->invoked_peratom = -1;
modify->compute[i]->invoked_local = -1;
// modify->compute[i]->eflag_global = -1;
}
update->ntimestep = ts;
if (update->laststep < ts)
update->laststep += ts;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::setTimeStep(Real ts) {
update->dt = ts;
for (int i = 0; i < modify->nfix; ++i)
modify->fix[i]->init();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> UInt DomainLammps<Dim>::getCurrentStep() {
return update->ntimestep;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getTimeStep() { return update->dt; }
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::updatePAD() {
// if (modify->n_pre_exchange)
modify->pre_exchange();
if (domain->triclinic)
domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
neighbor->setup_bins();
// if (domain->box_change) {
// domain->reset_box();
// comm->setup();
// if (neighbor->style)
// neighbor->setup_bins();
// }
timer->stamp();
comm->exchange();
comm->borders();
if (domain->triclinic)
domain->lamda2x(atom->nlocal + atom->nghost);
// if (modify->n_pre_neighbor)
modify->pre_neighbor();
neighbor->build();
// this->performStep2();
this->update_pad = false;
this->restart_start = false;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
int DomainLammps<Dim>::getOrientValue(UInt axis, UInt comp) {
if (axis == 0) {
return this->lattice_orient_x[comp];
} else if (axis == 1) {
return this->lattice_orient_y[comp];
} else if (axis == 2) {
return this->lattice_orient_z[comp];
} else {
LM_FATAL("wrong axis value: " << axis);
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> Real DomainLammps<Dim>::getLatticeConstant() {
return this->lattice_size;
}
/* -------------------------------------------------------------------------- */
/* LMDESC LAMMPS_BASE
This plugin make the interface with lammps
*/
/* LMHERITANCE domain_md */
/* LM EXAMPLE
Section MultiScale AtomsUnits\\
...\\
GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\
MODEL LAMMPS md\\
...\\
endSection\\ \\
Section LAMMPS:md AtomsUnits \\
LAMMPS_FILE lammps.in \\
DOMAIN_GEOMETRY 1 \\
endSection
*/
template <UInt Dim> void DomainLammps<Dim>::declareParams() {
DomainAtomic<ContainerLammps<Dim>, Dim>::declareParams();
/*LMKEYWORD LAMMPS_FILE
This specifies the configuration file to be used for lammps.
If the create header keyword is used then only the potential
definition should be into that file. Otherwise, a complete
lammps configuration file can be specified.
*/
this->parseKeyword("LAMMPS_FILE", lammpsfile);
/* LMKEYWORD CREATE_HEADER
A header is automatically generated based on the information
passed by keywords LATTICE, LATTICE\_SIZE, BOUNDARY and REPLICA
*/
this->parseTag("CREATE_HEADER", create_header_flag, false);
/* LMKEYWORD LATTICE
Specifies the lattice to be used to generate the lammps header.
Admissible values are lammps keywords
(see \url{http://lammps.sandia.gov/doc/lattice.html})
*/
this->parseKeyword("LATTICE", lattice, "");
/* LMKEYWORD LATTICE_SIZE
Lattice length parameter to generate the primitive lattice
*/
this->parseKeyword("LATTICE_SIZE", lattice_size, 0.);
/* LMKEYWORD LATTICE_ORIGIN
Lattice length parameter to generate the primitive lattice
*/
- this->parseVectorKeyword("LATTICE_ORIGIN", 3, lattice_origin,
- VEC_DEFAULTS(0., 0., 0.));
+ this->parseVectorKeyword("LATTICE_ORIGIN", 3, lattice_origin, VEC_DEFAULTS(0., 0., 0.));
/* LMKEYWORD LATTICE_ORIENTX
Lattice length parameter to generate the primitive lattice
*/
- this->parseVectorKeyword("LATTICE_ORIENTX", 3, lattice_orient_x,
- VEC_DEFAULTS(1, 0, 0));
+ this->parseVectorKeyword("LATTICE_ORIENTX", 3, lattice_orient_x, VEC_DEFAULTS(1, 0, 0));
/* LMKEYWORD LATTICE_ORIENTY
Lattice length parameter to generate the primitive lattice
*/
- this->parseVectorKeyword("LATTICE_ORIENTY", 3, lattice_orient_y,
- VEC_DEFAULTS(0, 1, 0));
+ this->parseVectorKeyword("LATTICE_ORIENTY", 3, lattice_orient_y, VEC_DEFAULTS(0, 1, 0));
/* LMKEYWORD LATTICE_ORIENTZ
Lattice length parameter to generate the primitive lattice
*/
- this->parseVectorKeyword("LATTICE_ORIENTZ", 3, lattice_orient_z,
- VEC_DEFAULTS(0, 0, 1));
+ this->parseVectorKeyword("LATTICE_ORIENTZ", 3, lattice_orient_z, VEC_DEFAULTS(0, 0, 1));
/* LMKEYWORD SPACING
The minimal normalized lattice sizes
*/
- this->parseVectorKeyword("SPACING", 3, lattice_spacing,
- VEC_DEFAULTS(1., 1., 1.));
+ this->parseVectorKeyword("SPACING", 3, lattice_spacing, VEC_DEFAULTS(1., 1., 1.));
/* LMKEYWORD BOUNDARY
Sequence of lammps code for boundary.
(see \url{http://lammps.sandia.gov/doc/boundary.html})
*/
this->parseVectorKeyword("BOUNDARY", 3, boundaries, VEC_DEFAULTS("", "", ""));
/* LMKEYWORD REPLICA
For the creation of atoms a region, specified as the number of
replicas of the lattice must be provided, before calling the
create atoms command.
*/
- this->parseVectorKeyword("REPLICA", 6, replica,
- VEC_DEFAULTS(0, 0, 0, 0, 0, 0));
+ this->parseVectorKeyword("REPLICA", 6, replica, VEC_DEFAULTS(0, 0, 0, 0, 0, 0));
/* LMKEYWORD CHANGE_DOMAIN_BOX
This allows to specify live reshaping of the bounding box.
Useful to enforce stresses states or dislocation displacement
field with periodic boundary conditions.
*/
this->parseKeyword("CHANGE_DOMAIN_BOX", newGeomBox, invalidGeom);
/* LMKEYWORD CHANGE_DOMAIN_BOX_TIME
Gives a timestep where the reshaping provided by CHANGE\_DOMAIN\_BOX
should occur.
*/
this->parseKeyword("CHANGE_DOMAIN_BOX_TIME", when_change_box, 0u);
/* LMKEYWORD TRICLINIC
This changes the orthogonal simulation box to triclinic box
*/
this->parseKeyword("TRICLINIC", triclinic, 0);
/* LMKEYWORD TILT
This tils the simulation box by given amount value
0: xy-direction; 1: xz-direction; 2: yz-direction
(TILT 0 0.25: tilt simulation box in xy-direction by 0.25)
*/
this->parseVectorKeyword("TILT", 2, tilt, VEC_DEFAULTS(0., 0.));
/* LMKEYWORD UNITS_CONVERSION
This is for debug purpose: force libmultiscale not to touch the
units conversion factors of lammps.
*/
this->parseTag("UNITS_CONVERSION", flag_units, true);
/* LMKEYWORD GEOM_BY_TYPE
Assign a geometry for the atom type provided
*/
this->parseKeyword("GEOM_BY_TYPE", geom_by_type);
/* LMKEYWORD MINCG_ETOL
obsolete
*/
this->parseKeyword("MINCG_ETOL", mincg_etol, 1e-8);
/* LMKEYWORD MINCG_FTOL
osolete
*/
this->parseKeyword("MINCG_FTOL", mincg_ftol, 1e-8);
/* LMKEYWORD MINCG_MAXEVAL
obsolete
*/
this->parseKeyword("MINCG_MAXEVAL", mincg_maxeval, 1000000u);
// // /* LMK EYWORD MINCG_DMAX
// // obsolete
// // */
// this->parseKeyword("MINCG_DMAX",mincg_dmax,.1);
}
/* -------------------------------------------------------------------------- */
template class DomainLammps<2>;
template class DomainLammps<3>;
__END_LIBMULTISCALE__
diff --git a/src/units/quantity.hh b/src/units/quantity.hh
index 7463106..a0f4640 100644
--- a/src/units/quantity.hh
+++ b/src/units/quantity.hh
@@ -1,203 +1,221 @@
/**
* @file quantity.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Nov 25 12:23:57 2013
*
* @brief This is a templated wrapper over float/doubles enabling units
* conversions
*
* @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/>.
*
*/
#ifndef __LIBMULTISCALE_QUANTITY_HH__
#define __LIBMULTISCALE_QUANTITY_HH__
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <PhysicalQuantity q, UInt nb = 1, typename T = Real,
bool wrapped = false>
class QuantityBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
QuantityBase() {
for (UInt i = 0; i < nb; ++i)
this->value[i] = T(0.0);
}
QuantityBase(T value) {
for (UInt i = 0; i < nb; ++i)
this->value[i] = value;
}
template <UInt n> QuantityBase(const QuantityBase<q, n, T> &_q) {
for (UInt i = 0; i < nb; ++i)
this->value[i] = _q[i];
}
+ void operator-=(Real val){
+ for (UInt i = 0; i < nb; ++i)
+ this->value[i] -= val;
+ }
+ void operator+=(Real val){
+ for (UInt i = 0; i < nb; ++i)
+ this->value[i] += val;
+ }
+ template <typename T2>
+ Vector<nb> operator-(T2 v) const{
+ return this->value - v;
+ }
+ Vector<nb> operator-(Real val) const{
+ Vector<nb> res(this->value);
+ for (UInt i = 0; i < nb; ++i)
+ res[i] -= val;
+ return res;
+ }
virtual T &operator[](UInt i) = 0;
virtual const T &operator[](UInt i) const = 0;
protected:
Vector<nb> value;
};
/* -------------------------------------------------------------------------- */
template <PhysicalQuantity q, UInt nb = 1, typename T = Real,
bool wrapped = false>
class Quantity : public QuantityBase<q, nb, T, wrapped> {
public:
Quantity() : QuantityBase<q, nb, T, wrapped>(){};
Quantity(T value) : QuantityBase<q, nb, T, wrapped>(value){};
template <UInt n>
Quantity(const QuantityBase<q, n, T, wrapped> &_q)
: QuantityBase<q, nb, T>(_q){};
template <typename T1> T &operator=(T1 v) {
for (UInt i = 0; i < nb; ++i)
this->value[i] = v;
return (*this);
}
T operator*(T v) {
Quantity<q, nb, T> res;
for (UInt i = 0; i < nb; ++i)
res.value[i] = this->value[i] * v;
return res;
}
-
+
T &operator[](UInt i) { return this->value[i]; }
const T &operator[](UInt i) const {
if (i >= nb)
LM_FATAL("internal error");
return this->value[i];
}
// operator const T * () const {
// return value;
// }
const T *getValues() const { return this->value.data(); }
// template <PhysicalQuantity q2>
// T operator / (Quantity<q2,1,T> v) {
// if (nb > 1)
// LM_FATAL("this quantity is a vector: "
// << "cannot implicitely cast to scalar");
// return value[0]/v[0];
// }
T operator/(T v) {
Quantity<q, nb, T> res;
for (UInt i = 0; i < nb; ++i)
res.value[i] = this->value[i] / v;
return res;
}
operator T &() {
if (nb > 1)
LM_FATAL("this quantity is a vector: "
<< "cannot implicitely cast to scalar");
return this->value[0];
}
operator Vector<nb> &() { return this->value; }
template <typename T1, PhysicalQuantity q1, UInt nb2, typename T2>
friend T1 operator/(T1 &a, Quantity<q1, nb2, T2> &b);
};
/* -------------------------------------------------------------------------- */
template <PhysicalQuantity q1, PhysicalQuantity q2, typename T>
inline T operator*(Quantity<q1, 1, T> &a, Quantity<q2, 1, T> &b) {
return b[0] * a[0];
}
template <typename T1, PhysicalQuantity q, typename T>
inline T1 operator*(T1 &a, Quantity<q, 1, T> &b) {
return b * a;
}
template <typename T1, PhysicalQuantity q, typename T>
inline T1 operator/(T1 &a, Quantity<q, 1, T> &b) {
return a / b[0];
}
inline std::ostream &operator<<(std::ostream &os, PhysicalQuantity q) {
switch (q) {
case Length:
os << "[Length]";
break;
case Mass:
os << "[Mass]";
break;
case Energy:
os << "[Energy]";
break;
case Time:
os << "[Time]";
break;
case MassDensity:
os << "[MassDensity]";
break;
case Force:
os << "[Force]";
break;
case Pressure:
os << "[Pressure]";
break;
case Temperature:
os << "[Temperature]";
break;
default:
LM_FATAL("unknown physical quantity " << q);
}
return os;
}
template <PhysicalQuantity q, UInt nb, typename T>
inline std::ostream &operator<<(std::ostream &os,
const Quantity<q, nb, T> &value) {
for (UInt i = 0; i < nb; ++i) {
const T tmp = value[i];
os << tmp << " ";
}
os << q;
return os;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_QUANTITY_HH__ */
diff --git a/test b/test
index 734570b..553bb7a 160000
--- a/test
+++ b/test
@@ -1 +1 @@
-Subproject commit 734570b66bb3961390e00eff1538ca9f9d39bd8f
+Subproject commit 553bb7ac8e721663a6cd8da21b05b4c4ccfa4c50