diff --git a/.gitignore b/.gitignore
index 1d1ae5ca..b4b35554 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,5 +7,6 @@ log.*
*~
**/linux*Gcc*/
+**/.vscode
lnInclude
diff --git a/doc/dataExchangeModel_twoWayOne2One.txt b/doc/dataExchangeModel_twoWayOne2One.txt
new file mode 100644
index 00000000..bbd609af
--- /dev/null
+++ b/doc/dataExchangeModel_twoWayOne2One.txt
@@ -0,0 +1,47 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+dataExchangeModel_twoWayOne2One command :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+dataExchangeModel twoWayOne2One;
+twoWayOne2OneProps
+\{
+ liggghtsPath "path";
+ useStaticProcMap useMap;
+\}; :pre
+
+{path} = path to the DEM simulation input file :ulb,l
+{useMap} = yes or no, determines if the map is build once (yes) or every coupling step (no) (default: no):l
+:ule
+
+[Examples:]
+
+dataExchangeModel twoWayOne2One;
+twoWayOne2OneProps
+\{
+ liggghtsPath "../DEM/in.liggghts_init";
+ useStaticProcMap yes;
+\} :pre
+
+[Description:]
+
+The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayOne2One model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using a more sophisticated mapping scheme than twoWayMPI / all2all and scales much better for large systems and many cores. The DEM run is executed by the coupling model, via a liggghtsCommandModel object. Only use staticProcMap yes if no load balancing is employed.
+
+[Restrictions:]
+
+Must be used in combination with the engineSearchMany2Many locate model! Use the "one2one" cfd datacoupling option in fix couple/cfd in LIGGGHTS!
+
+Some warnings may be given for particles that have not been located - this is due to LIGGGHTS' treatment of domain crossers.
+
+[Related commands:]
+
+"dataExchangeModel"_dataExchangeModel.html
+
diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files
index e4eb3b00..2d409055 100644
--- a/src/lagrangian/cfdemParticle/Make/files
+++ b/src/lagrangian/cfdemParticle/Make/files
@@ -149,6 +149,8 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
$(dataExchangeModels)/twoWayFiles/twoWayFiles.C
$(dataExchangeModels)/noDataExchange/noDataExchange.C
$(dataExchangeModels)/twoWayMPI/twoWayMPI.C
+$(dataExchangeModels)/twoWayOne2One/twoWayOne2One.C
+$(dataExchangeModels)/twoWayOne2One/one2one.C
$(averagingModels)/averagingModel/averagingModel.C
$(averagingModels)/averagingModel/newAveragingModel.C
diff --git a/src/lagrangian/cfdemParticle/Make/options b/src/lagrangian/cfdemParticle/Make/options
index 9c152540..8a84bb35 100644
--- a/src/lagrangian/cfdemParticle/Make/options
+++ b/src/lagrangian/cfdemParticle/Make/options
@@ -34,6 +34,4 @@ LIB_LIBS = \
-lmpi_cxx \
-Wl,-rpath,$(CFDEM_LIGGGHTS_BIN_DIR) \
-L$(CFDEM_LIGGGHTS_BIN_DIR) \
- -lliggghts \
- -L$(CFDEM_Many2ManyLIB_PATH) \
- -lcoupleMany2Many
+ -lliggghts
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H
index 6153b593..672f1138 100755
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H
@@ -245,17 +245,25 @@ public:
virtual int getNumberOfTypes() const;
virtual double* getTypeVol() const;
- inline void setPositions(label n,double* pos) const
+ inline void setPositions(label n,double* pos)
{
for (int i=0;i.
+
+Copyright 2018- Paul Kieckhefen, TUHH
+
+\*---------------------------------------------------------------------------*/
+
+//#define O2ODEBUG
+
+#ifdef O2ODEBUG
+#include
+#endif
+
+#include
+#include "one2one.H"
+
+template
+struct mpi_type_wrapper {
+ MPI_Datatype mpi_type;
+ mpi_type_wrapper();
+};
+template <> mpi_type_wrapper::mpi_type_wrapper()
+: mpi_type(MPI_FLOAT) {}
+template <> mpi_type_wrapper::mpi_type_wrapper()
+: mpi_type(MPI_DOUBLE) {}
+template <> mpi_type_wrapper::mpi_type_wrapper()
+: mpi_type(MPI_INT) {}
+
+/* ---------------------------------------------------------------------- */
+
+One2One::One2One(MPI_Comm caller)
+:
+ncollected_(-1),
+comm_(caller),
+nsrc_procs_(-1),
+src_procs_(nullptr),
+ndst_procs_(-1),
+dst_procs_(nullptr),
+nlocal_(-1),
+natoms_(nullptr),
+request_(nullptr),
+status_(nullptr)
+{
+ MPI_Comm_rank(comm_,&me_);
+ MPI_Comm_size(comm_,&nprocs_);
+}
+
+/* ---------------------------------------------------------------------- */
+
+One2One::~One2One()
+{
+ deallocate();
+}
+
+/* ----------------------------------------------------------------------
+ communicate particle ids based on processor communication pattern
+------------------------------------------------------------------------- */
+
+void One2One::setup
+(
+ int nsrc_procs,
+ int *src_procs,
+ int ndst_procs,
+ int *dst_procs,
+ int nlocal
+)
+{
+ // free any previous one2one info
+ deallocate();
+
+ src_procs_ = src_procs;
+ nsrc_procs_ = nsrc_procs;
+ dst_procs_ = dst_procs;
+ ndst_procs_ = ndst_procs;
+ nlocal_ = nlocal;
+
+ // gather number of ids for reserving memory
+ natoms_ = new int[nprocs_];
+ MPI_Allgather // may be replaced by send/irecv
+ (
+ &nlocal_,
+ 1,
+ MPI_INT,
+ natoms_,
+ 1,
+ MPI_INT,
+ comm_
+ );
+
+ ncollected_ = 0;
+ int nrequests = 0;
+ for (int i = 0; i < nsrc_procs_; i++)
+ {
+ if (natoms_[src_procs_[i]] > 0)
+ {
+ ncollected_ += natoms_[src_procs_[i]];
+
+ if (src_procs_[i] != me_) // no receive for on-proc info
+ {
+ nrequests++;
+ }
+ }
+ }
+
+ if (nrequests > 0)
+ {
+ request_ = new MPI_Request[nrequests];
+ status_ = new MPI_Status[nrequests];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ src: what is present on this proc
+ dst: what is received from other procs
+ all comm according to map set up in setup(...)
+------------------------------------------------------------------------- */
+template
+void One2One::exchange(T *&src, T *&dst, int data_length)
+{
+ mpi_type_wrapper wrap;
+
+ // post receives
+ int offset_local = -1;
+ int offset = 0;
+ int requesti = 0;
+ for (int i = 0; i < nsrc_procs_; i++)
+ {
+ // do post a receives for procs who own particles
+ if (natoms_[src_procs_[i]] > 0)
+ {
+ if (src_procs_[i] != me_)
+ {
+ #ifdef O2ODEBUG
+ std::cout<< "[" << me_ << "]"
+ << " RCV " << i
+ << " of " << nsrc_procs_
+ << " from: " << src_procs_[i]
+ << " natoms_[src_procs_[i]] " << natoms_[src_procs_[i]]
+ << " datalength " << data_length
+ << " offset " << offset
+ << std::endl;
+ #endif
+ MPI_Irecv
+ (
+ &dst[offset],
+ natoms_[src_procs_[i]]*data_length,
+ wrap.mpi_type,
+ src_procs_[i],
+ MPI_ANY_TAG,
+ comm_,
+ &request_[requesti]
+ );
+ requesti++;
+ }
+ else // data is available on-proc
+ {
+ offset_local = offset;
+ }
+ }
+ offset += natoms_[src_procs_[i]]*data_length;
+ }
+
+ // make sure all receives are posted
+ MPI_Barrier(comm_);
+
+ // blocking sends - do nonblocking instead
+ // since doing many-2-many here?
+ // only do sends if I have particles
+ if (nlocal_ > 0)
+ {
+ for (int i = 0; i < ndst_procs_; i++)
+ {
+ if (dst_procs_[i] != me_)
+ {
+ #ifdef O2ODEBUG
+ std::cout<< "[" << me_ << "]"
+ << " SEND to: " << dst_procs_[i]
+ << " nlocal_ " << nlocal_
+ << " data_length " << data_length
+ << std::endl;
+ #endif
+ MPI_Send
+ (
+ src,
+ nlocal_*data_length,
+ wrap.mpi_type,
+ dst_procs_[i],
+ 0,
+ comm_
+ );
+ }
+ }
+ }
+
+ // only wait if requests were actually posted
+ if (requesti > 0)
+ MPI_Waitall(requesti, request_, status_);
+
+ // copy on-proc data
+ if (offset_local > -1)
+ {
+ const int max_locali = nlocal_ * data_length;
+ for
+ (
+ int locali = 0;
+ locali < max_locali;
+ locali++
+ )
+ {
+ dst[locali+offset_local] = src[locali];
+ }
+ }
+}
+
+template void One2One::exchange(int*&, int*&, int);
+template void One2One::exchange(double*&, double*&, int);
+
+// there should be a way to do this without copying data
+template
+void One2One::exchange(T **&src, T **&dst, int data_length)
+{
+ mpi_type_wrapper wrap;
+
+ T* tmp_dst = new T[ncollected_*data_length];
+ T* tmp_src = new T[nlocal_*data_length];
+
+ for (int i = 0; i < nlocal_; i++)
+ for (int j = 0; j < data_length; j++)
+ tmp_src[data_length*i+j] = src[i][j];
+
+ exchange(tmp_src, tmp_dst, data_length);
+
+ for (int i = 0; i < ncollected_; i++)
+ for (int j = 0; j < data_length; j++)
+ dst[i][j] = tmp_dst[data_length*i+j];
+
+ delete [] tmp_src;
+ delete [] tmp_dst;
+}
+template void One2One::exchange(int**&, int**&, int);
+template void One2One::exchange(double**&, double**&, int);
+
+
+template
+void One2One::exchange(T **&src, T *&dst, int data_length)
+{
+ mpi_type_wrapper wrap;
+
+ T* tmp_src = new T[nlocal_*data_length];
+
+ for (int i = 0; i < nlocal_; i++)
+ for (int j = 0; j < data_length; j++)
+ tmp_src[data_length*i+j] = src[i][j];
+
+ exchange(tmp_src, dst, data_length);
+
+
+ delete [] tmp_src;
+}
+template void One2One::exchange(int**&, int*&, int);
+template void One2One::exchange(double**&, double*&, int);
+
+/* ---------------------------------------------------------------------- */
+
+void One2One::deallocate()
+{
+ delete [] src_procs_;
+ delete [] dst_procs_;
+ delete [] natoms_;
+
+ delete [] request_;
+ delete [] status_;
+}
\ No newline at end of file
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H
new file mode 100644
index 00000000..4205e11c
--- /dev/null
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/one2one.H
@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+License
+
+ This is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This code 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 General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this code. If not, see .
+
+Copyright 2018- Paul Kieckhefen, TUHH
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ONE2ONE_H
+#define ONE2ONE_H
+
+#include
+
+class One2One
+{
+ public:
+ One2One(MPI_Comm);
+ ~One2One();
+
+ void setup
+ (
+ int nsrc_procs,
+ int *src_procs,
+ int ndst_procs,
+ int* dst_procs,
+ int nlocal
+ );
+
+ template
+ void exchange(T *&, T *&, int data_length=1);
+ template
+ void exchange(T **&, T **&, int data_length=1);
+ template
+ void exchange(T **&, T *&, int data_length=1);
+
+ int ncollected_; // # of ids in from group
+
+ protected:
+
+ MPI_Comm comm_;
+ int me_, nprocs_; // rank and size
+
+ // communication partners
+ int nsrc_procs_; // # of off-processor IDs
+ int* src_procs_; // procs I receive data from
+ int ndst_procs_; // # of off-processor IDs
+ int* dst_procs_; // procs I receive data from
+
+ int nlocal_; // # particle ids I own
+ int* natoms_;
+
+ MPI_Request* request_;
+ MPI_Status* status_;
+
+ void deallocate();
+};
+
+#endif
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C
new file mode 100644
index 00000000..f921dde4
--- /dev/null
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayOne2One/twoWayOne2One.C
@@ -0,0 +1,941 @@
+/*---------------------------------------------------------------------------*\
+ CFDEMcoupling - Open Source CFD-DEM coupling
+
+ CFDEMcoupling is part of the CFDEMproject
+ www.cfdem.com
+ Christoph Goniva, christoph.goniva@cfdem.com
+ Copyright 2009-2012 JKU Linz
+ Copyright 2012- DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+ This file is part of CFDEMcoupling.
+
+ CFDEMcoupling is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 3 of the License, or (at your
+ option) any later version.
+
+ CFDEMcoupling 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 General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with CFDEMcoupling; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+ This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+ and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+Contributing authors
+ Paul Kieckhefen (TUHH) 2018-
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "twoWayOne2One.H"
+#include "addToRunTimeSelectionTable.H"
+#include "clockModel.H"
+#include "pair.h"
+#include "force.h"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(twoWayOne2One, 0);
+
+addToRunTimeSelectionTable
+(
+ dataExchangeModel,
+ twoWayOne2One,
+ dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+// Construct from components
+twoWayOne2One::twoWayOne2One
+(
+ const dictionary& dict,
+ cfdemCloud& sm
+)
+:
+ dataExchangeModel(dict,sm),
+ propsDict_(dict.subDict(typeName + "Props")),
+ thisLigPartner_(0),
+ thisFoamPartner_(0),
+ lig2foam_(nullptr),
+ foam2lig_(nullptr),
+ lig2foam_mask_(nullptr),
+ lig2foam_ids_(nullptr),
+ foam2lig_ids_(nullptr),
+ lig2foam_vec_tmp_(nullptr),
+ lig2foam_scl_tmp_(nullptr),
+ foam2lig_vec_tmp_(nullptr),
+ foam2lig_scl_tmp_(nullptr),
+ staticProcMap_(propsDict_.lookupOrDefault("useStaticProcMap", false)),
+ cellIdComm_(propsDict_.lookupOrDefault("useCellIdComm", false)),
+ my_prev_cell_ids_fix_(nullptr),
+ verbose_(propsDict_.lookupOrDefault("verbose", false)),
+ lmp(nullptr)
+{
+ Info<<"Starting up LIGGGHTS for first time execution"<input->file(liggghtsPath.c_str());
+
+ // get DEM time step size
+ DEMts_ = lmp->update->dt;
+ checkTSsize();
+
+ // calculate boundingBox of FOAM subdomain
+ primitivePatch tmpBoundaryFaces
+ (
+ SubList
+ (
+ sm.mesh().faces(),
+ sm.mesh().nFaces() - sm.mesh().nInternalFaces(),
+ sm.mesh().nInternalFaces()
+ ),
+ sm.mesh().points()
+ );
+ typedef PrimitivePatch bPatch;
+ bPatch boundaryFaces
+ (
+ tmpBoundaryFaces.localFaces(),
+ tmpBoundaryFaces.localPoints()
+ );
+ thisFoamBox_ = treeBoundBox(boundaryFaces.localPoints());
+ if (staticProcMap_)
+ {
+ createProcMap();
+ }
+
+ if (cellIdComm_)
+ {
+ my_prev_cell_ids_fix_ = static_cast
+ ( lmp->modify->find_fix_property
+ (
+ "prev_cell_ids",
+ "property/atom",
+ "scalar",
+ 0,
+ 0,
+ "cfd coupling",
+ true
+ )
+ );
+ }
+}
+
+void twoWayOne2One::createProcMap()
+{
+ List foamBoxes(Pstream::nProcs());
+ foamBoxes[Pstream::myProcNo()] = thisFoamBox_;
+ Pstream::gatherList(foamBoxes);
+ Pstream::scatterList(foamBoxes);
+
+ // calculate bounding box of LIG subdomain
+ // this may have to move to couple when dynamic LB occurs
+ List ligBoxes(Pstream::nProcs());
+ double** ligbb = o2o_liggghts_get_boundingbox(lmp);
+ boundBox thisLigBox
+ (
+ point(ligbb[0][0], ligbb[0][1], ligbb[0][2]),
+ point(ligbb[1][0], ligbb[1][1], ligbb[1][2])
+ );
+ ligBoxes[Pstream::myProcNo()] = thisLigBox;
+ Pstream::gatherList(ligBoxes);
+ Pstream::scatterList(ligBoxes);
+
+ thisLigPartner_.clear();
+ thisFoamPartner_.clear();
+
+ // detect LIG subdomains which this FOAM has to interact with
+ forAll(ligBoxes, ligproci)
+ {
+ if (thisFoamBox_.overlaps(ligBoxes[ligproci]))
+ {
+ thisLigPartner_.append(ligproci);
+ }
+ }
+ // detect FOAM subdomains this LIG has to interact with
+ // TODO: refactor to invert this list here
+ forAll(foamBoxes, foamproci)
+ {
+ if (thisLigBox.overlaps(foamBoxes[foamproci]))
+ {
+ thisFoamPartner_.append(foamproci);
+ }
+ }
+
+ if (verbose_)
+ {
+ Pout<< "FOAM bounding box: " << thisFoamBox_
+ << " LIG bounding box: " << thisLigBox
+ << nl
+ << "FOAM comm partners: " << thisFoamPartner_
+ << " LIG comm partners: " << thisLigPartner_
+ << endl;
+ }
+}
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+twoWayOne2One::~twoWayOne2One()
+{
+ delete foam2lig_;
+ delete lig2foam_;
+
+ destroy(lig2foam_ids_);
+ destroy(foam2lig_ids_);
+
+ destroy(lig2foam_vec_tmp_);
+ destroy(lig2foam_scl_tmp_);
+ destroy(foam2lig_vec_tmp_);
+ destroy(foam2lig_scl_tmp_);
+
+ delete lmp;
+}
+
+// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * //
+void twoWayOne2One::getData
+(
+ word name,
+ word type,
+ double ** const& field,
+ label /*step*/
+) const
+{
+ if (name == "x") // the location is transferred by couple()
+ {
+ return;
+ }
+ if (type == "vector-atom")
+ {
+ double **tmp_= static_cast(lammps_extract_atom(lmp,name.c_str()));
+ if (!tmp_)
+ {
+ LAMMPS_NS::Fix *fix = nullptr;
+ fix = lmp->modify->find_fix_property
+ (
+ name.c_str(),
+ "property/atom",
+ "vector",
+ 0,
+ 0,
+ "cfd coupling",
+ false
+ );
+ if (fix)
+ {
+ tmp_ = static_cast(fix)->array_atom;
+ }
+ else
+ {
+ Warning<< "coupling fix not found!" << endl;
+ }
+
+ if (!tmp_)
+ {
+ FatalError<< "find_fix_property " << name
+ << " array_atom not found."
+ << abort(FatalError);
+ }
+ }
+
+ lig2foam_->exchange
+ (
+ tmp_,
+ lig2foam_vec_tmp_,
+ 3
+ );
+ extractCollected
+ (
+ lig2foam_vec_tmp_,
+ const_cast(field),
+ 3
+ );
+ }
+ else if (type == "scalar-atom")
+ {
+ double *tmp_ = static_cast(lammps_extract_atom(lmp,name.c_str()));
+ if (!tmp_)
+ {
+ LAMMPS_NS::Fix *fix = nullptr;
+ fix = lmp->modify->find_fix_property
+ (
+ name.c_str(),
+ "property/atom",
+ "scalar",
+ 0,
+ 0,
+ "cfd coupling",
+ true
+ );
+
+ if (fix)
+ {
+ tmp_ = static_cast(fix)->vector_atom;
+ }
+ else
+ {
+ FatalError<< "coupling fix not found!" << abort(FatalError);
+ }
+
+ if (!tmp_)
+ {
+ FatalError<< "find_fix_property " << name
+ << " vector_atom not found."
+ << abort(FatalError);
+ }
+ }
+ lig2foam_->exchange
+ (
+ tmp_,
+ lig2foam_scl_tmp_
+ );
+ extractCollected
+ (
+ lig2foam_scl_tmp_,
+ const_cast(field)
+ );
+ }
+ else
+ {
+ FatalError << "requesting type " << type << " and name " << name << abort(FatalError);
+ }
+}
+
+void twoWayOne2One::getData
+(
+ word name,
+ word type,
+ int ** const& field,
+ label /*step*/
+) const
+{
+ FatalError << "do not use this getData!!!" << abort(FatalError);
+/*
+ o2o_data_liggghts_to_of
+ (
+ name.c_str(),
+ type.c_str(),
+ lmp,
+ (void*&) field,
+ "int",
+ comm_liggghts_
+ );
+*/
+}
+
+void twoWayOne2One::giveData
+(
+ word name,
+ word type,
+ double ** const& field,
+ const char* datatype
+) const
+{
+ if (type == "vector-atom")
+ {
+ foam2lig_->exchange
+ (
+ const_cast(field),
+ foam2lig_vec_tmp_,
+ 3
+ );
+ o2o_data_of_to_liggghts
+ (
+ name.c_str(),
+ type.c_str(),
+ lmp,
+ foam2lig_vec_tmp_,
+ datatype,
+ foam2lig_ids_,
+ foam2lig_->ncollected_
+ );
+ }
+ else if (type == "scalar-atom")
+ {
+ foam2lig_->exchange
+ (
+ const_cast(field),
+ foam2lig_scl_tmp_,
+ 1
+ );
+ o2o_data_of_to_liggghts
+ (
+ name.c_str(),
+ type.c_str(),
+ lmp,
+ foam2lig_scl_tmp_,
+ datatype,
+ foam2lig_ids_,
+ foam2lig_->ncollected_
+ );
+ }
+ else
+ {
+ FatalError<< "twoWayMany2Many::giveData requested type " << type
+ << " not implemented!"
+ << abort(FatalError);
+ }
+}
+
+//============
+// double **
+void twoWayOne2One::allocateArray
+(
+ double**& array,
+ double initVal,
+ int width,
+ int length
+) const
+{
+ int len = max(length,1);
+ lmp->memory->grow(array, len, width, "o2o:dbl**");
+ for (int i = 0; i < len; i++)
+ for (int j = 0; j < width; j++)
+ array[i][j] = initVal;
+}
+
+void twoWayOne2One::allocateArray
+(
+ double**& array,
+ double initVal,
+ int width,
+ const char* length
+) const
+{
+ int len = max(particleCloud_.numberOfParticles(),1);
+ lmp->memory->grow(array, len, width, "o2o:dbl**:autolen");
+ for (int i = 0; i < len; i++)
+ for (int j = 0; j < width; j++)
+ array[i][j] = initVal;
+}
+
+void inline twoWayOne2One::destroy(double** array,int len) const
+{
+ lmp->memory->destroy(array);
+}
+
+//============
+// int **
+void twoWayOne2One::allocateArray
+(
+ int**& array,
+ int initVal,
+ int width,
+ int length
+) const
+{
+ int len = max(length,1);
+ lmp->memory->grow(array, len, width, "o2o:int**");
+ for (int i = 0; i < len; i++)
+ for (int j = 0; j < width; j++)
+ array[i][j] = initVal;
+}
+
+void twoWayOne2One::allocateArray
+(
+ int**& array,
+ int initVal,
+ int width,
+ const char* length
+) const
+{
+ int len = max(particleCloud_.numberOfParticles(),1);
+ lmp->memory->grow(array, len, width, "o2o:int**:autolen");
+ for (int i = 0; i < len; i++)
+ for (int j = 0; j < width; j++)
+ array[i][j] = initVal;
+}
+
+void inline twoWayOne2One::destroy(int** array,int len) const
+{
+ lmp->memory->destroy(array);
+}
+
+//============
+// double *
+void twoWayOne2One::allocateArray(double*& array, double initVal, int length) const
+{
+ int len = max(length,1);
+ lmp->memory->grow(array, len, "o2o:dbl*");
+ for (int i = 0; i < len; i++)
+ array[i] = initVal;
+}
+
+void inline twoWayOne2One::destroy(double* array) const
+{
+ lmp->memory->destroy(array);
+}
+
+//==============
+// int *
+void twoWayOne2One::allocateArray(int*& array, int initVal, int length) const
+{
+ int len = max(length,1);
+ lmp->memory->grow(array, len, "o2o:int*");
+ for (int i = 0; i < len; i++)
+ array[i] = initVal;
+}
+
+void inline twoWayOne2One::destroy(int* array) const
+{
+ lmp->memory->destroy(array);
+}
+//==============
+
+bool twoWayOne2One::couple(int i)
+{
+ bool coupleNow = false;
+ if (i==0)
+ {
+ couplingStep_++;
+ coupleNow = true;
+
+
+ // run commands from liggghtsCommands dict
+ Info<< "Starting up LIGGGHTS" << endl;
+ particleCloud_.clockM().start(3,"LIGGGHTS");
+
+ // check if liggghtsCommandModels with exaxt timing are being run
+ bool exactTiming(false);
+ int runComNr = -10;
+ DynamicList interruptTimes(0);
+ DynamicList DEMstepsToInterrupt(0);
+ DynamicList lcModel(0);
+
+ forAll(particleCloud_.liggghtsCommandModelList(),i)
+ {
+ // Check if exact timing is needed
+ // get time for execution
+ // store time for execution in list
+ if(particleCloud_.liggghtsCommand(i).exactTiming())
+ {
+ exactTiming = true;
+ DynamicList h
+ = particleCloud_.liggghtsCommand(i).executionsWithinPeriod
+ (
+ TSstart(),
+ TSend()
+ );
+
+ forAll(h,j)
+ {
+ // save interrupt times (is this necessary)
+ interruptTimes.append(h[j]);
+
+ // calc stepsToInterrupt
+ DEMstepsToInterrupt.append(DEMstepsTillT(h[j]));
+
+ // remember which liggghtsCommandModel to run
+ lcModel.append(i);
+ }
+
+ // make cumulative
+ label len = DEMstepsToInterrupt.size();
+ label ind(0);
+ forAll(DEMstepsToInterrupt,i)
+ {
+ ind = len - i - 1;
+ if(ind > 0)
+ {
+ DEMstepsToInterrupt[ind] -= DEMstepsToInterrupt[ind-1];
+ }
+ }
+
+ Info<< "Foam::twoWayOne2One::couple(i): interruptTimes=" << interruptTimes << nl
+ << "Foam::twoWayOne2One::couple(i): DEMstepsToInterrupt=" << DEMstepsToInterrupt << nl
+ << "Foam::twoWayOne2One::couple(i): lcModel=" << lcModel
+ << endl;
+ }
+
+ if(particleCloud_.liggghtsCommand(i).type() == "runLiggghts")
+ {
+ runComNr = i;
+ }
+ }
+
+ // models with exact timing exists
+ label commandLines(0);
+ if(exactTiming)
+ {
+ // extension for more liggghtsCommands active the same time:
+ // sort interrupt list within this run period
+ // keep track of corresponding liggghtsCommand
+ int DEMstepsRun(0);
+
+ forAll(interruptTimes,j)
+ {
+ // set run command till interrupt
+ DEMstepsRun += DEMstepsToInterrupt[j];
+ particleCloud_.liggghtsCommand(runComNr).set(DEMstepsToInterrupt[j]);
+ const char* command = particleCloud_.liggghtsCommand(runComNr).command(0);
+ Info<< "Executing run command: '"<< command <<"'"<< endl;
+ lmp->input->one(command);
+
+ // run liggghts command with exact timing
+ command = particleCloud_.liggghtsCommand(lcModel[j]).command(0);
+ Info << "Executing command: '"<< command <<"'"<< endl;
+ lmp->input->one(command);
+ }
+
+ // do the run
+ if(particleCloud_.liggghtsCommand(runComNr).runCommand(couplingStep()))
+ {
+ particleCloud_.liggghtsCommand(runComNr).set(couplingInterval() - DEMstepsRun);
+ const char* command = particleCloud_.liggghtsCommand(runComNr).command(0);
+ Info<< "Executing run command: '"<< command <<"'"<< endl;
+ lmp->input->one(command);
+ }
+
+ // do the other non exact timing models
+ forAll(particleCloud_.liggghtsCommandModelList(),i)
+ {
+ if
+ (
+ ! particleCloud_.liggghtsCommand(i).exactTiming() &&
+ particleCloud_.liggghtsCommand(i).runCommand(couplingStep())
+ )
+ {
+ commandLines=particleCloud_.liggghtsCommand(i).commandLines();
+ for(int j=0;jinput->one(command);
+ }
+ }
+ }
+ }
+ // no model with exact timing exists
+ else
+ {
+ forAll(particleCloud_.liggghtsCommandModelList(),i)
+ {
+ if(particleCloud_.liggghtsCommand(i).runCommand(couplingStep()))
+ {
+ commandLines=particleCloud_.liggghtsCommand(i).commandLines();
+ for(int j=0;jinput->one(command);
+ }
+ }
+ }
+ }
+
+ particleCloud_.clockM().stop("LIGGGHTS");
+ Info<< "LIGGGHTS finished" << endl;
+
+ if (!staticProcMap_)
+ {
+ createProcMap();
+ }
+
+ setupLig2FoamCommunication();
+
+ locateParticles();
+
+ setupFoam2LigCommunication();
+
+ if (verbose_)
+ {
+ Pout<< "FOAM owns " << getNumberOfParticles()
+ << " LIG owns " << lmp->atom->nlocal
+ << nl
+ << "FOAM collects " << lig2foam_->ncollected_
+ << " LIG collects " << foam2lig_->ncollected_
+ << endl;
+ }
+ }
+
+ return coupleNow;
+}
+
+void twoWayOne2One::setupLig2FoamCommunication()
+{
+ int* src_procs = new int[thisLigPartner_.size()];
+ for (int proci = 0; proci < thisLigPartner_.size(); proci++)
+ {
+ src_procs[proci] = thisLigPartner_[proci];
+ }
+ int* dst_procs = new int[thisFoamPartner_.size()];
+ for (int proci = 0; proci < thisFoamPartner_.size(); proci++)
+ {
+ dst_procs[proci] = thisFoamPartner_[proci];
+ }
+
+ delete lig2foam_;
+ lig2foam_ = new One2One(comm_liggghts_);
+ lig2foam_->setup
+ (
+ thisLigPartner_.size(),
+ src_procs,
+ thisFoamPartner_.size(),
+ dst_procs,
+ lmp->atom->nlocal
+ );
+ allocateArray
+ (
+ lig2foam_vec_tmp_,
+ 0.,
+ 3 * lig2foam_->ncollected_
+ );
+ allocateArray
+ (
+ lig2foam_scl_tmp_,
+ 0.,
+ lig2foam_->ncollected_
+ );
+}
+
+
+void twoWayOne2One::locateParticles()
+{
+ // get positions for locate
+ double** my_positions = static_cast(lmp->atom->x);
+ double* my_flattened_positions = nullptr;
+ allocateArray(my_flattened_positions, 0., 3*lmp->atom->nlocal);
+ for (int atomi = 0; atomi < lmp->atom->nlocal; atomi++)
+ {
+ for (int coordi = 0; coordi < 3; coordi++)
+ {
+ my_flattened_positions[atomi*3+coordi] = my_positions[atomi][coordi];
+ }
+ }
+
+ double* collected_flattened_positions = nullptr;
+ allocateArray(collected_flattened_positions, 0., 3*lig2foam_->ncollected_);
+
+ lig2foam_->exchange(my_flattened_positions, collected_flattened_positions, 3);
+ destroy(my_flattened_positions);
+
+ double* my_prev_cell_ids = nullptr;
+ double* prev_cell_ids = nullptr;
+ if (cellIdComm_)
+ {
+ my_prev_cell_ids = my_prev_cell_ids_fix_->vector_atom;
+ allocateArray(prev_cell_ids, -1, lig2foam_->ncollected_);
+ lig2foam_->exchange(my_prev_cell_ids, prev_cell_ids);
+ }
+
+ if (lig2foam_mask_)
+ {
+ delete [] lig2foam_mask_;
+ }
+ lig2foam_mask_ = new bool[lig2foam_->ncollected_];
+
+ DynamicList