Contains.h 5.56 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//======================================================================================================================
//
//  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/>.
//
16
//! \file
17
18
19
20
21
22
23
24
25
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================

#pragma once

#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/Flags.h>
#include <mesa_pd/data/IAccessor.h>
26
27
28
29
30
31

#include <mesa_pd/common/ParticleFunctions.h>

#include <mesa_pd/data/shape/Box.h>
#include <mesa_pd/data/shape/CylindricalBoundary.h>
#include <mesa_pd/data/shape/Ellipsoid.h>
32
33
34
35
36
37
38
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>

namespace walberla {
namespace mesa_pd {

/*
39
40
41
 * "contains point" functionality
 * can either be formulated in world frame coordinates (then the rotation of the geometry is not taken into account)
 * or in body frame coordinates (BF) which requires the point to be first transformed
42
43
 */

44
45
bool isPointInsideSphere(const Vec3& point,
                         const Vec3& spherePosition, const real_t sphereRadius )
46
47
48
49
{
   return !((point - spherePosition).sqrLength() > sphereRadius * sphereRadius);
}

50
51
bool isPointInsideHalfSpace(const Vec3& point,
                            const Vec3& halfSpacePosition, const Vec3& halfSpaceNormal )
52
{
53
   return !((point - halfSpacePosition) * halfSpaceNormal > real_t(0));
54
55
}

56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
bool isPointInsideBoxBF(const Vec3& pointBF,
                        const Vec3& edgeLengths )
{
   return std::fabs(pointBF[0]) <= real_t(0.5)*edgeLengths[0] &&
          std::fabs(pointBF[1]) <= real_t(0.5)*edgeLengths[1] &&
          std::fabs(pointBF[2]) <= real_t(0.5)*edgeLengths[2];
}

bool isPointInsideEllipsoidBF(const Vec3& pointBF,
                              const Vec3& semiAxes )
{
   return ( (pointBF[0] * pointBF[0])/(semiAxes[0] * semiAxes[0]) + (pointBF[1] * pointBF[1])/(semiAxes[1] * semiAxes[1])
            + (pointBF[2] * pointBF[2])/(semiAxes[2] * semiAxes[2]) <= 1_r );
}

bool isPointInsideCylindricalBoundary(const Vec3& point,
                                      const Vec3& cylindricalBoundaryPosition, const real_t radius, const Vec3& axis  )
{
   const Vec3 distanceFromCylinderCenterLine = (point - cylindricalBoundaryPosition) - ((point - cylindricalBoundaryPosition) * axis) * axis;
   return distanceFromCylinderCenterLine.sqrLength() >= radius*radius;
}
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

struct ContainsPointFunctor
{
   template<typename ParticleAccessor_T, typename Shape_T>
   bool operator()(const size_t /*particleIdx*/, const Shape_T& /*shape*/, const ParticleAccessor_T& /*ac*/, const Vector3<real_t>& /*point*/ )
   {
      WALBERLA_ABORT("ContainsPoint not implemented!");
   }

   template<typename ParticleAccessor_T>
   bool operator()(const size_t particleIdx, const mesa_pd::data::Sphere& sphere, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
   {
      static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");

      return isPointInsideSphere(point, ac.getPosition(particleIdx), sphere.getRadius() );
   }

   template<typename ParticleAccessor_T>
   bool operator()(const size_t particleIdx, const mesa_pd::data::HalfSpace& halfSpace, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
   {
      static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");

99
      return isPointInsideHalfSpace(point, ac.getPosition(particleIdx), halfSpace.getNormal() );
100
101
   }

102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
   template<typename ParticleAccessor_T>
   bool operator()(const size_t particleIdx, const mesa_pd::data::Box& box, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
   {
      static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");

      return isPointInsideBoxBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, point), box.getEdgeLength());
   }

   template<typename ParticleAccessor_T>
   bool operator()(const size_t particleIdx, const mesa_pd::data::Ellipsoid& ellipsoid, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
   {
      static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");

      return isPointInsideEllipsoidBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, point), ellipsoid.getSemiAxes());
   }

   template<typename ParticleAccessor_T>
   bool operator()(const size_t particleIdx, const mesa_pd::data::CylindricalBoundary& cylindricalBoundary, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
   {
      static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");

      return isPointInsideCylindricalBoundary(point, ac.getPosition(particleIdx), cylindricalBoundary.getRadius(), cylindricalBoundary.getAxis());
   }

126
127
128
129
130
};


} //namespace mesa_pd
} //namespace walberla