MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TriangulatedWall.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 "TriangulatedWall.h"
27 #include "InteractionHandler.h"
28 #include "Particles/BaseParticle.h"
29 #include <fstream>
30 #include <stdio.h>
31 #include <stdlib.h>
32 
34 {
35  logger(DEBUG, "TriangulatedWall() constructed.");
36 }
37 
42  : BaseWall(other)
43 {
44  face_ = other.face_;
45  vertex_ = other.vertex_;
46  //now reset the pointers
47  for (unsigned f = 0; f < face_.size(); f++)
48  {
49  for (unsigned i = 0; i < 3; i++)
50  {
51  if (other.face_[f].neighbor[i])
52  {
53  face_[f].neighbor[i] = &face_[other.face_[f].neighbor[i] - &other.face_[0]];
54  } //else nullptr
55  face_[f].vertex[i] = &vertex_[other.face_[f].vertex[i] - &other.vertex_[0]];
56  }
57  }
58  logger(DEBUG, "TriangulatedWall(TriangulatedWall&) constructed.");
59 }
60 
61 TriangulatedWall::TriangulatedWall(std::string filename, const ParticleSpecies* species)
62 {
63  setSpecies(species);
64  readVTK(filename);
65 }
66 
70 void TriangulatedWall::readVTK(std::string filename)
71 {
72  std::fstream file;
73  file.open(filename.c_str(), std::ios::in);
74  logger.assert_always(file.is_open(), "File opening failed: %", filename);
75 
76  std::string dummy;
77  getline(file, dummy);
78  getline(file, dummy);
79  getline(file, dummy);
80  getline(file, dummy);
81  //read points
82  unsigned num;
83  file >> dummy >> num >> dummy;
84  vertex_.reserve(num);
85  Vec3D v;
86  for (unsigned i = 0; i < num; i++)
87  {
88  file >> v.X >> v.Y >> v.Z;
89  vertex_.push_back(v);
90  }
91  //read faces
92  file >> dummy >> num >> dummy;
93  face_.reserve(num);
94  Face f;
95  unsigned id0, id1, id2;
96  for (unsigned i = 0; i < num; i++)
97  {
98  file >> dummy >> id0 >> id1 >> id2;
99  f.vertex[0] = &vertex_[id0];
100  f.vertex[1] = &vertex_[id1];
101  f.vertex[2] = &vertex_[id2];
102  face_.push_back(f);
103  }
104  file >> dummy;
105  //set normals and positions
106  for (auto& face: face_)
107  {
108  face.normal = Vec3D::getUnitVector(
109  Vec3D::cross(*face.vertex[1] - *face.vertex[0], *face.vertex[2] - *face.vertex[0]));
110  }
111  //set neighbours
112  for (auto face0 = face_.begin(); face0 + 1 != face_.end(); face0++)
113  {
114  for (auto face1 = face0 + 1; face1 != face_.end(); face1++)
115  {
116  if (face0->vertex[0] == face1->vertex[0])
117  {
118  if (face0->vertex[1] == face1->vertex[2])
119  { //edge 0=2
120  face0->neighbor[0] = &*face1;
121  face1->neighbor[2] = &*face0;
122  }
123  else if (face0->vertex[2] == face1->vertex[1])
124  { //edge 2=0
125  face0->neighbor[2] = &*face1;
126  face1->neighbor[0] = &*face0;
127  }
128  }
129  else if (face0->vertex[0] == face1->vertex[1])
130  {
131  if (face0->vertex[1] == face1->vertex[0])
132  { //edge 0=0
133  face0->neighbor[0] = &*face1;
134  face1->neighbor[0] = &*face0;
135  }
136  else if (face0->vertex[2] == face1->vertex[2])
137  { //edge 2=1
138  face0->neighbor[2] = &*face1;
139  face1->neighbor[1] = &*face0;
140  }
141  }
142  else if (face0->vertex[0] == face1->vertex[2])
143  {
144  if (face0->vertex[1] == face1->vertex[1])
145  { //edge 0=1
146  face0->neighbor[0] = &*face1;
147  face1->neighbor[1] = &*face0;
148  }
149  else if (face0->vertex[2] == face1->vertex[0])
150  { //edge 2=2
151  face0->neighbor[2] = &*face1;
152  face1->neighbor[2] = &*face0;
153  }
154  }
155  else if (face0->vertex[1] == face1->vertex[0])
156  {
157  if (face0->vertex[2] == face1->vertex[2])
158  { //edge 1=2
159  face0->neighbor[1] = &*face1;
160  face1->neighbor[2] = &*face0;
161  }
162  }
163  else if (face0->vertex[1] == face1->vertex[1])
164  {
165  if (face0->vertex[2] == face1->vertex[0])
166  { //edge 1=0
167  face0->neighbor[1] = &*face1;
168  face1->neighbor[0] = &*face0;
169  }
170  }
171  else if (face0->vertex[1] == face1->vertex[2])
172  {
173  if (face0->vertex[2] == face1->vertex[1])
174  { //edge 1=1
175  face0->neighbor[1] = &*face1;
176  face1->neighbor[1] = &*face0;
177  }
178  }
179 
180  }
181  }
182  //set edge normals (inwards facing)
183  for (auto& face: face_)
184  {
185  face.edgeNormal[0] = Vec3D::getUnitVector(Vec3D::cross(face.normal, *face.vertex[1] - *face.vertex[0]));
186  face.edgeNormal[1] = Vec3D::getUnitVector(Vec3D::cross(face.normal, *face.vertex[2] - *face.vertex[1]));
187  face.edgeNormal[2] = Vec3D::getUnitVector(Vec3D::cross(face.normal, *face.vertex[0] - *face.vertex[2]));
188  }
189  file.close();
190 }
191 
193 {
194  logger(DEBUG, "~TriangulatedWall() has been called.");
195 }
196 
201 {
202  logger(DEBUG, "TriangulatedWall::operator= called.");
203  if (this == &other)
204  {
205  return *this;
206  }
207  return *(other.copy());
208 }
209 
214 {
215  return new TriangulatedWall(*this);
216 }
217 
234 bool TriangulatedWall::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return) const
235 {
236  Mdouble interactionRadius = p.getWallInteractionRadius(this);
237  //it's important to use a reference here, as the faces use pointers
238  for (const auto& face : face_)
239  {
240  if (face.getDistanceAndNormal(p, distance, normal_return, interactionRadius))
241  return true;
242  }
243  return false;
244 }
245 
257 void TriangulatedWall::move(const Vec3D& move)
258 {
260  for (auto& v : vertex_)
261  {
262  v += move;
263  }
264 }
265 
269 void TriangulatedWall::read(std::istream& is)
270 {
272 }
273 
278 void TriangulatedWall::write(std::ostream& os) const
279 {
280  os << "Vertices " << vertex_.size();
281  for (const auto& vertex: vertex_)
282  os << " " << vertex;
283  os << std::endl;
284  //os << "Faces " << face_.size() << std::endl;
285  unsigned counter = 0;
286  for (const auto& face: face_)
287  {
288  os << "Face " << counter++
289  << " vertex " << face.vertex[0] - &vertex_[0]
290  << " " << face.vertex[1] - &vertex_[0]
291  << " " << face.vertex[2] - &vertex_[0]
292  << " neighbor " << (face.neighbor[0] ? (face.neighbor[0] - &face_[0]) : -1)
293  << " " << (face.neighbor[1] ? (face.neighbor[1] - &face_[0]) : -1)
294  << " " << (face.neighbor[2] ? (face.neighbor[2] - &face_[0]) : -1)
295  << " normal " << face.normal
296  << " edgeNormal " << face.edgeNormal[0]
297  << " " << face.edgeNormal[1]
298  << " " << face.edgeNormal[2]
299  << std::endl;
300  }
301 }
302 
306 std::string TriangulatedWall::getName() const
307 {
308  return "TriangulatedWall";
309 }
310 
319  InteractionHandler* interactionHandler)
320 {
321  Mdouble distance;
322  Vec3D normal;
323  Mdouble interactionRadius = p->getWallInteractionRadius(this);
324  for (const auto& face : face_)
325  {
326  if (face.getDistanceAndNormal(*p, distance, normal,interactionRadius))
327  {
328  normal = -normal;
329  BaseInteraction* const c = interactionHandler->getInteraction(p, this, timeStamp, normal);
330  c->setNormal(normal);
331  c->setDistance(distance);
332  c->setOverlap(p->getRadius() - distance);
333  c->setContactPoint(p->getPosition() - distance * normal);
334  return c;
335  }
336  }
337  return nullptr;
338 }
339 
341 {
342  return Vec3D::dot(*vertex[0] - otherPosition, normal);
343 }
344 
349 bool TriangulatedWall::Face::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return, Mdouble interactionRadius) const
350 {
351  //check if particle overlaps
352  distance = getDistance(p.getPosition());
353  //return if distance is large, or particle has too much overlap
354  if (fabs(distance) >= interactionRadius)
355  return false;
356 
357  //Contact radius squared
358  Mdouble allowedDistanceSquared = interactionRadius * interactionRadius - distance * distance;
359  int touchingEdge = -1;
360  Vec3D* touchingVertex = nullptr;
361  //loop through edges to find out if there is a contact with face, edge, vertex, or neither.
362  for (unsigned i = 0; i < 3; i++)
363  {
364  //distance from the edge
365  Mdouble edgeDistance = Vec3D::dot(edgeNormal[i], *vertex[i] - p.getPosition());
366  if (edgeDistance <= 0) //if in contact with the face
367  continue;
368  if (edgeDistance * edgeDistance >= allowedDistanceSquared) //if not in contact with anything
369  return false;
370  //this point is only reached if in contact, but not with the face
371  //if edge is convex, treat as face, not edge contact
372  bool convex = (neighbor[i] && ((neighbor[i]->getDistance(p.getPosition()) < 0) ? (
373  Vec3D::dot(neighbor[i]->normal, edgeNormal[i]) > 1e-12) : (
374  Vec3D::dot(neighbor[i]->normal, edgeNormal[i]) < -1e-12)));
375  if (touchingEdge == -1)
376  {
377  //edge or vertex contact (depending if more edges are found)
378  touchingEdge = i;
379  }
380  else
381  {
382  //vertex contact
383  touchingVertex = vertex[(i == 2 && touchingEdge == 0) ? 0 : i];
384  }
385  //if (!convex) continue;
386  //do not compute if neighbor exists and has lower index or face contact.
387  if (neighbor[i] > this)
388  { //if neighbor has higher index
389  if (neighbor[i]->neighbor[0] == this)
390  {
391  edgeDistance = Vec3D::dot(neighbor[i]->edgeNormal[0], *vertex[i] - p.getPosition());
392  }
393  else if (neighbor[i]->neighbor[1] == this)
394  {
395  edgeDistance = Vec3D::dot(neighbor[i]->edgeNormal[1], *vertex[i] - p.getPosition());
396  }
397  else
398  { //if (neighbor[i]->neighbor[2]==this)
399  edgeDistance = Vec3D::dot(neighbor[i]->edgeNormal[2], *vertex[i] - p.getPosition());
400  }
401  if (edgeDistance <= 0 && !convex) //if neighbor has face contact, ignore the edge contact
402  return false;
403  }
404  else if (neighbor[i] && !convex)
405  {
406  return false;
407  }
408  }
409 
410 
411  //check if convex neighbours are overlapping as well
412  if (touchingVertex)
413  { //vertex contact
414  normal_return = *touchingVertex - p.getPosition();
415  distance = Vec3D::getLength(normal_return);
416  normal_return /= distance;
417  //return false;
418  //logger(INFO,"vertex contact");
419  }
420  else if (touchingEdge == -1)
421  { //face contact
422  if (distance >= 0)
423  { // front face contact
424  normal_return = normal;
425  }
426  else
427  { // back face contact
428  normal_return = -normal;
429  distance = -distance;
430  }
431  //return false;
432  //logger(INFO,"face contact");
433  }
434  else
435  { //edge contact
436  Vec3D VP = *vertex[touchingEdge] - p.getPosition();
437  Vec3D VW = *vertex[touchingEdge] - *vertex[(touchingEdge == 2) ? 0 : (touchingEdge + 1)];
438  normal_return = VP - Vec3D::dot(VP, VW) / Vec3D::getLengthSquared(VW) * VW;
439  distance = Vec3D::getLength(normal_return);
440  normal_return /= distance;
441  //return false;
442  //logger(INFO,"edge contact at c=% p=% d=% n=%",p.getPosition() + distance * normal_return,p.getPosition(), distance, normal_return);
443  }
444  return true;
445 }
446 
448 {
449  const int s = vtk.points.size();
450  for (auto v : vertex_)
451  {
452  vtk.points.push_back(v);
453  }
454  for (auto f : face_)
455  {
456  std::vector<double> cell;
457  cell.reserve(3);
458  cell.push_back(s + f.vertex[0] - &vertex_[0]);
459  cell.push_back(s + f.vertex[1] - &vertex_[0]);
460  cell.push_back(s + f.vertex[2] - &vertex_[0]);
461  vtk.triangleStrips.push_back(cell);
462  }
463 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
void setNormal(Vec3D normal)
Sets the normal vector between the two interacting objects.
Mdouble X
the vector components
Definition: Vector.h:65
TriangulatedWall()
Default constructor.
static Vec3D getUnitVector(const Vec3D &a)
Returns a unit Vec3D based on a.
Definition: Vector.cc:345
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
void setOverlap(Mdouble overlap)
Set the overlap between the two interacting object.
double Mdouble
Definition: GeneralDefine.h:34
A TriangulatedWall is a triangulation created from a set of vertices and a n-by-3 connectivity matrix...
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
BaseInteraction * getInteractionWith(BaseParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler) override
Get the interaction between this TriangulatedWall and given BaseParticle at a given time...
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...
void setContactPoint(Vec3D contactPoint)
Set the location of the contact point between the two interacting objects.
void write(std::ostream &os) const override
Writes an TriangulatedWall to an output stream, for example a restart file.
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
void setDistance(Mdouble distance)
Sets the interaction distance between the two interacting objects.
Mdouble getLengthSquared() const
Calculates the squared length of this Vec3D: .
Definition: Vector.cc:184
Stores information about interactions between two interactable objects; often particles but could be ...
std::vector< Vec3D > vertex_
Struct used to store the properties of a face needed for contact detection.
void read(std::istream &is) override
Reads an TriangulatedWall from an input stream, for example a restart file.
TriangulatedWall & operator=(const TriangulatedWall &other)
void readVTK(std::string filename)
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
BaseInteraction * getInteraction(BaseInteractable *P, BaseInteractable *I, unsigned timeStamp)
Returns the Interaction between the BaseInteractable's P and I.
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
Container to store Interaction objects.
std::string getName() const override
Returns the name of the object, here the string "TriangulatedWall".
Vec3D normal
face normal (vertices are ordered anticlockwise direction around the normal)
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
void writeVTK(VTKContainer &vtk) const override
Basic class for walls.
Definition: BaseWall.h:47
std::array< Vec3D *, 3 > vertex
defines the three vertices (anticlockwise direction around the normal)
Mdouble Y
Definition: Vector.h:65
void move(const Vec3D &move) override
Move the TriangulatedWall to a new position, which is a Vec3D from the old position.
std::vector< Vec3D > points
Definition: BaseWall.h:38
Mdouble getLength() const
Calculates the length of this Vec3D: .
Definition: Vector.cc:320
std::vector< Face > face_
~TriangulatedWall() override
Destructor.
TriangulatedWall * copy() const override
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
Mdouble getDistance(const Vec3D &otherPosition) const
computes the signed distance to the face in normal direction
bool getDistanceAndNormal(const BaseParticle &p, Mdouble &distance, Vec3D &normal_return, Mdouble interactionRadius) const
Returns true if contact with the face exists, false if not. If contact exists, then the distance and ...
Definition: Vector.h:49
virtual void move(const Vec3D &move)
Moves this BaseInteractable by adding an amount to the position.
Mdouble Z
Definition: Vector.h:65
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:171