MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SuperQuadricParticle.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 "DPMBase.h"
27 #include "BaseParticle.h"
28 #include <cmath>
29 #include "SuperQuadricParticle.h"
30 #include "InteractionHandler.h"
31 #include "ParticleHandler.h"
32 #include "SpeciesHandler.h"
33 
39  : BaseParticle()
40 {
41  axes_ = Vec3D(1.0, 1.0, 1.0);
42  eps1_ = 1.0;
43  eps2_ = 1.0;
44  logger(DEBUG, "SuperQuadricParticle::SuperQuadricParticle() finished");
45 }
46 
56  : BaseParticle(p)
57 {
58  axes_ = p.axes_;
59  eps1_ = p.eps1_;
60  eps2_ = p.eps2_;
61 }
62 
63 
65 {
66  Mdouble radius = p.getRadius();
67  axes_ = Vec3D(radius, radius, radius);
68  eps1_ = 1.0;
69  eps2_ = 1.0;
70 }
71 
77 {
78  if (getHandler() != nullptr)
79  {
81  }
82  logger(DEBUG, "SuperQuadricParticle::SuperQuadricParticleParticle() of particle % finished.", getId());
83 
84 }
85 
91 {
92  return new SuperQuadricParticle(*this);
93 }
94 
100 void SuperQuadricParticle::write(std::ostream& os) const
101 {
103  os << " axes " << axes_
104  << " exp1 " << eps1_
105  << " exp2 " << eps2_;
106 }
107 
112 std::string SuperQuadricParticle::getName() const
113 {
114  return "SuperQuadricParticle";
115 }
116 
117 
122 void SuperQuadricParticle::read(std::istream& is)
123 {
124  BaseParticle::read(is);
125  std::string dummy;
126  is >> dummy >> axes_ >> dummy >> eps1_ >> dummy >> eps2_;
127  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
128  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
129 }
130 
131 // Set and Get Functions
132 
133 void SuperQuadricParticle::setAxesAndExponents(const Mdouble& a1, const Mdouble& a2, const Mdouble& a3, const Mdouble& eps1,
134  const Mdouble& eps2)
135 {
136  eps1_ = eps1;
137  eps2_ = eps2;
138  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
139  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
140 
141  setAxes(a1,a2,a3);
142 }
143 
144 void SuperQuadricParticle::setAxesAndExponents(const Vec3D& axes, const Mdouble& eps1, const Mdouble& eps2)
145 {
146  eps1_ = eps1;
147  eps2_ = eps2;
148  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
149  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
150 
151  setAxes(axes);
152 }
153 
155 {
156  axes_ = axes;
158  if (getSpecies() != nullptr)
159  {
160  getSpecies()->computeMass(this);
161  }
162 }
163 
164 void SuperQuadricParticle::setAxes(const Mdouble& a1, const Mdouble& a2, const Mdouble& a3)
165 {
166  setAxes({a1,a2,a3});
167 }
168 
169 void SuperQuadricParticle::setExponents(const Mdouble& eps1, const Mdouble& eps2)
170 {
171  eps1_ = eps1;
172  eps2_ = eps2;
174  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
175  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
176  if (getSpecies() != nullptr)
177  {
178  getSpecies()->computeMass(this);
179  }
180 }
181 
183 {
184  return axes_;
185 }
186 
188 {
189  return eps1_;
190 }
191 
193 {
194  return eps2_;
195 }
196 
205 {
206  logger.assert(getHandler() != nullptr, "[SuperQuadricParticle::getVolume] no particle handler specified");
207 
208  return (2.0 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) * mathsFunc::beta(0.5 * eps1_ + 1.0, eps1_) *
209  mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_);
210 }
211 
217 {
218  MatrixSymmetric3D inertia;
219 
220  inertia.XX = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
221  (axes_.Y * axes_.Y * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
222  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0)
223  + 4.0 * axes_.Z * axes_.Z * mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1) *
224  mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0));
225 
226  inertia.YY = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
227  (axes_.X * axes_.X * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
228  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0)
229  + 4.0 * axes_.Z * axes_.Z * mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1) *
230  mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0));
231 
232  inertia.ZZ = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
233  (axes_.X * axes_.X + axes_.Y * axes_.Y) * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
234  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0);
235 
236  BaseParticle::setInertia(inertia);
237 }
238 
240 {
241  logger(ERROR,"This function should not be used");
242 }
243 
245 {
246 
247  if(mathsFunc::isEqual(eps2_,1,std::numeric_limits<Mdouble>::epsilon()) && mathsFunc::isEqual(eps1_,1,std::numeric_limits<Mdouble>::epsilon()))
248  {
249  BaseParticle::setRadius(std::max(std::max(axes_.Y,axes_.X),axes_.Z));
250  return;
251  }else
252  {
253 
254  const Mdouble axesX = std::max(axes_.X, axes_.Y);
255  const Mdouble axesY = std::min(axes_.X, axes_.Y);
256 
257  Mdouble alpha;
258 
259  const Mdouble eps1 = std::min(.96, eps1_);
260  const Mdouble eps2 = std::min(.96, eps2_);
261 
262  alpha = std::pow(axesY / axesX, 2.0 / (2.0 / eps2 - 2.0));
263 
264  const Mdouble help1 = std::pow(alpha, 2.0 / eps2);
265  const Mdouble gamma = std::pow(1.0 + help1, eps2 / eps1 - 1.0);
266  const Mdouble beta = std::pow(gamma * axes_.Z * axes_.Z / (axesX * axesX), 1.0 / (2.0 / eps1 - 2.0));
267  const Mdouble xTilde = std::pow(std::pow(1 + help1, eps2 / eps1) + std::pow(beta, 2.0 / eps1),
268  -eps1 / 2.0);
269  BaseParticle::setRadius(std::sqrt(mathsFunc::square(axesX * xTilde)
270  + mathsFunc::square(alpha * axesY * xTilde)
271  + mathsFunc::square(beta * axes_.Z * xTilde)));
272 
273  }
274 }
275 
287  InteractionHandler* const interactionHandler)
288 {
289  //get the normal (from P away from the contact)
290  const LabFixedCoordinates branchVector = p->getPosition() - getPosition();
291  //Get the square of the distance between particle i and particle j
292  const Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
293  auto mixedSpecies = getSpecies()->getHandler()->getMixedObject(getSpecies(),p->getSpecies());
294 
295  const Mdouble sumOfInteractionRadii = getMaxInteractionRadius()+p->getMaxInteractionRadius();
296  if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
297  {
298  if (p->isSphericalParticle())
299  {
301  pQuad->setAxes(p->getRadius(), p->getRadius(), p->getRadius());
302  BaseInteraction* contacts = getInteractionWithSuperQuad(pQuad, timeStamp, interactionHandler);
303  delete pQuad;
304  return contacts;
305  } else {
306  SuperQuadricParticle* pQuad = static_cast<SuperQuadricParticle*>(p);
307  return getInteractionWithSuperQuad(pQuad, timeStamp, interactionHandler);
308  }
309  }
310  return nullptr;
311 }
312 
324  InteractionHandler* const interactionHandler)
325 {
326  BaseInteraction* const C = interactionHandler->getInteraction(p, this, timeStamp);
327  const SmallVector<4> approximateContactPoint = getContactPoint(p, C);
328 
329  //set the new contact point:
330  LabFixedCoordinates contactPoint;
331  contactPoint.X = approximateContactPoint[0];
332  contactPoint.Y = approximateContactPoint[1];
333  contactPoint.Z = approximateContactPoint[2];
334 
335  //if the minimum is positive, there is no contact: return empty vector
336  if (computeShape(contactPoint) > -1e-10)
337  {
338  return nullptr;
339  }
340 
341  C->setContactPoint(contactPoint);
342  C->setLagrangeMultiplier(approximateContactPoint[3]);
343 
344  SmallVector<3> gradThis = computeShapeGradientLabFixed(contactPoint).getNormalised();
345  LabFixedCoordinates normal;
346  normal.X = gradThis[0];
347  normal.Y = gradThis[1];
348  normal.Z = gradThis[2];
349  C->setNormal(normal);
350 
351  const Mdouble alphaI = overlapFromContactPoint(contactPoint, normal);
352  const Mdouble alphaJ = p->overlapFromContactPoint(contactPoint, -normal);
353  C->setOverlap(alphaI + alphaJ);
355  C->setDistance((getPosition() - p->getPosition()).getLength() - C->getOverlap());
356 
357  return C;
358 }
359 
372 {
373  //Assuming contact between two spheres
374  SmallVector<4> approximateContactPoint = getInitialGuessForContact(p, C);
375  if (computeContactPoint(approximateContactPoint, this, p))
376  {
377  return approximateContactPoint;
378  }
379  else
380  {
381  return getContactPointPlanB(p, 4);
382  }
383 }
384 
393  const LabFixedCoordinates& normal) const
394 {
395  Mdouble alphaI = 0;
396  LabFixedCoordinates xEdge = contactPoint + alphaI * normal;
397  Mdouble dampingCoefficient = 1;
398  while (std::abs(computeShape(xEdge)) > 1e-10)
399  {
400  SmallVector<3> gradientShape = computeShapeGradientLabFixed(xEdge);
401  LabFixedCoordinates gradient;
402  gradient.X = gradientShape(0);
403  gradient.Y = gradientShape(1);
404  gradient.Z = gradientShape(2);
405  Mdouble newValue = alphaI - dampingCoefficient * computeShape(xEdge) / Vec3D::dot(gradient, normal);
406  const LabFixedCoordinates xEdgeNew = contactPoint + newValue * normal;
407 
408  if (std::abs(computeShape(xEdgeNew)) < std::abs(computeShape(xEdge)))
409  {
410  alphaI = newValue;
411  xEdge = xEdgeNew;
412  dampingCoefficient = 1;
413  }
414  else
415  {
416  dampingCoefficient /= 2;
417  if (dampingCoefficient < 1e-10)
418  {
419  logger(ERROR, "overlap detection does not converge");
420  }
421  }
422  }
423  return alphaI;
424 }
425 
436 {
437  SmallVector<4> approximateContactPoint;
438  if (C == nullptr || C->getOverlap() < 1e-15)
439  {
440  //this is a new interaction
441  const LabFixedCoordinates branchVector = pQuad->getPosition() - getPosition();
442  //Get the square of the distance between particle i and particle j
443  const Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
444  const Mdouble distance = sqrt(distanceSquared);
445  const LabFixedCoordinates normal = (branchVector / distance);
446  const Mdouble overlap = getSumOfInteractionRadii(pQuad) - distance;
447  const LabFixedCoordinates contactPoint = (pQuad->getPosition() -
448  (pQuad->getInteractionRadius(this) - 0.5 * overlap) * normal);
449 
450  approximateContactPoint[0] = contactPoint.X;
451  approximateContactPoint[1] = contactPoint.Y;
452  approximateContactPoint[2] = contactPoint.Z;
453  approximateContactPoint[3] = 1;
454  }
455  else
456  {
457  //this is an existing interaction
458  const LabFixedCoordinates contactPoint = C->getContactPoint();
459  const Mdouble multiplier = C->getLagrangeMultiplier();
460  approximateContactPoint[0] = contactPoint.X;
461  approximateContactPoint[1] = contactPoint.Y;
462  approximateContactPoint[2] = contactPoint.Z;
463  approximateContactPoint[3] = multiplier;
464  }
465  return approximateContactPoint;
466 }
467 
474 {
475  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
476  getOrientation().rotateBack(bodyFixedCoordinates);
477 
478  const Mdouble n1 = 2. / eps1_;
479  const Mdouble n2 = 2. / eps2_;
480 
481  return std::pow(std::pow(std::abs(bodyFixedCoordinates.X / axes_.X), n2)
482  + std::pow(std::abs(bodyFixedCoordinates.Y / axes_.Y), n2), n1 / n2)
483  + std::pow(std::abs(bodyFixedCoordinates.Z / axes_.Z), n1) - 1.0;
484 }
485 
494 {
495 
496  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
497  getOrientation().rotateBack(bodyFixedCoordinates);
498 
499  const Mdouble n1 = 2. / eps1_;
500  const Mdouble n2 = 2. / eps2_;
501 
502  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
503  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
504  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
505 
506  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
507 
508  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
509 
510  Vec3D gradientBodyFixed = {
511  (n1 / axes_.X) * std::pow(absXoa, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.X),
512  (n1 / axes_.Y) * std::pow(absYob, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.Y),
513  (n1 / axes_.Z) * std::pow(absZoc, n1 - 1.0) * mathsFunc::sign(bodyFixedCoordinates.Z)};
514  //we get the gradient in terms of the lab-fixed coordinates by \nabla_X F(X) = A \nabla_x f(x).
515  getOrientation().rotate(gradientBodyFixed);
516  return {gradientBodyFixed.X, gradientBodyFixed.Y, gradientBodyFixed.Z};
517 }
518 
527 {
528  SmallMatrix<3, 3> hessian;
529  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
530  getOrientation().rotateBack(bodyFixedCoordinates);
531 
532  const Mdouble n1 = 2.0 / eps1_;
533  const Mdouble n2 = 2.0 / eps2_;
534 
535  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
536  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
537  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
538 
539  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
540  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
541  const Mdouble help2 = std::pow(nu, n1 / n2 - 2.0);
542 
543  hessian(0, 0) = n1 * (n2 - 1) / axes_.X / axes_.X * std::pow(absXoa, n2 - 2.0) * help1
544  + ((n1 - n2) * n1 / axes_.X / axes_.X) * std::pow(absXoa, 2. * n2 - 2.0) * help2;
545  hessian(1, 1) = n1 * (n2 - 1) / axes_.Y / axes_.Y * std::pow(absYob, n2 - 2.0) * help1
546  + ((n1 - n2) * n1 / axes_.Y / axes_.Y) * std::pow(absYob, 2. * n2 - 2.0) * help2;
547  hessian(2, 2) = n1 * (n1 - 1) / axes_.Z / axes_.Z * std::pow(absZoc, n1 - 2.0);
548  hessian(1, 0) = n1 * (n1 - n2) / axes_.X / axes_.Y * std::pow(absXoa, n2 - 1.0) * std::pow(absYob, n2 - 1.0) *
549  help2 * mathsFunc::sign(bodyFixedCoordinates.X * bodyFixedCoordinates.Y);
550  hessian(0, 1) = hessian(1, 0);
553  return A * hessian * (A.transpose());
554 }
555 
567  const SuperQuadricParticle* const p1,
568  const SuperQuadricParticle* const p2) const
569 {
570  LabFixedCoordinates labFixedCoordinates;
571  labFixedCoordinates.X = position[0];
572  labFixedCoordinates.Y = position[1];
573  labFixedCoordinates.Z = position[2];
574  const Mdouble lagrangeMultiplier = position[3];
575 
576  // First compute the gradient of the superquadric shape function and then rotate it to the
577  // global frame of reference
578  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
579 
580  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
581 
582  // Eqn. (23) in Podhlozhnyuk et al, Comp. Part. Mech. 2017.
583  // where contactPoint[3] = mu
584  SmallVector<4> toBecomeZero;
585  for (unsigned int i = 0; i < 3; ++i)
586  {
587  toBecomeZero[i] = gradThis[i] + lagrangeMultiplier * lagrangeMultiplier * gradOther[i];
588  }
589  toBecomeZero[3] = p1->computeShape(labFixedCoordinates) - p2->computeShape(labFixedCoordinates);
590 
591  return toBecomeZero;
592 }
593 
624  const SuperQuadricParticle* const p1,
625  const SuperQuadricParticle* const p2) const
626 {
627  LabFixedCoordinates labFixedCoordinates;
628  labFixedCoordinates.X = contactPoint[0];
629  labFixedCoordinates.Y = contactPoint[1];
630  labFixedCoordinates.Z = contactPoint[2];
631  const Mdouble lagrangeMultiplier = contactPoint[3];
632 
633  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
634  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
635 
636  SmallMatrix<3, 3> hessianThis = p1->computeHessianLabFixed(labFixedCoordinates);
637  SmallMatrix<3, 3> hessianOther = p2->computeHessianLabFixed(labFixedCoordinates);
638 
639  SmallMatrix<3, 3> upperLeftCorner = hessianThis + lagrangeMultiplier * lagrangeMultiplier * hessianOther;
640 
641  SmallMatrix<4, 4> jacobian;
642  for (unsigned int i = 0; i < 3; ++i)
643  {
644  for (unsigned int j = 0; j < 3; ++j)
645  {
646  jacobian(i, j) = upperLeftCorner(i, j);
647  }
648  jacobian(i, 3) = 2 * lagrangeMultiplier * gradOther[i];
649  }
650 
651  for (unsigned int j = 0; j < 3; ++j)
652  {
653  jacobian(3, j) = gradThis[j] - gradOther[j];
654  }
655 
656  return jacobian;
657 }
658 
665 {
666  SmallVector<3> gradientVec = computeShapeGradientLabFixed(labFixedCoordinates);
667  SmallMatrix<3, 1> gradient = gradientVec;
668  SmallMatrix<3, 3> hessian = computeHessianLabFixed(labFixedCoordinates);
669  Mdouble helper = ((gradient.transpose() * hessian) * gradient)(0, 0);
670  Mdouble gradientLength = gradientVec.length();
671  Mdouble helper2 = gradientLength * gradientLength * (hessian(0, 0) + hessian(1, 1) + hessian(2, 2));
672  Mdouble denominator = 2 * gradientLength * gradientLength * gradientLength;
673 
675  return (helper2 - helper) / denominator;
676 }
677 
686 {
687  SmallVector<4> approximateContactPoint;
688 
689  if (p->isSphericalParticle())
690  {
692  pQuad->setAxes(p->getRadius(), p->getRadius(), p->getRadius());
693  approximateContactPoint = getContactPoint(pQuad, nullptr);
694  delete pQuad;
695  } else {
696  const SuperQuadricParticle* pQuad = static_cast<const SuperQuadricParticle*>(p);
697  approximateContactPoint = getContactPoint(pQuad, nullptr);
698  }
699 
700  //set the new contact point:
701  LabFixedCoordinates contactPoint;
702  contactPoint.X = approximateContactPoint[0];
703  contactPoint.Y = approximateContactPoint[1];
704  contactPoint.Z = approximateContactPoint[2];
705 
706  //if the minimum is positive, there is no contact: return false
707  return (computeShape(contactPoint) < 0);
708 
709 }
710 
720 SmallVector<4> SuperQuadricParticle::getContactPointPlanB(const SuperQuadricParticle* const pOther, unsigned numberOfSteps) const
721 {
722  logger(VERBOSE, "Number of iterations: %", numberOfSteps);
723  // Step 1: Compute contact point for two volume equivalent spheres
724  const Mdouble interactionRadiusThis = getInteractionRadius(pOther);
725  const Mdouble interactionRadiusOther = pOther->getInteractionRadius(this);
726  const Vec3D positionThis = getPosition();
727  const Vec3D positionOther = pOther->getPosition();
728 
729  //analytical solution for contact point between bounding spheres
730  const Vec3D initialPoint = (interactionRadiusOther * positionThis + interactionRadiusThis * positionOther) /
731  (interactionRadiusThis + interactionRadiusOther);
732 
733  SmallVector<4> contactPointPlanB;
734  // Analytical solution
735  contactPointPlanB[0] = initialPoint.X;
736  contactPointPlanB[1] = initialPoint.Y;
737  contactPointPlanB[2] = initialPoint.Z;
738  contactPointPlanB[3] = 1.0; //Need to check.
739 
740  writeDebugMessageStep1(pOther, contactPointPlanB);
741 
742  //Step 2: Compute the deltas
743  const Vec3D dAxesThis = (getAxes() - interactionRadiusThis * Vec3D(1, 1, 1)) / numberOfSteps;
744  const Mdouble dn11 = (2.0 / getExponentEps1() - 2.0) / numberOfSteps;
745  const Mdouble dn12 = (2.0 / getExponentEps2() - 2.0) / numberOfSteps;
746 
747  const Vec3D dAxesOther = (pOther->getAxes() - interactionRadiusOther * Vec3D(1, 1, 1)) / numberOfSteps;
748  const Mdouble dn21 = (2.0 / pOther->getExponentEps1() - 2.0) / numberOfSteps;
749  const Mdouble dn22 = (2.0 / pOther->getExponentEps2() - 2.0) / numberOfSteps;
750 
751  writeDebugMessageStep2(pOther, dAxesThis, dn11, dn12, dAxesOther, dn21, dn22);
752 
753  // Create two superquadrics with the above parameters for incrementally computing the contact point
756 
757  p1.setPosition(getPosition());
758  p2.setPosition(pOther->getPosition());
759 
761  p2.setOrientation(pOther->getOrientation());
762  bool success = true;
763  const unsigned maxIterations = 20;
764  for (unsigned int counter = 1; counter <= numberOfSteps; ++counter)
765  {
766 
767  // Step 3: Increment the shape and blockiness parameters
768  const Vec3D axesThis = interactionRadiusThis * Vec3D(1, 1, 1) + counter * dAxesThis;
769  const Mdouble n11 = 2.0 + counter * dn11;
770  const Mdouble n12 = 2.0 + counter * dn12;
771 
772  Vec3D axesOther = interactionRadiusOther * Vec3D(1, 1, 1) + counter * dAxesOther;
773  const Mdouble n21 = 2.0 + counter * dn21;
774  const Mdouble n22 = 2.0 + counter * dn22;
775 
776  p1.setAxesAndExponents(axesThis, 2.0 / n11, 2.0 / n12);
777  p2.setAxesAndExponents(axesOther, 2.0 / n21, 2.0 / n22);
778 
779  writeDebugMessageStep3(axesThis, n11, n12, axesOther, n21, n22);
780 
781  writeDebugMessageMiddleOfLoop(p1, p2, contactPointPlanB, counter);
782 
783  //compute the contact point of the new particles
784  success = computeContactPoint(contactPointPlanB, &p1, &p2);
785  if (!success)
786  {
787  if (numberOfSteps > maxIterations)
788  {
789  write(std::cout);
790  pOther->write(std::cout);
791  logger(ERROR, "Plan B fails even with more than 20 intermediate steps");
792  }
793  return getContactPointPlanB(pOther, 2 * numberOfSteps);
794  }
795  }
796  logger.assert(p1.getAxes().X == getAxes().X, "In getContactPointPlanB, final particle for contact-detection not "
797  "the same as original particle");
798 
799  logger(VERBOSE, "Plan B contact point: %", contactPointPlanB);
800 
801  return contactPointPlanB;
802 }
803 
817  const SuperQuadricParticle* const p2) const
818 {
819  // Damped Newton's method: (dampingFactor 1 is undamped)
820  Mdouble dampingFactor = 1;
821  int iteration = 1;
822  const Mdouble tolerance = 1e-5;
823  const unsigned maxIterations = 100;
824 
825  logger(VERBOSE, "func to be zero value: % \n",
826  computeResidualContactDetection(contactPoint, p1, p2));
827 
828  while (computeResidualContactDetection(contactPoint, p1, p2).length() > tolerance && iteration < maxIterations)
829  {
830  logger(VERBOSE, "Iteration: %\n", iteration);
831 
832  SmallMatrix<4, 4> jacobian = getJacobianOfContactDetectionObjective(contactPoint, p1, p2);
833  SmallVector<4> func = -computeResidualContactDetection(contactPoint, p1, p2);
834  jacobian.solve(func); //note: solve is in place
835 
836  SmallVector<4> newPoint = contactPoint + dampingFactor * func;
837  const auto residueOld = computeResidualContactDetection(contactPoint, p1, p2);
838  const auto residueNew = computeResidualContactDetection(newPoint, p1, p2);
839  logger(VERBOSE, "Old value: % (%) new value: % (%)",
840  computeResidualContactDetection(contactPoint, p1, p2),
841  computeResidualContactDetection(contactPoint, p1, p2).length(),
842  computeResidualContactDetection(newPoint, p1, p2),
843  computeResidualContactDetection(newPoint, p1, p2).length());
844 
845  logger(VERBOSE, "current contact point: %, new contact point: %\n", contactPoint, newPoint);
846  logger(VERBOSE, "damping factor: %", dampingFactor);
847 
848  if (residueNew.length()
849  < residueOld.length())
850  {
851  contactPoint = newPoint;
852  dampingFactor = 1;
853  }
854  else
855  {
856  dampingFactor /= 2;
857 
858  if (dampingFactor < 1e-10)
859  {
860  return false;
861  }
862  }
863 
864  iteration++;
865  }
866  return (iteration != maxIterations);
867 }
868 
870  SmallVector<4>& contactPointPlanB, const unsigned int& counter) const
871 {
872  logger(VERBOSE, "Position of particle 1 (p1): % \nPosition of particle 2 (p2): % \n",
873  p1.getPosition(), p2.getPosition());
874  logger(VERBOSE, "Orientation of particle 1 (p1): % \nOrientation of particle 2 (p2): %\n",
875  p1.getOrientation(), p2.getOrientation());
876 
877  // Step 4: Calculate the contact point for the updated shape parameters
878 
879  logger(VERBOSE, "----------------------------------------");
880  logger(VERBOSE, "STEP 4: Compute contact point");
881  logger(VERBOSE, "----------------------------------------");
882  logger(VERBOSE, " ");
883  logger(VERBOSE, "Counter: %", counter);
884  logger(VERBOSE, " ");
885 
886  logger(VERBOSE, "Analytical solution Contact Point: % ", contactPointPlanB);
887  logger(VERBOSE, " ");
888 }
889 
890 void
891 SuperQuadricParticle::writeDebugMessageStep3(const Vec3D& axesThis, const Mdouble& n11, const Mdouble& n12,
892  const Vec3D& axesOther, const Mdouble& n21, const Mdouble& n22) const
893 {
894  logger(VERBOSE, "-----------------------------------");
895  logger(VERBOSE, "STEP 3: First increment");
896  logger(VERBOSE, "-----------------------------------");
897  logger(VERBOSE, " ");
898 
899  logger(VERBOSE, "Particle 1 x-scale after increment: %", axesThis.X);
900  logger(VERBOSE, "Particle 1 y-scale after increment: %", axesThis.Y);
901  logger(VERBOSE, "Particle 1 z-scale after increment: %", axesThis.Z);
902  logger(VERBOSE, "Particle 1 n1 after increment: %", n11);
903  logger(VERBOSE, "Particle 1 n2 after increment: %", n12);
904  logger(VERBOSE, " ");
905 
906  logger(VERBOSE, "Particle 2 x-scale after increment: %", axesOther.X);
907  logger(VERBOSE, "Particle 2 y-scale after increment: %", axesOther.Y);
908  logger(VERBOSE, "Particle 2 z-scale after increment: %", axesOther.Z);
909  logger(VERBOSE, "Particle 2 n1 after increment: %", n21);
910  logger(VERBOSE, "Particle 2 n2 after increment: %", n22);
911  logger(VERBOSE, " ");
912 }
913 
914 void
916  const Mdouble& dn12, const Vec3D& dAxesOther, const Mdouble& dn21,
917  const Mdouble& dn22) const
918 {
919  logger(VERBOSE, "---------------------");
920  logger(VERBOSE, "STEP 2");
921  logger(VERBOSE, "---------------------");
922  logger(VERBOSE, " ");
923 
924  logger(VERBOSE, "Particle 1 x-scale: %", pQuad->getAxes().X);
925  logger(VERBOSE, "Particle 1 y-scale: %", pQuad->getAxes().Y);
926  logger(VERBOSE, "Particle 1 z-scale: %", pQuad->getAxes().Z);
927  logger(VERBOSE, "Particle 1 n1: %", 2. / pQuad->getExponentEps1());
928  logger(VERBOSE, "Particle 1 n2: %", 2. / pQuad->getExponentEps2());
929  logger(VERBOSE, " ");
930 
931  logger(VERBOSE, "Particle 1 x-scale increment: %", dAxesThis.X);
932  logger(VERBOSE, "Particle 1 y-scale increment: %", dAxesThis.Y);
933  logger(VERBOSE, "Particle 1 z-scale increment: %", dAxesThis.Z);
934  logger(VERBOSE, "Particle 1 n1 increment: %", dn11);
935  logger(VERBOSE, "Particle 1 n2 increment: %", dn12);
936  logger(VERBOSE, " ");
937 
938  logger(VERBOSE, "Particle 2 x-scale: %", getAxes().X);
939  logger(VERBOSE, "Particle 2 y-scale: %", getAxes().Y);
940  logger(VERBOSE, "Particle 2 z-scale: %", getAxes().Z);
941  logger(VERBOSE, "Particle 2 n1: %", 2. / getExponentEps1());
942  logger(VERBOSE, "Particle 2 n2: %", 2. / getExponentEps2());
943  logger(VERBOSE, " ");
944 
945  logger(VERBOSE, "Particle 2 x-scale increment: %", dAxesOther.X);
946  logger(VERBOSE, "Particle 2 y-scale increment: %", dAxesOther.Y);
947  logger(VERBOSE, "Particle 2 z-scale increment: %", dAxesOther.Z);
948  logger(VERBOSE, "Particle 2 n1 increment: %", dn21);
949  logger(VERBOSE, "Particle 2 n2 increment: %", dn22);
950  logger(VERBOSE, " ");
951 }
952 
953 void SuperQuadricParticle::writeDebugMessageStep1(const SuperQuadricParticle* pQuad, const SmallVector<4>& contactPointPlanB) const
954 {
955  logger(VERBOSE, "---------------------");
956  logger(VERBOSE, "STEP 1");
957  logger(VERBOSE, "---------------------\n");
958 
959  logger(VERBOSE, "Position of particle 1: % ", getPosition());
960  logger(VERBOSE, "Position of particle 2 (pQuad): % ", pQuad->getPosition());
961  logger(VERBOSE, " ");
962 
963  logger(VERBOSE, "Radius of particle 1: % ", getInteractionRadius(pQuad));
964  logger(VERBOSE, "Radius of particle 2 (pQuad): % \n", pQuad->getInteractionRadius(this));
965 
966  logger(VERBOSE, "Orientation of particle 1: % ", getOrientation());
967  logger(VERBOSE, "Orientation of particle 2 (pQuad): % ", pQuad->getOrientation());
968  logger(VERBOSE, " ");
969 
970  logger(VERBOSE, "Particle 1 axes: %", getAxes());
971  logger(VERBOSE, "Particle 2 axes (pQuad): % \n", pQuad->getAxes());
972 
973  logger(VERBOSE, "Analytical solution for two equivalent spheres in contact: % \n", contactPointPlanB);
974 }
975 
977 {
978  const auto mixedSpecies = getSpecies()->getHandler()->getMixedObject(getSpecies(),particle->getSpecies());
979  return getRadius() + 0.5 * mixedSpecies->getInteractionDistance();
980 }
981 
984  if (isFixed()) return;
985  if (getParticleDimensions()==3) {
986  Mdouble volume = getVolume();
987 
988  Mdouble help1 = mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_);
989  Mdouble help2 = mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0);
990  Mdouble help3 = mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1);
991  Mdouble help4 = mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0);
992 
993  invMass_ = 1.0 / (volume * s.getDensity());
994  invInertia_.XX = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
995  (axes_.Y * axes_.Y * help1 * help2
996  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
997  invInertia_.XY = 0.0;
998  invInertia_.XZ = 0.0;
999  invInertia_.YY = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
1000  (axes_.X * axes_.X * help1 * help2
1001  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
1002  invInertia_.YZ = 0.0;
1003  invInertia_.ZZ = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
1004  (axes_.X * axes_.X + axes_.Y * axes_.Y) * help1 * help2);
1005  } else {
1006  logger(ERROR, "[SuperQuadricParticle::computeMass] SuperQuadricParticle cannot be two-dimensional (yet)");
1007  }
1008 };
BaseInteraction * getInteractionWith(BaseParticle *P, unsigned timeStamp, InteractionHandler *interactionHandler) override
Checks if this superquadric is in interaction with the given particle, and if so, returns vector of p...
Mdouble getCurvature(const LabFixedCoordinates &labFixedCoordinates) const override
Get the mean curvature of this superquadric at the given (lab-fixed) position, see Podlozhyuk et al...
void read(std::istream &is) override
Particle read function, which accepts an std::istream as input.
void computeMass(const ParticleSpecies &s) override
Computes the particle's (inverse) mass and inertia.
SuperQuadricParticle * copy() const override
Copy method. It calls to copy constructor of this superquadric, useful for polymorphism.
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
std::string getName() const override
Returns the name of the class, here "SuperQuadricParticle".
bool computeContactPoint(SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Perform the actual Newton-iterations to find the contact point. Note, that it is given back as a para...
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
Mdouble getInteractionRadius(const BaseParticle *particle) const
returns the radius plus half the interactionDistance of the mixed species
void checkExtremaOnDelete(BaseParticle *P)
Checks if the extrema of this ParticleHandler needs updating when a particle is deleted.
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
void setBoundingRadius()
Get the radius of the sphere that fits precisely around the particle.
void read(std::istream &is) override
Read function: read in the information for this superquadric from the given input-stream, for example a restart file.
virtual void setInertia()
virtual bool isSphericalParticle() const
Definition: BaseParticle.h:642
void setRadius(const Mdouble radius) override
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
void setOverlap(Mdouble overlap)
Set the overlap between the two interacting object.
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
SuperQuadricParticle()
Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2...
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
Definition: BaseParticle.h:359
BaseInteraction * getInteractionWithSuperQuad(SuperQuadricParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler)
Checks if this superquadric is in interaction with the given superquadric, and if so...
SmallMatrix< numberOfColumns, numberOfRows > transpose() const
Definition: SmallMatrix.h:330
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
bool isInContactWith(const BaseParticle *p) const override
Get whether or not this superquadric is in contact with the given particle.
~SuperQuadricParticle() override
Destructor, needs to be implemented and checked to see if it is the largest or smallest particle curr...
Mdouble invMass_
Particle radius_.
Definition: BaseParticle.h:653
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:376
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
void setContactPoint(Vec3D contactPoint)
Set the location of the contact point between the two interacting objects.
Mdouble length() const
Definition: SmallVector.h:253
void computeMass(BaseParticle *p) const
Compute Particle mass function, which required a reference to the Species vector. It computes the Par...
MatrixSymmetric3D invInertia_
Inverse Particle mass (for computation optimization)
Definition: BaseParticle.h:654
void setLagrangeMultiplier(Mdouble multiplier)
const Vec3D & getContactPoint() const
Gets constant reference to contact point (vector).
int sign(T val)
This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0...
Definition: ExtendedMath.h:95
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
ParticleHandler * getHandler() const
Returns pointer to the particle's ParticleHandler.
Vec3D getAxes() const override
Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 ...
void setDistance(Mdouble distance)
Sets the interaction distance between the two interacting objects.
SmallVector< 4 > getContactPoint(const SuperQuadricParticle *p, BaseInteraction *C) const
Compute the contact point between this and the given superquadric particle.
bool isEqual(Mdouble v1, Mdouble v2, Mdouble absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
SmallVector< 4 > computeResidualContactDetection(const SmallVector< 4 > &position, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Objective function for contact detection between the two given superquadrics. See Podlozhyuk et al...
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 ...
void setExponents(const Mdouble &eps1, const Mdouble &eps2) override
Set the exponents to eps1 and eps2 for this superquadric. We use the super-ellipsoid definition state...
SpeciesHandler * getHandler() const
Returns the pointer to the handler to which this species belongs.
Definition: BaseSpecies.cc:99
SmallVector< 4 > getContactPointPlanB(const SuperQuadricParticle *pOther, unsigned numberOfSteps) const
If the "normal" procedure fails to find a contact point, use an alternative approach that involves st...
Mdouble getVolume() const override
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
BaseInteraction * getInteraction(BaseInteractable *P, BaseInteractable *I, unsigned timeStamp)
Returns the Interaction between the BaseInteractable's P and I.
SmallVector< 4 > getInitialGuessForContact(const SuperQuadricParticle *pQuad, BaseInteraction *C) const
Get an initial guess for the contact-point between this particle and the given particle.
void write(std::ostream &os) const override
Particle print function, which accepts an std::ostream as input.
Container to store Interaction objects.
Mdouble eps1_
Blockiness parameters.
void solve(SmallMatrix< numberOfRows, numberOfRightHandSideColumns > &B) const
solves Ax=B where A is the current matrix and B is passed in. The result is returned in B...
Mdouble getLagrangeMultiplier()
Mdouble computeShape(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the shape-functiion at the given (lab-fixed) position.
void setInertia() override
Compute and set the inertia-tensor for this superquadric. For internal use only.
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
bool isFixed() const override
Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Ma...
Definition: BaseParticle.h:93
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
Mdouble getExponentEps1() const override
Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter ...
void writeDebugMessageStep1(const SuperQuadricParticle *pQuad, const SmallVector< 4 > &contactPointPlanB) const
SmallMatrix< 4, 4 > getJacobianOfContactDetectionObjective(const SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Compute and return the derivative of functionThatShouldBecomeZeroForContactDetection, both to the position and the Lagrange multiplier, and evaluated at the contact point.
void setAxes(const Mdouble &a1, const Mdouble &a2, const Mdouble &a3)
Set the axes-lengths to a1, a2 and a3 for this superquadric. We use the super-ellipsoid definition st...
Mdouble Y
Definition: Vector.h:65
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
void getRotationMatrix(SmallMatrix< 3, 3 > &A) const
Retrieves the rotation matrix.
Definition: Quaternion.cc:508
void writeDebugMessageStep3(const Vec3D &axesThis, const Mdouble &n11, const Mdouble &n12, const Vec3D &axesOther, const Mdouble &n21, const Mdouble &n22) const
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
void setAxesAndExponents(const Mdouble &a1, const Mdouble &a2, const Mdouble &a3, const Mdouble &eps1, const Mdouble &eps2)
Set the geometrical properties of the superquadrics, namely the axes-lengths a1, a2 and a3...
Mdouble getDensity() const
Allows density_ to be accessed.
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
Definition: Vector.h:49
Mdouble overlapFromContactPoint(const LabFixedCoordinates &contactPoint, const LabFixedCoordinates &normal) const
Compute the distance between the contact-point and surface of this superquadric particle.
Data type for small dense matrix.
Definition: SmallMatrix.h:67
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
unsigned int getParticleDimensions() const
Returns the particle's dimensions (either 2 or 3).
Mdouble getExponentEps2() const override
Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter...
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
void writeDebugMessageStep2(const SuperQuadricParticle *pQuad, const Vec3D &dAxesThis, const Mdouble &dn11, const Mdouble &dn12, const Vec3D &dAxesOther, const Mdouble &dn21, const Mdouble &dn22) const
void writeDebugMessageMiddleOfLoop(const SuperQuadricParticle &p1, const SuperQuadricParticle &p2, SmallVector< 4 > &contactPointPlanB, const unsigned int &counter) const
Implementation of a 3D symmetric matrix.
void write(std::ostream &os) const override
Write function: write this superquadric to the given output-stream, for example a restart-file...
Mdouble XX
The six distinctive matrix elements.
SmallMatrix< 3, 3 > computeHessianLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the hessian ("second derivative") of the shape-function at the given (lab-fixed) posi...