From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 24 Jun 2019 14:43:36 +0200
Subject: Initialize at openlb-1-3
---
src/particles/hlbm/hlbmDynamics2D.h | 110 ++++++++++++
src/particles/hlbm/hlbmDynamics2D.hh | 261 ++++++++++++++++++++++++++++
src/particles/hlbm/hlbmDynamics3D.h | 108 ++++++++++++
src/particles/hlbm/hlbmDynamics3D.hh | 317 +++++++++++++++++++++++++++++++++++
4 files changed, 796 insertions(+)
create mode 100644 src/particles/hlbm/hlbmDynamics2D.h
create mode 100644 src/particles/hlbm/hlbmDynamics2D.hh
create mode 100644 src/particles/hlbm/hlbmDynamics3D.h
create mode 100644 src/particles/hlbm/hlbmDynamics3D.hh
(limited to 'src/particles/hlbm')
diff --git a/src/particles/hlbm/hlbmDynamics2D.h b/src/particles/hlbm/hlbmDynamics2D.h
new file mode 100644
index 0000000..afc7bb5
--- /dev/null
+++ b/src/particles/hlbm/hlbmDynamics2D.h
@@ -0,0 +1,110 @@
+/* DESCRIPTOR Boltzmann sample, written in C++, using the OpenLB
+ * library
+ *
+ * Copyright (C) 2006-2016 Thomas Henn, Fabian Klemens, Robin Trunk, Davide Dapelo
+ * E-mail contact: info@openlb.net
+ * The most recent release of OpenLB can be downloaded at
+ *
+ *
+ * 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.
+ */
+
+#ifndef PARTICLEDYNAMICS_2D_H
+#define PARTICLEDYNAMICS_2D_H
+
+#include "core/superLattice2D.h"
+#include "functors/analytical/indicator/smoothIndicatorBaseF2D.h"
+
+
+namespace olb {
+
+template
+class ParticleDynamics2D {
+private:
+ SuperLattice2D& _sLattice;
+ UnitConverter const& _converter;
+ SuperGeometry2D& _superGeometry;
+ std::shared_ptr > _indicatorF;
+ std::vector* > _vectorOfIndicator;
+ T _lengthX;
+ T _lengthY;
+ Vector _accExt;
+ bool _escapeFromDomain=false;
+ bool _oldWallCollision=true;
+
+public:
+ ParticleDynamics2D(SuperLattice2D& sLattice,
+ UnitConverter const& converter,
+ SuperGeometry2D& superGeometry,
+ T lengthX, T lengthY, Vector accExt = Vector (0.,0.))
+ : _sLattice(sLattice),
+ _converter(converter),
+ _superGeometry(superGeometry),
+ _lengthX(lengthX),
+ _lengthY(lengthY),
+ _accExt(accExt)
+ {}
+
+ ParticleDynamics2D(SuperLattice2D& sLattice,
+ UnitConverter const& converter,
+ SuperGeometry2D& superGeometry, std::shared_ptr > indicatorF,
+ bool escapeFromDomain, bool oldWallCollision=false,
+ Vector accExt = Vector (0.,0.))
+ : _sLattice(sLattice),
+ _converter(converter),
+ _superGeometry(superGeometry),
+ _indicatorF(indicatorF),
+ _accExt(accExt),
+ _escapeFromDomain(escapeFromDomain),
+ _oldWallCollision(oldWallCollision)
+ {}
+
+ void addCircle(Vector< T, 2> center, T radius, T density, T epsilon,
+ Vector vel = Vector (0.,0.));
+
+ void addCuboid(Vector< T, 2> center, T xLength, T yLength, T density,
+ T epsilon, T theta=0, Vector vel = Vector (0.,0.));
+
+ void addTriangle(Vector< T, 2> center, T radius, T density, T epsilon,
+ T theta, Vector vel = Vector (0.,0.));
+
+ void addParticle(SmoothIndicatorF2D& indicator);
+
+ void computeBoundaryForce(std::vector* >& indicator);
+
+ void addWallColl(SmoothIndicatorF2D& indicator, T delta);
+
+ void verletIntegration(SmoothIndicatorF2D& indicator);
+
+ void updateParticleDynamics(std::string name, SmoothIndicatorF2D& indicator);
+
+ void checkAndRemoveEscaped();
+
+ void addParticleField(SmoothIndicatorF2D& indicator);
+
+ void simulateTimestep(std::string name);
+
+ void print();
+
+ void load(std::string filename, T epsilon);
+
+ void save(std::string filename);
+
+};
+
+}
+
+#endif
diff --git a/src/particles/hlbm/hlbmDynamics2D.hh b/src/particles/hlbm/hlbmDynamics2D.hh
new file mode 100644
index 0000000..95e996f
--- /dev/null
+++ b/src/particles/hlbm/hlbmDynamics2D.hh
@@ -0,0 +1,261 @@
+/* DESCRIPTOR Boltzmann sample, written in C++, using the OpenLB
+ * library
+ *
+ * Copyright (C) 2006-2016 Thomas Henn, Fabian Klemens, Robin Trunk, Davide Dapelo
+ * E-mail contact: info@openlb.net
+ * The most recent release of OpenLB can be downloaded at
+ *
+ *
+ * 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.
+ */
+
+
+
+#ifndef PARTICLEDYNAMICS_2D_HH
+#define PARTICLEDYNAMICS_2D_HH
+
+#include "functors/analytical/indicator/indicCalc2D.h"
+#include "functors/lattice/superLatticeLocalF2D.h"
+#include "hlbmDynamics2D.h"
+
+namespace olb {
+
+template
+void ParticleDynamics2D::addCircle(Vector< T, 2> center, T radius, T density, T epsilon, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorCircle2D(center, radius, epsilon, density, vel) );
+}
+
+template
+void ParticleDynamics2D::addCuboid(Vector< T, 2> center, T xLength, T yLength, T density, T epsilon, T theta, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorCuboid2D(center, xLength, yLength, epsilon, theta, density, vel) );
+}
+
+template
+void ParticleDynamics2D::addTriangle(Vector< T, 2> center, T radius, T density, T epsilon, T theta, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorTriangle2D(center, radius, epsilon, theta, density, vel) );
+}
+
+template
+void ParticleDynamics2D::addParticle(SmoothIndicatorF2D& indicator)
+{
+ _vectorOfIndicator.push_back(&indicator);
+}
+
+template
+void ParticleDynamics2D::computeBoundaryForce(std::vector* >& indicator)
+{
+ SuperLatticePorousMomentumLossForce2D force(_sLattice, _superGeometry, indicator,_converter);
+ T sumF[force.getTargetDim()];
+ for (int i=0; i* >::size_type iInd=0; iInd!=indicator.size(); iInd++) {
+ /// get particle acceleration through boundary force and gravity (and buoyancy)
+ Vector acceleration2;
+ Vector force;
+ force[0] = sumF[0+4*iInd];
+ force[1] = sumF[1+4*iInd];
+ acceleration2[0] = sumF[0+4*iInd] / indicator[iInd]->getMass() + _accExt[0];
+ acceleration2[1] = sumF[1+4*iInd] / indicator[iInd]->getMass() + _accExt[1];
+ indicator[iInd]->setForce( force );
+ indicator[iInd]->setAcc2( acceleration2 );
+ indicator[iInd]->setAlpha2( sumF[2+4*iInd] / indicator[iInd]->getMofi() );
+ }
+}
+
+template
+void ParticleDynamics2D::addWallColl(SmoothIndicatorF2D& indicator, T delta)
+{
+ T w1 = 1e-7 / 2.;
+ T w = 1e-7 / 2.;
+
+ T rad = indicator.getRadius();
+ T massInv = 1. / indicator.getMass();
+
+ std::vector dx(2, T());
+ dx[0] = _lengthX - indicator.getPos()[0];
+ dx[1] = _lengthY - indicator.getPos()[1];
+
+ for (int i = 0; i < 2; i++) {
+ if (dx[i] <= rad) {
+ indicator.getAcc2()[i] += massInv * -dx[i] * (rad - dx[i]) / w1;
+ }
+ if (indicator.getPos()[i] <= rad) {
+ indicator.getAcc2()[i] += massInv * indicator.getPos()[i] * (rad - indicator.getPos()[i]) / w1;
+ }
+ if (dx[i] > rad && dx[i] <= rad + delta) {
+ indicator.getAcc2()[i] += massInv * -dx[i] * std::pow((rad + delta - dx[i]), 2) / w;
+ }
+ if (indicator.getPos()[i] > rad && indicator.getPos()[i] <= rad + delta) {
+ indicator.getAcc2()[i] += massInv * indicator.getPos()[i] * std::pow((rad + delta - indicator.getPos()[i]), 2) / w;
+ }
+ }
+}
+
+template
+void ParticleDynamics2D::verletIntegration(SmoothIndicatorF2D& indicator)
+{
+ T time = _converter.getConversionFactorTime();
+ T time2 = time*time;
+
+ Vector position, velocity;
+ Vector rotationMatrix;
+ for (int i=0; i<2; i++) {
+ position[i] = indicator.getPos()[i] + indicator.getVel()[i] * time + (0.5 * indicator.getAcc()[i] * time2);
+ T avgAcc = (indicator.getAcc()[i] + indicator.getAcc2()[i]) * 0.5;
+ velocity[i] = indicator.getVel()[i] + avgAcc * time;
+ }
+ indicator.setPos(position);
+ indicator.setVel(velocity);
+ indicator.setAcc(indicator.getAcc2());
+
+ indicator.setTheta( std::fmod( indicator.getTheta() + indicator.getOmega() * time + (0.5 * indicator.getAlpha() * time2), 2.*M_PI ) );
+ T avgAlpha = (indicator.getAlpha() + indicator.getAlpha2()) * 0.5;
+ indicator.setOmega( indicator.getOmega() + avgAlpha * time);
+ indicator.setAlpha( indicator.getAlpha2() );
+
+ T cos = std::cos(indicator.getTheta());
+ T sin = std::sin(indicator.getTheta());
+
+ rotationMatrix[0] = cos;
+ rotationMatrix[1] = -sin;
+ rotationMatrix[2] = sin;
+ rotationMatrix[3] = cos;
+ indicator.setRotationMatrix(rotationMatrix);
+}
+
+template
+void ParticleDynamics2D::updateParticleDynamics(std::string name, SmoothIndicatorF2D& indicator)
+{
+ this->verletIntegration(indicator);
+}
+
+template
+void ParticleDynamics2D::checkAndRemoveEscaped()
+{
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ auto internal = std::make_shared > (*_indicatorF.get(), -_converter.getConversionFactorLength()+(**i).getEpsilon() );
+ IndicMinus2D layer(_indicatorF, internal);
+ SuperLatticeIndicatorSmoothIndicatorIntersection2D intersection(
+ _sLattice, _superGeometry, layer, **i );
+ int input[1];
+ T output[1] = {0.};
+ intersection(output, input);
+ if (output[0] == 1) {
+ _vectorOfIndicator.erase(i);
+ if (i==_vectorOfIndicator.end() ) {
+ break;
+ }
+ }
+ }
+}
+
+
+template
+void ParticleDynamics2D::addParticleField(SmoothIndicatorF2D& indicator)
+{
+ /// Analytical2D functor for particle motion (trans+rot)
+ ParticleU2D velocity(indicator, _converter);
+ _sLattice.setExternalParticleField(_superGeometry, velocity, indicator);
+}
+
+template
+void ParticleDynamics2D::simulateTimestep(std::string name)
+{
+ // Remove particles from domain if _escapeFromDomain is toggled
+ if (_escapeFromDomain) checkAndRemoveEscaped();
+
+ // Compute force acting on particles boundary
+ computeBoundaryForce(_vectorOfIndicator);
+
+ // Update particle dynamics and porous particle field
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ updateParticleDynamics(name, **i);
+ addParticleField(**i);
+ }
+}
+
+template
+void ParticleDynamics2D::print()
+{
+ OstreamManager clout(std::cout, "ParticleDynamics2D");
+ clout << "Number of particles=" << _vectorOfIndicator.size() << std::endl;
+ int count = 0;
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ clout << "Particle " << ++count << " (" << (**i).name() << "):" << std::endl;
+ clout << " |Circum radius(m)= " << std::setw(13) << (**i).getCircumRadius() <> name;
+
+ if (name == "circle") {
+ T radius = T(), mass = T();
+ Vector center = {T(), T()};
+ Vector vel = {T(), T()};
+ iout >> center[0] >> center[1] >> radius >> mass >> vel[0] >> vel[1];
+ addCircle(center, radius, mass, epsilon, vel);
+ }
+ }
+
+ iout.close();
+}
+
+// WORKS ONLY FOR CIRCLES FOR NOW!!
+template
+void ParticleDynamics2D::save(std::string filename)
+{
+ std::ofstream fout( (singleton::directories().getLogOutDir() + filename).c_str() );
+
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ if ( (**i).name() == "circle") {
+ fout << (**i).name() << " "
+ << (**i).getPos()[0] << " " << (**i).getPos()[1] << " "
+ << (**i).getCircumRadius() << " " << (**i).getMass() << " "
+ << (**i).getVel()[0] << " " << (**i).getVel()[1]
+ << std::endl;
+ }
+ }
+
+ fout.close();
+}
+
+}
+
+#endif
diff --git a/src/particles/hlbm/hlbmDynamics3D.h b/src/particles/hlbm/hlbmDynamics3D.h
new file mode 100644
index 0000000..b39219b
--- /dev/null
+++ b/src/particles/hlbm/hlbmDynamics3D.h
@@ -0,0 +1,108 @@
+/* DESCRIPTOR Boltzmann sample, written in C++, using the OpenLB
+ * library
+ *
+ * Copyright (C) 2006-2016 Thomas Henn, Fabian Klemens, Robn Trunk, Davide Dapelo
+ * E-mail contact: info@openlb.net
+ * The most recent release of OpenLB can be downloaded at
+ *
+ *
+ * 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.
+ */
+
+#ifndef PARTICLEDYNAMICS_3D_H
+#define PARTICLEDYNAMICS_3D_H
+
+#include "core/superLattice3D.h"
+#include "functors/analytical/indicator/smoothIndicatorBaseF3D.h"
+
+namespace olb {
+
+template
+class ParticleDynamics3D {
+private:
+ SuperLattice3D& _sLattice;
+ UnitConverter const& _converter;
+ SuperGeometry3D& _superGeometry;
+ std::shared_ptr > _indicatorF;
+ std::vector* > _vectorOfIndicator;
+ T _lengthX;
+ T _lengthY;
+ T _lengthZ;
+ Vector _accExt;
+ bool _escapeFromDomain=false;
+ bool _oldWallCollision=true;
+
+public:
+ ParticleDynamics3D(SuperLattice3D& sLattice,
+ UnitConverter const& converter,
+ SuperGeometry3D& superGeometry,
+ T lengthX, T lengthY, T lengthZ,
+ Vector accExt = Vector (0.,0.,0.))
+ : _sLattice(sLattice),
+ _converter(converter),
+ _superGeometry(superGeometry),
+ _lengthX(lengthX),
+ _lengthY(lengthY),
+ _lengthZ(lengthZ),
+ _accExt(accExt)
+ {}
+
+ ParticleDynamics3D(SuperLattice3D& sLattice,
+ UnitConverter const& converter,
+ SuperGeometry3D& superGeometry,
+ std::shared_ptr > indicatorF,
+ bool escapeFromDomain, bool oldWallCollision=false,
+ Vector accExt = Vector (0.,0.,0.))
+ : _sLattice(sLattice),
+ _converter(converter),
+ _superGeometry(superGeometry),
+ _indicatorF(indicatorF),
+ _accExt(accExt),
+ _escapeFromDomain(escapeFromDomain),
+ _oldWallCollision(oldWallCollision)
+ {}
+
+ void addCuboid(Vector< T, 3> center, T xLength, T yLength, T zLength, T mass, T epsilon, Vector< T, 3 > theta, Vector vel = Vector (0.,0.,0.));
+ void addSphere(Vector< T, 3> center, T radius, T epsilon, T density, Vector vel = Vector (0.,0.,0.));
+ void addParticle(SmoothIndicatorF3D& indicator);
+
+ void computeBoundaryForce(std::vector* >& indicator);
+
+ void addWallColl(SmoothIndicatorF3D& indicator, T delta);
+
+ void verletIntegration(SmoothIndicatorF3D& indicator);
+
+ void updateParticleDynamics(std::string name, SmoothIndicatorF3D& indicator);
+
+ void checkAndRemoveEscaped();
+
+ void addParticleField(SmoothIndicatorF3D& indicator);
+
+ void simulateTimestep(std::string name);
+
+ void print();
+
+ void load(std::string filename, T epsilon);
+
+ void save(std::string filename);
+
+ void eulerIntegration(SmoothIndicatorF3D& indicator);
+
+};
+
+}
+
+#endif
diff --git a/src/particles/hlbm/hlbmDynamics3D.hh b/src/particles/hlbm/hlbmDynamics3D.hh
new file mode 100644
index 0000000..bc86b90
--- /dev/null
+++ b/src/particles/hlbm/hlbmDynamics3D.hh
@@ -0,0 +1,317 @@
+/* DESCRIPTOR Boltzmann sample, written in C++, using the OpenLB
+ * library
+ *
+ * Copyright (C) 2006-2016 Thomas Henn, Fabian Klemens, Robin Trunk, Davide Dapelo
+ * E-mail contact: info@openlb.net
+ * The most recent release of OpenLB can be downloaded at
+ *
+ *
+ * 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.
+ */
+
+
+#ifndef PARTICLEDYNAMICS_3D_HH
+#define PARTICLEDYNAMICS_3D_HH
+
+#include "hlbmDynamics3D.h"
+
+namespace olb {
+
+template
+void ParticleDynamics3D::addCuboid(Vector< T, 3> center, T xLength, T yLength, T zLength, T density, T epsilon, Vector< T, 3 > theta, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorCuboid3D(center, xLength, yLength, zLength, density, epsilon, theta, vel) );
+}
+
+template
+void ParticleDynamics3D::addSphere(Vector< T, 3> center, T radius, T epsilon, T density, Vector vel)
+{
+ _vectorOfIndicator.push_back (
+ new SmoothIndicatorSphere3D(center, radius, density, epsilon, vel) );
+}
+
+template
+void ParticleDynamics3D::addParticle(SmoothIndicatorF3D& indicator)
+{
+ _vectorOfIndicator.push_back(&indicator);
+}
+
+template
+void ParticleDynamics3D::computeBoundaryForce(std::vector* >& indicator)
+{
+ SuperLatticePorousMomentumLossForce3D force(_sLattice, _superGeometry, indicator, _converter);
+ T sumF[force.getTargetDim()];
+ for (int i=0; i* >::size_type iInd=0; iInd!=indicator.size(); iInd++) {
+ /// get particle acceleration through boundary force and gravity (and buoyancy)
+ Vector acceleration2;
+ Vector force;
+ Vector alpha2;
+ force[0] = sumF[0+7*iInd];
+ force[1] = sumF[1+7*iInd];
+ force[2] = sumF[2+7*iInd];
+ acceleration2[0] = sumF[0+7*iInd] / indicator[iInd]->getMass() + _accExt[0];
+ acceleration2[1] = sumF[1+7*iInd] / indicator[iInd]->getMass() + _accExt[1];
+ acceleration2[2] = sumF[2+7*iInd] / indicator[iInd]->getMass() + _accExt[2];
+ alpha2[0] = sumF[3+7*iInd] / indicator[iInd]->getMofi()[0];
+ alpha2[1] = sumF[4+7*iInd] / indicator[iInd]->getMofi()[1];
+ alpha2[2] = sumF[5+7*iInd] / indicator[iInd]->getMofi()[2];
+ indicator[iInd]->setForce( force );
+ indicator[iInd]->setAcc2( acceleration2 );
+ indicator[iInd]->setAlpha2( alpha2 );
+ }
+}
+
+template
+void ParticleDynamics3D::addWallColl(SmoothIndicatorF3D& indicator, T delta)
+{
+ T w1 = delta * .5;
+ T w = delta * .5;
+
+ T rad = indicator.getRadius();
+ T massInv = 1. / indicator.getMass();
+
+ std::vector dx(3, T());
+ dx[0] = _lengthX - _converter.getPhysDeltaX() - indicator.getPos()[0];
+ dx[1] = _lengthY - _converter.getPhysDeltaX() - indicator.getPos()[1];
+ dx[2] = indicator.getPos()[2] - _converter.getPhysDeltaX();
+
+ for (int i = 0; i < 3; i++) {
+ if (dx[i] <= rad) {
+ indicator.getAcc2()[i] += massInv * -dx[i] * (rad - dx[i]) / w1;
+ }
+ if (indicator.getPos()[i] <= rad) {
+ indicator.getAcc2()[i] += massInv * indicator.getPos()[i] * (rad - indicator.getPos()[i]) / w1;
+ }
+ if (dx[i] > rad && dx[i] <= rad + delta) {
+ indicator.getAcc2()[i] += massInv * -dx[i] * std::pow((rad + delta - dx[i]), 2) / w;
+ }
+ if (indicator.getPos()[i] > rad && indicator.getPos()[i] <= rad + delta) {
+ indicator.getAcc2()[i] += massInv * indicator.getPos()[i] * std::pow((rad + delta - indicator.getPos()[i]), 2) / w;
+ }
+ }
+}
+
+template
+void ParticleDynamics3D::verletIntegration(SmoothIndicatorF3D& indicator)
+{
+ T time = _converter.getConversionFactorTime();
+ T time2 = time*time;
+
+ Vector position, velocity, theta, omega, alpha;
+ Vector rotationMatrix;
+ for (int i=0; i<3; i++) {
+ position[i] = indicator.getPos()[i] + indicator.getVel()[i] * time + (0.5 * indicator.getAcc()[i] * time2);
+ T avgAcc = (indicator.getAcc()[i] + indicator.getAcc2()[i]) * 0.5;
+ velocity[i] = indicator.getVel()[i] + avgAcc * time;
+ theta[i] = std::fmod( indicator.getTheta()[i] + indicator.getOmega()[i] * time + (0.5 * indicator.getAlpha()[i] * time2), 2.*M_PI );
+ T avgAlpha = (indicator.getAlpha()[i] + indicator.getAlpha2()[i]) * 0.5;
+ omega[i] = indicator.getOmega()[i] + avgAlpha * time;
+ }
+ indicator.setPos( position );
+ indicator.setVel( velocity );
+ indicator.setAcc( indicator.getAcc2() );
+ indicator.setTheta( theta );
+ indicator.setOmega( omega );
+ indicator.setAlpha( indicator.getAlpha2() );
+
+ T cos0 = std::cos(indicator.getTheta()[0]);
+ T cos1 = std::cos(indicator.getTheta()[1]);
+ T cos2 = std::cos(indicator.getTheta()[2]);
+ T sin0 = std::sin(indicator.getTheta()[0]);
+ T sin1 = std::sin(indicator.getTheta()[1]);
+ T sin2 = std::sin(indicator.getTheta()[2]);
+
+ rotationMatrix[0] = cos1 * cos2;
+ rotationMatrix[1] = sin0*sin1*cos2 - cos0*sin2;
+ rotationMatrix[2] = cos0*sin1*cos2 + sin0*sin2;
+ rotationMatrix[3] = cos1*sin2;
+ rotationMatrix[4] = sin0*sin1*sin2 + cos0*cos2;
+ rotationMatrix[5] = cos0*sin1*sin2 - sin0*cos2;
+ rotationMatrix[6] = -sin1;
+ rotationMatrix[7] = sin0*cos1;
+ rotationMatrix[8] = cos0*cos1;
+ indicator.setRotationMatrix( rotationMatrix );
+}
+
+template
+void ParticleDynamics3D::updateParticleDynamics(std::string name, SmoothIndicatorF3D& indicator)
+{
+ if (name == "euler") {
+ this->eulerIntegration(indicator);
+ } else if (name == "verlet") {
+ this->verletIntegration(indicator);
+ } else {
+ std::cout << "ERROR: no valid integration...use 'euler' or 'verlet'"
+ << std::endl;
+ }
+}
+
+template
+void ParticleDynamics3D::checkAndRemoveEscaped()
+{
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ auto internal = std::make_shared > (*_indicatorF.get(), -_converter.getConversionFactorLength()+(**i).getEpsilon() );
+ IndicMinus3D layer(_indicatorF, internal);
+ SuperLatticeIndicatorSmoothIndicatorIntersection3D intersection(
+ _sLattice, _superGeometry, layer, **i );
+ int input[1];
+ T output[1] = {0.};
+ intersection(output, input);
+ if (output[0] == 1) {
+ _vectorOfIndicator.erase(i);
+ if (i==_vectorOfIndicator.end() ) {
+ break;
+ }
+ }
+ }
+}
+
+template
+void ParticleDynamics3D::addParticleField(SmoothIndicatorF3D& indicator)
+{
+ /// Analytical3D functor for particle motion (trans+rot)
+ ParticleU3D velocity(indicator, _converter);
+ _sLattice.setExternalParticleField(_superGeometry, velocity, indicator);
+}
+
+template
+void ParticleDynamics3D::simulateTimestep(std::string name)
+{
+ // Remove particles from domain if _escapeFromDomain is toggled
+ if (_escapeFromDomain) checkAndRemoveEscaped();
+
+ // Compute force acting on particles boundary
+ computeBoundaryForce(_vectorOfIndicator);
+
+ // Update particle dynamics and porous particle field
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ updateParticleDynamics(name, **i);
+ addParticleField(**i);
+ }
+}
+
+
+template
+void ParticleDynamics3D::print()
+{
+ OstreamManager clout(std::cout, "ParticleDynamics3D");
+ clout << "Number of particles=" << _vectorOfIndicator.size() << std::endl;
+ int count = 0;
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ clout << "Particle " << ++count << " (" << (**i).name() << "):" << std::endl;
+ clout << " |Circum radius(m)= " << setw(13) << (**i).getCircumRadius() << std::endl;
+ clout << " |Mass(kg)= " << setw(13) << (**i).getMass() << std::endl;
+ clout << " |Position(m)= (" << setw(13) << (**i).getPos()[0] << ", " << setw(13) << (**i).getPos()[1] << ", " << setw(13) << (**i).getPos()[2] << ")" << std::endl;
+ clout << " |Angle(°)= (" << setw(13) << (**i).getTheta()[0]*(180/M_PI) << ", " << setw(13) << (**i).getTheta()[1]*(180/M_PI) << ", " << setw(13) << (**i).getTheta()[2]*(180/M_PI) << ")" << std::endl;
+ clout << " |Velocity(m/s)= (" << setw(13) << (**i).getVel()[0] << ", " << setw(13) << (**i).getVel()[1] << ", " << setw(13) << (**i).getVel()[2] << ")" << std::endl;
+ clout << " |Ang. velocity(°/s)= (" << setw(13) << (**i).getOmega()[0]*(180/M_PI) << ", " << setw(13) << (**i).getOmega()[1]*(180/M_PI) << ", " << setw(13) << (**i).getOmega()[2]*(180/M_PI) << ")" << std::endl;
+ clout << " |Hydro. Force(N)= (" << setw(13) << (**i).getForce()[0] << ", " << setw(13) << (**i).getForce()[1] << ", " << setw(13) << (**i).getForce()[2] << ")" << std::endl;
+ clout << " |Acceleration(m/s^2)= (" << setw(13) << (**i).getAcc()[0] << ", " << setw(13) << (**i).getAcc()[1] << ", " << setw(13) << (**i).getAcc()[2] << ")" << std::endl;
+ clout << " |Ang. acc.(°/s^2)= (" << setw(13) << (**i).getAlpha()[0]*(180/M_PI) << ", " << setw(13) << (**i).getAlpha()[1]*(180/M_PI) << ", " << setw(13) << (**i).getAlpha()[1]*(180/M_PI) << ")" << std::endl;
+ }
+}
+
+// WORKS ONLY FOR SPHERES FOR NOW!!
+template
+void ParticleDynamics3D::load(std::string filename, T epsilon)
+{
+ std::ifstream iout( (singleton::directories().getLogOutDir() + filename).c_str() );
+
+ while (iout) {
+ std::string name = "";
+ iout >> name;
+
+ if (name == "sphere") {
+ T radius = T(), mass = T();
+ Vector center = {T(), T(), T()};
+ Vector vel = {T(), T(), T()};
+ iout >> center[0] >> center[1] >> center[2]
+ >> radius >> mass
+ >> vel[0] >> vel[1] >> vel[2];
+ addSphere(center, radius, epsilon, mass, vel);
+ }
+ }
+
+ iout.close();
+}
+
+// WORKS ONLY FOR SPHERES FOR NOW!!
+template
+void ParticleDynamics3D::save(std::string filename)
+{
+ std::ofstream fout( (singleton::directories().getLogOutDir() + filename).c_str() );
+
+ for (auto i=_vectorOfIndicator.begin(); i!=_vectorOfIndicator.end(); i++) {
+ if ((**i).name() == "sphere") {
+ fout << (**i).name() << " "
+ << (**i).getPos()[0] << " " << (**i).getPos()[1] << " " << (**i).getPos()[2] << " "
+ << (**i).getCircumRadius() << " " << (**i).getMass() << " "
+ << (**i).getVel()[0] << " " << (**i).getVel()[1] << " " << (**i).getVel()[2]
+ << std::endl;
+ }
+ }
+ fout.close();
+}
+
+template
+void ParticleDynamics3D::eulerIntegration(SmoothIndicatorF3D& indicator) {
+ T time = _converter.getConversionFactorTime();
+
+ Vector position, velocity, theta, omega, alpha;
+ Vector rotationMatrix;
+ for (int i=0; i<3; i++) {
+ velocity[i] = indicator.getVel()[i] + indicator.getAcc2()[i] * time;
+ position[i] = indicator.getPos()[i] + indicator.getVel()[i] * time;
+ omega[i] = indicator.getOmega()[i] + indicator.getAlpha2()[i] * time;
+ theta[i] = indicator.getTheta()[i] + indicator.getOmega()[i] * time;
+ }
+
+ indicator.setPos( position );
+ indicator.setVel( velocity );
+ indicator.setAcc( indicator.getAcc2() );
+ indicator.setTheta( theta );
+ indicator.setOmega( omega );
+ indicator.setAlpha( indicator.getAlpha2() );
+
+ T cos0 = std::cos(indicator.getTheta()[0]);
+ T cos1 = std::cos(indicator.getTheta()[1]);
+ T cos2 = std::cos(indicator.getTheta()[2]);
+ T sin0 = std::sin(indicator.getTheta()[0]);
+ T sin1 = std::sin(indicator.getTheta()[1]);
+ T sin2 = std::sin(indicator.getTheta()[2]);
+
+ rotationMatrix[0] = cos1 * cos2;
+ rotationMatrix[1] = sin0*sin1*cos2 - cos0*sin2;
+ rotationMatrix[2] = cos0*sin1*cos2 + sin0*sin2;
+ rotationMatrix[3] = cos1*sin2;
+ rotationMatrix[4] = sin0*sin1*sin2 + cos0*cos2;
+ rotationMatrix[5] = cos0*sin1*sin2 - sin0*cos2;
+ rotationMatrix[6] = -sin1;
+ rotationMatrix[7] = sin0*cos1;
+ rotationMatrix[8] = cos0*cos1;
+ indicator.setRotationMatrix( rotationMatrix );
+}
+
+}
+
+#endif
--
cgit v1.2.3