Commit b88a8ec8 authored by Sebastian Eibl's avatar Sebastian Eibl

Merge branch 'DetectAndStoreContacts' into 'master'

Added DetectAndStoreContacts-Kernel and its testcase, performed other...

See merge request walberla/walberla!211
parents 0fe3bd47 4d96d96a
Pipeline #16245 failed with stages
in 350 minutes and 16 seconds
......@@ -74,19 +74,15 @@ class ParticleAccessorWithShape : public data::ParticleAccessor
{
public:
ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
: ParticleAccessor(ps)
, ss_(ss)
: ParticleAccessor(ps)
, ss_(ss)
{}
const walberla::real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
walberla::real_t& getInvMassRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
void setInvMass(const size_t p_idx, const walberla::real_t& v) { ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass() = v;}
const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
auto& getInvInertiaBFRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
void setInvInertiaBF(const size_t p_idx, const Mat3& v) { ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF() = v;}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
private:
std::shared_ptr<data::ShapeStorage> ss_;
};
......
......@@ -99,19 +99,15 @@ class ParticleAccessorWithShape : public data::ParticleAccessor
{
public:
ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
: ParticleAccessor(ps)
, ss_(ss)
: ParticleAccessor(ps)
, ss_(ss)
{}
const walberla::real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
walberla::real_t& getInvMassRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
void setInvMass(const size_t p_idx, const walberla::real_t& v) { ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass() = v;}
const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
auto& getInvInertiaBFRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
void setInvInertiaBF(const size_t p_idx, const Mat3& v) { ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF() = v;}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
private:
std::shared_ptr<data::ShapeStorage> ss_;
};
......
......@@ -33,6 +33,7 @@ if __name__ == '__main__':
ch = data.ContactHistory()
lc = data.LinkedCells()
ss = data.ShapeStorage(ps, shapes)
cs = data.ContactStorage()
ps.addProperty("position", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
ps.addProperty("linearVelocity", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
......@@ -54,14 +55,37 @@ if __name__ == '__main__':
ps.addProperty("oldContactHistory", "std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>", defValue="", syncMode="ALWAYS")
ps.addProperty("newContactHistory", "std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>", defValue="", syncMode="NEVER")
ps.addProperty("temperature", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS")
ps.addProperty("heatFlux", "walberla::real_t", defValue="real_t(0)", syncMode="NEVER")
ps.addProperty("temperature", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS")
ps.addProperty("heatFlux", "walberla::real_t", defValue="real_t(0)", syncMode="NEVER")
# Properties for HCSITS
ps.addProperty("dv", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.addProperty("dw", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ch.addProperty("tangentialSpringDisplacement", "walberla::mesa_pd::Vec3", defValue="real_t(0)")
ch.addProperty("isSticking", "bool", defValue="false")
ch.addProperty("impactVelocityMagnitude", "real_t", defValue="real_t(0)")
cs.addProperty("id1", "walberla::id_t", defValue = "walberla::id_t(-1)", syncMode="NEVER")
cs.addProperty("id2", "walberla::id_t", defValue = "walberla::id_t(-1)", syncMode="NEVER")
cs.addProperty("distance", "real_t", defValue = "real_t(1)", syncMode="NEVER")
cs.addProperty("normal", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("position", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("t", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("o", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("r1", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("r2", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("mu", "real_t", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("p", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("diag_nto", "walberla::mesa_pd::Mat3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("diag_nto_inv", "walberla::mesa_pd::Mat3", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("diag_to_inv", "walberla::mesa_pd::Mat2", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("diag_n_inv", "real_t", defValue = "real_t(0)", syncMode="NEVER")
cs.addProperty("p", "walberla::mesa_pd::Vec3", defValue = "real_t(0)", syncMode="NEVER")
kernels = []
kernels.append( kernel.DetectAndStoreContacts() )
kernels.append( kernel.DoubleCast(shapes) )
kernels.append( kernel.ExplicitEuler() )
kernels.append( kernel.ExplicitEulerWithShape() )
......@@ -76,6 +100,7 @@ if __name__ == '__main__':
kernels.append( kernel.VelocityVerlet() )
kernels.append( kernel.VelocityVerletWithShape() )
ac = Accessor()
for k in kernels:
ac.mergeRequirements(k.getRequirements())
......@@ -93,6 +118,7 @@ if __name__ == '__main__':
ch.generate(args.path + "/src/mesa_pd/")
lc.generate(args.path + "/src/mesa_pd/")
ss.generate(args.path + "/src/mesa_pd/")
cs.generate(args.path + "/src/mesa_pd/")
for k in kernels:
k.generate(args.path + "/src/mesa_pd/")
......
# -*- coding: utf-8 -*-
import numpy as np
from ..Container import Container
from ..utility import generateFile
class ContactStorage(Container):
def __init__(self):
super().__init__()
self.addProperty("uid", "walberla::id_t", defValue = "walberla::id_t(-1)", syncMode="NEVER")
def generate(self, path):
self.unrollDimension()
print("="*90)
print("Creating ContactStorage Datastructure:")
print("")
print("{0: <20}{1: <30}{2: <20}{3: <10}".format("Type", "Name", "Def. Value", "SyncMode"))
print("="*90)
for prop in self.properties:
print("{0: <20.19}{1: <30.29}{2: <20.19}{3: <10.9}".format(prop.type, prop.name, prop.defValue, prop.syncMode))
print("="*90)
context = dict()
context["includes"] = self.includes
context["properties"] = self.properties
generateFile(path, 'data/ContactStorage.templ.h', context, filename='data/ContactStorage.h')
generateFile(path, 'data/ContactAccessor.templ.h', context, filename='data/ContactAccessor.h')
# -*- coding: utf-8 -*-
from .ContactHistory import ContactHistory
from .ContactStorage import ContactStorage
from .LinkedCells import LinkedCells
from .ParticleStorage import ParticleStorage
from .ShapeStorage import ShapeStorage
__all__ = ['ContactHistory',
'ContactStorage',
'GeometryStorage',
'LinkedCells',
'ParticleStorage']
# -*- coding: utf-8 -*-
from mesa_pd.accessor import Accessor
from mesa_pd.utility import generateFile
class DetectAndStoreContacts:
def __init__(self):
self.accessor = Accessor()
self.accessor.require("uid", "walberla::id_t", access="g")
self.accessor.require("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g")
self.accessor.require("position", "walberla::mesa_pd::Vec3", access="g")
self.accessor.require("rotation", "walberla::mesa_pd::Rot3", access="g")
self.accessor.require("shape", "BaseShape*", access="g")
def getRequirements(self):
return self.accessor
def generate(self, path):
context = dict()
context["interface"] = self.accessor.properties
generateFile(path, 'kernel/DetectAndStoreContacts.templ.h', context)
# -*- coding: utf-8 -*-
from .DetectAndStoreContacts import DetectAndStoreContacts
from .DoubleCast import DoubleCast
from .ExplicitEuler import ExplicitEuler
from .ExplicitEulerWithShape import ExplicitEulerWithShape
......@@ -15,6 +15,7 @@ from .VelocityVerlet import VelocityVerlet
from .VelocityVerletWithShape import VelocityVerletWithShape
__all__ = ['DoubleCast',
'DetectAndStoreContacts',
'ExplicitEuler',
'ExplicitEulerWithShape',
'ForceLJ',
......
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla 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.
//
// waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file ContactAccessor.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
//======================================================================================================================
//
// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
//
//======================================================================================================================
#pragma once
#include <mesa_pd/data/IAccessor.h>
#include <mesa_pd/data/ContactStorage.h>
#include <core/UniqueID.h>
#include <limits>
namespace walberla {
namespace mesa_pd {
namespace data {
/**
* @brief Basic ContactAccessor for the ContactStorage
*
* Provides get, set and getRef for all members of the ContactStorage.
* Can be used as a basis class for a more advanced ContactAccessor.
*/
class ContactAccessor : public IAccessor
{
public:
ContactAccessor(const std::shared_ptr<data::ContactStorage>& ps) : ps_(ps) {}
virtual ~ContactAccessor() = default;
{%- for prop in properties %}
const {{prop.type}}& get{{prop.name | capFirst}}(const size_t p_idx) const {return ps_->get{{prop.name | capFirst}}(p_idx);}
{{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t p_idx) {return ps_->get{{prop.name | capFirst}}Ref(p_idx);}
void set{{prop.name | capFirst}}(const size_t p_idx, const {{prop.type}}& v) { ps_->set{{prop.name | capFirst}}(p_idx, v);}
{% endfor %}
id_t getInvalidUid() const {return UniqueID<data::Contact>::invalidID();}
size_t getInvalidIdx() const {return std::numeric_limits<size_t>::max();}
/**
* @brief Returns the index of Contact specified by uid.
* @param uid unique id of the Contact to be looked up
* @return the index of the Contact or std::numeric_limits<size_t>::max() if the Contact is not found
*/
size_t uidToIdx(const id_t& uid) const {auto it = ps_->find(uid); return it != ps_->end() ? it.getIdx() : std::numeric_limits<size_t>::max();}
size_t size() const { return ps_->size(); }
inline size_t create(const id_t& uid);
inline size_t erase(const size_t& idx);
inline size_t find(const id_t& uid);
protected:
std::shared_ptr<data::ContactStorage> ps_;
};
inline size_t ContactAccessor::create(const id_t& uid)
{
auto it = ps_->create(uid);
return it.getIdx();
}
inline size_t ContactAccessor::erase(const size_t& idx)
{
data::ContactStorage::iterator it(ps_.get(), idx);
it = ps_->erase(it);
return it.getIdx();
}
inline size_t ContactAccessor::find(const id_t& uid)
{
auto it = ps_->find(uid);
return it.getIdx();
}
/**
* @brief Basic ContactAccessor which emulates a single Contact in a ContactStorage
* without actually needing a ContactStorage. This class is used mainly for testing purposes.
*
* Provides get, set and getRef.
*/
class SingleContactAccessor : public IAccessor
{
public:
virtual ~SingleContactAccessor() = default;
{%- for prop in properties %}
const {{prop.type}}& get{{prop.name | capFirst}}(const size_t /*p_idx*/) const {return {{prop.name}}_;}
void set{{prop.name | capFirst}}(const size_t /*p_idx*/, const {{prop.type}}& v) { {{prop.name}}_ = v;}
{{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t /*p_idx*/) {return {{prop.name}}_;}
{% endfor %}
id_t getInvalidUid() const {return UniqueID<data::Contact>::invalidID();}
size_t getInvalidIdx() const {return std::numeric_limits<size_t>::max();}
/**
* @brief Returns the index of Contact specified by uid.
* @param uid unique id of the Contact to be looked up
* @return the index of the Contact or std::numeric_limits<size_t>::max() if the Contact is not found
*/
size_t uidToIdx(const id_t& uid) const {return uid == uid_ ? 0 : std::numeric_limits<size_t>::max();}
size_t size() const { return 1; }
private:
{%- for prop in properties %}
{{prop.type}} {{prop.name}}_;
{%- endfor %}
};
} //namespace data
} //namespace mesa_pd
} //namespace walberla
This diff is collapsed.
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla 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.
//
// waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file DetectAndStoreContacts.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//! \author Tobias Leemann <tobias.leemann@fau.de>
//
//======================================================================================================================
//======================================================================================================================
//
// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
//
//======================================================================================================================
#pragma once
#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/mpi/ContactFilter.h>
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/ContactStorage.h>
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/data/Flags.h>
namespace walberla {
namespace mesa_pd {
namespace kernel {
/**
* Kernel which performes collision detection on a pair of two particles
* and inserts the contact (if existent) into the contact storage passed in the constructor.
* Call this kernel on each particle pair to perform contact detection and insert each contact in the contact
* storage.
* \ingroup mesa_pd_kernel
*/
class DetectAndStoreContacts
{
public:
explicit DetectAndStoreContacts(data::ContactStorage &cs) : cs_(cs) {}
/**
* Compute if particle 1 and particle 2 collide. If so, insert the contact into the contact storage.
* \param idx1 The index of particle 1
* \param idx2 The index of particle 2
* \param ac The accessor used to access the values of the particles.
* \param domain The domain of the block (used to decide if the contact has to be treated by this process)
* \param acd The collision detection to be used. Default parameter: AnalyticContactDetection()
*/
template <typename Accessor>
void operator()(size_t idx1, size_t idx2, Accessor &ac, const domain::IDomain& domain, collision_detection::AnalyticContactDetection acd = collision_detection::AnalyticContactDetection());
private:
data::ContactStorage& cs_;
};
template <typename Accessor>
inline void DetectAndStoreContacts::operator()(size_t idx1, size_t idx2, Accessor &ac, const domain::IDomain& domain, collision_detection::AnalyticContactDetection acd)
{
using namespace data::particle_flags;
kernel::DoubleCast double_cast;
mpi::ContactFilter contact_filter;
if (double_cast(idx1, idx2, ac, acd, ac ))
{
if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), domain))
{
auto c = cs_.create();
c->setId1(acd.getIdx1());
c->setId2(acd.getIdx2());
c->setDistance(acd.getPenetrationDepth());
c->setNormal(acd.getContactNormal());
c->setPosition(acd.getContactPoint());
}
}
}
} //namespace kernel
} //namespace mesa_pd
} //namespace walberla
This diff is collapsed.
This diff is collapsed.
......@@ -21,9 +21,11 @@
#pragma once
#include <core/DataTypes.h>
#include <core/math/Matrix2.h>
#include <core/math/Matrix3.h>
#include <core/math/Quaternion.h>
#include <core/math/Rot3.h>
#include <core/math/Vector2.h>
#include <core/math/Vector3.h>
#include <mesa_pd/data/Flags.h>
......@@ -32,9 +34,10 @@ namespace walberla {
namespace mesa_pd {
using Mat3 = math::Matrix3<real_t>;
using Mat2 = math::Matrix2<real_t>;
using Rot3 = math::Rot3<real_t>;
using Quat = math::Quaternion<real_t>;
using Vec3 = math::Vector3<real_t>;
using Vec2 = math::Vector2<real_t>;
}
}
......@@ -132,6 +132,14 @@ public:
walberla::real_t& getHeatFluxRef(const size_t p_idx) {return ps_->getHeatFluxRef(p_idx);}
void setHeatFlux(const size_t p_idx, const walberla::real_t& v) { ps_->setHeatFlux(p_idx, v);}
const walberla::mesa_pd::Vec3& getDv(const size_t p_idx) const {return ps_->getDv(p_idx);}
walberla::mesa_pd::Vec3& getDvRef(const size_t p_idx) {return ps_->getDvRef(p_idx);}
void setDv(const size_t p_idx, const walberla::mesa_pd::Vec3& v) { ps_->setDv(p_idx, v);}
const walberla::mesa_pd::Vec3& getDw(const size_t p_idx) const {return ps_->getDw(p_idx);}
walberla::mesa_pd::Vec3& getDwRef(const size_t p_idx) {return ps_->getDwRef(p_idx);}
void setDw(const size_t p_idx, const walberla::mesa_pd::Vec3& v) { ps_->setDw(p_idx, v);}
id_t getInvalidUid() const {return UniqueID<data::Particle>::invalidID();}
size_t getInvalidIdx() const {return std::numeric_limits<size_t>::max();}
......@@ -261,6 +269,14 @@ public:
void setHeatFlux(const size_t /*p_idx*/, const walberla::real_t& v) { heatFlux_ = v;}
walberla::real_t& getHeatFluxRef(const size_t /*p_idx*/) {return heatFlux_;}
const walberla::mesa_pd::Vec3& getDv(const size_t /*p_idx*/) const {return dv_;}
void setDv(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { dv_ = v;}
walberla::mesa_pd::Vec3& getDvRef(const size_t /*p_idx*/) {return dv_;}
const walberla::mesa_pd::Vec3& getDw(const size_t /*p_idx*/) const {return dw_;}
void setDw(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { dw_ = v;}
walberla::mesa_pd::Vec3& getDwRef(const size_t /*p_idx*/) {return dw_;}
id_t getInvalidUid() const {return UniqueID<data::Particle>::invalidID();}
size_t getInvalidIdx() const {return std::numeric_limits<size_t>::max();}
......@@ -293,6 +309,8 @@ private:
std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> newContactHistory_;
walberla::real_t temperature_;
walberla::real_t heatFlux_;
walberla::mesa_pd::Vec3 dv_;
walberla::mesa_pd::Vec3 dw_;
};
} //namespace data
......
......@@ -152,6 +152,14 @@ public:
walberla::real_t& getHeatFluxRef() {return storage_.getHeatFluxRef(i_);}
void setHeatFlux(const walberla::real_t& v) { storage_.setHeatFlux(i_, v);}
const walberla::mesa_pd::Vec3& getDv() const {return storage_.getDv(i_);}
walberla::mesa_pd::Vec3& getDvRef() {return storage_.getDvRef(i_);}
void setDv(const walberla::mesa_pd::Vec3& v) { storage_.setDv(i_, v);}
const walberla::mesa_pd::Vec3& getDw() const {return storage_.getDw(i_);}
walberla::mesa_pd::Vec3& getDwRef() {return storage_.getDwRef(i_);}
void setDw(const walberla::mesa_pd::Vec3& v) { storage_.setDw(i_, v);}
size_t getIdx() const {return i_;}
public:
......@@ -296,6 +304,14 @@ public:
walberla::real_t& getHeatFluxRef(const size_t idx) {return heatFlux_[idx];}
void setHeatFlux(const size_t idx, const walberla::real_t& v) { heatFlux_[idx] = v; }
const walberla::mesa_pd::Vec3& getDv(const size_t idx) const {return dv_[idx];}
walberla::mesa_pd::Vec3& getDvRef(const size_t idx) {return dv_[idx];}
void setDv(const size_t idx, const walberla::mesa_pd::Vec3& v) { dv_[idx] = v; }
const walberla::mesa_pd::Vec3& getDw(const size_t idx) const {return dw_[idx];}
walberla::mesa_pd::Vec3& getDwRef(const size_t idx) {return dw_[idx];}
void setDw(const size_t idx, const walberla::mesa_pd::Vec3& v) { dw_[idx] = v; }
/**
* @brief creates a new particle and returns an iterator pointing to it
......@@ -406,6 +422,8 @@ public:
std::vector<std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>> newContactHistory_ {};
std::vector<walberla::real_t> temperature_ {};
std::vector<walberla::real_t> heatFlux_ {};
std::vector<walberla::mesa_pd::Vec3> dv_ {};
std::vector<walberla::mesa_pd::Vec3> dw_ {};
std::unordered_map<id_t, size_t> uidToIdx_;
static_assert(std::is_same<decltype(uid_)::value_type, id_t>::value,
"Property uid of type id_t is missing. This property is required!");
......@@ -436,6 +454,8 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(const ParticleSt
getNewContactHistoryRef() = rhs.getNewContactHistory();
getTemperatureRef() = rhs.getTemperature();
getHeatFluxRef() = rhs.getHeatFlux();
getDvRef() = rhs.getDv();
getDwRef() = rhs.getDw();
return *this;
}
......@@ -463,6 +483,8 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
getNewContactHistoryRef() = std::move(rhs.getNewContactHistoryRef());
getTemperatureRef() = std::move(rhs.getTemperatureRef());
getHeatFluxRef() = std::move(rhs.getHeatFluxRef());
getDvRef() = std::move(rhs.getDvRef());
getDwRef() = std::move(rhs.getDwRef());
return *this;
}
......@@ -492,6 +514,8 @@ std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
"newContactHistory : " << p.getNewContactHistory() << "\n" <<
"temperature : " << p.getTemperature() << "\n" <<
"heatFlux : " << p.getHeatFlux() << "\n" <<
"dv : " << p.getDv() << "\n" <<
"dw : " << p.getDw() << "\n" <<
"================================" << std::endl;
return os;
}
......@@ -590,6 +614,8 @@ inline ParticleStorage::iterator ParticleStorage::create(const id_t& uid)
newContactHistory_.emplace_back();
temperature_.emplace_back(real_t(0));
heatFlux_.emplace_back(real_t(0));
dv_.emplace_back(real_t(0));
dw_.emplace_back(real_t(0));
uid_.back() = uid;
uidToIdx_[uid] = uid_.size() - 1;
return iterator(this, size() - 1);
......@@ -643,6 +669,8 @@ inline ParticleStorage::iterator ParticleStorage::erase(iterator& it)
newContactHistory_.pop_back();
temperature_.pop_back();
heatFlux_.pop_back();
dv_.pop_back();
dw_.pop_back();
return it;
}
......@@ -683,6 +711,8 @@ inline void ParticleStorage::reserve(const size_t size)
newContactHistory_.reserve(size);
temperature_.reserve(size);
heatFlux_.reserve(size);
dv_.reserve(size);
dw_.reserve(size);
}
inline void ParticleStorage::clear()
......@@ -708,6 +738,8 @@ inline void ParticleStorage::clear()
newContactHistory_.clear();
temperature_.clear();
heatFlux_.clear();
dv_.clear();
dw_.clear();
uidToIdx_.clear();
}
......@@ -734,6 +766,8 @@ inline size_t ParticleStorage::size() const
//WALBERLA_ASSERT_EQUAL( uid_.size(), newContactHistory.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), temperature.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), heatFlux.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), dv.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), dw.size() );
return uid_.size();
}
template <typename Selector, typename Accessor, typename Func, typename... Args>
......@@ -1065,6 +1099,24 @@ public:
walberla::real_t& operator()(data::Particle&& p) const {return p.getHeatFluxRef();}
const walberla::real_t& operator()(const data::Particle& p) const {return p.getHeatFlux();}
};
///Predicate that selects a certain property from a Particle
class SelectParticleDv
{
public:
using return_type = walberla::mesa_pd::Vec3;
walberla::mesa_pd::Vec3& operator()(data::Particle& p) const {return p.getDvRef();}
walberla::mesa_pd::Vec3& operator()(data::Particle&& p) const {return p.getDvRef();}
const walberla::mesa_pd::Vec3& operator()(const data::Particle& p) const {return p.getDv();}
};
///Predicate that selects a certain property from a Particle
class SelectParticleDw
{
public:
using return_type = walberla::mesa_pd::Vec3;
walberla::mesa_pd::Vec3& operator()(data::Particle& p) const {return p.getDwRef();}
walberla::mesa_pd::Vec3& operator()(data::Particle&& p) const {return p.getDwRef();}
const walberla::mesa_pd::Vec3& operator()(const data::Particle& p) const {return p.getDw();}
};
} //namespace data
} //namespace mesa_pd
......
......@@ -43,17 +43,19 @@ public:
virtual real_t getVolume() const = 0;
real_t& getInvMass() {return invMass_;}
const real_t& getMass() const {return mass_;}
const real_t& getInvMass() const {return invMass_;}
Mat3& getInvInertiaBF() {return invInertiaBF_;}
const Mat3& getInertiaBF() const {return inertiaBF_;}
const Mat3& getInvInertiaBF() const {return invInertiaBF_;}
const int& getShapeType() const {return shapeType_;}
static const int INVALID_SHAPE = -1; ///< Unique *invalid* shape type identifier.\ingroup mesa_pd_shape
private:
protected:
real_t mass_ = real_t(0); ///< mass
real_t invMass_ = real_t(0); ///< inverse mass
Mat3 inertiaBF_ = Mat3(real_t(0)); ///< inertia matrix in the body frame
Mat3 invInertiaBF_ = Mat3(real_t(0)); ///< inverse inertia matrix in the body frame
int shapeType_ = INVALID_SHAPE; ///< \ingroup mesa_pd_shape
};
......
......@@ -56,8 +56,11 @@ private:
inline
void HalfSpace::updateMassAndInertia(const real_t /*density*/)
{
getInvMass() = real_t(0);
getInvInertiaBF() = Mat3(real_t(0));
mass_ = std::numeric_limits<real_t>::infinity();
invMass_ = real_t(0);
inertiaBF_ = Mat3(std::numeric_limits<real_t>::infinity());
invInertiaBF_ = Mat3(real_t(0));
}
} //namespace data
......
......@@ -51,8 +51,11 @@ void Sphere::updateMassAndInertia(const real_t density)
const real_t m = (real_c(4.0)/real_c(3.0) * math::M_PI) * getRadius() * getRadius() * getRadius() * density;
const Mat3 I = Mat3::makeDiagonalMatrix( real_c(0.4) * m * getRadius() * getRadius() );
getInvMass() = real_t(1.0) / m;
getInvInertiaBF() = I.getInverse();
mass_ = m;
invMass_ = real_t(1.0) / m;
inertiaBF_ = I;
invInertiaBF_ = I.getInverse();
}
} //namespace data
......
......@@ -30,11 +30,11 @@ class InfiniteDomain : public IDomain
{