MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AxisymmetricIntersectionOfWalls.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 
27 #include "Particles/BaseParticle.h"
28 #include "WallHandler.h"
29 #include "DPMBase.h"
30 
32 {
33  logger(DEBUG, "AxisymmetricIntersectionOfWalls() finished");
34 }
35 
40  : IntersectionOfWalls(other)
41 {
42  logger(DEBUG, "AxisymmetricIntersectionOfWalls(const AxisymmetricIntersectionOfWalls &p) finished");
43 }
44 
46  std::vector<AxisymmetricIntersectionOfWalls::normalAndPosition> walls,
47  const ParticleSpecies* species)
48  : IntersectionOfWalls(walls, species)
49 {
50  setPosition(position);
51  setOrientationViaNormal(orientation);
52 }
53 
55 {
56  logger(DEBUG, "~AxisymmetricIntersectionOfWalls() finished.");
57 }
58 
64 {
65  if (this == &other)
66  {
67  return *this;
68  }
69  else
70  {
71  return *(other.copy());
72  }
73 }
74 
79 {
80  return new AxisymmetricIntersectionOfWalls(*this);
81 }
82 
94  Vec3D& normalReturn) const
95 {
96  //transform to axisymmetric coordinates
97  //move the coordinate system to the axis origin, so pOrigin=(xhat,yhat,zhat)
98  Vec3D pOrigin = p.getPosition() - getPosition();
100  Vec3D axis = getOrientation().getAxis();
101  Mdouble a = Vec3D::dot(pOrigin, axis);
102  //Vec3D(r,a) is the projection into the radial-axial plane
103  Vec3D radialDirection = pOrigin - a * axis;
104  Mdouble r = radialDirection.getLength();
105  Vec3D normal;
106  //determine wall distance, normal and contact in axissymmetric coordinates
107  //and transform from axisymmetric coordinates
108  if (!IntersectionOfWalls::getDistanceAndNormal(Vec3D(r, 0, a), p.getWallInteractionRadius(this), distance, normal))
109  {
110  //if not in contact
111  return false;
112  }
113  else
114  {
115  //if in contact
116  if (r != 0)
117  radialDirection /= r;
118  else //in this case the tangential vector is irrelevant
119  logger(WARN, "Warning: Particle % is exactly on the symmetry axis of wall %", p.getIndex(), getIndex());
120  normalReturn = normal.Z * axis + normal.X * radialDirection;
121  //logger.assert(normalReturn.Y==0,"Error");
122  return true;
123  }
124 }
125 
131 {
133 }
134 
139 void AxisymmetricIntersectionOfWalls::write(std::ostream& os) const
140 {
142 }
143 
148 {
149  return "AxisymmetricIntersectionOfWalls";
150 }
151 
153 {
155 }
156 
158 {
160  Vec3D rMin = min - getPosition();
161  q.rotateBack(rMin); //set min/max initial values to values of first corner point
162  Vec3D rMax = max - getPosition();
163  q.rotateBack(rMax); //set min/max initial values to values of first corner point
164 
165  Mdouble r = std::sqrt(std::max(rMax.Y * rMax.Y + rMax.Z * rMax.Z, rMin.Y * rMin.Y + rMin.Z * rMin.Z));
166  max = Vec3D(r, 0.001, std::max(rMin.X,rMax.X));
167  min = Vec3D(0, 0, std::min(rMin.X,rMax.X));
168  //std::cout << "r=" << r << std::endl;
169 }
170 
172 {
173  for (auto wall = wallObjects_.begin(); wall != wallObjects_.end(); wall++)
174  {
175  //plot each of the intersecting walls
176  std::vector<Vec3D> myPoints;
177 
178  //first create a slice of non-rotated wall in the xz plane, 0<y<1
179  Vec3D min = getHandler()->getDPMBase()->getMin();
180  Vec3D max = getHandler()->getDPMBase()->getMax();
181  convertLimits(min, max);
182 
183  //create the basic slice for the first wall using the InfiniteWall routine
184  wall->createVTK(myPoints, min, max);
185 
186  //create intersections with the other walls, similar to the IntersectionOfWalls routine
187  for (auto other = wallObjects_.begin(); other != wallObjects_.end(); other++)
188  {
189  if (other != wall)
190  {
191  intersectVTK(myPoints, -other->getNormal(), other->getPosition());
192  }
193  }
194 
195  //only keep the y=0 values
196  std::vector<Vec3D> rzVec;
197  for (auto& p: myPoints)
198  {
199  if (p.Y == 0)
200  {
201  rzVec.push_back(p);
202  }
203  }
204  if (rzVec.empty())
205  return;
206 
207  //create points on the unit circle
208  unsigned nr = 180;
209  struct XY
210  {
211  double X;
212  double Y;
213  };
214  std::vector<XY> xyVec;
215  for (unsigned ir = 0; ir < nr; ir++)
216  {
217  Mdouble angle = 2.0 * constants::pi * ir / nr;
218  xyVec.push_back({cos(angle), sin(angle)});
219  }
220 
221  //now create rings of points on the axisym. shape
223  unsigned long nPoints = vtk.points.size();
224  Vec3D p;
225  //Vec3D o = getOrientation().getAxis();
226  for (auto rz : rzVec)
227  {
228  for (auto xy : xyVec)
229  {
230  p = Vec3D(rz.Z, rz.X * xy.X, rz.X * xy.Y);
231  getOrientation().rotate(p);
232  p += getPosition();
233  vtk.points.push_back(p);
234  }
235  }
236 
237  //finally create the connectivity matri to plot shell-like triangle strips.
238  unsigned long nz = rzVec.size();
239  unsigned long nCells = vtk.triangleStrips.size();
240  for (unsigned iz = 0; iz < nz - 1; iz++)
241  {
242  std::vector<double> cell;
243  cell.reserve(2 * nr + 2);
244  for (unsigned ir = 0; ir < nr; ir++)
245  {
246  cell.push_back(nPoints + ir + iz * nr);
247  cell.push_back(nPoints + ir + (iz + 1) * nr);
248  }
249  cell.push_back(nPoints + iz * nr);
250  cell.push_back(nPoints + (iz + 1) * nr);
251  vtk.triangleStrips.push_back(cell);
252  }
253  }
254 
255 }
Implementation of a 3D quaternion (by Vitaliy).
Definition: Quaternion.h:62
unsigned int getIndex() const
Returns the index of the object in the handler.
Definition: BaseObject.h:118
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
A IntersectionOfWalls is convex polygon defined as an intersection of InfiniteWall's.
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
void intersectVTK(std::vector< Vec3D > &points, Vec3D normal, Vec3D position) const
Definition: BaseWall.cc:243
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
void convertLimits(Vec3D &min, Vec3D &max) const
void read(std::istream &is) override
Move the IntersectionOfWalls to a new position, which is a Vec3D from the old position.
Vec3D getMin() const
Definition: DPMBase.h:623
void write(std::ostream &os) const final
outputs wall
bool getDistanceAndNormal(const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const final
Computes the distance from the wall for a given BaseParticle and returns true if there is a collision...
std::string getName() const final
Returns the name of the object.
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...
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
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
Definition: Vector.cc:331
std::vector< InfiniteWall > wallObjects_
The wall "segments"/directions that together make up the finite wall.
Vec3D getAxis() const
Converts the quaternions into a normal vector by rotating the vector x=(1,0,0); see See Wiki for deta...
Definition: Quaternion.cc:501
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:136
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
void writeVTK(VTKContainer &vtk) const override
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
const Mdouble pi
Definition: ExtendedMath.h:45
AxisymmetricIntersectionOfWalls * copy() const final
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
AxisymmetricIntersectionOfWalls & operator=(const AxisymmetricIntersectionOfWalls &other)
Copy assignment operator.
void read(std::istream &is) final
reads wall
Vec3D getMax() const
Definition: DPMBase.h:629
Use AxisymmetricIntersectionOfWalls to Screw Screw::read Screw::read Screw::read define axisymmetric ...
Mdouble Y
Definition: Vector.h:65
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
std::vector< Vec3D > points
Definition: BaseWall.h:38
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
Definition: Vector.h:49
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
void write(std::ostream &os) const override
Writes an IntersectionOfWalls to an output stream, for example a restart file.