Initialize at openlb-1-3

5 files changed, 855 insertions, 0 deletions

diff --git a/examples/particles/magneticParticles3d/Makefile b/examples/particles/magneticParticles3d/Makefile
new file mode 100644
index 0000000..a953954
--- /dev/null
+++ b/examples/particles/magneticParticles3d/Makefile
@@ -0,0 +1,105 @@
+#  This file is part of the OpenLB library
+#
+#  Copyright (C) 2007 Mathias Krause
+#  E-mail contact: info@openlb.net
+#  The most recent release of OpenLB can be downloaded at
+#  <http://www.openlb.net/>
+#
+#  This program 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 2
+#  of the License, or (at your option) any later version.
+#
+#  This program 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 program; if not, write to the Free
+#  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+#  Boston, MA  02110-1301, USA.
+
+###########################################################################
+## definitions
+
+include definitions.mk
+
+include $(ROOT)/global.mk
+
+OBJECTS    := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.o))
+DEPS       := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.d))
+
+###########################################################################
+## all
+
+all : depend compile link
+
+
+###########################################################################
+## dependencies
+
+depend : $(DEPS)
+
+$(PWD)/%.d : %.cpp
+	@echo Create dependencies for $<
+	@$(SHELL) -ec '$(CXX) -M $(CXXFLAGS) $(IDIR) $< \
+                       | sed -e "s!$*\.o!$(PWD)\/$*\.o!1" > .tmpfile; \
+                       cp -f .tmpfile $@;'
+
+###########################################################################
+## compile
+
+compile : $(OBJECTS)
+
+$(PWD)/%.o: %.cpp
+	@echo Compile $<
+	$(CXX) $(CXXFLAGS) $(IDIR) -c $< -o $@
+
+###########################################################################
+## clean
+
+clean : cleanrub cleanobj cleandep
+
+cleanrub:
+	@echo Clean rubbish files
+	@rm -f *~ core .tmpfile tmp/*.* $(OUTPUT)
+	@rm -f tmp/vtkData/*.* tmp/vtkData/data/*.* tmp/imageData/*.* tmp/gnuplotData/*.* tmp/gnuplotData/data/*.*
+
+cleanobj:
+	@echo Clean object files
+	@rm -f $(OBJECTS)
+
+cleandep:
+	@echo Clean dependencies files
+	@rm -f $(DEPS)
+
+cleanbuild:
+	@cd $(ROOT); \
+	 $(MAKE) cleanlib;
+
+###########################################################################
+## update lib
+
+$(ROOT)/$(LIBDIR)/lib$(LIB).a :
+	@cd $(ROOT); \
+	$(MAKE) all
+
+###########################################################################
+## link
+
+link: $(OUTPUT)
+
+$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a
+	@echo Link $@
+	$(CXX) $(foreach file, $(SRC), $(file:.cpp=.o)) $(LDFLAGS) -L$(ROOT)/$(LIBDIR) -l$(LIB) -lz -o $@
+
+###########################################################################
+## include dependencies
+
+ifneq "$(strip $(wildcard *.d))" ""
+   include $(foreach file,$(DEPS),$(file))
+endif
+
+###########################################################################
+###########################################################################
diff --git a/examples/particles/magneticParticles3d/definitions.mk b/examples/particles/magneticParticles3d/definitions.mk
new file mode 100644
index 0000000..02bf7e1
--- /dev/null
+++ b/examples/particles/magneticParticles3d/definitions.mk
@@ -0,0 +1,30 @@
+#  This file is part of the OpenLB library
+#
+#  Copyright (C) 2007 Mathias Krause
+#  E-mail contact: info@openlb.net
+#  The most recent release of OpenLB can be downloaded at
+#  <http://www.openlb.net/>
+#
+#  This program 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 2
+#  of the License, or (at your option) any later version.
+#
+#  This program 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 program; if not, write to the Free
+#  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+#  Boston, MA  02110-1301, USA.
+
+
+###########################################################################
+###########################################################################
+## DEFINITIONS TO BE CHANGED
+
+ROOT            := ../../..
+SRC             := magneticParticles3d.cpp
+OUTPUT          := magneticParticles3d
diff --git a/examples/particles/magneticParticles3d/magneticParticles3d.cpp b/examples/particles/magneticParticles3d/magneticParticles3d.cpp
new file mode 100644
index 0000000..673733d
--- /dev/null
+++ b/examples/particles/magneticParticles3d/magneticParticles3d.cpp
@@ -0,0 +1,644 @@
+/*  Lattice Boltzmann sample, written in C++, using the OpenLB
+ *  library
+ *
+ *  Copyright (C) 2018 Sascha Janz, Marie-Luise Maier
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program 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 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program 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 program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+ */
+
+/* magneticParticles3d.cpp:
+ * High-gradient magnetic separation is a method
+ * to separate ferromagnetic particles from a suspension.
+ * The simulation shows the deposition of magnetic particles
+ * on a single magnetized wire and models
+ * the magnetic separation step of the complete process.
+ */
+
+
+#include "olb3D.h"
+#ifndef OLB_PRECOMPILED // Unless precompiled version is used
+#include "olb3D.hh"     // Include full template code
+#endif
+
+using namespace olb;
+using namespace olb::descriptors;
+using namespace olb::graphics;
+using namespace olb::util;
+using namespace std;
+
+#define DESCRIPTOR D3Q19<>
+
+#define PARTICLE MagneticParticle3D
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+typedef double T;
+
+int iT = 0;
+
+// simulation geometry dimensions
+T geometryLengthX = 3.e-3 ; // length x-dircetion [m]
+T geometryLengthY = 2.e-3 ; // length y-dircetion [m]
+T geometryLengthZ = 1.e-3 ; // length z-dircetion [m]
+
+// magnetic particle paramters
+T pRadius = 1.e-5 ;  // particle radius [m]
+T pDensity = 3250.; // particle density [kg / m^3]
+T pVolume = 4. / 3. * M_PI * std::pow(pRadius, 3.); // particle volume [m^3]
+T pMass = pVolume * pDensity; // particle mass [kg]
+T pMag = 67. * pDensity * 0.06; // particle sat. magnetization [A / m]
+T elastModulus = 2.e4; // elastic modulus [Pa]
+T poissonRatio = 0.3; // poisson's ratio [-]
+T shearModulus = elastModulus / (2. * (1. + poissonRatio)); // shear modulus [Pa]
+
+// magnetic wire paramters
+T wRadius = 5.e-4;  // wire radius [m]
+T wLength = geometryLengthZ; // wire length [m]
+T wDensity = 7874.; // wire (iron) density [kg / m^3]
+T wMag = 215. * wDensity * 8.5e-3; // wire sat. magnetization [A / m]
+// mass-specific sat. magnetization iron = 215. [A m^2 / kg]
+
+// Rayleigh time step
+// can be used as an orientation value for particle time step sizes
+T eta = 0.8766 + 0.1631 * poissonRatio; // [-]
+T t_Rayleigh = (M_PI * pRadius / eta) * std::sqrt(pDensity / shearModulus); // [s]
+
+void prepareGeometry(UnitConverter<T, DESCRIPTOR> const& converter, IndicatorF3D<T>& wire,
+                     IndicatorF3D<T>& indicator, IndicatorF3D<T>& extendedDomain,
+                     SuperGeometry3D<T>& superGeometry)
+{
+
+  OstreamManager clout(std::cout, "prepareGeometry");
+  clout << "Prepare Geometry ..." << std::endl;
+
+  superGeometry.rename(0, 2, extendedDomain);
+  superGeometry.rename(2, 1, indicator);
+  superGeometry.rename(1, 5, wire);
+
+  T minR = superGeometry.getStatistics().getMinPhysR(2)[0];
+  minR -= 0.5 * converter.getConversionFactorLength();
+  std::vector < T > center = superGeometry.getStatistics().getCenterPhysR(2);
+  std::vector<T> physExtend = superGeometry.getStatistics().getPhysExtend(1);
+
+  Vector<T, 3> origin = { minR, center[1], center[2] };
+  Vector<T, 3> extend = { converter.getConversionFactorLength(), physExtend[1], physExtend[2] };
+
+  IndicatorCuboid3D<T> inlet(extend[0], extend[1], extend[2], origin);
+  origin[0] += superGeometry.getStatistics().getPhysExtend(2)[0];
+  IndicatorCuboid3D<T> outlet(extend[0], extend[1], extend[2], origin);
+
+  // rename the material at the inlet
+  superGeometry.rename(2, 3, inlet);
+  // rename the material at the outlet
+  superGeometry.rename(2, 4, outlet);
+
+  superGeometry.clean();
+  superGeometry.innerClean();
+  superGeometry.checkForErrors();
+  superGeometry.print();
+
+  clout << "Prepare Geometry ... OK" << std::endl;
+}
+
+/// Set up the geometry of the simulation
+void prepareLattice(SuperLattice3D<T, DESCRIPTOR>& sLattice,
+                    UnitConverter<T, DESCRIPTOR> const& converter, IndicatorF3D<T>& wire,
+                    Dynamics<T, DESCRIPTOR>& bulkDynamics,
+                    sOnLatticeBoundaryCondition3D<T, DESCRIPTOR>& sOnBC,
+                    sOffLatticeBoundaryCondition3D<T, DESCRIPTOR>& offBc,
+                    SuperGeometry3D<T>& superGeometry)
+{
+  OstreamManager clout(std::cout, "prepareLattice");
+  clout << "Prepare Lattice ..." << std::endl;
+
+  const T omega = (1. / converter.getLatticeRelaxationTime());
+
+  /// Material=0 -->do nothing
+  sLattice.defineDynamics(superGeometry, 0,
+                          &instances::getNoDynamics<T, DESCRIPTOR>());
+
+  /// Material=1 -->bulk dynamics
+  sLattice.defineDynamics(superGeometry, 1, &bulkDynamics);
+
+  /// Material=2 -->bounce back
+  sLattice.defineDynamics(superGeometry, 2,
+                          &instances::getBounceBack<T, DESCRIPTOR>());
+
+  sLattice.defineDynamics(superGeometry, 3, &bulkDynamics);
+  sLattice.defineDynamics(superGeometry, 4, &bulkDynamics);
+
+  /// Material=5 -->do nothing
+  sLattice.defineDynamics(superGeometry, 5,
+                          &instances::getNoDynamics<T, DESCRIPTOR>());
+  offBc.addZeroVelocityBoundary(superGeometry, 5, wire);
+
+  // boundary conditions for fluid
+
+  // inlet
+  sOnBC.addVelocityBoundary(superGeometry, 3, omega);
+  // outlet
+  sOnBC.addPressureBoundary(superGeometry, 4, omega);
+
+  // initialisation
+  AnalyticalConst3D<T, T> roh(1.);
+  AnalyticalConst3D<T, T> u0(0., 0., 0.);
+
+  sLattice.defineRhoU(superGeometry, 1, roh, u0);
+  sLattice.iniEquilibrium(superGeometry, 1, roh, u0);
+  sLattice.defineRhoU(superGeometry, 3, roh, u0);
+  sLattice.iniEquilibrium(superGeometry, 3, roh, u0);
+  sLattice.defineRhoU(superGeometry, 4, roh, u0);
+  sLattice.iniEquilibrium(superGeometry, 4, roh, u0);
+
+  sLattice.initialize();
+
+  clout << "Prepare Lattice ... OK" << std::endl;
+}
+
+void setBoundaryValues(SuperLattice3D<T, DESCRIPTOR>& sLattice,
+                       UnitConverter<T, DESCRIPTOR> const& converter, SuperParticleSystem3D<T, PARTICLE>& spSys,
+                       int iT, int itStartScaleT, SuperGeometry3D<T>& superGeometry, bool outNS, bool outLP)
+{
+
+  OstreamManager clout(std::cout, "setBoundaryValues");
+  std::vector < T > maxVelocity(3, T());
+  T distanceToBoundary = converter.getConversionFactorLength() / 2.;
+
+  if (outNS && iT <= itStartScaleT) {
+
+    SinusStartScale<T, int> startScale(itStartScaleT, T(1));
+    int help[1] = { iT };
+    T frac[3] = { T() };
+    startScale(frac, help);
+
+    maxVelocity[0] = converter.getCharPhysVelocity() * frac[0];
+
+    // outlet
+    RectanglePoiseuille3D<T> u3(superGeometry, 3, maxVelocity,
+                                distanceToBoundary, distanceToBoundary, distanceToBoundary);
+
+    sLattice.defineU(superGeometry, 3, u3);
+  }
+
+  // to reset and adopt boundaries in case of loaded fluid data
+  if (outLP && iT == 0) {
+
+    maxVelocity[0] = converter.getCharPhysVelocity();
+
+    // inlet
+    RectanglePoiseuille3D<T> u3(superGeometry, 3, maxVelocity,
+                                distanceToBoundary, distanceToBoundary, distanceToBoundary);
+
+    clout << "BC for loaded fluid data is reset on mat 3 " << std::endl;
+
+    sLattice.defineU(superGeometry, 3, u3);
+  }
+}
+
+void getResults(SuperGeometry3D<T>& superGeometry,
+                SuperLattice3D<T, DESCRIPTOR>& sLattice,
+                UnitConverter<T, DESCRIPTOR> const& converter,
+                AnalyticalF3D<T, S>& magForce, AnalyticalF3D<T, S>& magField,
+                SuperParticleSystem3D<T, PARTICLE>& spSys,
+                SuperParticleSysVtuWriterMag<T>& particleOut, Timer<double>& timer, int& iT,
+                int& itConsoleOutputFluid, int& itVtkOutputFluid,
+                int& itConsoleOutputMagParticles, int& itVtkOutputMagParticles,
+                bool& outNS, bool& outLP)
+{
+
+  OstreamManager clout(std::cout, "getResults");
+
+  // start vtm objects data until fluid reaches equilibrium
+  SuperVTMwriter3D<T> vtmWriterFluidStart("fluidReachingEquilibrium");
+
+  SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
+  SuperLatticeDensity3D<T, DESCRIPTOR> density(sLattice);
+  SuperLatticeFfromAnalyticalF3D<T, DESCRIPTOR> superMagPForceOne(magForce,
+      sLattice);
+  SuperLatticeFfromAnalyticalF3D<T, DESCRIPTOR> superMagPFieldOne (magField, sLattice);
+
+  vtmWriterFluidStart.addFunctor(velocity);
+  vtmWriterFluidStart.addFunctor(density);
+  vtmWriterFluidStart.addFunctor(superMagPForceOne);
+
+  if (outNS && iT == 0) {
+    SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
+    SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLattice);
+    SuperLatticeRank3D<T, DESCRIPTOR> rank(sLattice);
+
+    vtmWriterFluidStart.write(geometry);
+    vtmWriterFluidStart.write(cuboid);
+    vtmWriterFluidStart.write(rank);
+
+    vtmWriterFluidStart.createMasterFile();
+
+    converter.write("magnetics 1. loop");
+
+    clout << "fluid computation" << std::endl;
+  }
+  if (outLP && iT == 0) {
+
+    converter.write("magnetics 2. loop");
+    clout << "magParticles computation" << std::endl;
+  }
+
+  // === fluid
+  // console output fluid
+  if (outNS && iT % itConsoleOutputFluid == 0) {
+    timer.update(iT);
+    timer.print( iT );
+    // output for latticeStatistics
+    sLattice.getStatistics().print(iT, converter.getPhysTime(iT));
+  }
+  // vtk output fluid
+  if (outNS && iT % itVtkOutputFluid == 0) {
+    vtmWriterFluidStart.write(iT);
+  }
+
+  // === particles
+  // Writes the console output files of the magnetic particles
+  if (outLP && iT % itConsoleOutputMagParticles == 0) {
+    /// Lattice statistics console output
+    timer.print( iT );
+  }
+  // Writes the pvd files of the magnetic particles
+  if (outLP && iT % itVtkOutputMagParticles == 0) {
+    particleOut.write(iT);
+  }
+}
+
+int main(int argc, char* argv[])
+{
+
+  /// === 1st Step: Initialization ===
+  olbInit(&argc, &argv);
+  singleton::directories().setOutputDir("./tmp/");
+  OstreamManager clout(std::cout, "main");
+
+  string fName("magneticParticles3d.xml");
+  XMLreader config(fName);
+
+  // converter contains parameters of fluid simulation
+  UnitConverter<T, DESCRIPTOR>* converter = createUnitConverter<T, DESCRIPTOR>(config);
+
+  clout << "fluid converter: ..." << std::endl;
+  converter->print(); // gives overview of converter into console
+
+  // particles time step size
+  T physDeltaTParticles = t_Rayleigh * 0.4; // [s]
+  // particles relaxation time
+  T tau_particles = (physDeltaTParticles * 3 * converter->getPhysViscosity() /
+                     std::pow(converter->getConversionFactorLength(), 2.) + 0.5);
+
+  // converter contains parameters of particle simulation
+  UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converterParticles(
+  int {converter->getResolution()},    // resolution: number of voxels per charPhysL
+  (T) tau_particles, // latticeRelaxationTime: relaxation time, has to be greater than 0.5!
+  (T) converter->getCharPhysLength(), // (hydraulic raius [m]) charPhysLength: reference length of simulation geometry
+  (T) converter->getCharPhysVelocity(), // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
+  (T) converter->getPhysViscosity(),  // physViscosity: physical kinematic viscosity in __m^2 / s__
+  (T) converter->getPhysDensity()   // physDensity: physical density in __kg / m^3__
+  );
+
+  clout << "particle converter: ..." << std::endl;
+  converterParticles.print();
+
+  // name of particle data output file
+  string filename = "particleData.txt";
+  std::ofstream outputFile;
+
+  bool multiOutput = false;
+  config["Output"]["MultiOutput"].read(multiOutput);
+  clout.setMultiOutput(multiOutput);
+
+  std::string olbdir, outputdir;
+  config["Application"]["OlbDir"].read(olbdir);
+  config["Output"]["OutputDir"].read(outputdir);
+  singleton::directories().setOlbDir(olbdir);
+  singleton::directories().setOutputDir(outputdir);
+
+  /// Instantiation of a cuboidGeometry with weights
+#ifdef PARALLEL_MODE_MPI
+  const int noOfCuboids = singleton::mpi().getSize();
+#else
+  const int noOfCuboids = 7;
+#endif
+
+  // read physical simulation times from xml-file
+  T physFluidNST;   // time for fluid simulation
+  config["Application"]["SimulationTimes"]["PhysFluidNST"].read(physFluidNST);
+
+  T physStartScT;    // factor * physFluidNST;
+  config["Application"]["SimulationTimes"]["PhysStartScT"].read(physStartScT);
+
+  T physParticleT;
+  config["Application"]["SimulationTimes"]["PhysParticleT"].read(physParticleT);
+
+  // convert physical times to lattice time steps
+  int itFluidNST = converter->getLatticeTime(physFluidNST);
+  int itStartScaleT = converter->getLatticeTime(physStartScT);
+  int itParticleT = converterParticles.getLatticeTime(physParticleT);
+
+  int itConsoleOutputFluid = itFluidNST / 10. , itVtkOutputFluid = itFluidNST / 4.;
+  int itConsoleOutputMagParticles = itParticleT / 75., itVtkOutputMagParticles = itParticleT / 249.;
+
+  clout << "Fluid: physFluidNST = " << physFluidNST << "s itFluidNST = "
+        << itFluidNST << std::endl;
+
+  clout << "Particles: physParticleT = " << physParticleT << "s itParticleT = "
+        << itParticleT << std::endl;
+
+  clout << "noOfCuboids = " << noOfCuboids << std::endl;
+
+  /// === 2nd Step: Prepare Geometry ===
+
+  // A: Cuboid as Channel
+  std::vector < T > center(3, T());
+  std::vector<T> length(3, T());
+  length[0] = geometryLengthX;
+  length[1] = geometryLengthY;
+  length[2] = geometryLengthZ;
+
+  IndicatorCuboid3D<T> cuboid(length, center);
+  IndicatorLayer3D<T> extendedDomainCuboid(cuboid, converter->getConversionFactorLength());
+
+  std::vector<T> centerW(3, T());
+  centerW[0] = geometryLengthX - 0.675e-3;
+  centerW[1] = geometryLengthY / 2;
+  centerW[2] = geometryLengthZ / 2.;
+
+  std::vector<T> normalW(3, T());
+  normalW[0] = T(0);
+  normalW[1] = T(0);
+  normalW[2] = T(1);
+  wLength += 2 * converter->getConversionFactorLength();
+  IndicatorCylinder3D<T> wire(centerW, normalW, wRadius, wLength);
+
+  CuboidGeometry3D<T> cuboidGeometry(extendedDomainCuboid,
+                                     converter->getConversionFactorLength(), noOfCuboids);
+
+  /// Instantiation of a loadBalancer
+  HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
+
+  /// Instantiation of a superGeometry
+  SuperGeometry3D<T> superGeometry(cuboidGeometry, loadBalancer, 2);
+
+  prepareGeometry(*converter, wire, cuboid, extendedDomainCuboid,
+                  superGeometry);
+
+  /// === 3rd Step: Prepare Lattice ===
+
+  SuperLattice3D<T, DESCRIPTOR> sLattice(superGeometry);
+
+  BGKdynamics<T, DESCRIPTOR> bulkDynamics((1. / converter->getLatticeRelaxationTime()),
+                                          instances::getBulkMomenta<T, DESCRIPTOR>());
+
+  // choose between local and non-local boundary condition
+  sOnLatticeBoundaryCondition3D<T, DESCRIPTOR> sOnBoundaryCondition(sLattice);
+  createInterpBoundaryCondition3D<T, DESCRIPTOR>(sOnBoundaryCondition);
+
+  // for the velocity field around the wire
+  sOffLatticeBoundaryCondition3D<T, DESCRIPTOR> sOffBoundaryCondition(sLattice);
+  createBouzidiBoundaryCondition3D<T, DESCRIPTOR>(sOffBoundaryCondition);
+
+  // gives dynamics to cells
+  prepareLattice(sLattice, *converter, wire, bulkDynamics, sOnBoundaryCondition,
+                 sOffBoundaryCondition, superGeometry);
+
+  /// === 4th Step: Prepare Lagrange Particles
+
+  SuperParticleSystem3D<T, PARTICLE> spSys(cuboidGeometry, loadBalancer,
+      superGeometry);
+
+  // contact detection
+  NanoflannContact<T, PARTICLE> nanoflannContact(
+    *(spSys.getParticleSystems()[0]), 6. * pRadius);
+  spSys.setContactDetection(nanoflannContact);
+
+  // magnetic force field from wire
+  std::vector < T > origin = superGeometry.getStatistics().getCenterPhysR(1);
+  origin[0] = geometryLengthX - 0.675e-3;
+  origin[2] = superGeometry.getStatistics().getMinPhysR(1)[2];
+
+  std::vector<T> orientation(3, T());
+  orientation[0] = 1.;
+  T degree = 0.;
+
+  // car2cyl helps to know which cylinder coordinates fit to the cartesion ones
+  CartesianToCylinder3D<T, T> car2cyl(origin, degree, orientation);
+
+  /// Forces of Lagrange particles
+
+  // Object of magnetic force to know magnetic value at any place in geometry
+  // but needs to be called by operator function.
+  // inhomogeneous magnetic field round cylinder
+  MagneticForceFromCylinder3D<T, T> magForceFromCylOnParticle(car2cyl, wLength,
+      wRadius, 4., wMag, pMag, pRadius);
+  MagneticFieldFromCylinder3D<T, T> magFieldFromCylOnParticle(car2cyl, wLength,
+      wRadius, 4., wMag);
+
+  // external magnetic force from the wire
+  auto magneticPForceOne = make_shared
+                           < MagneticForceForMagP3D<T, PARTICLE, DESCRIPTOR>
+                           > (magForceFromCylOnParticle, magFieldFromCylOnParticle, 1.e0);
+  spSys.addForce(magneticPForceOne);
+
+  // interparticular magnetic force
+  auto interpMagF = make_shared
+                    < InterpMagForceForMagP3D<T, PARTICLE, DESCRIPTOR>
+                    > (1.e0, 1.e0);
+  spSys.addForce(interpMagF);
+
+  // dynamic viscosity
+  T dynVisc = converter->getPhysViscosity() * converter->getPhysDensity();
+
+  // damping force for the particle rotation
+  auto rotDampingF = make_shared
+                     < LinearDampingForceForMagDipoleMoment3D<T, PARTICLE, DESCRIPTOR>
+                     > (dynVisc, 1.e0);
+  spSys.addForce(rotDampingF);
+
+  // mechanic contact force
+  auto hertzMindlinContF = make_shared
+                           < HertzMindlinDeresiewicz3D<T, PARTICLE, DESCRIPTOR>
+                           > (shearModulus, shearModulus, poissonRatio, poissonRatio, 1.e0, 1.e0, false);
+  spSys.addForce(hertzMindlinContF);
+
+  // stokes drag force
+  SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR> getVel( sLattice, *converter );
+  auto stokesDragF = make_shared
+                     < StokesDragForce3D<T, PARTICLE, DESCRIPTOR>
+                     > (getVel, converterParticles);
+  spSys.addForce(stokesDragF);
+
+  /// Boundarys
+
+  T dT = converter->getConversionFactorTime() ;
+  std::set<int> wireBMaterial = { 5, 4 };
+  std::set<int> reflBMat = { 2, 3 };
+
+  // particle boundary for deposited particles
+  auto wireBoundary = make_shared
+                      < WireBoundaryForMagP3D<T, PARTICLE>
+                      > (superGeometry, wireBMaterial);
+  spSys.addBoundary(wireBoundary);
+
+  // particle boundary
+  auto materialreflectBoundary = make_shared
+                                 < SimpleReflectBoundary3D<T, PARTICLE>
+                                 > (dT, superGeometry, reflBMat);
+  spSys.addBoundary(materialreflectBoundary);
+
+  // cuboid overlap
+  spSys.setOverlap(2.*converter->getConversionFactorLength());
+
+  // particles vtu output
+  std::string particleOutputName = "vtuOutMagP";
+  std::bitset < 9 > properties;
+  properties.set();
+
+  SuperParticleSysVtuWriterMag<T> particleOut(spSys, particleOutputName, properties);
+
+  clout << "Number of Forces: " << spSys.numOfForces()[0] << std::endl;
+
+  /// Paramters for magnetic particle generation
+  std::vector<T> pPos = { 0., 0., 0. }; // position
+  std::vector<T> pVel = { 0., 0., 0. }; // velocity
+  std::vector<T> pAVel = { 0., 0., 0. }; // angular velocity
+  std::vector<T> pTrq = { 0., 0., 0. }; // torque
+  std::vector<T> pDMoment = { -1., 0., 0. }; // orientation mag. dipole moment
+  if(!util::nearZero(util::norm(pDMoment))){util::normalize(pDMoment);}
+  else {clout << "Norm of pDMoment near zero!" << endl; exit(0);}
+
+  // magnetic particle magPartTemplate is used as copy template
+  PARTICLE<T> magPartTemplate(pPos, pVel, pMass, pRadius, 0, pDMoment, pAVel, pTrq, pMag, 1);
+
+  // intialization geometry particles
+  std::vector<T> particlesOrigin(3, T(0));
+  particlesOrigin = superGeometry.getStatistics().getMinPhysR(1);
+  std::vector<T> particlesExtend(3, T(0));
+
+  particlesOrigin[0] = 1.e-4;
+  particlesOrigin[1] = (superGeometry.getStatistics().getMaxPhysR(1)[1]) * (1./2. - 1./6.);
+  particlesOrigin[2] = (superGeometry.getStatistics().getMaxPhysR(1)[2]) * (1./2. - 1./6.);
+  particlesExtend[0] = (superGeometry.getStatistics().getMaxPhysR(1)[0]) * 1./15. ;
+  particlesExtend[1] = (superGeometry.getStatistics().getMaxPhysR(1)[1]) * 2./6. ;
+  particlesExtend[2] = (superGeometry.getStatistics().getMaxPhysR(1)[2]) * 2./6. ;
+  IndicatorCuboid3D<T> cuboidPart(particlesExtend, particlesOrigin);
+
+  clout << "starting simulation" << endl;
+
+  /// === 5th Step: Main Loop with Timer ===
+
+  // is set on true for the appropriate loop of fluid / particle simulation
+  bool outNS = false;
+  bool outLP = false;
+
+  // === Calculating Navier-Stokes Fluid
+
+  // if there is no computed fluid data in a external data files
+  // compute fluid data in data file, else load it and jump to next loop
+  if (!(sLattice.load("water"))) {
+
+    clout << std::endl;
+    clout << "Fluid data has to be computed " << std::endl;
+    clout << std::endl;
+
+    outNS = true;
+
+    Timer<double> timerFluid(itFluidNST, superGeometry.getStatistics().getNvoxel());
+    timerFluid.start();
+
+    for (int iT = 0; iT < itFluidNST; ++iT) {
+
+      setBoundaryValues(sLattice, *converter, spSys, iT, itStartScaleT,
+                        superGeometry, outNS, outLP);
+
+      sLattice.collideAndStream();
+
+      getResults(superGeometry, sLattice, *converter,
+                 magForceFromCylOnParticle, magFieldFromCylOnParticle, spSys,
+                 particleOut, timerFluid, iT, itConsoleOutputFluid, itVtkOutputFluid,
+                 itConsoleOutputMagParticles, itVtkOutputMagParticles, outNS, outLP);
+    }
+
+    outNS = false;
+
+    // save computed fluid data in external data files
+    clout << "save fluid data " << std::endl;
+
+    sLattice.communicate();
+    sLattice.save("water");
+
+    timerFluid.stop();
+    timerFluid.printSummary();
+
+  } else {
+    sLattice.communicate();
+    sLattice.collideAndStream();
+    sLattice.getStatistics().print(iT, converter->getPhysTime(iT));
+  }
+
+  // === Calculating Particles
+
+  clout << std::endl;
+  clout << "particles computation" << std::endl;
+  clout << std::endl;
+
+  outLP = true;
+
+  // number of Particles in initialization geometry
+  int noOfPartX = 4 ;
+  int noOfPartY = 12 ;
+  int noOfPartZ = 6 ;
+  int noOfPart = noOfPartX * noOfPartY * noOfPartZ;
+
+  Timer<double> timerParticles(itParticleT, noOfPart);
+  timerParticles.start();
+
+  setBoundaryValues(sLattice, converterParticles, spSys, iT, itStartScaleT,
+                    superGeometry, outNS, outLP);
+
+  for (int iT = 0; iT < itParticleT; ++iT) {
+
+    // particles initialization
+    if (iT == 0) {
+
+      spSys.addParticleEquallyDistributed(cuboidPart, noOfPartX, noOfPartY, noOfPartZ, magPartTemplate);
+      spSys.setParticlesPosRandom(0.75 * pRadius, 0.75 * pRadius, 0.55 * pRadius);
+    }
+
+    getResults(superGeometry, sLattice, converterParticles,
+               magForceFromCylOnParticle, magFieldFromCylOnParticle, spSys,
+               particleOut, timerParticles, iT, itConsoleOutputFluid, itVtkOutputFluid,
+               itConsoleOutputMagParticles, itVtkOutputMagParticles,
+               outNS, outLP);
+
+    // particles simulation
+    spSys.simulate(converterParticles.getConversionFactorTime());
+
+  }
+
+  timerParticles.stop();
+  timerParticles.printSummary();
+
+  clout << "End of simulation" << std::endl;
+}
diff --git a/examples/particles/magneticParticles3d/magneticParticles3d.xml b/examples/particles/magneticParticles3d/magneticParticles3d.xml
new file mode 100644
index 0000000..8235e9f
--- /dev/null
+++ b/examples/particles/magneticParticles3d/magneticParticles3d.xml
@@ -0,0 +1,47 @@
+<Param>
+  <Application>
+    <Name> sedimentation </Name>
+    <dim> 3 </dim>
+    <OlbDir> ../../../ </OlbDir>
+    <Discretization>
+      <LatticeRelaxationTime> 0.75 </LatticeRelaxationTime>
+      <Resolution> 50 </Resolution>
+    </Discretization>
+    <PhysParameters>
+      <CharPhysLength> .001 </CharPhysLength>
+      <CharPhysVelocity> 0.002 </CharPhysVelocity>
+      <viscosityType> Kinematic viscosity </viscosityType>
+      <PhysViscosity> 1.e-6  </PhysViscosity>
+      <PhysDensity> 1000. </PhysDensity>
+      <CharPhysPressure> 0. </CharPhysPressure>
+    </PhysParameters>
+    <SimulationTimes>  
+      <PhysFluidNST>   .12 </PhysFluidNST>
+      <PhysStartScT>   .1 </PhysStartScT>   
+      <PhysParticleT>  1.55 </PhysParticleT>
+    </SimulationTimes>
+  </Application>
+  <Output>
+    <OutputDir> ./tmp/ </OutputDir>
+    <MultiOutput> false </MultiOutput>
+    <Log>
+      <Filename> magneticParticles3d.log </Filename>
+      <SaveTime> 1.  </SaveTime>
+    </Log>
+    <VisualizationVTK>
+      <Filename> magneticParticles3dvtk</Filename>  
+      <SaveTime> 0.01 </SaveTime>
+    </VisualizationVTK>
+    <VisualizationImages>
+      <Filename> magneticParticles3dimage</Filename>
+      <SaveTime> 0.01 </SaveTime>
+    </VisualizationImages>    
+    <Timer>
+      <PrintMode>0</PrintMode>
+      <SkipType> 0 </SkipType>
+      <PhysTime> 0.1 </PhysTime>
+      <TimeSteps> 0.1 </TimeSteps>
+      <PrintSummary> true </PrintSummary>
+    </Timer>
+  </Output>
+</Param>
diff --git a/examples/particles/magneticParticles3d/module.mk b/examples/particles/magneticParticles3d/module.mk
new file mode 100644
index 0000000..1190482
--- /dev/null
+++ b/examples/particles/magneticParticles3d/module.mk
@@ -0,0 +1,29 @@
+#  This file is part of the OpenLB library
+#
+#  Copyright (C) 2017 Markus Mohrhard
+#  E-mail contact: info@openlb.net
+#  The most recent release of OpenLB can be downloaded at
+#  <http://www.openlb.net/>
+#
+#  This program 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 2
+#  of the License, or (at your option) any later version.
+#
+#  This program 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 program; if not, write to the Free
+#  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+#  Boston, MA  02110-1301, USA.
+
+current_dir := $(dir $(word $(words $(MAKEFILE_LIST)),$(MAKEFILE_LIST)))
+
+include global.mk
+include rules.mk
+include $(addsuffix definitions.mk, $(current_dir))
+
+$(eval $(call sample,$(current_dir)$(OUTPUT),$(addprefix $(current_dir), $(SRC))))