diff --git a/src/dd/container_dd.hh b/src/dd/container_dd.hh index 2c102d6..32c7625 100644 --- a/src/dd/container_dd.hh +++ b/src/dd/container_dd.hh @@ -1,149 +1,150 @@ /** * @file container_dd.hh * * @author Guillaume Anciaux * * @date Thu Jul 24 14:21:58 2014 * * @brief Generic container for DD 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 . * */ #ifndef __LIBMULTISCALE_CONTAINER_DD_HH__ #define __LIBMULTISCALE_CONTAINER_DD_HH__ /* -------------------------------------------------------------------------- */ #include "container.hh" #include "container_array.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ template class ContainerDD : public Container_base { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: typedef ContNodes ContainerNodes; typedef ContSources ContainerSources; typedef ContObstacles ContainerObstacles; typedef typename ContNodes::iterator iterator; typedef typename ContSources::iterator iteratorSources; typedef typename ContObstacles::iterator iteratorObstacles; typedef typename ContNodes::Ref Ref; typedef typename ContSources::Ref RefSource; typedef typename ContObstacles::Ref RefObstacle; typedef ContainerDD, ContainerArray, ContainerArray> ContainerSubset; static constexpr UInt Dim = Ref::Dim; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: iterator begin(DOFType dt = dt_local) { return container_nodes->begin(dt); }; iterator end(DOFType dt = dt_local){ return container_nodes->end(dt); }; ContainerDD(ContNodes &cN, ContSources &cS, ContObstacles &cO); ContainerDD(); virtual ~ContainerDD(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void incRelease() { - Container_base::incRelease(); - container_nodes->incRelease(); - container_sources->incRelease(); - container_obstacles->incRelease(); + LM_TOIMPLEMENT; + // Container_base::incRelease(); + // container_nodes->incRelease(); + // container_sources->incRelease(); + // container_obstacles->incRelease(); }; void setRelease(UInt rel) { Container_base::setRelease(rel); container_nodes->setRelease(rel); container_sources->setRelease(rel); container_obstacles->setRelease(rel); }; operator ContainerNodes &() { LM_FATAL( "silent conversion from ContainerDD to ContainerNode is not allowed"); return *container_nodes; } void empty() { container_nodes->empty(); container_sources->empty(); container_obstacles->empty(); } std::vector &getArray() { return getContainerNodes().getArray(); } UInt size(DOFType dt = dt_local){ return getContainerNodes().size(); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: ContNodes &getContainerNodes() { return *container_nodes; }; ContSources &getContainerSources() { return *container_sources; }; ContObstacles &getContainerObstacles() { return *container_obstacles; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: ContNodes *container_nodes; ContSources *container_sources; ContObstacles *container_obstacles; }; /* -------------------------------------------------------------------------- */ template ContainerDD::ContainerDD( ContNodes &cN, ContSources &cS, ContObstacles &cO) : container_nodes(&cN), container_sources(&cS), container_obstacles(&cO) {} /* -------------------------------------------------------------------------- */ template ContainerDD::ContainerDD() : container_nodes(new ContNodes()), container_sources(new ContSources()), container_obstacles(new ContObstacles()) {} /* -------------------------------------------------------------------------- */ template ContainerDD::~ContainerDD() {} /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_CONTAINER_DD_HH__ */ diff --git a/src/dd/domain_dd.hh b/src/dd/domain_dd.hh index 4b3b16b..58cf29d 100644 --- a/src/dd/domain_dd.hh +++ b/src/dd/domain_dd.hh @@ -1,163 +1,159 @@ /** * @file domain_dd.hh * * @author Guillaume Anciaux * * @date Fri Jul 11 15:47:44 2014 * * @brief Generic DD 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 . * */ #ifndef __LIBMULTISCALE_DOMAIN_DD_HH__ #define __LIBMULTISCALE_DOMAIN_DD_HH__ /* -------------------------------------------------------------------------- */ #include "component_libmultiscale.hh" #include "container_dd.hh" #include "domain_dd_interface.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /** * Class DomainDD * */ template class DomainDD : public DomainDDInterface { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DomainDD(LMID ID, CommGroup GID) - : DomainDDInterface(ID, GID), dd_container(nodes, sources, obstacles){}; + : DomainDDInterface(GID), dd_container(nodes, sources, obstacles){}; virtual ~DomainDD(){}; + void compute_make_call(){}; + /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ typedef ContNodes ContainerNodes; typedef typename ContainerNodes::iterator IteratorNodes; typedef typename ContainerNodes::Ref RefNode; typedef typename ContainerNodes::Ref RefPoint; typedef ContSources ContainerSources; typedef typename ContainerSources::iterator IteratorSources; typedef typename ContainerSources::Ref RefSource; typedef ContObstacles ContainerObstacles; typedef typename ContainerObstacles::iterator IteratorObstacles; typedef typename ContainerObstacles::Ref RefObstacle; typedef ContainerDD ContainerPoints; - typedef typename ContainerPoints::ContainerSubset ContainerSubset; - //**************************************************************** - // TEMPORARY - //**************************************************************** - typedef ContainerPoints Output; - static const COutputType output_type = OUTPUT_REF; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ ContainerPoints &getContainer(); ContainerPoints &getOutput() { return this->getContainer(); }; //! return the container of nodes ContainerNodes &getContainerNodes(); //! return the container of sources virtual ContainerSources &getContainerSources(); //! return the container of obstacles virtual ContainerObstacles &getContainerObstacles(); //! return a reference to the geometry Geometry &getGeom(); //! return the potential energy Real getEpot() { return 0; }; void build(){}; virtual void declareParams(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ //! restart from xml file function void readXMLFile(const std::string &filename) { LM_FATAL("not yet implemented"); } // //! connect and build the adapted object // void connect(FactoryMultiScale &f); // //! connect and build the adapted object // void connect(FactoryMultiScale &f); /* ------------------------------------------------------------------------ */ /* Functions fo AMELCG */ /* ------------------------------------------------------------------------ */ //! return product of force by descent direction (for AMELCG) Real getFdotDir() { return 0; }; //! return max of forces (for AMELCG) Real getFMax() { return 0; }; //! return max of direction vector (for AMELCG) Real getDirMax() { return 0; }; //! return norm 2 of forces (for AMELCG) Real getFNorm2() { return 0; }; //! return stuff for AMELCG Real getFdotOldF() { return 0; }; //! update direction AMELCG void updateDirection(Real beta){}; //! return stuff for AMELCG void saveForceVector(){}; //! displace in descent direction with given magnitude void displaceTowardsDirection(Real alpha){}; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! domain geometry LMID geom; //! constrained geometry LMID geom_constrained; ContainerNodes nodes; ContainerSources sources; ContainerObstacles obstacles; ContainerPoints dd_container; }; /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ #endif /* __LIBMULTISCALE_DOMAIN_DD_HH__ */ diff --git a/src/dd/paradis/domain_paradis.cc b/src/dd/paradis/domain_paradis.cc index 4df9176..6685747 100644 --- a/src/dd/paradis/domain_paradis.cc +++ b/src/dd/paradis/domain_paradis.cc @@ -1,913 +1,915 @@ /** * @file domain_paradis.cc * * @author Guillaume Anciaux * @author Till Junge * @author Moseley Philip Arthur * @author Jaehyun Cho * @author Max Hodapp * * @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 . * */ /* 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 #endif } #undef PARALLEL /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ void DomainPARADIS::remesh() { RemeshAfterCoupling(this->paradis_ptr); // //ParadisStep.c from line 336 to 418 // InitTopologyExemptions(this->paradis_ptr); // ClearOpList(this->paradis_ptr); // SortNodesForCollision(this->paradis_ptr); // SplitMultiNodes(this->paradis_ptr); // CrossSlip(this->paradis_ptr); // HandleCollisions(this->paradis_ptr); // ProximityCollisions(this->paradis_ptr); // #ifdef PARALLEL // #ifdef PARADIS_IN_LIBMULTISCALE // CommSendRemesh(this->paradis_ptr); // #endif // #endif // FixRemesh(this->paradis_ptr); // Remesh(this->paradis_ptr); // Migrate(this->paradis_ptr); // #ifdef PARALLEL // #ifdef PARADIS_IN_LIBMULTISCALE // RecycleGhostNodes(this->paradis_ptr); // #endif // #endif // SortNativeNodes(this->paradis_ptr); // #ifdef PARALLEL // #ifdef PARADIS_IN_LIBMULTISCALE // CommSendGhosts(this->paradis_ptr); // #endif // #endif } /* -------------------------------------------------------------------------- */ void DomainPARADIS::computeDisplacements(std::vector &positions, std::vector &displacements) { if (positions.empty()) { displacements.clear(); return; } if (!this->paradis_ptr) { return; } // Pass data to Paradis. unsigned int size = positions.size(); LM_ASSERT(size % 3 == 0, "paradis computeDisplacements requires 3d."); if (displacements.size() != size) displacements.resize(size); this->paradis_ptr->dispPointCount = size / 3; this->paradis_ptr->dispPointX = &positions[0]; this->paradis_ptr->dispPointU = &displacements[0]; // Calculate the displacements. PointDisplacements(this->paradis_ptr); } /* -------------------------------------------------------------------------- */ void DomainPARADIS::computeBarnett(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() { Communicator &comm = Communicator::getCommunicator(); if (!comm.amIinGroup(this->getGroupID())) return; const UInt nbprocs = comm.getNBprocsOnGroup(this->getGroupID()); // Mostly taken from ParaDiS source files Main.c and ParadisInit.c std::string dlb_flag = "-r"; std::string datafile_flag = "-d"; std::string ctrlfile_flag = "-c"; std::string fake_argv0 = "IGNORED"; std::string dlb_num; std::stringstream(dlb_num) << numDLBCycles; const char *argv[8] = {fake_argv0.c_str(), datafile_flag.c_str(), data_filename.c_str(), dlb_flag.c_str(), dlb_num.c_str(), ctrlfile_flag.c_str(), control_filename.c_str(), mobility_filename.c_str()}; unsigned int argc = 8; /* * On some systems, the getrusage() call made by Meminfo() to get * the memory resident set size does not work properly. In those * cases, the function will try to return the current heap size * instead. This initial call allows meminfo() to get a copy of * the original heap pointer so subsequent calls can calculate the * heap size by taking the difference of the original and current * heap pointers. */ Meminfo(&(this->memSize)); /* * On linux systems (e.g. MCR) if built to have floating point exceptions * turned on, invoke macro to do so */ #ifdef FPES_ON unmask_std_fpes(); #endif this->paradis_ptr = InitHome(); TimerInit(this->paradis_ptr); TimerStart(this->paradis_ptr, TOTAL_TIME); TimerStart(this->paradis_ptr, INITIALIZE); if (comm.groupRank(lm_my_proc_id, this->getGroupID()) == 0) std::cout << "BEGIN: Initialize ParaDiS on " << nbprocs << " processors." << std::endl; Initialize(this->paradis_ptr, argc, const_cast(argv), comm.getMpiGroup(this->getGroupID())); if (strcmp(this->paradis_ptr->param->timestepIntegrator, "trapezoidtimestep") != 0) { LM_FATAL("The timestep integrator for ParaDiS needs to be " << "\"trapezoidtimestep\" in order for ParaDis to work with " << "LibMultiscale. You specified \"timestepIntegrator = " << this->paradis_ptr->param->timestepIntegrator << "\" in the control file \"" << control_filename << "\"."); } TimerStop(this->paradis_ptr, INITIALIZE); if (nbprocs > 1) MPI_Barrier(this->paradis_ptr->MPI_COMM_PARADIS); this->paradis_ptr->param->timeMax = this->timestep; this->paradis_ptr->cycle = this->paradis_ptr->param->cycleStart; /* * Perform the needed number (if any) of load-balance-only * steps before doing the main processing loop. These steps * perform only the minimal amount of stuff needed to * estimate per-process load, move boundaries and migrate * nodes among processsors to get a good initial balance. */ time_t tp; int initialDLBCycles = this->paradis_ptr->param->numDLBCycles; TimerStart(this->paradis_ptr, INITIALIZE); if ((this->paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) { time(&tp); printf(" +++ Beginning %d load-balancing steps at %s", initialDLBCycles, asctime(localtime(&tp))); } while (this->paradis_ptr->param->numDLBCycles > 0) { ParadisStep(this->paradis_ptr); this->paradis_ptr->cycle++; this->paradis_ptr->param->numDLBCycles--; } if ((this->paradis_ptr->myDomain == 0) && (initialDLBCycles != 0)) { time(&tp); printf(" +++ Completed load-balancing steps at %s", asctime(localtime(&tp))); } TimerStop(this->paradis_ptr, INITIALIZE); /* * Any time spent doing the initial DLB-only steps should * just be attributed to initialization time, so be sure to * reset the other timers before going into the main * computational loop */ TimerInitDLBReset(this->paradis_ptr); /* * The cycle number has been incremented during the initial * load-balance steps, and previous steps so reset it to the * proper startingvalue before entering the main processing loop. * In paradis, cycle uses 1-based indexing. */ this->paradis_ptr->cycle = this->paradis_ptr->param->cycleStart; this->paradis_ptr->cycleLMS = 1; this->getContainerNodes().setGrain(this->paradis_ptr); // this->getContainerDisplacements().setGrain(this->paradis_ptr); #ifndef LM_OPTIMIZED // Debug tests. unsigned int counter = 0; DUMP("Printing all initial ParaDiS nodes...", DBG_INFO); for (auto &&node : this->getContainerNodes()) { DUMP("Processor " << lm_my_proc_id << ", Node " << ++counter << "(" << node.position() << ")", DBG_INFO); } #endif // ContainerPARADISDisplacement displs = this->getContainerDisplacements(); // DUMP("Number of ParaDiS displacement fieldpoints: "<paradis_ptr->myDomain == 0) std::cout << "END: Initialize ParaDiS on " << nbprocs << " processors." << std::endl; - this->dd_container.setRelease(INITIAL_MODEL_RELEASE); + + LM_TOIMPLEMENT; + // this->dd_container.setRelease(INITIAL_MODEL_RELEASE); } /* -------------------------------------------------------------------------- */ void DomainPARADIS::performStep1() { this->paradis_ptr->param->timeLMSCycle = 0.; while (this->paradis_ptr->param->timeLMSCycle < this->paradis_ptr->param->timeMax) { STARTTIMER("ParadisStep"); ParadisStep(this->paradis_ptr); TimerClearAll(this->paradis_ptr); STOPTIMER("ParadisStep"); } this->paradis_ptr->cycleLMS++; this->dd_container.incRelease(); if (this->boundary) { - Cube &c = + 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(*dd); // to mimic strong grain boundary... // Real xmin = this->paradis_ptr->param->xBoundMin + 10.0; // Real xmax = this->paradis_ptr->param->xBoundMax - 10.0; // Real ymin = this->paradis_ptr->param->yBoundMin + 10.0; // Real ymax = this->paradis_ptr->param->yBoundMax - 10.0; // Real zmin = this->paradis_ptr->param->zBoundMin + 10.0; // Real zmax = this->paradis_ptr->param->zBoundMax - 10.0; // Real xmin = this->paradis_ptr->param->xBoundMin; // Real xmax = this->paradis_ptr->param->xBoundMax; // Real ymin = this->paradis_ptr->param->yBoundMin; // Real ymax = this->paradis_ptr->param->yBoundMax; // Real zmin = this->paradis_ptr->param->zBoundMin; // Real zmax = this->paradis_ptr->param->zBoundMax; Real xmin = minval[0]; Real xmax = maxval[0]; Real ymin = minval[1]; Real ymax = maxval[1]; Real zmin = minval[2]; Real zmax = maxval[2]; for (auto &&node : this->getContainerNodes()) { if ((xmin > node.position()[0]) || (node.position()[0] > xmax) || (ymin > node.position()[1]) || (node.position()[1] > ymax) || (zmin > node.position()[2]) || (node.position()[2] > zmax)) { node.slaving(); node.constrain(); } } } } /* -------------------------------------------------------------------------- */ void DomainPARADIS::performStep2() { // this->paradis_ptr->param->timeLMSCycle = 0.; // while (this->paradis_ptr->param->timeLMSCycle < // this->paradis_ptr->param->timeMax){ // ParadisStep(this->paradis_ptr); // TimerClearAll(this->paradis_ptr); // } // this->paradis_ptr->cycleLMS++; // this->dd_container.incRelease(); } void DomainPARADIS::performStep3() {} void DomainPARADIS::performStep4() {} void DomainPARADIS::performStep5() {} /* -------------------------------------------------------------------------- */ void DomainPARADIS::addSegments(std::vector> &positions, std::vector> &conn, std::vector &burgers, std::vector &normals) { UInt nb_points = positions.size(); std::vector created_nodes; std::vector 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); 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> &positions, // DD nodes in an atomistic domain std::vector> &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 numNbrs; // ParaDis Node_t *nh, *nh_conInt1 = NULL, *nh_conInt2 = NULL; // ParaDis - Iterators std::vector::iterator nh_it; /* ------------------------------------------------------------------------ */ // Delete all existing DD segments in the atomistic domain // NOTE: the FreeNode() function doesn't remove the segments associated with // the node to delete // --> more precisely, they are not deallocated // --> this is of no particular issue since one can call // ReallocNodeArms() in step 3) for (nh_it = created_nodes.begin(); nh_it != created_nodes.end(); ++nh_it) { FreeNode(paradis_ptr, (*nh_it)->myTag.index); } created_nodes.clear(); /* ------------------------------------------------------------------------ */ /* 1) Import new nodes in ParaDis */ for (UInt i = 0; i < nb_points; ++i) { nh = GetNewNativeNode(this->paradis_ptr); auto &&pos = positions[i].position(); nh->x = pos[0]; nh->y = pos[1]; nh->z = pos[2]; nh->native = 1; nh->constraint = PINNED_NODE; // nodes in an atomistic domain remain always fixed nh->oldvX = 0.0; nh->oldvY = 0.0; nh->oldvZ = 0.0; AssignNodeToCell(this->paradis_ptr, nh); created_nodes.push_back(nh); } numNbrs.resize(created_nodes.size()); /* ------------------------------------------------------------------------ */ /* 2) Find nearest neighbor(s) in the continuum */ // 2.1) First node nh = created_nodes[0]; pos_ato[0] = nh->x; pos_ato[1] = nh->y; pos_ato[2] = nh->z; for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end(); ++nh_it) { // Position of "continuum" node pos_con[0] = (*nh_it)->x; pos_con[1] = (*nh_it)->y; pos_con[2] = (*nh_it)->z; // Compute the distance between the continuum and the dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) + (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) + (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2])); // If distance is smaller than a prescribed threshold // --> mark the node as a neighbor of the first "atomistic" DD node if (dist < threshold) { nh_conInt1 = *nh_it; numNbrs[0] = 1; break; } } // 2.2) Last node nh = created_nodes[nb_points - 1]; pos_ato[0] = nh->x; pos_ato[1] = nh->y; pos_ato[2] = nh->z; for (nh_it = this->nh_l_conInt.begin(); nh_it != this->nh_l_conInt.end(); ++nh_it) { // Position of "continuum" node pos_con[0] = (*nh_it)->x; pos_con[1] = (*nh_it)->y; pos_con[2] = (*nh_it)->z; // Compute the distance between the continuum and the dist = sqrt((pos_ato[0] - pos_con[0]) * (pos_ato[0] - pos_con[0]) + (pos_ato[1] - pos_con[1]) * (pos_ato[1] - pos_con[1]) + (pos_ato[2] - pos_con[2]) * (pos_ato[2] - pos_con[2])); // If distance is smaller than a prescribed threshold // --> mark the node as a neighbor of the last "atomistic" DD node if (dist < threshold) { nh_conInt2 = *nh_it; numNbrs[nb_points - 1] = 1; break; } } /* ------------------------------------------------------------------------ */ /* 3) Import connectivity and allocate memory for segment info */ UInt nb_elem = conn.size(); for (UInt i = 0; i < nb_elem; ++i) { RefGenericElem<3> &el = conn[i]; LM_ASSERT(el.conn.size() == 2, "we cannot add other things than segments to DD model"); for (UInt n = 0; n < 2; ++n) { UInt nd = el.conn[n]; ++numNbrs[nd]; } } for (UInt i = 0; i < nb_points; ++i) { ReallocNodeArms(created_nodes[i], numNbrs[i]); // AllocNodeArms(newNode, numNbrs[i]); numNbrs[i] = 0; } /* ------------------------------------------------------------------------ */ /* 4) Connect atomistic and continuum DD line */ // 4.1) First interface node in the continuum connects to the atomistic DD // line if ((nh_conInt1 != NULL) && (nh_conInt2 == NULL)) { // Get the Burgers vector of the dislocation burgersVec[0] = nh_conInt1->burgX[0]; burgersVec[1] = nh_conInt1->burgY[0]; burgersVec[2] = nh_conInt1->burgZ[0]; // Get the normal vector of the glide plane gpNormVec[0] = nh_conInt1->nx[0]; gpNormVec[1] = nh_conInt1->ny[0]; gpNormVec[2] = nh_conInt1->nz[0]; // Re-allocate memory for the hybrid segment if necessary (--> should only // be done once before the first iteration) if (nh_conInt1->numNbrs < 2) { ReallocNodeArms(nh_conInt1, 2); nh_conInt1->constraint = UNCONSTRAINED; } // Establish connectivity between nodes in the continuum and the atomistic // domain ... setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec); // // ... and set homogeneous Dirichlet conditions on the other node // (--> u=0 <=> pinned node) // created_nodes[nb_points-1]->constraint = PINNED_NODE; // Project atomistic DD node on the glide plane projectNodeOnGP(nh_conInt1, created_nodes[0], 1); ++numNbrs[0]; } // 4.2) Last interface node in the continuum connects to the atomistic DD line else if ((nh_conInt1 == NULL) && (nh_conInt2 != NULL)) { // Get the Burgers vector of the dislocation burgersVec[0] = nh_conInt2->burgX[0]; burgersVec[1] = nh_conInt2->burgY[0]; burgersVec[2] = nh_conInt2->burgZ[0]; // Get the normal vector of the glide plane gpNormVec[0] = nh_conInt2->nx[0]; gpNormVec[1] = nh_conInt2->ny[0]; gpNormVec[2] = nh_conInt2->nz[0]; // Re-allocate memory for the hybrid segment if necessary (--> should only // be done once before the first iteration) if (nh_conInt2->numNbrs < 2) { ReallocNodeArms(nh_conInt2, 2); nh_conInt2->constraint = UNCONSTRAINED; } // Establish connectivity between nodes in the continuum and the atomistic // domain ... setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec, gpNormVec); // // ... and set homogeneous Dirichlet conditions on the other node // (--> u=0 <=> pinned node) // created_nodes[0]->constraint = PINNED_NODE; // Project atomistic DD node on the glide plane projectNodeOnGP(nh_conInt2, created_nodes[nb_points - 1], 1); ++numNbrs[nb_points - 1]; } // 4.3) Both interfaces nodes in the continuum connects to the atomistic DD // line else if ((nh_conInt1 != NULL) && (nh_conInt2 != NULL)) { // Get the Burgers vector of the dislocation burgersVec[0] = nh_conInt1->burgX[0]; burgersVec[1] = nh_conInt1->burgY[0]; burgersVec[2] = nh_conInt1->burgZ[0]; // Get the normal vector of the glide plane gpNormVec[0] = nh_conInt1->nx[0]; gpNormVec[1] = nh_conInt1->ny[0]; gpNormVec[2] = nh_conInt1->nz[0]; // Re-allocate memory for the hybrid segment if necessary (--> should only // be done once before the first iteration) if (nh_conInt1->numNbrs < 2) { ReallocNodeArms(nh_conInt1, 2); nh_conInt1->constraint = UNCONSTRAINED; } if (nh_conInt2->numNbrs < 2) { ReallocNodeArms(nh_conInt2, 2); nh_conInt2->constraint = UNCONSTRAINED; } // Establish connectivity between nodes in the continuum and the atomistic // domain setNodeArm(nh_conInt1, created_nodes[0], 1, 0, burgersVec, gpNormVec); setNodeArm(nh_conInt2, created_nodes[nb_points - 1], 1, 0, burgersVec, gpNormVec); // Project atomistic DD nodes on the glide plane projectNodeOnGP(nh_conInt1, created_nodes[0], 1); projectNodeOnGP(nh_conInt1, created_nodes[nb_points - 1], 1); ++numNbrs[0]; ++numNbrs[nb_points - 1]; } else { LM_FATAL(std::endl << "Hybrid dislocation failed to connect!!" << std::endl); } /* ------------------------------------------------------------------------ */ /* 5) Now connect the remainder of the hybrid dislocation line in the * atomistic domain */ for (UInt i = 0; i < nb_elem; ++i) { // Get connectivity RefGenericElem<3> &el = conn[i]; UInt idxNhSrc = el.conn[0]; UInt idxNhDst = el.conn[1]; UInt &idxArmSrc = numNbrs[idxNhSrc]; UInt &idxArmDst = numNbrs[idxNhDst]; // Update arm information for each node of segment i setNodeArm(created_nodes[idxNhSrc], created_nodes[idxNhDst], idxArmSrc, idxArmDst, burgersVec, gpNormVec); // Project atomistic DD node on the glide plane projectNodeOnGP(created_nodes[idxNhSrc], created_nodes[idxNhDst], idxArmSrc); ++idxArmSrc; ++idxArmDst; } } /* -------------------------------------------------------------------------- */ UInt DomainPARADIS::getNBNodes() { return this->getContainerNodes().size(); } /* -------------------------------------------------------------------------- */ Real DomainPARADIS::getShearModulus() { return this->paradis_ptr->param->shearModulus; } /* -------------------------------------------------------------------------- */ VectorProxy<3u> DomainPARADIS::getXmin() { return VectorProxy<3>(this->paradis_ptr->param->xBoundMin, this->paradis_ptr->param->yBoundMin, this->paradis_ptr->param->zBoundMin); } /* -------------------------------------------------------------------------- */ VectorProxy<3u> DomainPARADIS::getXmax() { return VectorProxy<3>(this->paradis_ptr->param->xBoundMax, this->paradis_ptr->param->yBoundMax, this->paradis_ptr->param->zBoundMax); } /* -------------------------------------------------------------------------- */ void DomainPARADIS::freenodearm(UInt index) { Node_t node = *paradis_ptr->nodeKeys[index]; FreeNodeArms(&node); } /* -------------------------------------------------------------------------- */ void DomainPARADIS::freenode(UInt index) { Node_t node = *paradis_ptr->nodeKeys[index]; FreeNodeArms(&node); } /* -------------------------------------------------------------------------- */ bool DomainPARADIS::is_it_bound() { // ContainerPARADISNode::iterator dbgNodes = // this->getContainerNodes().begin(); // UInt counter = 0; // for(RefPARADISNode node=dbgNodes.getFirst(); !dbgNodes.end(); // node=dbgNodes.getNext()){ // ++counter; // } return this->boundary; } // CHANGE THIS? /* -------------------------------------------------------------------------- */ /* LMDESC PARADIS This section is used to configure the ParaDiS plugin */ /* LMEXAMPLE Section MultiScale AtomsUnits\\ ...\\ GEOMETRY 1 CUBE BBOX -10 10 -10 10 -10 10 MODEL PARADIS dislomodel \\ ...\\ endSection\\ \\ Section PARADIS:dislomodel RealUnits \\ NUMDLBSTEPS 0 \\ CONTROLFILE frank_read_src.ctrl \\ DATAFILE frank_read_src.data \\ TIMESTEP 1e-10 \\ DOMAIN_GEOMETRY 1 \\ endSection */ /* LMHERITATE domain_dd */ void DomainPARADIS::declareParams() { DomainDD::declareParams(); /* LMKEYWORD TIMESTEP Set the maximum time step used for ParaDiS */ this->parseKeyword("TIMESTEP", timestep); /* LMKEYWORD NUMDLBSTEPS Number of dynamic load balancing steps */ this->parseKeyword("NUMDLBSTEPS", numDLBCycles); /* LMKEYWORD CONTROLFILE Path to the ParaDiS control file */ this->parseKeyword("CONTROLFILE", control_filename); /* LMKEYWORD DATAFILE Path to the ParaDis data file */ this->parseKeyword("DATAFILE", data_filename); /* LMKEYWORD MOBILITYFILE Mobility input acquired from MD */ this->parseKeyword("MOBILITYFILE", mobility_filename); /* LMKEYWORD BOUNDARY Constraint boundary to mimic grain boundary */ this->parseTag("BOUNDARY", this->boundary, false); } /* -------------------------------------------------------------------------- */ __END_LIBMULTISCALE__ diff --git a/src/dd/paradis/domain_paradis.hh b/src/dd/paradis/domain_paradis.hh index 186ab21..38f93cc 100644 --- a/src/dd/paradis/domain_paradis.hh +++ b/src/dd/paradis/domain_paradis.hh @@ -1,220 +1,221 @@ /** * @file domain_paradis.hh * * @author Till Junge * @author Jaehyun Cho * @author Max Hodapp * * @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 . * */ #ifndef __LIBMULTISCALE_DOMAIN_PARADIS_HH__ #define __LIBMULTISCALE_DOMAIN_PARADIS_HH__ /* -------------------------------------------------------------------------- */ #include "container_array.hh" #include "container_paradis.hh" #include "domain_dd.hh" /* -------------------------------------------------------------------------- */ __BEGIN_LIBMULTISCALE__ /** * Class DummySource * DD dummy container (defect fictive container ?) */ typedef ContainerPARADISNode DummySource; typedef ContainerPARADISNode DummyObstacle; /* -------------------------------------------------------------------------- */ /** * Class DomainPARADIS * domain to implement the plugin interface to PARADIS dislocation dynamics */ class DomainPARADIS : public DomainDD { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - DomainPARADIS(DomainID ID, CommGroup GID) - : DomainDD(ID, GID) { + DomainPARADIS(LMID ID, CommGroup GID) + : LMObject(ID), + DomainDD(ID, GID) { paradis_ptr = nullptr; } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void performStep1(); void performStep2(); void performStep3(); void performStep4(); void performStep5(); void imageForce(const double *position, const double *normal, double *force); void getSubDomainDimensions(Real *xmin, Real *xmax) { LM_TOIMPLEMENT; }; void getDomainDimensions(Real *xmin, Real *xmax) { for (UInt dir = 0; dir < this->getDim(); ++dir) { xmin[dir] = this->paradis_ptr->param->minCoordinates[dir]; xmax[dir] = this->paradis_ptr->param->maxCoordinates[dir]; } }; //! compute displacement field for a set of positions void computeDisplacements(std::vector &positions, std::vector &displacements); //! 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(); //! intialisation stage void init(); //! parse the keywords void declareParams(); //! add segments to the dd model virtual void addSegments(std::vector> &pos, std::vector> &conn, std::vector &burgers, std::vector &glide_normals); 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> &positions, // DD nodes in an atomistic domain std::vector> &conn, // connectivity between DD nodes in the atomistic domain Real threshold); // max. distance between two adjacent DD nodes possibly // generating a hybrid segment /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ //! set current step to a given value void setCurrentStep(UInt ts) { LM_FATAL("to be implemented"); } // GUILLAUME? //! get current step UInt getCurrentStep() { LM_FATAL("to be implemented"); return 0; }; // CHANGE THIS //! set timestep to a given value void setTimeStep(Real ts){LM_FATAL("to be implemented")}; // CHANGE THIS //! GetTimeStep Real getTimeStep() { return this->timestep; }; // CHANGE THIS //! return the stress component (what is i and j ?) double &get_stress_component(int i, int j) { int indices[3][3] = {{0, 5, 4}, {5, 1, 3}, {4, 3, 2}}; return this->paradis_ptr->param->appliedStress[indices[i][j]]; } UInt getDim() { return 3; }; Real getPoisson() { return this->paradis_ptr->param->pois; }; Real getMaxSeg() { return this->paradis_ptr->param->maxSeg; }; Real getMinSeg() { return this->paradis_ptr->param->minSeg; }; Real getBurgersMagnitude() { return this->paradis_ptr->param->burgMag; }; UInt getNBpairsInfinite() { return this->paradis_ptr->param->numPairs; }; bool getPairedIndex(int index0, int *index0_infinite) { UInt num_pairs = this->paradis_ptr->param->numPairs; for (int i = 0; i < (int)num_pairs; ++i) { Pair_t pair_ind; pair_ind = *this->paradis_ptr->pairKeys[i]; if (pair_ind.index1 == index0) { *index0_infinite = pair_ind.index2; return true; } else if (pair_ind.index2 == index0) { *index0_infinite = pair_ind.index1; return true; } } return false; }; /* ------------------------------------------------------------------------ */ /* Functions for AMELCG */ /* ------------------------------------------------------------------------ */ double getExternalWork() { return 0; } void setExternalWork(double) {} /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: bool boundary; int memSize; Home_t *paradis_ptr; Quantity