diff --git a/packages/paradis_plugin.cmake b/packages/paradis_plugin.cmake
index 1c6c090..6d79b2b 100644
--- a/packages/paradis_plugin.cmake
+++ b/packages/paradis_plugin.cmake
@@ -1,59 +1,61 @@
#===============================================================================
# @file paradis_plugin.cmake
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Till Junge <till.junge@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
# @author Moseley Philip Arthur <philip.moseley@epfl.ch>
#
# @date Fri Jul 11 15:47:44 2014
#
# @brief ParaDiS 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(paradis_plugin
DESCRIPTION "Plugin ParaDiS for LM"
SYSTEM OFF
DEFAULT OFF
DEPENDS PARADIS
)
set(LIBMULTISCALE_DD_HDRS_LIST_PARADIS
dd/paradis/domain_paradis.hh
+ dd/paradis/paradis_home.hh
)
set(LIBMULTISCALE_PARADIS_FILES
${LIBMULTISCALE_DD_HDRS_LIST_PARADIS}
dd/paradis/domain_paradis.cc
dd/paradis/iterator_paradis.hh
dd/paradis/container_paradis.hh
dd/paradis/ref_node_paradis.hh
dd/paradis/ref_elem_paradis.hh
+ dd/paradis/paradis_home.cc
)
set(LIBMULTISCALE_DD_MODEL_LIST_PARADIS
"((DomainParaDiS,3,PARADIS))"
)
package_declare_sources(paradis_plugin ${LIBMULTISCALE_PARADIS_FILES})
set_module_headers_from_package(dd paradis_plugin ${LIBMULTISCALE_DD_HDRS_LIST_PARADIS})
declare_boost_list_from_package(dd_model paradis_plugin ${LIBMULTISCALE_DD_MODEL_LIST_PARADIS})
diff --git a/src/common/container_mesh.hh b/src/common/container_mesh.hh
index 10d51ad..bb49fbc 100644
--- a/src/common/container_mesh.hh
+++ b/src/common/container_mesh.hh
@@ -1,222 +1,217 @@
/**
* @file container_mesh.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Sep 08 23:40:22 2014
*
* @brief This is the abstract container for mesh data
*
* @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_CONTAINER_MESH_HH__
#define __LIBMULTISCALE_CONTAINER_MESH_HH__
/* -------------------------------------------------------------------------- */
#include "container.hh"
#include "container_array.hh"
#include "ref_element.hh"
#include "ref_point_data.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename ContNodes, typename ContElems>
class ContainerMesh : public virtual LMObject,
public Container_base<typename ContNodes::Ref> {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
typedef ContNodes ContainerNodes;
typedef ContElems ContainerElems;
typedef typename ContNodes::iterator iterator;
typedef typename ContNodes::Ref Ref;
typedef typename ContElems::Ref RefElem;
typedef ContainerMesh<ContainerArray<Ref>, ContainerArray<RefElem>>
ContainerSubset;
static constexpr UInt Dim = ContainerNodes::Dim;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ContainerMesh(const LMID &id);
virtual ~ContainerMesh();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
operator ContainerNodes &() {
LM_FATAL(
"silent conversion from ContainerMesh to ContainerNode is not allowed");
return *container_nodes;
}
iterator begin(DOFType dt = dt_local) { return container_nodes.begin(dt); };
iterator end(DOFType dt = dt_local) { return container_nodes.end(dt); };
- void add(const Ref &node) { container_nodes.add(node); }
-
Ref &get(UInt index) { return container_nodes.get(index); }
void clear() {
container_nodes.clear();
container_elems.clear();
}
UInt size(DOFType dt = dt_local) { return container_nodes.size(dt); };
template <typename ContNodesPristine>
void computeAlteredConnectivity(ContNodesPristine &contNodes);
void setCommGroup(CommGroup &group) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
ContNodes &getArray() { return this->getContainerNodes(); };
ContNodes &getContainerNodes() { return container_nodes; };
ContElems &getContainerElems() { return container_elems; };
UInt subIndex2Index(UInt index) {
LM_ASSERT(nodeIndexList.size() > 0,
"nodeIndexList was not initialized: did you call "
<< "computeAlteredConnectivity method ?");
return nodeIndexList.searchInSorted(index);
};
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
ContNodes container_nodes;
ContElems container_elems;
ContainerArray<UInt> nodeIndexList;
};
/* -------------------------------------------------------------------------- */
template <typename ContNodes, typename ContElems>
ContainerMesh<ContNodes, ContElems>::ContainerMesh(const LMID &id)
: LMObject(id), container_nodes(id + ":nodes"),
container_elems(id + ":elems"), nodeIndexList(id + ":node_index") {}
/* -------------------------------------------------------------------------- */
template <typename ContNodes, typename ContElems>
ContainerMesh<ContNodes, ContElems>::~ContainerMesh() {}
/* -------------------------------------------------------------------------- */
template <typename ContNodes, typename ContElems>
template <typename ContNodesPristine>
void ContainerMesh<ContNodes, ContElems>::computeAlteredConnectivity(
ContNodesPristine &contNodes) {
for (auto &&el : container_elems) {
auto &&connectivity = el.globalIndexes();
for (UInt n = 0; n < connectivity.size(); ++n)
nodeIndexList.addInSorted(connectivity[n]);
}
for (auto &&n_index : nodeIndexList) {
container_nodes.add(contNodes.get(n_index));
}
// recompute the connectivity
std::vector<UInt> altered_connectivity;
for (UInt i = 0; i < container_elems.size(); ++i) {
auto &el = container_elems.get(i);
auto &&connectivity = el.globalIndexes();
altered_connectivity.resize(connectivity.size());
for (UInt n = 0; n < connectivity.size(); ++n)
altered_connectivity[n] = nodeIndexList.searchInSorted(connectivity[n]);
el.setAlteredConnectivity(altered_connectivity);
}
this->copyContainerInfo(contNodes);
}
/* -------------------------------------------------------------------------- */
template <typename ContNodes, typename ContElems>
void ContainerMesh<ContNodes, ContElems>::setCommGroup(CommGroup &group) {
Container_base<typename ContNodes::Ref>::setCommGroup(group);
container_nodes.setCommGroup(group);
container_elems.setCommGroup(group);
};
/* -------------------------------------------------------------------------- */
template <UInt _Dim>
class RefGenericElem : public RefElement<RefGenericElem<_Dim>> {
public:
static constexpr UInt Dim = _Dim;
template <typename Vec> void setConnectivity(Vec &&conn) {
- this->conn = conn.template cast<UInt>();
+ auto n_nodes = conn.rows() * conn.cols();
+ this->conn.clear();
+ for (UInt i = 0; i < n_nodes; ++i)
+ this->conn.push_back(conn[i]);
}
bool operator==(RefGenericElem &e) {
LM_TOIMPLEMENT;
return true;
};
//! return the local indexes (in DOF vector) for connected nodes
- inline std::vector<UInt> localIndexes() {
- std::vector<UInt> res;
- for (UInt i = 0; i < Dim; ++i) {
- res.push_back(this->conn[i]);
- }
- return res;
- }
+ inline std::vector<UInt> localIndexes() { return this->conn; }
//! return the global indexes (in DOF vector) for connected nodes
virtual std::vector<UInt> globalIndexes() { return this->localIndexes(); };
inline void setAlteredConnectivity(std::vector<UInt> &altered_connectivity) {
LM_TOIMPLEMENT;
};
inline std::vector<UInt> &getAlteredConnectivity() { LM_TOIMPLEMENT; };
inline bool isAltered() { return false; };
Real stress(UInt coord) { LM_TOIMPLEMENT; };
Real strain(UInt coord) { LM_TOIMPLEMENT; };
- Vector<Dim, UInt> conn;
+ std::vector<UInt> conn;
};
template <UInt Dim>
using ContainerGenericMesh = ContainerMesh<ContainerArray<RefPointData<Dim>>,
ContainerArray<RefGenericElem<Dim>>>;
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_CONTAINER_MESH_HH__ */
diff --git a/src/common/id_manager.hh b/src/common/id_manager.hh
index f984a88..30c943f 100644
--- a/src/common/id_manager.hh
+++ b/src/common/id_manager.hh
@@ -1,282 +1,284 @@
/**
* @file id_manager.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Sep 08 23:40:22 2014
*
* @brief This is the generic manager of objects from their IDs (string)
*
* @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_ID_MANAGER_HH__
#define __LIBMULTISCALE_ID_MANAGER_HH__
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "lm_object.hh"
#include <algorithm>
#include <fstream>
#include <list>
#include <map>
#include <memory>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
-struct UnknownID : public std::runtime_error{
- UnknownID(const LMID & id, const std::string & manager_type): std::runtime_error(
- "Unknown id: " + id + " from manager type: " + manager_type), id(id){}
+struct UnknownID : public std::runtime_error {
+ UnknownID(const LMID &id, const std::string &manager_type)
+ : std::runtime_error("Unknown id: " + id +
+ " from manager type: " + manager_type),
+ id(id) {}
LMID id;
};
-
template <typename T> class IDManager {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
typedef std::string typeID;
typedef std::map<typeID, T *> Cont;
-
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
friend struct std::default_delete<IDManager<T>>;
protected:
IDManager();
virtual ~IDManager();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
static void destroy();
//! add a geometry as a shared pointer
typeID addObject(std::shared_ptr<T> g);
//! add a geometry as an object
typeID addObject(T &g);
//! generate a new id automatically
typeID generateNewID();
//! remove a managed object
void removeObject(const typeID &ID);
//! remove a managed object
void removeObject(const T &ID);
//! true if singleton already allocated
static bool allocated() { return static_pointer != nullptr; }
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! return a reference to a geometry from its ID
T &getObject(const typeID &ID);
//! return a shared_pointer from its ID
std::shared_ptr<T> getSharedObject(const typeID &ID);
//! return a casted object from its ID
template <typename Cast> inline Cast &getCastedObject(const typeID &ID);
//! search id of a given geometry
typeID getID(T &obj);
protected:
//! accessor to the singleton manager
template <typename IT> static IT &getManager();
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! mapping between IDs and pointers
std::map<typeID, T *> objects;
//! mapping between IDs and shared pointers
std::map<typeID, std::shared_ptr<T>> shared_pointer_objects;
//! file to which output the object infos
std::ofstream fout;
static std::unique_ptr<IDManager<T>> static_pointer;
std::list<typeID> insert_order;
};
/* -------------------------------------------------------------------------- */
template <typename T>
template <typename IT>
inline IT &IDManager<T>::getManager() {
if (!static_pointer)
static_pointer.reset(new IT());
IT &ref = dynamic_cast<IT &>(*static_pointer);
return ref;
}
/* -------------------------------------------------------------------------- */
template <typename T> inline IDManager<T>::IDManager() {
std::string filename = "log.";
- filename += typeid(*this).name();
+ filename += typeinfo<decltype(*this)>();
fout.open(filename.c_str(), std::ios_base::out);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline typename IDManager<T>::typeID IDManager<T>::generateNewID() {
UInt i = objects.size();
std::stringstream str;
str << "geom" << i;
return str.str();
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline typename IDManager<T>::typeID
IDManager<T>::addObject(std::shared_ptr<T> obj) {
if (obj == NULL)
LM_FATAL("attempt to register a null pointer as " << typeid(*this).name());
typeID id = obj->getID();
if (objects.count(id) > 0 && objects[id])
LM_FATAL("ID " << id << " has already been used for " << *objects[id]);
shared_pointer_objects[id] = obj;
objects[id] = obj.get();
DUMPFILE(fout,
- id << "(" << *obj << ") added to " << typeid(*this).name()
+ id << "(" << *obj << ") added to " << typeinfo<decltype(*this)>()
<< " manager");
insert_order.push_back(id);
return id;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline typename IDManager<T>::typeID IDManager<T>::addObject(T &obj) {
typeID id = obj.getID();
if (objects.count(id) > 0 && objects[id])
LM_FATAL("ID " << id << " has already been used for " << *objects[id]);
objects[id] = &obj;
DUMPFILE(fout,
- id << "(" << obj << ") added to " << typeid(*this).name()
+ id << "(" << obj << ") added to " << typeinfo<decltype(*this)>()
<< " manager");
insert_order.push_back(id);
return id;
}
/* -------------------------------------------------------------------------- */
template <typename T> inline void IDManager<T>::removeObject(const T &obj) {
typeID ID = obj.getID();
removeObject(ID);
- DUMPFILE(fout, "removed " << ID << " from manager " << typeid(*this).name());
+ DUMPFILE(fout,
+ "removed " << ID << " from manager " << typeinfo<decltype(*this)>());
}
/* -------------------------------------------------------------------------- */
template <typename T> inline void IDManager<T>::removeObject(const typeID &ID) {
if (objects.count(ID) == 0) {
DUMP("impossible to delete entry -> unknown id: " << ID, DBG_INFO);
return;
}
objects.erase(ID);
if (shared_pointer_objects.count(ID))
shared_pointer_objects.erase(ID);
insert_order.remove(ID);
}
/* -------------------------------------------------------------------------- */
template <typename T> inline T &IDManager<T>::getObject(const typeID &ID) {
if (objects.count(ID) == 0) {
- throw UnknownID(ID, typeid(*this).name());
+ throw UnknownID(ID, typeinfo<decltype(*this)>());
}
return *objects[ID];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline std::shared_ptr<T> IDManager<T>::getSharedObject(const typeID &ID) {
if (shared_pointer_objects.count(ID) == 0) {
- throw UnknownID(ID, typeid(*this).name());
+ throw UnknownID(ID, typeinfo<decltype(*this)>());
}
return shared_pointer_objects[ID];
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <typename Cast>
inline Cast &IDManager<T>::getCastedObject(const typeID &ID) {
if (objects.count(ID) == 0) {
- throw UnknownID(ID, typeid(*this).name());
+ throw UnknownID(ID, typeinfo<decltype(*this)>());
}
auto *ptr = dynamic_cast<Cast *>(objects[ID]);
if (!ptr) {
- LM_FATAL("bad cast for ID: " << ID << " from " << typeid(*this).name()
- << " to " << typeid(Cast).name());
+ LM_FATAL("bad cast for ID: " << ID << " from "
+ << typeinfo<decltype(*this)>() << " to "
+ << typeinfo<Cast>());
}
return *ptr;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline typename IDManager<T>::typeID IDManager<T>::getID(T &obj) {
for (auto &&o : objects) {
if (o.second == obj)
return o.first;
}
DUMPFILE(fout,
"Warning : auto registering of " << this->getID() << " {\n"
<< *obj << "\n}");
return (addObject(obj));
}
/* -------------------------------------------------------------------------- */
template <typename T> inline IDManager<T>::~IDManager() {
destroy();
- DUMPFILE(fout, "all " << typeid(*this).name() << " cleaned...");
+ DUMPFILE(fout, "all " << typeinfo<decltype(*this)>() << " cleaned...");
fout.close();
}
/* -------------------------------------------------------------------------- */
template <typename T> inline void IDManager<T>::destroy() {
if (static_pointer == nullptr)
return;
static_pointer->shared_pointer_objects.clear();
if (static_pointer)
static_pointer.release();
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_ID_MANAGER_HH__ */
diff --git a/src/common/lm_types.hh b/src/common/lm_types.hh
index cb41187..ac0db5a 100644
--- a/src/common/lm_types.hh
+++ b/src/common/lm_types.hh
@@ -1,447 +1,459 @@
/**
* @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 <exception>
#include <ostream>
#include <vector>
/* -------------------------------------------------------------------------- */
__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;
}
/* -------------------------------------------------------------------------- */
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,
// 1u>;
template <UInt Dim, typename T = Real>
struct Vector : public Eigen::Matrix<T, Dim, 1u> {
using Eigen::Matrix<T, Dim, 1u>::Matrix;
// Vector(const std::vector<T> & vec){
// if (vec.size() != Dim)
// throw std::runtime_error("not matching vector size");
// for (UInt i = 0 ; i < Dim ; ++i ){
// (*this)[i] = vec[i];
// }
// }
};
template <typename T> struct Vector<1, T> : public Eigen::Matrix<T, 1, 1u> {
using Eigen::Matrix<T, 1, 1u>::Matrix;
Vector<1, T>(const T &val) { (*this)[0] = val; }
// operator T() const { return (*this)[0]; }
};
// 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;
using Eigen::Map<Eigen::Array<Real, Dim, 1u>>::operator=;
};
// 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;
using Eigen::Matrix<Real, Dim, Dim>::operator=;
};
// 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>
struct MatrixView : public Eigen::Map<Eigen::Matrix<Real, Dim, Dim>> {
using Eigen::Map<Eigen::Matrix<Real, Dim, Dim>>::Map;
using Eigen::Map<Eigen::Matrix<Real, Dim, Dim>>::operator=;
};
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,
_id = 13,
_proc = 14,
_charge = 15,
ft_uninitialised = 32767
};
template <UInt Dim, FieldType ft> struct FieldTraits {
using TensorType = Vector<Dim>;
using TensorViewType = VectorView<Dim>;
};
template <UInt Dim> struct FieldTraits<Dim, _mass> {
using TensorType = Vector<1>;
using TensorViewType = VectorView<1>;
};
/* -------------------------------------------------------------------------- */
// 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, typename T1> void init_array(T *var, T1 &arg) {
*var = &arg;
}
template <typename T, typename T1, typename... Args>
void init_array(T *var, T1 &arg, Args &... args) {
*var = &arg;
init_array(var + 1, args...);
}
template <UInt Dim, typename T = Real> struct VectorProxy {
- VectorProxy(VectorProxy &v) {
+ explicit VectorProxy(VectorProxy &v) {
for (UInt i = 0; i < Dim; ++i)
tup[i] = v.tup[i];
};
- VectorProxy(VectorProxy &&v) {
+ VectorProxy(const VectorProxy &&v) {
for (UInt i = 0; i < Dim; ++i)
- tup[i] = v.tup[i];
+ tup[i] = std::move(v.tup[i]);
+ };
+
+ VectorProxy &operator=(VectorProxy v) {
+ for (UInt i = 0; i < Dim; ++i)
+ *this->tup[i] = *v.tup[i];
+ return *this;
+ };
+
+ VectorProxy &operator=(Vector<Dim, T> &v) {
+ for (UInt i = 0; i < Dim; ++i)
+ *this->tup[i] = v[i];
+ return *this;
};
template <typename... Args> VectorProxy(Args &... args) {
init_array(tup, args...);
};
auto &operator[](UInt i) { return *tup[i]; }
const auto &operator[](UInt i) const { return *tup[i]; }
template <typename Tscal> auto operator*(const Tscal &scalar) {
Vector<Dim> result;
for (UInt i = 0; i < Dim; ++i)
result[i] = *tup[i] * scalar;
return result;
}
template <typename Tscal> void operator*=(const Tscal &scalar) {
for (UInt i = 0; i < Dim; ++i)
*tup[i] *= scalar;
}
operator Eigen::Matrix<T, Dim, 1u>() {
Eigen::Matrix<T, Dim, 1u> res;
for (UInt i = 0; i < Dim; ++i)
res[i] = (*this)[i];
return res;
}
operator Vector<Dim, T>() const {
Vector<Dim, T> res;
for (UInt i = 0; i < Dim; ++i)
res[i] = (*this)[i];
return res;
}
UInt size() { return Dim; };
T *tup[Dim];
};
template <UInt Dim, typename T>
inline std::ostream &operator<<(std::ostream &os, VectorProxy<Dim, T> v) {
Vector<Dim, T> v2 = v;
os << v2;
return os;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
#endif /* __LIBMULTISCALE_LM_TYPES_HH__ */
diff --git a/src/dd/domain_dd_interface.hh b/src/dd/domain_dd_interface.hh
index f859f8d..2bac9b7 100644
--- a/src/dd/domain_dd_interface.hh
+++ b/src/dd/domain_dd_interface.hh
@@ -1,102 +1,75 @@
/**
* @file domain_dd_interface.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Max Hodapp <max.hodapp@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief Common interface to all DD domains
*
* @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_DD_INTERFACE_HH__
#define __LIBMULTISCALE_DOMAIN_DD_INTERFACE_HH__
/* -------------------------------------------------------------------------- */
#include "container_mesh.hh"
#include "domain_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
class Geometry;
/* -------------------------------------------------------------------------- */
/**
* Class DomainDD
* domain dd
*/
class DomainDDInterface : public DomainInterface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DomainDDInterface(CommGroup &group) : DomainInterface(group) {
type = DDTYPE;
};
virtual ~DomainDDInterface(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
//! return the domain geometry
virtual Geometry &getGeom() = 0;
//! initiation function
virtual void init() = 0;
-
- //! compute displacement field for a set of positions
- virtual void computeDisplacements(std::vector<Real> &positions,
- std::vector<Real> &displacements) {
- LM_TOIMPLEMENT;
- }
-
virtual void remesh() { LM_TOIMPLEMENT; }
-
- //! 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) {
- LM_TOIMPLEMENT;
- };
-
- //! 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
- {
- LM_TOIMPLEMENT;
- };
};
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_DOMAIN_DD_INTERFACE_HH__ */
diff --git a/src/dd/paradis/container_paradis.hh b/src/dd/paradis/container_paradis.hh
index 1c02b80..ab24834 100644
--- a/src/dd/paradis/container_paradis.hh
+++ b/src/dd/paradis/container_paradis.hh
@@ -1,237 +1,321 @@
/**
* @file container_paradis.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date Mon Jul 28 09:52:51 2014
*
* @brief ParaDiS container of DOFs
*
* @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_CONTAINER_PARADIS_HH__
#define __LIBMULTISCALE_CONTAINER_PARADIS_HH__
/* -------------------------------------------------------------------------- */
#include "container.hh"
+#include "paradis_home.hh"
#include "ref_elem_paradis.hh"
#include "ref_node_paradis.hh"
/* -------------------------------------------------------------------------- */
-#define PARALLEL
-extern "C" {
-#include "Home.h"
-#include "Typedefs.h"
-#undef X
-#undef Y
-}
-#undef PARALLEL
-/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class ReferenceManagerParaDis {
public:
template <typename Obj>
void attachObject(Obj &obj, ContainerArray<RefNodeParaDiS>) {
LM_TOIMPLEMENT;
}
};
/* -------------------------------------------------------------------------- */
/**
* Class ContainerParadis
*
*/
class IteratorNodesParaDiS;
class ContainerNodesParaDiS : public virtual LMObject,
public Container_base<RefNodeParaDiS> {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
using typename Container_base<RefNodeParaDiS>::iterator_base;
friend class IteratorNodesParaDiS;
public:
static constexpr UInt Dim = 3;
using iterator = IteratorNodesParaDiS;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ContainerNodesParaDiS(const LMID &id)
- : LMObject(id), Container_base<RefNodeParaDiS>(), paradis_ptr(nullptr){};
+ : LMObject(id), Container_base<RefNodeParaDiS>(), node(0){};
~ContainerNodesParaDiS(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
virtual ReferenceManagerParaDis &getRefManager() {
LM_TOIMPLEMENT;
return refmanager;
};
- void setGrain(Home_t *home) {
- node.setHome(home);
- paradis_ptr = home;
- // CHANGE THIS
- // node.setStressVector(paradis_ptr->getDisloNodes()->getStress());//CHANGE
- // THIS//PROBLEM, see ref_node_paradis.h!!!!!!
- };
-
IteratorNodesParaDiS begin(DOFType dt = dt_local);
IteratorNodesParaDiS end(DOFType dt = dt_local);
RefNodeParaDiS &get(UInt index);
inline UInt size(DOFType dt = dt_local);
void init(){};
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
RefNodeParaDiS node;
private:
- Home_t *paradis_ptr;
-
//! reference manager for ParaDis
ReferenceManagerParaDis refmanager;
};
/* -------------------------------------------------------------------------- */
inline UInt ContainerNodesParaDiS::size(DOFType dt) {
- return this->paradis_ptr->newNodeKeyPtr -
- this->paradis_ptr->recycledNodeHeapSize;
+ return ParaDiSHome::paradis_ptr->newNodeKeyPtr -
+ ParaDiSHome::paradis_ptr->recycledNodeHeapSize;
// if (this->paradis_ptr == NULL)
// return 0;
// unsigned int nb_elem(0);
// for (int i = 0; i < this->paradis_ptr->newNodeKeyPtr; ++i) {
// if (this->paradis_ptr->nodeKeys[i] != (Node_t *)NULL)
// ++nb_elem;
// }
// return nb_elem;
};
/* -------------------------------------------------------------------------- */
///
/// Doing elements now
///
/* -------------------------------------------------------------------------- */
class IteratorElemsParaDiS;
class ContainerElemsParaDiS : public virtual LMObject,
public Container_base<RefElemParaDiS> {
public:
static constexpr UInt Dim = 3;
public:
ContainerElemsParaDiS(const LMID &id) : LMObject(id), paradis_ptr(nullptr) {
nodes = nullptr;
};
~ContainerElemsParaDiS(){};
using iterator = IteratorElemsParaDiS;
friend class IteratorElemsParaDiS;
+ template <typename ContMesh, typename VecBurgers, typename VecNormals>
+ void add(ContMesh &segments, VecBurgers &&burgers, VecNormals &&normals);
+
void setNodes(ContainerNodesParaDiS &nodes) { this->nodes = &nodes; }
void setGrain(Home_t *home) { paradis_ptr = home; }
//! return an associated iterator
IteratorElemsParaDiS begin(DOFType dt = dt_local);
IteratorElemsParaDiS end(DOFType dt = dt_local);
//! get an item from its index
RefElemParaDiS &get(UInt index);
//! return the number of contained items
inline UInt size(DOFType dt = dt_local);
private:
Home_t *paradis_ptr;
ContainerNodesParaDiS *nodes;
};
__END_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
#include "iterator_paradis.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
inline IteratorNodesParaDiS ContainerNodesParaDiS::begin(DOFType dt) {
IteratorNodesParaDiS it(*this);
return it;
}
/* -------------------------------------------------------------------------- */
inline UInt ContainerElemsParaDiS::size(DOFType dt) {
UInt cpt = 0;
for (auto &n : *nodes) {
cpt += n.nbArms();
}
return cpt;
};
/* -------------------------------------------------------------------------- */
inline IteratorNodesParaDiS ContainerNodesParaDiS::end(DOFType dt) {
- IteratorNodesParaDiS it(*this, node.paradis_ptr->newNodeKeyPtr);
+ IteratorNodesParaDiS it(*this, ParaDiSHome::paradis_ptr->newNodeKeyPtr);
return it;
}
/* -------------------------------------------------------------------------- */
inline RefNodeParaDiS &ContainerNodesParaDiS::get(unsigned int index) {
node.setIndex(index);
return node;
}
/* -------------------------------------------------------------------------- */
inline IteratorElemsParaDiS ContainerElemsParaDiS::begin(DOFType dt) {
return IteratorElemsParaDiS(*this, nodes->begin());
}
/* -------------------------------------------------------------------------- */
inline IteratorElemsParaDiS ContainerElemsParaDiS::end(DOFType dt) {
return IteratorElemsParaDiS(*this, nodes->end());
}
/* -------------------------------------------------------------------------- */
inline RefElemParaDiS &ContainerElemsParaDiS::get(unsigned int index) {
LM_TOIMPLEMENT;
}
+/* -------------------------------------------------------------------------- */
+
+template <typename ContMesh, typename VecBurgers, typename VecNormals>
+inline void ContainerElemsParaDiS::add(ContMesh &mesh, VecBurgers &&burgers,
+ VecNormals &&normals) {
+
+ auto &&points = mesh.getContainerNodes();
+ UInt nb_points = points.size();
+
+ std::vector<Node_t *> created_nodes;
+ std::vector<UInt> numNbrs;
+
+ // adding the points
+ for (auto &&p : points) {
+ Node_t *newNode = GetNewNativeNode(this->paradis_ptr);
+
+ auto &&pos = p.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);
+ }
+
+ // adding the arms/elements
+ numNbrs.resize(nb_points);
+
+ auto &&elems = mesh.getContainerElems();
+ for (auto &&el : elems) {
+ 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 (auto &&el : elems) {
+ 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[0];
+ srcNode->burgY[current_arm_src] = burgers[1];
+ srcNode->burgZ[current_arm_src] = burgers[2];
+
+ srcNode->nx[current_arm_src] = normals[0];
+ srcNode->ny[current_arm_src] = normals[1];
+ srcNode->nz[current_arm_src] = normals[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[0];
+ dstNode->burgY[current_arm_dst] = -1. * burgers[1];
+ dstNode->burgZ[current_arm_dst] = -1. * burgers[2];
+
+ dstNode->nx[current_arm_dst] = normals[0];
+ dstNode->ny[current_arm_dst] = normals[1];
+ dstNode->nz[current_arm_dst] = normals[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]);
+ // }
+ // }
+}
+
+/* -------------------------------------------------------------------------- */
+
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_CONTAINER_PARADIS_HH__ */
diff --git a/src/dd/paradis/domain_paradis.cc b/src/dd/paradis/domain_paradis.cc
index 6f0d98c..15c7bbb 100644
--- a/src/dd/paradis/domain_paradis.cc
+++ b/src/dd/paradis/domain_paradis.cc
@@ -1,944 +1,865 @@
/**
* @file domain_paradis.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Moseley Philip Arthur <philip.moseley@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
* @author Max Hodapp <max.hodapp@epfl.ch>
*
* @date Fri Jul 11 15:47:44 2014
*
* @brief ParaDiS domain
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* TODO:
* - the function addHybSegments() is currently a rudimentary implementation
* for simulations containing only one dislocations
* ==> method has to be generalized and additional checks (e.g. if hybrid
* segments are actually self-consistent, etc...) have to be performed to
* improve robustness
*
*/
/* -------------------------------------------------------------------------- */
#define TIMER
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#define PARALLEL
#include "communicator.hh"
#include "domain_paradis.hh"
#include "geometry_manager.hh"
#include "stdio.h"
/* -------------------------------------------------------------------------- */
// ParaDiS includes
extern "C" {
// void TimerClearAll(Home_t *home);
#include "Home.h"
#include "Param.h"
#include "Typedefs.h"
/* -------------------------------------------------------------------------- */
Home_t *InitHome();
void ParadisStep(Home_t *home);
void RemeshAfterCoupling(Home_t *home);
void PointDisplacements(Home_t *home);
void SegmentDisplacementsOfAtom(libmultiscale::Vector<3u> &A,
libmultiscale::Vector<3u> &B,
libmultiscale::Vector<3u> &C,
libmultiscale::Vector<3u> &burg,
libmultiscale::Vector<3u> &dispPointU,
libmultiscale::Vector<3u> &dispPointX,
libmultiscale::Real pois);
void ComputeClosurePoint(libmultiscale::Vector<3u> &A,
libmultiscale::Vector<3u> &B,
libmultiscale::Vector<3u> &burg,
libmultiscale::Vector<3u> &C);
void ComputeProjectedPoint(libmultiscale::Vector<3u> &A,
libmultiscale::Vector<3u> &C,
libmultiscale::Vector<3u> &burg,
libmultiscale::Vector<3u> &P);
void GetFieldPointStress(Home_t *home, real8 x, real8 y, real8 z,
real8 totStress[3][3]);
#ifdef FPES_ON
#include <fpcontrol.h>
#endif
}
#undef PARALLEL
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
void DomainParaDiS::remesh() {
- RemeshAfterCoupling(this->paradis_ptr);
+ RemeshAfterCoupling(ParaDiSHome::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);
+ // InitTopologyExemptions(ParaDiSHome::paradis_ptr);
+ // ClearOpList(ParaDiSHome::paradis_ptr);
+ // SortNodesForCollision(ParaDiSHome::paradis_ptr);
+ // SplitMultiNodes(ParaDiSHome::paradis_ptr);
- // CrossSlip(this->paradis_ptr);
- // HandleCollisions(this->paradis_ptr);
- // ProximityCollisions(this->paradis_ptr);
+ // CrossSlip(ParaDiSHome::paradis_ptr);
+ // HandleCollisions(ParaDiSHome::paradis_ptr);
+ // ProximityCollisions(ParaDiSHome::paradis_ptr);
// #ifdef PARALLEL
// #ifdef PARADIS_IN_LIBMULTISCALE
- // CommSendRemesh(this->paradis_ptr);
+ // CommSendRemesh(ParaDiSHome::paradis_ptr);
// #endif
// #endif
- // FixRemesh(this->paradis_ptr);
+ // FixRemesh(ParaDiSHome::paradis_ptr);
- // Remesh(this->paradis_ptr);
- // Migrate(this->paradis_ptr);
+ // Remesh(ParaDiSHome::paradis_ptr);
+ // Migrate(ParaDiSHome::paradis_ptr);
// #ifdef PARALLEL
// #ifdef PARADIS_IN_LIBMULTISCALE
- // RecycleGhostNodes(this->paradis_ptr);
+ // RecycleGhostNodes(ParaDiSHome::paradis_ptr);
// #endif
// #endif
- // SortNativeNodes(this->paradis_ptr);
+ // SortNativeNodes(ParaDiSHome::paradis_ptr);
// #ifdef PARALLEL
// #ifdef PARADIS_IN_LIBMULTISCALE
- // CommSendGhosts(this->paradis_ptr);
+ // CommSendGhosts(ParaDiSHome::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::computeDisplacements(std::vector<Real> &positions,
+// std::vector<Real> &displacements) {
+// if (positions.empty()) {
+// displacements.clear();
+// return;
+// }
+// if (!ParaDiSHome::paradis_ptr) {
+// return;
+// }
+// // Pass data to Paradis.
+// unsigned int size = positions.size();
+// LM_ASSERT(size % 3 == 0, "paradis computeDisplacements requires 3d.");
+// if (displacements.size() != size)
+// displacements.resize(size);
+// ParaDiSHome::paradis_ptr->dispPointCount = size / 3;
+// ParaDiSHome::paradis_ptr->dispPointX = &positions[0];
+// ParaDiSHome::paradis_ptr->dispPointU = &displacements[0];
+// // Calculate the displacements.
+// PointDisplacements(ParaDiSHome::paradis_ptr);
+// }
/* -------------------------------------------------------------------------- */
-void DomainParaDiS::computeBarnett(Vector<3u> &A, Vector<3u> &B, Vector<3u> &C,
- Vector<3u> &burg, Vector<3u> &dispPointU,
- Vector<3u> &dispPointX, Real pois) {
- SegmentDisplacementsOfAtom(A, B, C, burg, dispPointU, dispPointX, pois);
-}
-
-/* -------------------------------------------------------------------------- */
-void DomainParaDiS::computeClosurePoint(Vector<3u> &A, Vector<3u> &B,
- Vector<3u> &burg, Vector<3u> &C) {
- ComputeClosurePoint(A, B, burg, C);
-}
-/* -------------------------------------------------------------------------- */
-void DomainParaDiS::computeProjectedPoint(Vector<3u> &A, Vector<3u> &C,
- Vector<3u> &burg, Vector<3u> &P) {
- ComputeProjectedPoint(A, C, burg, P);
-}
-/* -------------------------------------------------------------------------- */
+// void DomainParaDiS::computeBarnett(Vector<3u> &A, Vector<3u> &B, Vector<3u>
+// &C,
+// Vector<3u> &burg, Vector<3u> &dispPointU,
+// Vector<3u> &dispPointX, Real pois) {
+// SegmentDisplacementsOfAtom(A, B, C, burg, dispPointU, dispPointX, pois);
+// }
+
+// /* --------------------------------------------------------------------------
+// */
+// void DomainParaDiS::computeClosurePoint(Vector<3u> &A, Vector<3u> &B,
+// Vector<3u> &burg, Vector<3u> &C) {
+// ComputeClosurePoint(A, B, burg, C);
+// }
+// /* --------------------------------------------------------------------------
+// */
+// void DomainParaDiS::computeProjectedPoint(Vector<3u> &A, Vector<3u> &C,
+// Vector<3u> &burg, Vector<3u> &P) {
+// ComputeProjectedPoint(A, C, burg, P);
+// }
+// /* --------------------------------------------------------------------------
+// */
void DomainParaDiS::init() {
auto &group = this->getCommGroup();
if (!group.amIinGroup())
return;
const UInt nbprocs = group.size();
// Mostly taken from ParaDiS source files Main.c and ParadisInit.c
std::string dlb_flag = "-r";
std::string datafile_flag = "-d";
std::string ctrlfile_flag = "-c";
std::string fake_argv0 = "IGNORED";
std::string dlb_num;
std::stringstream(dlb_num) << numDLBCycles;
std::vector<const char *> argv({fake_argv0.c_str(), datafile_flag.c_str(),
data_filename.c_str(), dlb_flag.c_str(),
dlb_num.c_str(), ctrlfile_flag.c_str(),
control_filename.c_str()});
if (mobility_filename != "")
argv.push_back(mobility_filename.c_str());
unsigned int argc = argv.size();
/*
* On some systems, the getrusage() call made by Meminfo() to get
* the memory resident set size does not work properly. In those
* cases, the function will try to return the current heap size
* instead. This initial call allows meminfo() to get a copy of
* the original heap pointer so subsequent calls can calculate the
* heap size by taking the difference of the original and current
* heap pointers.
*/
Meminfo(&(this->memSize));
/*
* On linux systems (e.g. MCR) if built to have floating point exceptions
* turned on, invoke macro to do so
*/
#ifdef FPES_ON
unmask_std_fpes();
#endif
- this->paradis_ptr = InitHome();
- TimerInit(this->paradis_ptr);
- TimerStart(this->paradis_ptr, TOTAL_TIME);
- TimerStart(this->paradis_ptr, INITIALIZE);
+ ParaDiSHome::paradis_ptr = InitHome();
+ TimerInit(ParaDiSHome::paradis_ptr);
+ TimerStart(ParaDiSHome::paradis_ptr, TOTAL_TIME);
+ TimerStart(ParaDiSHome::paradis_ptr, INITIALIZE);
if (group.getMyRank() == 0)
std::cout << "BEGIN: Initialize ParaDiS on " << nbprocs << " processors."
<< std::endl;
- Initialize(this->paradis_ptr, argc, const_cast<char **>(&argv[0]),
+ Initialize(ParaDiSHome::paradis_ptr, argc, const_cast<char **>(&argv[0]),
group.getMPICommGroup());
- if (strcmp(this->paradis_ptr->param->timestepIntegrator,
+ if (strcmp(ParaDiSHome::paradis_ptr->param->timestepIntegrator,
"trapezoidtimestep") != 0) {
LM_FATAL("The timestep integrator for ParaDiS needs to be "
<< "\"trapezoidtimestep\" in order for ParaDis to work with "
<< "LibMultiscale. You specified \"timestepIntegrator = "
- << this->paradis_ptr->param->timestepIntegrator
+ << ParaDiSHome::paradis_ptr->param->timestepIntegrator
<< "\" in the control file \"" << control_filename << "\".");
}
- TimerStop(this->paradis_ptr, INITIALIZE);
+ TimerStop(ParaDiSHome::paradis_ptr, INITIALIZE);
if (nbprocs > 1)
- MPI_Barrier(this->paradis_ptr->MPI_COMM_PARADIS);
+ MPI_Barrier(ParaDiSHome::paradis_ptr->MPI_COMM_PARADIS);
- this->paradis_ptr->param->timeMax = this->timestep;
- this->paradis_ptr->cycle = this->paradis_ptr->param->cycleStart;
+ ParaDiSHome::paradis_ptr->param->timeMax = this->timestep;
+ ParaDiSHome::paradis_ptr->cycle = ParaDiSHome::paradis_ptr->param->cycleStart;
/*
* Perform the needed number (if any) of load-balance-only
* steps before doing the main processing loop. These steps
* perform only the minimal amount of stuff needed to
* estimate per-process load, move boundaries and migrate
* nodes among processsors to get a good initial balance.
*/
time_t tp;
- int initialDLBCycles = this->paradis_ptr->param->numDLBCycles;
- TimerStart(this->paradis_ptr, INITIALIZE);
- if ((this->paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
+ int initialDLBCycles = ParaDiSHome::paradis_ptr->param->numDLBCycles;
+ TimerStart(ParaDiSHome::paradis_ptr, INITIALIZE);
+ if ((ParaDiSHome::paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
time(&tp);
printf(" +++ Beginning %d load-balancing steps at %s", initialDLBCycles,
asctime(localtime(&tp)));
}
- while (this->paradis_ptr->param->numDLBCycles > 0) {
- ParadisStep(this->paradis_ptr);
- this->paradis_ptr->cycle++;
- this->paradis_ptr->param->numDLBCycles--;
+ while (ParaDiSHome::paradis_ptr->param->numDLBCycles > 0) {
+ ParadisStep(ParaDiSHome::paradis_ptr);
+ ParaDiSHome::paradis_ptr->cycle++;
+ ParaDiSHome::paradis_ptr->param->numDLBCycles--;
}
- if ((this->paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
+ if ((ParaDiSHome::paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) {
time(&tp);
printf(" +++ Completed load-balancing steps at %s",
asctime(localtime(&tp)));
}
- TimerStop(this->paradis_ptr, INITIALIZE);
+ TimerStop(ParaDiSHome::paradis_ptr, INITIALIZE);
/*
* Any time spent doing the initial DLB-only steps should
* just be attributed to initialization time, so be sure to
* reset the other timers before going into the main
* computational loop
*/
- TimerInitDLBReset(this->paradis_ptr);
+ TimerInitDLBReset(ParaDiSHome::paradis_ptr);
/*
* The cycle number has been incremented during the initial
* load-balance steps, and previous steps so reset it to the
* proper startingvalue before entering the main processing loop.
* In paradis, cycle uses 1-based indexing.
*/
- this->paradis_ptr->cycle = this->paradis_ptr->param->cycleStart;
- this->paradis_ptr->cycleLMS = 1;
+ ParaDiSHome::paradis_ptr->cycle = ParaDiSHome::paradis_ptr->param->cycleStart;
+ ParaDiSHome::paradis_ptr->cycleLMS = 1;
- this->container_dd.getContainerNodes().setGrain(this->paradis_ptr);
- this->container_dd.getContainerElems().setGrain(this->paradis_ptr);
this->container_dd.getContainerElems().setNodes(
this->container_dd.getContainerNodes());
#ifndef LM_OPTIMIZED
// Debug tests.
unsigned int counter = 0;
DUMP("Printing all initial ParaDiS nodes...", DBG_INFO);
for (auto &&node : this->container_dd.getContainerNodes()) {
DUMP("Processor " << lm_my_proc_id << ", Node " << ++counter << "("
<< node.position() << ")",
DBG_INFO);
}
#endif
// ContainerParaDiSDisplacement displs = this->getContainerDisplacements();
// DUMP("Number of ParaDiS displacement fieldpoints: "<<displs.size(),
// DBG_INFO);
// End debug tests.
- if (this->paradis_ptr->myDomain == 0)
+ if (ParaDiSHome::paradis_ptr->myDomain == 0)
std::cout << "END: Initialize ParaDiS on " << nbprocs << " processors."
<< std::endl;
}
/* -------------------------------------------------------------------------- */
void DomainParaDiS::performStep1() {
- this->paradis_ptr->param->timeLMSCycle = 0.;
- while (this->paradis_ptr->param->timeLMSCycle <
- this->paradis_ptr->param->timeMax) {
+ ParaDiSHome::paradis_ptr->param->timeLMSCycle = 0.;
+ while (ParaDiSHome::paradis_ptr->param->timeLMSCycle <
+ ParaDiSHome::paradis_ptr->param->timeMax) {
STARTTIMER("ParadisStep");
- ParadisStep(this->paradis_ptr);
- TimerClearAll(this->paradis_ptr);
+ ParadisStep(ParaDiSHome::paradis_ptr);
+ TimerClearAll(ParaDiSHome::paradis_ptr);
STOPTIMER("ParadisStep");
}
- this->paradis_ptr->cycleLMS++;
+ ParaDiSHome::paradis_ptr->cycleLMS++;
this->getOutput()->incRelease();
if (this->boundary) {
auto c =
GeometryManager::getManager().getGeometry(this->geom)->getBoundingBox();
auto &&maxval = c.getXmax();
auto &&minval = c.getXmin();
// Geometry* dd = GeometryManager::getManager().getGeometry(this->geom);
// Cube & dd_cube = dynamic_cast<Cube&>(*dd);
// to mimic strong grain boundary...
- // Real xmin = 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 = ParaDiSHome::paradis_ptr->param->xBoundMin + 10.0;
+ // Real xmax = ParaDiSHome::paradis_ptr->param->xBoundMax - 10.0;
+ // Real ymin = ParaDiSHome::paradis_ptr->param->yBoundMin + 10.0;
+ // Real ymax = ParaDiSHome::paradis_ptr->param->yBoundMax - 10.0;
+ // Real zmin = ParaDiSHome::paradis_ptr->param->zBoundMin + 10.0;
+ // Real zmax = ParaDiSHome::paradis_ptr->param->zBoundMax - 10.0;
+
+ // Real xmin = ParaDiSHome::paradis_ptr->param->xBoundMin;
+ // Real xmax = ParaDiSHome::paradis_ptr->param->xBoundMax;
+ // Real ymin = ParaDiSHome::paradis_ptr->param->yBoundMin;
+ // Real ymax = ParaDiSHome::paradis_ptr->param->yBoundMax;
+ // Real zmin = ParaDiSHome::paradis_ptr->param->zBoundMin;
+ // Real zmax = ParaDiSHome::paradis_ptr->param->zBoundMax;
Real xmin = minval[0];
Real xmax = maxval[0];
Real ymin = minval[1];
Real ymax = maxval[1];
Real zmin = minval[2];
Real zmax = maxval[2];
for (auto &&node : this->container_dd) {
if ((xmin > node.position()[0]) || (node.position()[0] > xmax) ||
(ymin > node.position()[1]) || (node.position()[1] > ymax) ||
(zmin > node.position()[2]) || (node.position()[2] > zmax)) {
node.slaving();
node.constrain();
}
}
}
}
/* -------------------------------------------------------------------------- */
void DomainParaDiS::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);
+ // ParaDiSHome::paradis_ptr->param->timeLMSCycle = 0.;
+ // while (ParaDiSHome::paradis_ptr->param->timeLMSCycle <
+ // ParaDiSHome::paradis_ptr->param->timeMax){
+ // ParadisStep(ParaDiSHome::paradis_ptr);
+ // TimerClearAll(ParaDiSHome::paradis_ptr);
// }
- // this->paradis_ptr->cycleLMS++;
+ // ParaDiSHome::paradis_ptr->cycleLMS++;
// this->dd_container.incRelease();
}
void DomainParaDiS::performStep3() {}
void DomainParaDiS::performStep4() {}
void DomainParaDiS::performStep5() {}
/* -------------------------------------------------------------------------- */
-void DomainParaDiS::addSegments(std::vector<RefPointData<3>> &positions,
- std::vector<RefGenericElem<3>> &conn,
- std::vector<Real> &burgers,
- std::vector<Real> &normals) {
-
- UInt nb_points = positions.size();
-
- std::vector<Node_t *> created_nodes;
- std::vector<UInt> numNbrs;
-
- for (UInt i = 0; i < nb_points; ++i) {
- Node_t *newNode = GetNewNativeNode(this->paradis_ptr);
-
- auto &&pos = positions[i].position();
- newNode->x = pos[0];
- newNode->y = pos[1];
- newNode->z = pos[2];
-
- newNode->native = 1;
- newNode->constraint = 7;
-
- newNode->oldvX = 0.0;
- newNode->oldvY = 0.0;
- newNode->oldvZ = 0.0;
-
- AssignNodeToCell(this->paradis_ptr, newNode);
-
- created_nodes.push_back(newNode);
- }
-
- numNbrs.resize(created_nodes.size());
-
- UInt nb_elem = conn.size();
- for (UInt i = 0; i < nb_elem; ++i) {
- RefGenericElem<3> &el = conn[i];
- LM_ASSERT(el.conn.size() == 2,
- "we cannot add other things than segments to DD model");
- for (UInt n = 0; n < 2; ++n) {
- UInt nd = el.conn[n];
- ++numNbrs[nd];
- }
- }
-
- for (UInt i = 0; i < nb_points; ++i) {
- Node_t *newNode = created_nodes[i];
- AllocNodeArms(newNode, numNbrs[i]);
- numNbrs[i] = 0;
- }
-
- for (UInt i = 0; i < nb_elem; ++i) {
- RefGenericElem<3> &el = conn[i];
- UInt nd_src = el.conn[0];
- UInt nd_dst = el.conn[1];
-
- Node_t *srcNode = created_nodes[nd_src];
- Node_t *dstNode = created_nodes[nd_dst];
-
- UInt ¤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);
+ GetFieldPointStress(ParaDiSHome::paradis_ptr, position[0], position[1],
+ position[2], totStress);
for (unsigned int i = 0; i < n_dim; ++i) {
force[i] = 0.;
for (unsigned int j = 0; j < n_dim; ++j) {
force[i] += totStress[i][j] * normal[j];
}
}
} // CHANGE THIS?
/* -------------------------------------------------------------------------- */
/* Configure segment connecting two nodes */
void setNodeArm(Node_t *nh_src, Node_t *nh_dst, const UInt idxArmSrc,
const UInt idxArmDst, const Real *burgersVec,
const Real *gpNormVec) {
/* Source node */
nh_src->nbrTag[idxArmSrc].domainID = nh_dst->myTag.domainID;
nh_src->nbrTag[idxArmSrc].index = nh_dst->myTag.index;
nh_src->burgX[idxArmSrc] = burgersVec[0];
nh_src->burgY[idxArmSrc] = burgersVec[1];
nh_src->burgZ[idxArmSrc] = burgersVec[2];
nh_src->nx[idxArmSrc] = gpNormVec[0];
nh_src->ny[idxArmSrc] = gpNormVec[1];
nh_src->nz[idxArmSrc] = gpNormVec[2];
/* Distributive node */
nh_dst->nbrTag[idxArmDst].domainID = nh_src->myTag.domainID;
nh_dst->nbrTag[idxArmDst].index = nh_src->myTag.index;
nh_dst->burgX[idxArmDst] = -burgersVec[0];
nh_dst->burgY[idxArmDst] = -burgersVec[1];
nh_dst->burgZ[idxArmDst] = -burgersVec[2];
nh_dst->nx[idxArmDst] = gpNormVec[0];
nh_dst->ny[idxArmDst] = gpNormVec[1];
nh_dst->nz[idxArmDst] = gpNormVec[2];
}
/* -------------------------------------------------------------------------- */
/* Orthogonal projection of a node onto the glide plane */
inline void projectNodeOnGP(Node_t *nh_src, // source node
Node_t *nh_dst, // inheriting node
const UInt idxArmSrc) // arm index of source node
// connecting to the
// inheriting node
{
Real alpha = nh_src->nx[idxArmSrc] * (nh_src->x - nh_dst->x) +
nh_src->ny[idxArmSrc] * (nh_src->y - nh_dst->y) +
nh_src->nz[idxArmSrc] * (nh_src->z - nh_dst->z);
nh_dst->x += alpha * nh_src->nx[idxArmSrc];
nh_dst->y += alpha * nh_src->ny[idxArmSrc];
nh_dst->z += alpha * nh_src->nz[idxArmSrc];
}
/* -------------------------------------------------------------------------- */
/* Add hybrid segments to an existing DD line crossing an interface */
// NOTE: test version --> the method assumes that the mesh struct carries only 1
// dislocation!!
// NOTE: it is further assumed that the node of hybrid DD line in the continuum
// connecting to a node in the atomistic domain remains the nearest node
// to the interface during a simulation!
// NOTE: there is no particular check performed if "continuum" DD nodes
// previously connected to a DD node in the atomistic domain retain the
// this connectivity!!!
// NOTE:
-void DomainParaDiS::addHybSegments(
- std::vector<RefPointData<3>> &positions, // DD nodes in an atomistic domain
- std::vector<RefGenericElem<3>>
- &conn, // connectivity between DD nodes in the atomistic domain
- Real threshold) // max. distance between two adjacent DD nodes possibly
- // generating a hybrid segment
-{
- UInt nb_points = positions.size();
- Real dist;
- Real pos_ato[3], pos_con[3];
- Real burgersVec[3], gpNormVec[3];
- // Containers/Lists
- std::vector<UInt> numNbrs;
- // ParaDis
- Node_t *nh, *nh_conInt1 = NULL, *nh_conInt2 = NULL;
- // ParaDis - Iterators
- std::vector<Node_t *>::iterator nh_it;
- /* ------------------------------------------------------------------------ */
- // Delete all existing DD segments in the atomistic domain
- // NOTE: the FreeNode() function doesn't remove the segments associated with
- // the node to delete
- // --> more precisely, they are not deallocated
- // --> this is of no particular issue since one can call
- // ReallocNodeArms() in step 3)
-
- for (nh_it = created_nodes.begin(); nh_it != created_nodes.end(); ++nh_it) {
- FreeNode(paradis_ptr, (*nh_it)->myTag.index);
- }
-
- created_nodes.clear();
- /* ------------------------------------------------------------------------ */
-
- /* 1) Import new nodes in ParaDis */
-
- for (UInt i = 0; i < nb_points; ++i) {
- nh = GetNewNativeNode(this->paradis_ptr);
-
- auto &&pos = positions[i].position();
- nh->x = pos[0];
- nh->y = pos[1];
- nh->z = pos[2];
-
- nh->native = 1;
- nh->constraint =
- PINNED_NODE; // nodes in an atomistic domain remain always fixed
-
- nh->oldvX = 0.0;
- nh->oldvY = 0.0;
- nh->oldvZ = 0.0;
-
- AssignNodeToCell(this->paradis_ptr, nh);
-
- created_nodes.push_back(nh);
- }
-
- numNbrs.resize(created_nodes.size());
- /* ------------------------------------------------------------------------ */
-
- /* 2) Find nearest neighbor(s) in the continuum */
-
- // 2.1) First node
- nh = created_nodes[0];
-
- pos_ato[0] = nh->x;
- pos_ato[1] = nh->y;
- pos_ato[2] = nh->z;
-
- for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
- ++nh_it) {
- // Position of "continuum" node
- pos_con[0] = (*nh_it)->x;
- pos_con[1] = (*nh_it)->y;
- pos_con[2] = (*nh_it)->z;
- // Compute the distance between the continuum and the
- dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
- (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
- (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
- // If distance is smaller than a prescribed threshold
- // --> mark the node as a neighbor of the first "atomistic" DD node
- if (dist < threshold) {
- nh_conInt1 = *nh_it;
- numNbrs[0] = 1;
- break;
- }
- }
-
- // 2.2) Last node
- nh = created_nodes[nb_points - 1];
-
- pos_ato[0] = nh->x;
- pos_ato[1] = nh->y;
- pos_ato[2] = nh->z;
-
- for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
- ++nh_it) {
- // Position of "continuum" node
- pos_con[0] = (*nh_it)->x;
- pos_con[1] = (*nh_it)->y;
- pos_con[2] = (*nh_it)->z;
- // Compute the distance between the continuum and the
- dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
- (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
- (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
- // If distance is smaller than a prescribed threshold
- // --> mark the node as a neighbor of the last "atomistic" DD node
- if (dist < threshold) {
- nh_conInt2 = *nh_it;
- numNbrs[nb_points - 1] = 1;
- break;
- }
- }
- /* ------------------------------------------------------------------------ */
-
- /* 3) Import connectivity and allocate memory for segment info */
-
- UInt nb_elem = conn.size();
- for (UInt i = 0; i < nb_elem; ++i) {
- RefGenericElem<3> &el = conn[i];
- LM_ASSERT(el.conn.size() == 2,
- "we cannot add other things than segments to DD model");
- for (UInt n = 0; n < 2; ++n) {
- UInt nd = el.conn[n];
- ++numNbrs[nd];
- }
- }
-
- for (UInt i = 0; i < nb_points; ++i) {
- ReallocNodeArms(created_nodes[i], numNbrs[i]);
- // AllocNodeArms(newNode, numNbrs[i]);
- numNbrs[i] = 0;
- }
- /* ------------------------------------------------------------------------ */
-
- /* 4) Connect atomistic and continuum DD line */
-
- // 4.1) First interface node in the continuum connects to the atomistic DD
- // line
- if ((nh_conInt1 != NULL) && (nh_conInt2 == NULL)) {
- // Get the Burgers vector of the dislocation
- burgersVec[0] = nh_conInt1->burgX[0];
- burgersVec[1] = nh_conInt1->burgY[0];
- burgersVec[2] = nh_conInt1->burgZ[0];
- // Get the normal vector of the glide plane
- gpNormVec[0] = nh_conInt1->nx[0];
- gpNormVec[1] = nh_conInt1->ny[0];
- gpNormVec[2] = nh_conInt1->nz[0];
- // Re-allocate memory for the hybrid segment if necessary (--> should only
- // be done once before the first iteration)
- if (nh_conInt1->numNbrs < 2) {
- ReallocNodeArms(nh_conInt1, 2);
- nh_conInt1->constraint = UNCONSTRAINED;
- }
- // Establish connectivity between nodes in the continuum and the atomistic
- // domain ...
- setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
- // // ... and set homogeneous Dirichlet conditions on the other node
- // (--> u=0 <=> pinned node)
- // created_nodes[nb_points-1]->constraint = PINNED_NODE;
- // Project atomistic DD node on the glide plane
- projectNodeOnGP(nh_conInt1, created_nodes[0], 1);
-
- ++numNbrs[0];
- }
- // 4.2) Last interface node in the continuum connects to the atomistic DD line
- else if ((nh_conInt1 == NULL) && (nh_conInt2 != NULL)) {
- // Get the Burgers vector of the dislocation
- burgersVec[0] = nh_conInt2->burgX[0];
- burgersVec[1] = nh_conInt2->burgY[0];
- burgersVec[2] = nh_conInt2->burgZ[0];
- // Get the normal vector of the glide plane
- gpNormVec[0] = nh_conInt2->nx[0];
- gpNormVec[1] = nh_conInt2->ny[0];
- gpNormVec[2] = nh_conInt2->nz[0];
- // Re-allocate memory for the hybrid segment if necessary (--> should only
- // be done once before the first iteration)
- if (nh_conInt2->numNbrs < 2) {
- ReallocNodeArms(nh_conInt2, 2);
- nh_conInt2->constraint = UNCONSTRAINED;
- }
- // Establish connectivity between nodes in the continuum and the atomistic
- // domain ...
- setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
- gpNormVec);
- // // ... and set homogeneous Dirichlet conditions on the other node
- // (--> u=0 <=> pinned node)
- // created_nodes[0]->constraint = PINNED_NODE;
- // Project atomistic DD node on the glide plane
- projectNodeOnGP(nh_conInt2, created_nodes[nb_points - 1], 1);
-
- ++numNbrs[nb_points - 1];
- }
- // 4.3) Both interfaces nodes in the continuum connects to the atomistic DD
- // line
- else if ((nh_conInt1 != NULL) && (nh_conInt2 != NULL)) {
- // Get the Burgers vector of the dislocation
- burgersVec[0] = nh_conInt1->burgX[0];
- burgersVec[1] = nh_conInt1->burgY[0];
- burgersVec[2] = nh_conInt1->burgZ[0];
- // Get the normal vector of the glide plane
- gpNormVec[0] = nh_conInt1->nx[0];
- gpNormVec[1] = nh_conInt1->ny[0];
- gpNormVec[2] = nh_conInt1->nz[0];
- // Re-allocate memory for the hybrid segment if necessary (--> should only
- // be done once before the first iteration)
- if (nh_conInt1->numNbrs < 2) {
- ReallocNodeArms(nh_conInt1, 2);
- nh_conInt1->constraint = UNCONSTRAINED;
- }
- if (nh_conInt2->numNbrs < 2) {
- ReallocNodeArms(nh_conInt2, 2);
- nh_conInt2->constraint = UNCONSTRAINED;
- }
- // Establish connectivity between nodes in the continuum and the atomistic
- // domain
- setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
- setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
- gpNormVec);
- // Project atomistic DD nodes on the glide plane
- projectNodeOnGP(nh_conInt1, created_nodes[0], 1);
- projectNodeOnGP(nh_conInt1, created_nodes[nb_points - 1], 1);
-
- ++numNbrs[0];
- ++numNbrs[nb_points - 1];
- } else {
- LM_FATAL(std::endl
- << "Hybrid dislocation failed to connect!!" << std::endl);
- }
- /* ------------------------------------------------------------------------ */
-
- /* 5) Now connect the remainder of the hybrid dislocation line in the
- * atomistic domain */
-
- for (UInt i = 0; i < nb_elem; ++i) {
- // Get connectivity
- RefGenericElem<3> &el = conn[i];
- UInt idxNhSrc = el.conn[0];
- UInt idxNhDst = el.conn[1];
- UInt &idxArmSrc = numNbrs[idxNhSrc];
- UInt &idxArmDst = numNbrs[idxNhDst];
- // Update arm information for each node of segment i
- setNodeArm(created_nodes[idxNhSrc], created_nodes[idxNhDst], idxArmSrc,
- idxArmDst, burgersVec, gpNormVec);
- // Project atomistic DD node on the glide plane
- projectNodeOnGP(created_nodes[idxNhSrc], created_nodes[idxNhDst],
- idxArmSrc);
-
- ++idxArmSrc;
- ++idxArmDst;
- }
-}
+// void DomainParaDiS::addHybSegments(
+// std::vector<RefPointData<3>> &positions, // DD nodes in an atomistic
+// domain
+// std::vector<RefGenericElem<3>>
+// &conn, // connectivity between DD nodes in the atomistic domain
+// Real threshold) // max. distance between two adjacent DD nodes possibly
+// // generating a hybrid segment
+// {
+// UInt nb_points = positions.size();
+// Real dist;
+// Real pos_ato[3], pos_con[3];
+// Real burgersVec[3], gpNormVec[3];
+// // Containers/Lists
+// std::vector<UInt> numNbrs;
+// // ParaDis
+// Node_t *nh, *nh_conInt1 = NULL, *nh_conInt2 = NULL;
+// // ParaDis - Iterators
+// std::vector<Node_t *>::iterator nh_it;
+// /* ------------------------------------------------------------------------
+// */
+// // Delete all existing DD segments in the atomistic domain
+// // NOTE: the FreeNode() function doesn't remove the segments associated
+// with
+// // the node to delete
+// // --> more precisely, they are not deallocated
+// // --> this is of no particular issue since one can call
+// // ReallocNodeArms() in step 3)
+
+// for (nh_it = created_nodes.begin(); nh_it != created_nodes.end(); ++nh_it)
+// {
+// FreeNode(ParaDiSHome::paradis_ptr, (*nh_it)->myTag.index);
+// }
+
+// created_nodes.clear();
+// /* ------------------------------------------------------------------------
+// */
+
+// /* 1) Import new nodes in ParaDis */
+
+// for (UInt i = 0; i < nb_points; ++i) {
+// nh = GetNewNativeNode(ParaDiSHome::paradis_ptr);
+
+// auto &&pos = positions[i].position();
+// nh->x = pos[0];
+// nh->y = pos[1];
+// nh->z = pos[2];
+
+// nh->native = 1;
+// nh->constraint =
+// PINNED_NODE; // nodes in an atomistic domain remain always fixed
+
+// nh->oldvX = 0.0;
+// nh->oldvY = 0.0;
+// nh->oldvZ = 0.0;
+
+// AssignNodeToCell(ParaDiSHome::paradis_ptr, nh);
+
+// created_nodes.push_back(nh);
+// }
+
+// numNbrs.resize(created_nodes.size());
+// /* ------------------------------------------------------------------------
+// */
+
+// /* 2) Find nearest neighbor(s) in the continuum */
+
+// // 2.1) First node
+// nh = created_nodes[0];
+
+// pos_ato[0] = nh->x;
+// pos_ato[1] = nh->y;
+// pos_ato[2] = nh->z;
+
+// for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
+// ++nh_it) {
+// // Position of "continuum" node
+// pos_con[0] = (*nh_it)->x;
+// pos_con[1] = (*nh_it)->y;
+// pos_con[2] = (*nh_it)->z;
+// // Compute the distance between the continuum and the
+// dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
+// (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
+// (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
+// // If distance is smaller than a prescribed threshold
+// // --> mark the node as a neighbor of the first "atomistic" DD node
+// if (dist < threshold) {
+// nh_conInt1 = *nh_it;
+// numNbrs[0] = 1;
+// break;
+// }
+// }
+
+// // 2.2) Last node
+// nh = created_nodes[nb_points - 1];
+
+// pos_ato[0] = nh->x;
+// pos_ato[1] = nh->y;
+// pos_ato[2] = nh->z;
+
+// for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end();
+// ++nh_it) {
+// // Position of "continuum" node
+// pos_con[0] = (*nh_it)->x;
+// pos_con[1] = (*nh_it)->y;
+// pos_con[2] = (*nh_it)->z;
+// // Compute the distance between the continuum and the
+// dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) +
+// (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) +
+// (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2]));
+// // If distance is smaller than a prescribed threshold
+// // --> mark the node as a neighbor of the last "atomistic" DD node
+// if (dist < threshold) {
+// nh_conInt2 = *nh_it;
+// numNbrs[nb_points - 1] = 1;
+// break;
+// }
+// }
+// /* ------------------------------------------------------------------------
+// */
+
+// /* 3) Import connectivity and allocate memory for segment info */
+
+// UInt nb_elem = conn.size();
+// for (UInt i = 0; i < nb_elem; ++i) {
+// RefGenericElem<3> &el = conn[i];
+// LM_ASSERT(el.conn.size() == 2,
+// "we cannot add other things than segments to DD model");
+// for (UInt n = 0; n < 2; ++n) {
+// UInt nd = el.conn[n];
+// ++numNbrs[nd];
+// }
+// }
+
+// for (UInt i = 0; i < nb_points; ++i) {
+// ReallocNodeArms(created_nodes[i], numNbrs[i]);
+// // AllocNodeArms(newNode, numNbrs[i]);
+// numNbrs[i] = 0;
+// }
+// /* ------------------------------------------------------------------------
+// */
+
+// /* 4) Connect atomistic and continuum DD line */
+
+// // 4.1) First interface node in the continuum connects to the atomistic DD
+// // line
+// if ((nh_conInt1 != NULL) && (nh_conInt2 == NULL)) {
+// // Get the Burgers vector of the dislocation
+// burgersVec[0] = nh_conInt1->burgX[0];
+// burgersVec[1] = nh_conInt1->burgY[0];
+// burgersVec[2] = nh_conInt1->burgZ[0];
+// // Get the normal vector of the glide plane
+// gpNormVec[0] = nh_conInt1->nx[0];
+// gpNormVec[1] = nh_conInt1->ny[0];
+// gpNormVec[2] = nh_conInt1->nz[0];
+// // Re-allocate memory for the hybrid segment if necessary (--> should
+// only
+// // be done once before the first iteration)
+// if (nh_conInt1->numNbrs < 2) {
+// ReallocNodeArms(nh_conInt1, 2);
+// nh_conInt1->constraint = UNCONSTRAINED;
+// }
+// // Establish connectivity between nodes in the continuum and the
+// atomistic
+// // domain ...
+// setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
+// // // ... and set homogeneous Dirichlet conditions on the other node
+// // (--> u=0 <=> pinned node)
+// // created_nodes[nb_points-1]->constraint = PINNED_NODE;
+// // Project atomistic DD node on the glide plane
+// projectNodeOnGP(nh_conInt1, created_nodes[0], 1);
+
+// ++numNbrs[0];
+// }
+// // 4.2) Last interface node in the continuum connects to the atomistic DD
+// line
+// else if ((nh_conInt1 == NULL) && (nh_conInt2 != NULL)) {
+// // Get the Burgers vector of the dislocation
+// burgersVec[0] = nh_conInt2->burgX[0];
+// burgersVec[1] = nh_conInt2->burgY[0];
+// burgersVec[2] = nh_conInt2->burgZ[0];
+// // Get the normal vector of the glide plane
+// gpNormVec[0] = nh_conInt2->nx[0];
+// gpNormVec[1] = nh_conInt2->ny[0];
+// gpNormVec[2] = nh_conInt2->nz[0];
+// // Re-allocate memory for the hybrid segment if necessary (--> should
+// only
+// // be done once before the first iteration)
+// if (nh_conInt2->numNbrs < 2) {
+// ReallocNodeArms(nh_conInt2, 2);
+// nh_conInt2->constraint = UNCONSTRAINED;
+// }
+// // Establish connectivity between nodes in the continuum and the
+// atomistic
+// // domain ...
+// setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
+// gpNormVec);
+// // // ... and set homogeneous Dirichlet conditions on the other node
+// // (--> u=0 <=> pinned node)
+// // created_nodes[0]->constraint = PINNED_NODE;
+// // Project atomistic DD node on the glide plane
+// projectNodeOnGP(nh_conInt2, created_nodes[nb_points - 1], 1);
+
+// ++numNbrs[nb_points - 1];
+// }
+// // 4.3) Both interfaces nodes in the continuum connects to the atomistic DD
+// // line
+// else if ((nh_conInt1 != NULL) && (nh_conInt2 != NULL)) {
+// // Get the Burgers vector of the dislocation
+// burgersVec[0] = nh_conInt1->burgX[0];
+// burgersVec[1] = nh_conInt1->burgY[0];
+// burgersVec[2] = nh_conInt1->burgZ[0];
+// // Get the normal vector of the glide plane
+// gpNormVec[0] = nh_conInt1->nx[0];
+// gpNormVec[1] = nh_conInt1->ny[0];
+// gpNormVec[2] = nh_conInt1->nz[0];
+// // Re-allocate memory for the hybrid segment if necessary (--> should
+// only
+// // be done once before the first iteration)
+// if (nh_conInt1->numNbrs < 2) {
+// ReallocNodeArms(nh_conInt1, 2);
+// nh_conInt1->constraint = UNCONSTRAINED;
+// }
+// if (nh_conInt2->numNbrs < 2) {
+// ReallocNodeArms(nh_conInt2, 2);
+// nh_conInt2->constraint = UNCONSTRAINED;
+// }
+// // Establish connectivity between nodes in the continuum and the
+// atomistic
+// // domain
+// setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec);
+// setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec,
+// gpNormVec);
+// // Project atomistic DD nodes on the glide plane
+// projectNodeOnGP(nh_conInt1, created_nodes[0], 1);
+// projectNodeOnGP(nh_conInt1, created_nodes[nb_points - 1], 1);
+
+// ++numNbrs[0];
+// ++numNbrs[nb_points - 1];
+// } else {
+// LM_FATAL(std::endl
+// << "Hybrid dislocation failed to connect!!" << std::endl);
+// }
+// /* ------------------------------------------------------------------------
+// */
+
+// /* 5) Now connect the remainder of the hybrid dislocation line in the
+// * atomistic domain */
+
+// for (UInt i = 0; i < nb_elem; ++i) {
+// // Get connectivity
+// RefGenericElem<3> &el = conn[i];
+// UInt idxNhSrc = el.conn[0];
+// UInt idxNhDst = el.conn[1];
+// UInt &idxArmSrc = numNbrs[idxNhSrc];
+// UInt &idxArmDst = numNbrs[idxNhDst];
+// // Update arm information for each node of segment i
+// setNodeArm(created_nodes[idxNhSrc], created_nodes[idxNhDst], idxArmSrc,
+// idxArmDst, burgersVec, gpNormVec);
+// // Project atomistic DD node on the glide plane
+// projectNodeOnGP(created_nodes[idxNhSrc], created_nodes[idxNhDst],
+// idxArmSrc);
+
+// ++idxArmSrc;
+// ++idxArmDst;
+// }
+// }
/* -------------------------------------------------------------------------- */
UInt DomainParaDiS::getNBNodes() {
return this->container_dd.getContainerNodes().size();
}
/* -------------------------------------------------------------------------- */
Real DomainParaDiS::getShearModulus() {
- return this->paradis_ptr->param->shearModulus;
+ return ParaDiSHome::paradis_ptr->param->shearModulus;
}
/* -------------------------------------------------------------------------- */
VectorProxy<3u> DomainParaDiS::getXmin() {
- return VectorProxy<3u>(this->paradis_ptr->param->xBoundMin,
- this->paradis_ptr->param->yBoundMin,
- this->paradis_ptr->param->zBoundMin);
+ return VectorProxy<3u>(ParaDiSHome::paradis_ptr->param->xBoundMin,
+ ParaDiSHome::paradis_ptr->param->yBoundMin,
+ ParaDiSHome::paradis_ptr->param->zBoundMin);
}
/* -------------------------------------------------------------------------- */
VectorProxy<3u> DomainParaDiS::getXmax() {
- return VectorProxy<3u>(this->paradis_ptr->param->xBoundMax,
- this->paradis_ptr->param->yBoundMax,
- this->paradis_ptr->param->zBoundMax);
+ return VectorProxy<3u>(ParaDiSHome::paradis_ptr->param->xBoundMax,
+ ParaDiSHome::paradis_ptr->param->yBoundMax,
+ ParaDiSHome::paradis_ptr->param->zBoundMax);
}
/* -------------------------------------------------------------------------- */
void DomainParaDiS::freenodearm(UInt index) {
- Node_t node = *paradis_ptr->nodeKeys[index];
+ Node_t node = *ParaDiSHome::paradis_ptr->nodeKeys[index];
FreeNodeArms(&node);
}
/* -------------------------------------------------------------------------- */
void DomainParaDiS::freenode(UInt index) {
- Node_t node = *paradis_ptr->nodeKeys[index];
+ Node_t node = *ParaDiSHome::paradis_ptr->nodeKeys[index];
FreeNodeArms(&node);
}
/* -------------------------------------------------------------------------- */
bool DomainParaDiS::is_it_bound() {
// ContainerParaDiSNode::iterator dbgNodes =
// this->getContainerNodes().begin();
// UInt counter = 0;
// for(RefParaDiSNode node=dbgNodes.getFirst(); !dbgNodes.end();
// node=dbgNodes.getNext()){
// ++counter;
// }
return this->boundary;
} // CHANGE THIS?
/* -------------------------------------------------------------------------- */
Cube DomainParaDiS::getDomainBoundingBox() {
Cube c;
- c.getXmin() = VectorView<3>(this->paradis_ptr->param->minCoordinates);
- c.getXmax() = VectorView<3>(this->paradis_ptr->param->maxCoordinates);
+ c.getXmin() = VectorView<3>(ParaDiSHome::paradis_ptr->param->minCoordinates);
+ c.getXmax() = VectorView<3>(ParaDiSHome::paradis_ptr->param->maxCoordinates);
return c;
}
/* -------------------------------------------------------------------------- */
Cube DomainParaDiS::getSubDomainBoundingBox() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
void DomainParaDiS::setDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
void DomainParaDiS::setSubDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }
// void getDomainDimensions(Real *xmin, Real *xmax) {
// };
/* -------------------------------------------------------------------------- */
/* LMDESC PARADIS
This section is used to configure the ParaDiS plugin
*/
/* LMEXAMPLE
Section MultiScale AtomsUnits\\
...\\
GEOMETRY 1 CUBE BBOX -10 10 -10 10 -10 10
MODEL PARADIS dislomodel \\
...\\
endSection\\ \\
Section PARADIS:dislomodel RealUnits \\
NUMDLBSTEPS 0 \\
CONTROLFILE frank_read_src.ctrl \\
DATAFILE frank_read_src.data \\
TIMESTEP 1e-10 \\
DOMAIN_GEOMETRY 1 \\
endSection
*/
/* LMHERITATE domain_dd */
void DomainParaDiS::declareParams() {
DomainDD<ContainerNodesParaDiS, ContainerElemsParaDiS, 3>::declareParams();
/* LMKEYWORD TIMESTEP
Set the maximum time step used for ParaDiS
*/
this->parseKeyword("TIMESTEP", timestep);
/* LMKEYWORD NUMDLBSTEPS
Number of dynamic load balancing steps
*/
this->parseKeyword("NUMDLBSTEPS", numDLBCycles);
/* LMKEYWORD CONTROLFILE
Path to the ParaDiS control file
*/
this->parseKeyword("CONTROLFILE", control_filename);
/* LMKEYWORD DATAFILE
Path to the ParaDis data file
*/
this->parseKeyword("DATAFILE", data_filename);
/* LMKEYWORD MOBILITYFILE
Mobility input acquired from MD
*/
this->parseKeyword("MOBILITYFILE", mobility_filename, "");
/* LMKEYWORD BOUNDARY
Constraint boundary to mimic grain boundary
*/
this->parseTag("BOUNDARY", this->boundary, false);
this->makeItOptional("DOMAIN_GEOMETRY");
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
diff --git a/src/dd/paradis/domain_paradis.hh b/src/dd/paradis/domain_paradis.hh
index 2d82283..4aa0b3c 100644
--- a/src/dd/paradis/domain_paradis.hh
+++ b/src/dd/paradis/domain_paradis.hh
@@ -1,218 +1,189 @@
/**
* @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 ?)
- */
-
-/* -------------------------------------------------------------------------- */
-
/**
* Class DomainParaDiS
* domain to implement the plugin interface to PARADIS dislocation dynamics
*/
class DomainParaDiS
: public DomainDD<ContainerNodesParaDiS, ContainerElemsParaDiS, 3> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DomainParaDiS(LMID ID, CommGroup &group)
: LMObject(ID), DomainDD<ContainerNodesParaDiS, ContainerElemsParaDiS, 3>(
- ID, group) {
- paradis_ptr = nullptr;
- }
+ ID, group) {}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
void performStep1() override;
void performStep2() override;
void performStep3() override;
void performStep4() override;
void performStep5() override;
void imageForce(const double *position, const double *normal, double *force);
Cube getDomainBoundingBox() override;
Cube getSubDomainBoundingBox() override;
void setDomainBoundingBox(const Cube &) override;
void setSubDomainBoundingBox(const Cube &) override;
- //! compute displacement field for a set of positions
- void computeDisplacements(std::vector<Real> &positions,
- std::vector<Real> &displacements) override;
-
- //! compute displacement field for a position with a given segment
- void computeBarnett(Vector<3> &A, Vector<3> &B, Vector<3> &C, Vector<3> &burg,
- Vector<3> &dispPointU, Vector<3> &dispPointX, Real pois);
-
- void computeClosurePoint(Vector<3> &A, Vector<3> &B, Vector<3> &burg,
- Vector<3> &C);
- void computeProjectedPoint(Vector<3> &A, Vector<3> &C, Vector<3> &burg,
- Vector<3> &P);
-
//! remesh
void remesh() override;
//! intialisation stage
void init() override;
//! parse the keywords
void declareParams() override;
-
- //! 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) override;
-
UInt getNBNodes();
bool is_it_bound();
VectorProxy<3u> getXmin();
VectorProxy<3u> 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)
- override; // max. distance between two adjacent DD nodes possibly
- // generating a hybrid segment
+ // 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)
+ // override; // max. distance between two adjacent DD nodes possibly
+ // // generating a hybrid segment
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
//! set current step to a given value
void setCurrentStep(UInt ts) override {
LM_FATAL("to be implemented");
} // GUILLAUME?
//! get current step
UInt getCurrentStep() override {
LM_FATAL("to be implemented");
return 0;
}; // CHANGE THIS
//! set timestep to a given value
void setTimeStep(Real ts) override{
LM_FATAL("to be implemented")}; // CHANGE THIS
//! GetTimeStep
Real getTimeStep() override { 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]];
+ return ParaDiSHome::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; };
+ Real getPoisson() { return ParaDiSHome::paradis_ptr->param->pois; };
+ Real getMaxSeg() { return ParaDiSHome::paradis_ptr->param->maxSeg; };
+ Real getMinSeg() { return ParaDiSHome::paradis_ptr->param->minSeg; };
+ Real getBurgersMagnitude() {
+ return ParaDiSHome::paradis_ptr->param->burgMag;
+ };
+ UInt getNBpairsInfinite() {
+ return ParaDiSHome::paradis_ptr->param->numPairs;
+ };
bool getPairedIndex(int index0, int *index0_infinite) {
- UInt num_pairs = this->paradis_ptr->param->numPairs;
+ UInt num_pairs = ParaDiSHome::paradis_ptr->param->numPairs;
for (int i = 0; i < (int)num_pairs; ++i) {
Pair_t pair_ind;
- pair_ind = *this->paradis_ptr->pairKeys[i];
+ pair_ind = *ParaDiSHome::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/iterator_paradis.hh b/src/dd/paradis/iterator_paradis.hh
index 617f42a..8c9957a 100644
--- a/src/dd/paradis/iterator_paradis.hh
+++ b/src/dd/paradis/iterator_paradis.hh
@@ -1,164 +1,159 @@
/**
* @file iterator_paradis.hh
*
* @author Till Junge <till.junge@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief ParaDiS iterator over DOFs
*
* @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_ITERATOR_PARADIS_HH__
#define __LIBMULTISCALE_ITERATOR_PARADIS_HH__
/* -------------------------------------------------------------------------- */
// ParaDiS includes
#define PARALLEL
extern "C" {
#include "Home.h"
#include "Typedefs.h"
}
#undef PARALLEL
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
/**
* Class IteratorNodeParaDiS
*
*/
class IteratorNodesParaDiS : public ContainerNodesParaDiS::iterator_base {
public:
IteratorNodesParaDiS(ContainerNodesParaDiS &c, UInt index_node = 0)
- : ContainerNodesParaDiS::iterator_base(c), index_node(index_node) {
- node.setHome(c.paradis_ptr);
- node.index_node = index_node;
- };
+ : ContainerNodesParaDiS::iterator_base(c), node(index_node),
+ index_node(index_node){};
void operator++();
RefNodeParaDiS &operator*();
const RefNodeParaDiS &operator*() const;
bool operator!=(const iterator_base &) const;
bool operator==(const iterator_base &) const;
RefNodeParaDiS node;
private:
UInt index_node;
};
/* -------------------------------------------------------------------------- */
inline RefNodeParaDiS &IteratorNodesParaDiS::operator*() { return node; }
/* -------------------------------------------------------------------------- */
inline const RefNodeParaDiS &IteratorNodesParaDiS::operator*() const {
return node;
}
/* -------------------------------------------------------------------------- */
inline void IteratorNodesParaDiS::operator++() {
++index_node;
node.index_node = index_node;
}
/* -------------------------------------------------------------------------- */
inline bool IteratorNodesParaDiS::operator!=(const iterator_base &it) const {
return !(*this == it);
}
/* -------------------------------------------------------------------------- */
inline bool IteratorNodesParaDiS::operator==(const iterator_base &it) const {
return index_node == static_cast<const IteratorNodesParaDiS &>(it).index_node;
}
/* -------------------------------------------------------------------------- */
///
/// Doing elements now
///
/* -------------------------------------------------------------------------- */
class IteratorElemsParaDiS : public ContainerElemsParaDiS::iterator_base {
public:
IteratorElemsParaDiS(ContainerElemsParaDiS &c,
IteratorNodesParaDiS iterator_nodes)
- : it_nodes(iterator_nodes), current_arm(0), paradis_ptr(c.paradis_ptr) {
- elem.setHome(c.paradis_ptr);
- }
+ : it_nodes(iterator_nodes), current_arm(0) {}
void operator++() {
++this->current_arm;
if (this->current_arm == (*it_nodes).nbArms()) {
current_arm = 0;
++it_nodes;
}
auto &node = *it_nodes;
int tag1 = node.index_node;
- if (tag1 < this->paradis_ptr->newNodeKeyPtr) {
- int tag2 = paradis_ptr->nodeKeys[node.getIndex()]
+ if (tag1 < ParaDiSHome::paradis_ptr->newNodeKeyPtr) {
+ int tag2 = ParaDiSHome::paradis_ptr->nodeKeys[node.getIndex()]
->nbrTag[this->current_arm]
.index;
if (tag1 > tag2)
++(*this);
}
};
RefElemParaDiS &operator*() {
elem.setIndex((*it_nodes).getIndex(), current_arm);
return elem;
};
bool operator!=(const iterator_base &other) const {
return !((*this) == other);
};
bool operator==(const iterator_base &other) const {
if (static_cast<const IteratorElemsParaDiS &>(other).it_nodes !=
this->it_nodes)
return false;
return current_arm ==
static_cast<const IteratorElemsParaDiS &>(other).current_arm;
};
private:
IteratorNodesParaDiS it_nodes;
UInt current_arm;
RefElemParaDiS elem;
- Home_t *paradis_ptr;
};
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_ITERATOR_PARADIS_HH__ */
diff --git a/src/dd/paradis/paradis_home.cc b/src/dd/paradis/paradis_home.cc
new file mode 100644
index 0000000..76ff14b
--- /dev/null
+++ b/src/dd/paradis/paradis_home.cc
@@ -0,0 +1,9 @@
+#include "lm_common.hh"
+#include <mpi.h>
+#include "paradis_home.hh"
+
+__BEGIN_LIBMULTISCALE__
+
+Home_t *ParaDiSHome::paradis_ptr = nullptr;
+
+__END_LIBMULTISCALE__
diff --git a/src/dd/paradis/paradis_home.hh b/src/dd/paradis/paradis_home.hh
new file mode 100644
index 0000000..82b2748
--- /dev/null
+++ b/src/dd/paradis/paradis_home.hh
@@ -0,0 +1,25 @@
+#ifndef __LIBMULTISCALE_PARADIS_HOME_HH__
+#define __LIBMULTISCALE_PARADIS_HOME_HH__
+/* -------------------------------------------------------------------------- */
+#define PARALLEL
+extern "C" {
+#include "Home.h"
+#include "Typedefs.h"
+#undef X
+#undef Y
+}
+#undef PARALLEL
+
+/* -------------------------------------------------------------------------- */
+
+__BEGIN_LIBMULTISCALE__
+
+struct ParaDiSHome {
+ static Home_t *paradis_ptr;
+};
+
+__END_LIBMULTISCALE__
+
+/* -------------------------------------------------------------------------- */
+
+#endif // __LIBMULTISCALE_PARADIS_HOME_HH__
diff --git a/src/dd/paradis/ref_elem_paradis.hh b/src/dd/paradis/ref_elem_paradis.hh
index 25508ae..7eb0576 100644
--- a/src/dd/paradis/ref_elem_paradis.hh
+++ b/src/dd/paradis/ref_elem_paradis.hh
@@ -1,112 +1,113 @@
/**
* @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_ELEM_PARADIS_HH__
#define __LIBMULTISCALE_REF_ELEM_PARADIS_HH__
/* -------------------------------------------------------------------------- */
// ParaDis includes
#define PARALLEL
#include <mpi.h>
extern "C" {
#include "Home.h"
#include "Typedefs.h"
}
#undef PARALLEL
/* -------------------------------------------------------------------------- */
#include "ref_elem_dd.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
/**
* Class RefElemParaDiS
*
*/
class RefElemParaDiS : public RefDDElem<RefElemParaDiS> {
public:
static const UInt Dim = 3;
- RefElemParaDiS() : paradis_ptr(nullptr) {}
+ RefElemParaDiS() {}
~RefElemParaDiS() {}
bool operator==(RefElemParaDiS &ref) { LM_TOIMPLEMENT; };
- void setHome(Home_t *paradis_ptr) { this->paradis_ptr = paradis_ptr; };
VectorProxy<3> burgers() {
return VectorProxy<3>(
- paradis_ptr->nodeKeys[this->index_node]->burgX[arm_number],
- paradis_ptr->nodeKeys[this->index_node]->burgY[arm_number],
- paradis_ptr->nodeKeys[this->index_node]->burgZ[arm_number]);
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->burgX[arm_number],
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->burgY[arm_number],
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]
+ ->burgZ[arm_number]);
}
VectorProxy<3> normal() {
return VectorProxy<3>(
- paradis_ptr->nodeKeys[this->index_node]->nx[arm_number],
- paradis_ptr->nodeKeys[this->index_node]->ny[arm_number],
- paradis_ptr->nodeKeys[this->index_node]->nz[arm_number]);
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->nx[arm_number],
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->ny[arm_number],
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->nz[arm_number]);
}
void setIndex(UInt index, UInt arm) {
index_node = index;
arm_number = arm;
}
std::vector<UInt> localIndexes() override {
std::vector<UInt> nodes;
- auto tag1 = paradis_ptr->nodeKeys[this->index_node]->myTag.index;
- auto tag2 =
- paradis_ptr->nodeKeys[this->index_node]->nbrTag[arm_number].index;
+ auto tag1 =
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->myTag.index;
+ auto tag2 = ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]
+ ->nbrTag[arm_number]
+ .index;
nodes.push_back(tag1);
nodes.push_back(tag2);
return nodes;
};
private:
- Home_t *paradis_ptr;
UInt index_node;
UInt arm_number;
};
/* -------------------------------------------------------------------------- */
inline std::ostream &operator<<(std::ostream &os, RefElemParaDiS &ref) {
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_REF_ELEM_PARADIS_HH__ */
diff --git a/src/dd/paradis/ref_node_paradis.hh b/src/dd/paradis/ref_node_paradis.hh
index e1966c5..6753bed 100644
--- a/src/dd/paradis/ref_node_paradis.hh
+++ b/src/dd/paradis/ref_node_paradis.hh
@@ -1,329 +1,354 @@
/**
* @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 ComparatorRefNodeParaDiS;
/* -------------------------------------------------------------------------- */
/**
* Class RefNodeParaDiS
*
*/
class RefNodeParaDiS : public RefDDNode<3, RefNodeParaDiS> {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
//! generic container of the original model container
using Domain = DomainParaDiS;
//! node comparator class to be usable in maps
using RefComparator = ComparatorRefNodeParaDiS;
static const UInt Dim = 3;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- RefNodeParaDiS() : fakestress(0) {}
+
+ RefNodeParaDiS() : fakestress(0) {
+ // Request creation of a new node
+ Node_t *newNode = GetNewNativeNode(ParaDiSHome::paradis_ptr);
+ AssignNodeToCell(ParaDiSHome::paradis_ptr, newNode);
+ this->index_node = newNode->myTag.index;
+ }
RefNodeParaDiS(int index) : fakestress(0) { index_node = index; }
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
bool operator==(RefNodeParaDiS &ref) {
return (ref.index_node == index_node);
};
VectorProxy<3> position0() { return this->position(); }
VectorProxy<3> position() {
- return VectorProxy<3>(paradis_ptr->nodeKeys[this->index_node]->x,
- paradis_ptr->nodeKeys[this->index_node]->y,
- paradis_ptr->nodeKeys[this->index_node]->z);
+ return VectorProxy<3>(
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->x,
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->y,
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->z);
}
VectorProxy<3> velocity() {
- return VectorProxy<3>(paradis_ptr->nodeKeys[this->index_node]->vX,
- paradis_ptr->nodeKeys[this->index_node]->vY,
- paradis_ptr->nodeKeys[this->index_node]->vZ);
+ return VectorProxy<3>(
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->vX,
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->vY,
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->vZ);
}
VectorProxy<3> force() {
- return VectorProxy<3>(paradis_ptr->nodeKeys[this->index_node]->fX,
- paradis_ptr->nodeKeys[this->index_node]->fY,
- paradis_ptr->nodeKeys[this->index_node]->fZ);
+ return VectorProxy<3>(
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->fX,
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->fY,
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->fZ);
}
- UInt nbArms() const { return paradis_ptr->nodeKeys[index_node]->numNbrs; }
+ UInt nbArms() const {
+ return ParaDiSHome::paradis_ptr->nodeKeys[index_node]->numNbrs;
+ }
// Real &prev_position(unsigned int i) {
// switch (i) {
// case 0:
- // return paradis_ptr->nodeKeys[this->index_node]->oldx;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldx;
// break;
// case 1:
- // return paradis_ptr->nodeKeys[this->index_node]->oldy;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldy;
// break;
// case 2:
- // return paradis_ptr->nodeKeys[this->index_node]->oldz;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldz;
// break;
// }
// Real *p = NULL;
// return *p;
// }
// Real &getNeighborNodeCoord(unsigned int n, unsigned int i) {
- // UInt my_domain = paradis_ptr->myDomain;
- // Node_t current = *paradis_ptr->nodeKeys[this->index_node];
+ // UInt my_domain = ParaDiSHome::paradis_ptr->myDomain;
+ // Node_t current = *ParaDiSHome::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;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[idx_for_neighbor]->x;
// break;
// case 1:
- // return paradis_ptr->nodeKeys[idx_for_neighbor]->y;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[idx_for_neighbor]->y;
// break;
// case 2:
- // return paradis_ptr->nodeKeys[idx_for_neighbor]->z;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[idx_for_neighbor]->z;
// break;
// }
// } else {
// switch (i) {
// case 0:
- // return paradis_ptr->remoteDomainKeys[tag_domain]
+ // return ParaDiSHome::paradis_ptr->remoteDomainKeys[tag_domain]
// ->nodeKeys[idx_for_neighbor]
// ->x;
// break;
// case 1:
- // return paradis_ptr->remoteDomainKeys[tag_domain]
+ // return ParaDiSHome::paradis_ptr->remoteDomainKeys[tag_domain]
// ->nodeKeys[idx_for_neighbor]
// ->y;
// break;
// case 2:
- // return paradis_ptr->remoteDomainKeys[tag_domain]
+ // return ParaDiSHome::paradis_ptr->remoteDomainKeys[tag_domain]
// ->nodeKeys[idx_for_neighbor]
// ->z;
// break;
// }
// }
// Real *p = NULL;
// return *p;
// }
// int &getNeighborNodeIndex0(unsigned int n) {
- // Node_t current = *paradis_ptr->nodeKeys[this->index_node];
+ // Node_t current = *ParaDiSHome::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];
+ // Node_t current = *ParaDiSHome::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];
+ // UInt my_domain = ParaDiSHome::paradis_ptr->myDomain;
+ // Node_t current = *ParaDiSHome::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;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[idx_for_neighbor]->slave;
// else
- // return paradis_ptr->remoteDomainKeys[tag_domain]
+ // return ParaDiSHome::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];
+ // UInt my_domain = ParaDiSHome::paradis_ptr->myDomain;
+ // Node_t current = *ParaDiSHome::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;
+ // return ParaDiSHome::paradis_ptr->nodeKeys[idx_for_neighbor]->numNbrs;
// else
- // return paradis_ptr->remoteDomainKeys[tag_domain]
+ // return ParaDiSHome::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;
+ return ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->myTag.index0;
}
int &getIndex() {
- return paradis_ptr->nodeKeys[this->index_node]->myTag.index;
+ return ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->myTag.index;
}
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);
+ // return
+ // UInt(ParaDiSHome::ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->numNbrs);
// }
int getSlave() {
- return this->paradis_ptr->nodeKeys[this->index_node]->slave;
+ return ParaDiSHome::ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]
+ ->slave;
}
//! constraint tag:
int getConstraint() {
- return this->paradis_ptr->nodeKeys[this->index_node]->constraint;
+ return ParaDiSHome::ParaDiSHome::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 slaving() {
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->slave = 1;
+ }
+ void unslaving() {
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->slave = 0;
+ }
- void constrain() { paradis_ptr->nodeKeys[this->index_node]->constraint = 7; }
+ void constrain() {
+ ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->constraint = 7;
+ }
void unconstrain() {
- paradis_ptr->nodeKeys[this->index_node]->constraint = 0;
+ ParaDiSHome::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;
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldvX = oldv;
// return;
// case 1:
- // paradis_ptr->nodeKeys[this->index_node]->oldvY = oldv;
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldvY = oldv;
// return;
// case 2:
- // paradis_ptr->nodeKeys[this->index_node]->oldvZ = oldv;
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldvZ = oldv;
// return;
// }
// return;
// }
// void release(UInt nbr_index) {
- // paradis_ptr->nodeKeys[this->index_node]->constraint = 0;
+ // ParaDiSHome::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;
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->nbrTag[0].index;
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldvX =
+ // 0.5 * (ParaDiSHome::paradis_ptr->nodeKeys[index_for_neighbor]->oldvX
+ // +
+ // ParaDiSHome::paradis_ptr->nodeKeys[nbr_index]->oldvX);
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldvY =
+ // 0.5 * (ParaDiSHome::paradis_ptr->nodeKeys[index_for_neighbor]->oldvY
+ // +
+ // ParaDiSHome::paradis_ptr->nodeKeys[nbr_index]->oldvY);
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->oldvZ =
+ // 0.5 * (ParaDiSHome::paradis_ptr->nodeKeys[index_for_neighbor]->oldvZ
+ // +
+ // ParaDiSHome::paradis_ptr->nodeKeys[nbr_index]->oldvZ);
+
+ // InsertArm(ParaDiSHome::paradis_ptr,
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node],
+ // &ParaDiSHome::paradis_ptr->nodeKeys[nbr_index]->myTag,
+ // -1.0 *
+ // (ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->burgX[0]),
+ // -1.0 *
+ // (ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->burgY[0]),
+ // -1.0 *
+ // (ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->burgZ[0]),
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->nx[0],
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->ny[0],
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->nz[0], 0,
+ // 1);
+ // ParaDiSHome::paradis_ptr->nodeKeys[this->index_node]->numNbrs = 2;
// }
friend class ContainerNodesParaDiS;
friend class IteratorNodesParaDiS;
friend class IteratorElemsParaDiS;
private:
- Home_t *paradis_ptr;
unsigned int index_node;
unsigned int index0_node;
Real fakestress; // This is a temprary fix to make the dumpers work: FIXMETILL
};
/* -------------------------------------------------------------------------- */
class ComparatorRefNodeParaDiS {
public:
bool operator()(RefNodeParaDiS &r1, RefNodeParaDiS &r2) const {
// return r1.index_atome < r2.index_atome;
LM_TOIMPLEMENT;
}
bool operator()(const RefNodeParaDiS &r1, const RefNodeParaDiS &r2) const {
LM_TOIMPLEMENT;
}
};
/* -------------------------------------------------------------------------- */
inline std::ostream &operator<<(std::ostream &os, RefNodeParaDiS &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_elem_dd.hh b/src/dd/ref_elem_dd.hh
index b8d1eac..df47000 100644
--- a/src/dd/ref_elem_dd.hh
+++ b/src/dd/ref_elem_dd.hh
@@ -1,67 +1,69 @@
/**
* @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_ELEM_DD_HH__
#define __LIBMULTISCALE_REF_ELEM_DD_HH__
/* ----------------------------------------------------------------- */
#include "lm_common.hh"
#include "ref_element.hh"
/* ----------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* ----------------------------------------------------------------- */
/**
* Class RefDDElem
*
*/
template <typename DaughterClass>
class RefDDElem : public RefElement<DaughterClass> {
public:
auto &&burgers() { return static_cast<DaughterClass *>(this)->burgers(); };
auto &&normal() { return static_cast<DaughterClass *>(this)->normal(); };
std::vector<UInt> globalIndexes() override { LM_TOIMPLEMENT; };
template <FieldType ftype> Real field() { LM_TOIMPLEMENT; };
template <FieldType ftype> Real field(UInt i) { LM_TOIMPLEMENT; };
+
+ void setConnectivity(std::vector<UInt> &conn) { LM_TOIMPLEMENT; };
};
/* --------------------------------------------------------------------------
*/
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_REF_ELEM_DD_HH__ */
diff --git a/src/filter/compute_add_ddloop.cc b/src/filter/compute_add_ddloop.cc
index 63d0f5a..ed5e967 100644
--- a/src/filter/compute_add_ddloop.cc
+++ b/src/filter/compute_add_ddloop.cc
@@ -1,147 +1,145 @@
/**
* @file compute_add_ddloop.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Jul 30 09:55:56 2014
*
* @brief This allows to add a loop to a DD model from a set of points
*
* @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 "compute_add_ddloop.hh"
#include "communicator.hh"
#include "compute_point_set.hh"
#include "domain_dd_interface.hh"
#include "domain_multiscale.hh"
#include "factory_multiscale.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "ref_point_data.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename ContIn, typename ContOut>
-void ComputeAddDDLoop::build(ContIn &cont_in, ContOut &cont_out) {
-
- auto dom_ptr = DomainMultiScale::getManager().getSharedObject(this->dom);
-
- std::vector<Real> v_burgers;
- std::vector<Real> v_glide_normals;
-
- ComputePointSet v_nodes(this->getID() + ":loop-nodes");
- LM_TOIMPLEMENT;
-
- // v_nodes.build(cont);
- // UInt npoints =
- // v_nodes.getCastedOutput<ContainerArray<RefPointData<Dim>>>().size();
- LM_TOIMPLEMENT;
- // std::vector<RefGenericElem<Cont::Dim>> v_connectivity;
- LM_TOIMPLEMENT;
- //for (UInt i = 0; i < npoints; ++i) {
- // RefGenericElem<Cont::Dim> tmp;
- // tmp.conn.push_back(i);
- // tmp.conn.push_back(i + 1);
- // v_connectivity.push_back(tmp);
- LM_TOIMPLEMENT;
-
- for (UInt k = 0; k < 3; ++k) {
- v_burgers.push_back(burgers[k]);
- v_glide_normals.push_back(glide_normal[k]);
- }
-// }
-
- // v_connectivity[npoints - 1].conn[1] = 0;
- LM_TOIMPLEMENT;
-
- // if (DomainDDInterface *ptr = dynamic_cast<DomainDDInterface *>(dom_ptr)) {
- // LM_TOIMPLEMENT;
- // // ptr->addSegments(v_nodes, v_connectivity, v_burgers,v_glide_normals);
- // } else {
- // LM_FATAL("Invalid domain input (" << this->dom << ")"
- // << " for the compute " <<
- // this->getID());
- // }
+std::enable_if_t<ContIn::Dim == 3 and ContOut::Dim == 3>
+ComputeAddDDLoop::build(ContIn &cont_in, ContOut &cont_out) {
+
+ auto &&nodes = cont_in.getContainerNodes();
+ auto &&elems = cont_in.getContainerElems();
+ // UInt nb_points = cont_in.size();
+
+ using RefNode =
+ std::remove_reference_t<decltype(cont_out.getContainerNodes().get(0))>;
+ using RefElem =
+ std::remove_reference_t<decltype(cont_out.getContainerElems().get(0))>;
+
+ auto new_mesh =
+ ContainerMesh<ContainerArray<RefNode>, ContainerArray<RefElem>>(
+ this->getID() + ":created_mesh");
+ auto &&new_nodes = new_mesh.getContainerNodes();
+ auto &&new_elems = new_mesh.getContainerElems();
+
+ // adding the points
+ for (auto &&n : nodes) {
+ RefNode new_node;
+ new_node.position() = n.position();
+ new_nodes.add(new_node);
+ }
+
+ // adding the points
+ for (auto &&el : elems) {
+ LM_ASSERT(el.conn.size() == 2,
+ "we cannot add other things than segments to DD model");
+
+ RefElem new_elem;
+ auto new_connectivity = el.conn;
+ new_connectivity[0] = new_nodes[el.conn[0]].getIndex();
+ new_connectivity[1] = new_nodes[el.conn[1]].getIndex();
+ new_elem.setConnectivity(new_connectivity);
+ new_elems.add(new_elem);
+ }
}
/* -------------------------------------------------------------------------- */
/* LMDESC ADD_DDLOOP
Add a loop into a dislocation dynamics domain from a set of points \\
*/
/* LMEXAMPLE
- COMPUTE ddLoop ADD_DDLOOP INPUT input=points INPUT output=dd_part BURGERS 0 0 1
+ COMPUTE ddLoop ADD_DDLOOP INPUT input=points INPUT output=dd_part BURGERS 0 0
+ 1
GLIDE_PLANE 0 0 1\\
*/
void ComputeAddDDLoop::declareParams() {
ComputeInterface::declareParams();
// /* LMKEYWORD DOMAIN
// Specify the domain ID to be used to generate the displacement data
// */
// this->parseKeyword("DOMAIN", this->dom);
/* LMKEYWORD BURGERS
Specify the burgers to use for the loop segments
*/
- this->parseVectorKeyword("BURGERS", spatial_dimension, this->burgers);
+ // this->parseVectorKeyword("BURGERS", spatial_dimension, this->burgers);
/* LMKEYWORD GLIDE_PLANE
Specify the glide plane normal
*/
- this->parseVectorKeyword("GLIDE_PLANE", spatial_dimension,
- this->glide_normal);
+ // this->parseVectorKeyword("GLIDE_PLANE", spatial_dimension,
+ // this->glide_normal);
};
/* -------------------------------------------------------------------------- */
ComputeAddDDLoop::ComputeAddDDLoop(const std::string &name) : LMObject(name) {
ActionManager::getManager().addObject(*this);
this->removeInput("input1");
this->createInput("input");
this->createInput("output");
}
/* -------------------------------------------------------------------------- */
ComputeAddDDLoop::~ComputeAddDDLoop() {}
/* -------------------------------------------------------------------------- */
void ComputeAddDDLoop::compute_make_call() {
call_compute(*this, [&](auto &... args) { this->build(args...); },
this->getInput("input"), this->getInput("output"));
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
diff --git a/src/filter/compute_add_ddloop.hh b/src/filter/compute_add_ddloop.hh
index 2b5a08c..323e643 100644
--- a/src/filter/compute_add_ddloop.hh
+++ b/src/filter/compute_add_ddloop.hh
@@ -1,61 +1,68 @@
/**
* @file compute_add_ddloop.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Jul 11 15:53:17 2014
*
* @brief This allows to add a loop to a DD model from a set of points
*
* @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_COMPUTE_ADD_DDLOOP_HH__
#define __LIBMULTISCALE_COMPUTE_ADD_DDLOOP_HH__
/* -------------------------------------------------------------------------- */
#include "compute_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class ComputeAddDDLoop : public ComputeInterface {
public:
DECLARE_COMPUTE(ComputeAddDDLoop);
DECLARE_MULTIPLE_INPUT(std::tuple<MESH>, std::tuple<DD>);
//! most generic build function
template <typename ContIn, typename ContOut>
- void build(ContIn &cont_in, ContOut &cont_out);
+ std::enable_if_t<ContIn::Dim == 3 and ContOut::Dim == 3>
+ build(ContIn &cont_in, ContOut &cont_out);
+
+ template <typename ContIn, typename ContOut>
+ std::enable_if_t<ContIn::Dim != 3 or ContOut::Dim != 3>
+ build(ContIn &cont_in, ContOut &cont_out) {
+ LM_TOIMPLEMENT;
+ }
private:
//! the domain ID to be used to insert new points
- LMID dom;
+ // LMID dom;
//! the burgers vectors
Quantity<Length, 3> burgers;
//! the glide plane normals
Real glide_normal[3];
};
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_COMPUTE_ADD_DDLOOP_HH__ */
diff --git a/src/filter/compute_mesh.cc b/src/filter/compute_mesh.cc
index 096dace..85f32d6 100644
--- a/src/filter/compute_mesh.cc
+++ b/src/filter/compute_mesh.cc
@@ -1,113 +1,117 @@
/**
* @file compute_mesh.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief This computes generate an independent mesh from points and
* connectivities provided as computes
*
* @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 "communicator.hh"
#include "compute_interface.hh"
#include "compute_point_set.hh"
#include "container_mesh.hh"
#include "factory_multiscale.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_dumper.hh"
#include "lib_filter.hh"
#include "lib_md.hh"
#include "lib_stimulation.hh"
#include "lm_common.hh"
/* ----------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* ----------------------------------------------------------------------- */
ComputeMesh::ComputeMesh(const std::string &name)
- : LMObject(name), FilterInterface(){};
+ : LMObject(name), FilterInterface(){
+ this->createInput("nodes");
+ this->createInput("connectivity");
+ this->createOutput("output");
+};
/* ----------------------------------------------------------------------- */
ComputeMesh::~ComputeMesh(){};
/* ----------------------------------------------------------------------- */
template <typename Nodes, typename Connectivity>
void ComputeMesh::build(Nodes &nodes, Connectivity &connectivity) {
constexpr UInt Dim = Nodes::Dim;
using MeshType = ContainerGenericMesh<Dim>;
auto &output = allocAndGetCastedOutput<MeshType>(this->getID());
// copy the nodes
auto &nds = output.getContainerNodes();
nds = nodes;
auto &els = output.getContainerElems();
els.clear();
auto nodes_per_elem = connectivity.getDim();
auto size = connectivity.size();
auto nb_elem = size / nodes_per_elem;
RefGenericElem<Dim> tmp_el;
- Eigen::VectorXd conn(nodes_per_elem);
for (UInt i = 0; i < nb_elem; ++i) {
+ Eigen::VectorXd conn(nodes_per_elem);
for (UInt j = 0; j < nodes_per_elem; ++j) {
- conn[i] = connectivity[i * nodes_per_elem + j];
+ conn[j] = UInt(connectivity[i * nodes_per_elem + j]);
+ tmp_el.setConnectivity(conn);
}
- tmp_el.setConnectivity(conn);
els.add(tmp_el);
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC MESH
This compute generates an independent mesh from an input
containing points and a compute which contains connectivities.
*/
/* LMEXAMPLE COMPUTE pset MESH INPUT coords INPUT conn */
-void ComputeMesh::declareParams() { this->declareParams(); };
+void ComputeMesh::declareParams() { FilterInterface::declareParams(); };
/* -------------------------------------------------------------------------- */
void ComputeMesh::compute_make_call() {
call_compute(*this, [&](auto &... args) { this->build(args...); },
this->getInput("nodes"), this->getInput("connectivity"));
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__