MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TriangleWall.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include <MercuryBase.h>
27 #include "TriangleWall.h"
28 #include "InteractionHandler.h"
29 #include "WallHandler.h"
30 #include "DPMBase.h"
31 //#include "Species/BaseSpecies.h"
32 #include "Particles/BaseParticle.h"
33 
34 bool setDistance(Mdouble& distance, const Vec3D& branch, Mdouble distanceMax)
35 {
36  const Mdouble distance2 = branch.getLengthSquared();
37  if (distance2 > distanceMax * distanceMax) return false;
38  distance = sqrt(distance2);
39  return true;
40 }
41 
42 
57 bool TriangleWall::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal) const
58 {
59  const Vec3D position = p.getPosition(); // note: no transfer to lab coordinates
60  const Mdouble distanceMax = p.getWallInteractionRadius(this);
61 
62  // compute distance from face
63  //get distance from particle position to the face
64  const Mdouble signedDistance = Vec3D::dot(position-vertex_[0], faceNormal_);
65  distance = fabs(signedDistance);
66 
67  // check if any contact is possible
68  if (distance >= distanceMax) return false;
69 
70  // compute distance from edges
71  const std::array<Vec3D,3> edgeBranch {position - vertex_[0], position - vertex_[1], position - vertex_[2]};
72  const std::array<double,3> edgeDistance {Vec3D::dot(edgeBranch[0], edgeNormal_[0]), Vec3D::dot(edgeBranch[1], edgeNormal_[1]), Vec3D::dot(edgeBranch[2], edgeNormal_[2])};
73 
74  // find edge with largest distance (this will be the edge if there is a edge contact)
75  const size_t id = (edgeDistance[1] > edgeDistance[2]) ?
76  (edgeDistance[0] > edgeDistance[1] ? 0 : 1) : (edgeDistance[0] > edgeDistance[2] ? 0 : 2);
77 
78  // check if there will be contact
79  if (edgeDistance[id] > distanceMax) return false;
80 
81  // determine if it is a face contact
82  const Vec3D posProjected = position - signedDistance * faceNormal_;
83  if (edgeDistance[id] <= 0 && isInsideTriangle(posProjected)){
84  normal = (signedDistance >= 0) ? -faceNormal_ : faceNormal_;
85  return true;
86  }
87 
88  // determine if it is an edge or vertex contact
89  const double positionAlongEdge = Vec3D::dot(edgeBranch[id], edge_[id]);
90  if (positionAlongEdge < 0) {
91  //possible contact with left vertex
92  distance = edgeBranch[id].getLength();
93  if (distance > distanceMax) return false;
94  normal = edgeBranch[id] / -distance;
95  } else if (positionAlongEdge > edgeLength_[id]) {
96  //contact with right vertex
97  const size_t idRight = (id + 1) % 3;
98  distance = edgeBranch[idRight].getLength();
99  if (distance > distanceMax) return false;
100  normal = edgeBranch[idRight] / -distance;
101  } else {
102  // edge contact
103  normal = edge_[id] * positionAlongEdge - edgeBranch[id];
104  distance = normal.getLength();
105  if (distance > distanceMax) return false;
106  normal /= distance;
107  }
108  return true;
109 }
110 
111 void TriangleWall::rotate(const Vec3D& angularVelocityDt)
112 {
113  if (!angularVelocityDt.isZero())
114  {
115  BaseInteractable::rotate(angularVelocityDt);
117  }
118 }
119 
123 void TriangleWall::read(std::istream& is)
124 {
125  BaseWall::read(is);
126  std::string dummy;
127  is >> dummy;
128  for (int i = 0; i < 3; i++)
129  {
130  is >> vertexInLabFrame_[i];
131  }
133 }
134 
139 void TriangleWall::write(std::ostream& os) const
140 {
141  BaseWall::write(os);
142  os << " vertexInLabFrame ";
143  for (int i = 0; i < 3; i++)
144  {
145  os << ' ' << vertexInLabFrame_[i];
146  }
147 }
148 
150 {
151  const unsigned long s = vtk.points.size();
152  for (auto v : vertex_)
153  {
154  vtk.points.push_back(v);
155  }
156  std::vector<double> cell;
157  cell.reserve(3);
158  cell.push_back(s);
159  cell.push_back(s + 1);
160  cell.push_back(s + 2);
161  //cell.push_back(s);
162  vtk.triangleStrips.push_back(cell);
163 }
164 
165 void TriangleWall::setVertices(const Vec3D A, const Vec3D B, const Vec3D C)
166 {
167  setPosition((A + B + C) / 3);
168  setOrientation({1, 0, 0, 0});
169  vertexInLabFrame_[0] = A - getPosition();
170  vertexInLabFrame_[1] = B - getPosition();
171  vertexInLabFrame_[2] = C - getPosition();
172  //vertexInLabFrame_[0] = A;
173  //vertexInLabFrame_[1] = B;
174  //vertexInLabFrame_[2] = C;
176 }
177 
184 void TriangleWall::move(const Vec3D& move)
185 {
188 }
189 
190 void TriangleWall::setVertices(const Vec3D A, const Vec3D B, const Vec3D C, const Vec3D position)
191 {
192  setPosition(position);
193  setOrientation({1, 0, 0, 0});
194  vertexInLabFrame_[0] = A - getPosition();
195  vertexInLabFrame_[1] = B - getPosition();
196  vertexInLabFrame_[2] = C - getPosition();
198 }
199 
207 {
208  for (int i = 0; i < 3; i++)
209  {
212  vertex_[i] += getPosition();
213  }
216 
217  edge_ = {vertex_[1] - vertex_[0], vertex_[2] - vertex_[1], vertex_[0] - vertex_[2]};
220 
221  for (int i = 0; i < 3; i++)
222  {
224  edgeLength_[i] = edge_[i].getLength();
225  edge_[i] /= edgeLength_[i];
226  edgeNormal_[i].normalise();
227  }
228  //logger(INFO,"vertex %,%,% edge %,%,% face %",vertex_[0],vertex_[1],vertex_[2],edgeNormal_[0],edgeNormal_[1],edgeNormal_[2],faceNormal_);
229 }
230 
231 bool TriangleWall::isLocal(Vec3D& min, Vec3D& max) const
232 {
233  min = vertexMin_;
234  max = vertexMax_;
235  return true;
236 }
237 
238 bool TriangleWall::isInsideTriangle(const Vec3D &point) const
239 {
240  const Vec3D branch0 = point - vertex_[0];
241  const Vec3D branch1 = point - vertex_[1];
242  const Vec3D branch2 = point - vertex_[2];
243 
244  //compute total area
245  Mdouble s = sqrt(Vec3D::cross(vertex_[1] - vertex_[0], vertex_[2] - vertex_[1]).getLengthSquared());
246  //compute barycentric coordinates
247  Mdouble s0 = sqrt(Vec3D::cross(branch0, branch1).getLengthSquared())/s;
248  Mdouble s1 = sqrt(Vec3D::cross(branch1, branch2).getLengthSquared())/s;
249  Mdouble s2 = sqrt(Vec3D::cross(branch2, branch0).getLengthSquared())/s;
250  return (1 > s0 > 0 && 1 > s1 > 0 && 1 > s2 > 0);
251 }
Vec3D faceNormal_
Definition: TriangleWall.h:164
bool setDistance(Mdouble &distance, const Vec3D &branch, Mdouble distanceMax)
Definition: TriangleWall.cc:34
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
double Mdouble
Definition: GeneralDefine.h:34
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
std::array< Vec3D, 3 > edge_
Definition: TriangleWall.h:158
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void updateVertexAndNormal()
This function should be called after setting either position_ or vertexInLabFrame_.
Vec3D vertexMax_
Definition: TriangleWall.h:152
void write(std::ostream &os) const override
Function that writes a BaseWall to an output stream, usually a restart file.
Definition: BaseWall.cc:102
void writeVTK(VTKContainer &vtk) const override
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
void rotate(Vec3D &position) const
Applies the rotation to a position.
Definition: Quaternion.cc:563
bool isZero() const
Checks if ALL elements are zero.
Definition: Vector.h:102
void read(std::istream &is) override
Reads an TriangleWall from an input stream, for example a restart file.
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
Definition: Vector.cc:331
std::array< double, 3 > edgeLength_
Definition: TriangleWall.h:159
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
bool isInsideTriangle(const Vec3D &point) const
static Vec3D min(const Vec3D &a, const Vec3D &b)
Calculates the pointwise minimum of two Vec3D.
Definition: Vector.cc:102
bool getDistanceAndNormal(const BaseParticle &p, Mdouble &distance, Vec3D &normal_return) const override
Compute the distance from the wall for a given BaseParticle and return if there is a collision...
Definition: TriangleWall.cc:57
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
Vec3D vertexMin_
Definition: TriangleWall.h:151
void setVertices(Vec3D A, Vec3D B, Vec3D C)
Sets member variables such that the wall represents a triangle with vertices A, B, C.
std::array< Vec3D, 3 > vertex_
Definition: TriangleWall.h:146
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
void write(std::ostream &os) const override
Writes an TriangleWall to an output stream, for example a restart file.
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:80
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
std::vector< Vec3D > points
Definition: BaseWall.h:38
virtual void rotate(const Vec3D &angularVelocityDt)
Rotates this BaseInteractable.
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:316
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
Definition: Vector.h:49
virtual void move(const Vec3D &move)
Moves this BaseInteractable by adding an amount to the position.
std::array< Vec3D, 3 > edgeNormal_
Definition: TriangleWall.h:157
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
void move(const Vec3D &move) override
void rotate(const Vec3D &angularVelocity) override
Rotates this BaseInteractable.
std::array< Vec3D, 3 > vertexInLabFrame_
Definition: TriangleWall.h:141
bool isLocal(Vec3D &min, Vec3D &max) const override
static Vec3D max(const Vec3D &a, const Vec3D &b)
Calculates the pointwise maximum of two Vec3D.
Definition: Vector.cc:89