Commit 4d96d96a authored by Tobias Leemann's avatar Tobias Leemann Committed by Sebastian Eibl
Browse files

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

Added DetectAndStoreContacts-Kernel and its testcase, performed other adaptions in preperation of HCSITS-Solver integreation
parent 0fe3bd47
......@@ -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 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 ContactStorage.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//! \author Tobias Leemann <tobias.leemann@fau.de>
//
//======================================================================================================================
/**
* Storage for detected contacts which can be used to perform actions
* for all contacts, e.g. resolving, relaxation schemes...
* The InsertIntoContactStorage-Kernel can be used to insert a contact.
*/
//======================================================================================================================
//
// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
//
//======================================================================================================================
#pragma once
#include <mesa_pd/data/DataTypes.h>
{%- for include in includes %}
#include <{{include}}>
{%- endfor %}
#include <mesa_pd/data/STLOverloads.h>
#include <core/Abort.h>
#include <core/debug/Debug.h>
#include <core/math/AABB.h>
#include <core/OpenMP.h>
#include <core/STLIO.h>
#include <core/UniqueID.h>
#include <atomic>
#include <limits>
#include <map>
#include <type_traits>
#include <unordered_map>
#include <vector>
namespace walberla {
namespace mesa_pd {
namespace data {
class ContactStorage;
class ContactStorage
{
public:
class Contact
{
public:
constexpr Contact(ContactStorage& storage, const size_t i) : storage_(storage), i_(i) {}
constexpr Contact(const Contact&) = default;
constexpr Contact(Contact&&) = default;
Contact& operator=(const Contact& rhs);
Contact& operator=(Contact&& rhs);
Contact* operator->(){return this;}
{% for prop in properties %}
const {{prop.type}}& get{{prop.name | capFirst}}() const {return storage_.get{{prop.name | capFirst}}(i_);}
{{prop.type}}& get{{prop.name | capFirst}}Ref() {return storage_.get{{prop.name | capFirst}}Ref(i_);}
void set{{prop.name | capFirst}}(const {{prop.type}}& v) { storage_.set{{prop.name | capFirst}}(i_, v);}
{% endfor %}
size_t getIdx() const {return i_;}
public:
ContactStorage& storage_;
const size_t i_;
};
class iterator
{
public:
using iterator_category = std::random_access_iterator_tag;
using value_type = Contact;
using pointer = Contact*;
using reference = Contact&;
using difference_type = std::ptrdiff_t;
explicit iterator(ContactStorage* storage, const size_t i) : storage_(storage), i_(i) {}
iterator(const iterator& it) = default;
iterator(iterator&& it) = default;
iterator& operator=(const iterator&) = default;
iterator& operator=(iterator&&) = default;
Contact operator*(){return Contact{*storage_, i_};}
Contact operator->(){return Contact{*storage_, i_};}
iterator& operator++(){ ++i_; return *this; }
iterator operator++(int){ iterator tmp(*this); ++(*this); return tmp; }
iterator& operator--(){ --i_; return *this; }
iterator operator--(int){ iterator tmp(*this); --(*this); return tmp; }
iterator& operator+=(const size_t n){ i_+=n; return *this; }
iterator& operator-=(const size_t n){ i_-=n; return *this; }
friend iterator operator+(const iterator& it, const size_t n);
friend iterator operator+(const size_t n, const iterator& it);
friend iterator operator-(const iterator& it, const size_t n);
friend difference_type operator-(const iterator& lhs, const iterator& rhs);
friend bool operator==(const iterator& lhs, const iterator& rhs);
friend bool operator!=(const iterator& lhs, const iterator& rhs);
friend bool operator<(const iterator& lhs, const iterator& rhs);
friend bool operator>(const iterator& lhs, const iterator& rhs);
friend bool operator<=(const iterator& lhs, const iterator& rhs);
friend bool operator>=(const iterator& lhs, const iterator& rhs);
friend void swap(iterator& lhs, iterator& rhs);
size_t getIdx() const {return i_;}
private:
ContactStorage* storage_;
size_t i_;
};
explicit ContactStorage(const size_t size);
iterator begin() { return iterator(this, 0); }
iterator end() { return iterator(this, size()); }
iterator operator[](const size_t n) { return iterator(this, n); }
{% for prop in properties %}
const {{prop.type}}& get{{prop.name | capFirst}}(const size_t idx) const {return {{prop.name}}_[idx];}
{{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t idx) {return {{prop.name}}_[idx];}
void set{{prop.name | capFirst}}(const size_t idx, const {{prop.type}}& v) { {{prop.name}}_[idx] = v; }
{% endfor %}
/**
* @brief creates a new Contact and returns an iterator pointing to it
*
* \attention Use this function only if you know what you are doing!
* Messing with the uid might break the simulation!
* If you are unsure use create(bool) instead.
* @param uid unique id of the Contact to be created
* @return iterator to the newly created Contact
*/
inline iterator create(const id_t& uid);
inline iterator create();
inline iterator erase(iterator& it);
/// Finds the entry corresponding to \p uid.
/// \return iterator to the object or end iterator
inline iterator find(const id_t& uid);
inline void reserve(const size_t size);
inline void clear();
inline size_t size() const;
/**
* Calls the provided functor \p func for all Contacts selected by the selector.
*
* Additional arguments can be provided.
* Call syntax for the provided functor
* \code
* func( *this, i, std::forward<Args>(args)... );
* \endcode
* \param openmp enables/disables OpenMP parallelization of the kernel call
*/
template <typename Selector, typename Accessor, typename Func, typename... Args>
inline void forEachContact(const bool openmp, const Selector& selector,
Accessor& acForPS, Func&& func, Args&&... args);
template <typename Selector, typename Accessor, typename Func, typename... Args>
inline void forEachContact(const bool openmp, const Selector& selector,
Accessor& acForPS, Func&& func, Args&&... args) const;
private:
{%- for prop in properties %}
std::vector<{{prop.type}}> {{prop.name}}_ {};
{%- endfor %}
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!");
};
using Contact = ContactStorage::Contact;
inline
ContactStorage::Contact& ContactStorage::Contact::operator=(const ContactStorage::Contact& rhs)
{
{%- for prop in properties %}
get{{prop.name | capFirst}}Ref() = rhs.get{{prop.name | capFirst}}();
{%- endfor %}
return *this;
}
inline
ContactStorage::Contact& ContactStorage::Contact::operator=(ContactStorage::Contact&& rhs)
{
{%- for prop in properties %}
get{{prop.name | capFirst}}Ref() = std::move(rhs.get{{prop.name | capFirst}}Ref());
{%- endfor %}
return *this;
}
inline
std::ostream& operator<<( std::ostream& os, const ContactStorage::Contact& p )
{
os << "========== {{StorageType | upper}} ==========" << "\n" <<
"idx : " << p.getIdx() << "\n" <<
{%- for prop in properties %}
"{{'%-20s'|format(prop.name)}}: " << p.get{{prop.name | capFirst}}() << "\n" <<
{%- endfor %}
"================================" << std::endl;
return os;
}
inline
ContactStorage::iterator operator+(const ContactStorage::iterator& it, const size_t n)
{
return ContactStorage::iterator(it.storage_, it.i_+n);
}
inline
ContactStorage::iterator operator+(const size_t n, const ContactStorage::iterator& it)
{
return it + n;
}
inline
ContactStorage::iterator operator-(const ContactStorage::iterator& it, const size_t n)
{
return ContactStorage::iterator(it.storage_, it.i_-n);
}
inline
ContactStorage::iterator::difference_type operator-(const ContactStorage::iterator& lhs, const ContactStorage::iterator& rhs)
{
return int64_c(lhs.i_) - int64_c(rhs.i_);
}