MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HorizontalBaseScrew.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 "HorizontalBaseScrew.h"
27 #include "Particles/BaseParticle.h"
28 #include "WallHandler.h"
29 #include "DPMBase.h"
30 
32 {
33  //BEGIN CHANGE
34  sinA2Max_ = 0;
35  timeMin_ = 0;
36  //END CHANGE
37  logger(DEBUG, "HorizontalBaseScrew() finished");
38 }
39 
44  : IntersectionOfWalls(other)
45 {
46  //BEGIN CHANGE
47  sinA2Max_ = other.sinA2Max_;
48  timeMin_ = other.timeMin_;
49  //END CHANGE
50  logger(DEBUG, "HorizontalBaseScrew(const HorizontalBaseScrew &p) finished");
51 }
52 
53 HorizontalBaseScrew::HorizontalBaseScrew(Vec3D position, Vec3D orientation, std::vector<HorizontalBaseScrew::normalAndPosition> walls, const ParticleSpecies* species, Mdouble sinA2Max, Mdouble timeMin)
54  : IntersectionOfWalls(walls, species), sinA2Max_(sinA2Max), timeMin_(timeMin)
55 {
56  setPosition(position);
57  setOrientationViaNormal(orientation);
58 }
59 
61 {
62  logger(DEBUG, "~HorizontalBaseScrew() finished.");
63 }
64 
69 {
70  if (this == &other)
71  {
72  return *this;
73  }
74  else
75  {
76  return *(other.copy());
77  }
78 }
79 
84 {
85  return new HorizontalBaseScrew(*this);
86 }
87 
98 bool HorizontalBaseScrew::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normalReturn) const
99 {
100  //transform to axisymmetric coordinates
101  //move the coordinate system to the axis origin, so pOrigin=(xhat,yhat,zhat)
102  Vec3D pOrigin = p.getPosition() -getPosition();
104  Vec3D axis = getOrientation().getAxis();
105  Mdouble a = Vec3D::dot(pOrigin, axis);
106  //Vec3D(r,a) is the projection into the radial-axial plane
107  Vec3D radialDirection = pOrigin - a * axis;
108  Mdouble r = radialDirection.getLength();
109 
110  //BEGIN CHANGE
111  //if in contact
112  if (r!=0)
113  radialDirection /= r;
114  else //in this case the tangential vector is irrelevant
115  logger(WARN, "Warning: Particle % is exactly on the symmetry axis of wall %", p.getIndex(), getIndex());
116 
117  //tangential is the projection into the (xhat,yhat) plane
119  Vec3D orientation = Vec3D(cos(angle),sin(angle),0);
120  Mdouble sinA2 = Vec3D::getLengthSquared(Vec3D::cross(radialDirection,orientation));
121  if (sinA2>sinA2Max_) {
122  a+=sinA2-sinA2Max_;
123  }
124  //END CHANGE
125 
126  //determine wall distance, normal and contact in axissymmetric coordinates
127  //and transform from axisymmetric coordinates
128  Vec3D normal;
129  if (!IntersectionOfWalls::getDistanceAndNormal(Vec3D(r, 0, a), p.getWallInteractionRadius(this), distance, normal))
130  {
131  //if not in contact
132  return false;
133  }
134  else
135  {
136  //if in contact
137  if (r!=0)
138  radialDirection /= r;
139  else //in this case the tangential vector is irrelevant
140  logger(WARN, "Warning: Particle % is exactly on the symmetry axis of wall %", p.getIndex(), getIndex());
141  normalReturn = normal.Z * axis + normal.X * radialDirection;
142  //logger.assert(normalReturn.Y==0,"Error");
143  return true;
144  }
145 }
146 
151 void HorizontalBaseScrew::read(std::istream& is)
152 {
154 }
155 
160 void HorizontalBaseScrew::write(std::ostream& os) const
161 {
163 }
164 
168 std::string HorizontalBaseScrew::getName() const
169 {
170  return "HorizontalBaseScrew";
171 }
172 
174 {
176 }
177 
179  //define the box in the rz plane
180  double x[2] = {min.X, max.X};
181  double y[2] = {min.Y, max.Y};
182  double z[2] = {min.Z, max.Z};
184  Vec3D o = q.getAxis();
185 
186  //look at the corner points of the domain and see that the wall is displayed in the domain (possibly a bit outside as well)
187  min -= getPosition();
188  q.rotate(min); //set min/max initial values to values of first corner point
189  max = min;
190  Mdouble maxXSquared = 0;
191  //
192  for (unsigned i=0; i<2; i++)
193  for (unsigned j=0; j<2; j++)
194  for (unsigned k=0; k<2; k++) {
195  Vec3D p = Vec3D(x[i],y[j],z[k])-getPosition();
196  Mdouble Z = Vec3D::dot(p,o);
197  Mdouble XSquared = Vec3D::getLengthSquared(p-Z*o);
198  if (XSquared>maxXSquared) maxXSquared=XSquared;
199  if (Z<min.Z) min.Z=Z;
200  if (Z>max.Z) max.Z=Z;
201  }
202  min.X = 0;
203  max.X = sqrt(maxXSquared);
204  min.Y = 0;
205  max.Y = 0.001; //unused
206  //logger(ERROR,"min % max %",min,max);
207 }
208 
210 {
211  for (auto wall=wallObjects_.begin(); wall!=wallObjects_.end(); wall++)
212  {
213  //plot each of the intersecting walls
214  std::vector<Vec3D> myPoints;
215 
216  //first create a slice of non-rotated wall in the xz plane, 0<y<1
217  Vec3D min = getHandler()->getDPMBase()->getMin();
218  Vec3D max = getHandler()->getDPMBase()->getMax();
219  convertLimits(min, max);
220 
221  //create the basic slice for the first wall using the InfiniteWall routine
222  wall->createVTK (myPoints,min,max);
223 
224  //create intersections with the other walls, similar to the IntersectionOfWalls routine
225  for (auto other=wallObjects_.begin(); other!=wallObjects_.end(); other++)
226  {
227  if (other!=wall)
228  {
229  intersectVTK(myPoints, -other->getNormal(), other->getPosition());
230  }
231  }
232 
233  //only keep the y=0 values
234  std::vector<Vec3D> rzVec;
235  for (auto& p: myPoints)
236  {
237  if (p.Y==0)
238  {
239  rzVec.push_back(p);
240  }
241  }
242  if (rzVec.size()==0)
243  return;
244 
245  //create points on the unit circle
246  unsigned nr = 180;
247  //BEGIN CHANGE
248  std::vector<Vec3D> xyVec;
250  for (unsigned ir=0; ir<nr; ir++) {
251  Mdouble angle = 2.0 * constants::pi * ir/nr;
252 
253  Mdouble sinA2 = sin(a1 - angle); sinA2 *= sinA2;
254  Mdouble dSin = sinA2-sinA2Max_;
255 
256  xyVec.push_back(Vec3D(cos(angle),sin(angle),fmax(dSin,0)));
257  }
258  //END CHANGE
259 
260  //now create rings of points on the axisym. shape
262  unsigned nPoints = vtk.points.size();
263  Vec3D p;
264  Vec3D o = getOrientation().getAxis();
265  for (auto rz : rzVec)
266  {
267  for (auto xy : xyVec)
268  {
269  //BEGIN CHANGE
270  p = Vec3D(rz.X * xy.X, rz.X * xy.Y,fmax(rz.Z - xy.Z,0.0));
271  //getOrientation().rotate(p);
272  //END CHANGE
273  p += getPosition();
274  vtk.points.push_back(p);
275  }
276  }
277 
278  //finally create the connectivity matri to plot shell-like triangle strips.
279  unsigned nz = rzVec.size();
280  unsigned nCells = vtk.triangleStrips.size();
281  for (unsigned iz=0; iz<nz-1; iz++) {
282  std::vector<double> cell;
283  cell.reserve(2*nr+2);
284  for (unsigned ir=0; ir<nr; ir++) {
285  cell.push_back(nPoints+ir+iz*nr);
286  cell.push_back(nPoints+ir+(iz+1)*nr);
287  }
288  cell.push_back(nPoints+iz*nr);
289  cell.push_back(nPoints+(iz+1)*nr);
290  vtk.triangleStrips.push_back(cell);
291  }
292  }
293 
294 }
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
void write(std::ostream &os) const final
outputs wall
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 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 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
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void convertLimits(Vec3D &min, Vec3D &max) const
HorizontalBaseScrew()
Default constructor.
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...
virtual const Vec3D & getAngularVelocity() const
Returns the angular velocity of this interactable.
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
Mdouble getLengthSquared() const
Calculates the squared length of this Vec3D: .
Definition: Vector.cc:184
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
void read(std::istream &is) final
reads wall
A HorizontalBaseScrew is a copy of AxisymmetricIntersectionOfWalls, with an additional, angle-dependent component.
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
const Mdouble pi
Definition: ExtendedMath.h:45
void writeVTK(VTKContainer &vtk) const override
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
HorizontalBaseScrew & operator=(const HorizontalBaseScrew &other)
Copy assignment operator.
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
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...
Vec3D getMax() const
Definition: DPMBase.h:629
~HorizontalBaseScrew()
Destructor.
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
std::string getName() const final
Returns the name of the object.
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
HorizontalBaseScrew * copy() const final
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
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
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:797
void write(std::ostream &os) const override
Writes an IntersectionOfWalls to an output stream, for example a restart file.