v0.13.0
Public Member Functions | Public Attributes | List of all members
CylinderRigidBody Struct Reference

#include <users_modules/multifield_plasticity/src/RigidBodies.hpp>

Inheritance diagram for CylinderRigidBody:
[legend]
Collaboration diagram for CylinderRigidBody:
[legend]

Public Member Functions

 CylinderRigidBody (VectorDouble c_coords, VectorDouble roller_disp, int id)
 
MoFEMErrorCode getRollerDataForTag ()
 
Tensor1< double, 3 > getNormal (Tensor1< TPack3, 3 > &t_coords, Tensor1< TPack1, 3 > &t_disp)
 
Tensor1< double, 3 > getNormal (Tensor1< TPack3, 3 > &t_coords, Tensor1< double, 3 > &t_disp)
 
Tensor1< double, 3 > getNormal (Tensor1< double, 3 > &t_coords, Tensor1< double, 3 > &t_disp)
 
Tensor2< double, 3, 3 > getDiffNormal (Tensor1< TPack3, 3 > &t_coords, Tensor1< TPack1, 3 > &t_disp, Tensor1< TPack1, 3 > &t_normal)
 
double getGap (Tensor1< TPack3, 3 > &t_coords)
 
Tensor1< double, 3 > getdGap (Tensor1< TPack3, 3 > &t_coords, Tensor1< TPack1, 3 > &t_normal)
 
template<typename T1 , typename T2 >
Tensor1< double, 3 > & getNormalImpl (Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp)
 
template<typename T1 , typename T2 , typename T3 >
Tensor2< double, 3, 3 > getDiffNormalImpl (Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp, Tensor1< T3, 3 > &t_normal)
 
template<typename T1 >
double getGapImpl (Tensor1< T1, 3 > &t_coords)
 
template<typename T1 , typename T2 >
Tensor1< double, 3 > getdGapImpl (Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_normal)
 
MoFEMErrorCode getBodyOptions ()
 
- Public Member Functions inherited from RigidBodyData
 RigidBodyData (VectorDouble c_coords, VectorDouble roller_disp, int id)
 
 RigidBodyData ()=delete
 
virtual ~RigidBodyData ()
 
MoFEMErrorCode computeRotationMatrix ()
 
Tensor1< double, 3 > getBodyOffset ()
 
template<typename T1 , typename T2 >
Tensor1< double, 3 > & getPointCoords (Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp)
 
MoFEMErrorCode saveBasicDataOnTag (moab::Interface &moab_debug, EntityHandle &vertex)
 

Public Attributes

double rAdius
 
double hEight
 
- Public Attributes inherited from RigidBodyData
const int iD
 
Index< 'i', 3 > i
 
Index< 'j', 3 > j
 
Index< 'k', 3 > k
 
VectorDouble originCoords
 
VectorDouble rollerDisp
 
VectorDouble BodyDispScaled
 
double gAp
 
Tensor1< double, 3 > tNormal
 
Tensor1< double, 3 > dGap
 
Tensor1< double, 3 > pointCoords
 
Tensor1< double, 3 > closestPoint
 
Tensor1< double, 3 > defaultOrientation
 
Tensor1< double, 3 > oRientation
 
Tensor2< double, 3, 3 > rotationMat
 
Tensor2< double, 3, 3 > diffNormal
 
string positionDataParamName
 
int bodyType
 
array< double, 9 > dataForTags
 

Additional Inherited Members

- Public Types inherited from RigidBodyData
using TPack3 = PackPtr< double *, 3 >
 
using TPack1 = PackPtr< double *, 1 >
 

Detailed Description

Definition at line 353 of file RigidBodies.hpp.

Constructor & Destructor Documentation

◆ CylinderRigidBody()

CylinderRigidBody::CylinderRigidBody ( VectorDouble  c_coords,
VectorDouble  roller_disp,
int  id 
)

Definition at line 356 of file RigidBodies.hpp.

357  : RigidBodyData(c_coords, roller_disp, id), rAdius(1), hEight(INFINITY) {}
RigidBodyData()=delete

Member Function Documentation

◆ getBodyOptions()

MoFEMErrorCode CylinderRigidBody::getBodyOptions ( )
virtual

Implements RigidBodyData.

Definition at line 540 of file RigidBodies.hpp.

540  {
542  PetscBool rflg;
543  PetscInt nb_dirs = 3;
544  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
545  string param = "-direction" + to_string(iD);
546  CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(), &oRientation(0),
547  &nb_dirs, &rflg);
548  param = "-radius" + to_string(iD);
549  CHKERR PetscOptionsScalar(param.c_str(), "set roller radius", "", rAdius,
550  &rAdius, PETSC_NULL);
551  param = "-height" + to_string(iD);
552  CHKERR PetscOptionsScalar(param.c_str(), "set cylinder height", "", hEight,
553  &hEight, PETSC_NULL);
554  if (rflg)
556  ierr = PetscOptionsEnd();
557  CHKERRG(ierr);
559  }
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Tensor1< double, 3 > oRientation
Definition: RigidBodies.hpp:43
MoFEMErrorCode computeRotationMatrix()
Definition: RigidBodies.hpp:85
const int iD
Definition: RigidBodies.hpp:30

◆ getdGap()

Tensor1<double, 3> CylinderRigidBody::getdGap ( Tensor1< TPack3, 3 > &  t_coords,
Tensor1< TPack1, 3 > &  t_normal 
)
virtual

Implements RigidBodyData.

Definition at line 390 of file RigidBodies.hpp.

391  {
392  return getdGapImpl(t_coords, t_normal);
393  }
Tensor1< double, 3 > getdGapImpl(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_normal)

◆ getdGapImpl()

template<typename T1 , typename T2 >
Tensor1<double, 3> CylinderRigidBody::getdGapImpl ( Tensor1< T1, 3 > &  t_coords,
Tensor1< T2, 3 > &  t_normal 
)

Definition at line 531 of file RigidBodies.hpp.

532  {
533  // dGap(i) =
534  // pointCoords(j) * diffNormal(i, j) +
535  // tNormal(j) * (kronecker_delta(i, j) - 2. * rAdius * diffNormal(i,
536  // j));
537  return dGap;
538  };
Tensor1< double, 3 > dGap
Definition: RigidBodies.hpp:39

◆ getDiffNormal()

Tensor2<double, 3, 3> CylinderRigidBody::getDiffNormal ( Tensor1< TPack3, 3 > &  t_coords,
Tensor1< TPack1, 3 > &  t_disp,
Tensor1< TPack1, 3 > &  t_normal 
)
virtual

Implements RigidBodyData.

Definition at line 380 of file RigidBodies.hpp.

382  {
383  return getDiffNormalImpl(t_coords, t_disp, t_normal);
384  }
Tensor2< double, 3, 3 > getDiffNormalImpl(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp, Tensor1< T3, 3 > &t_normal)

◆ getDiffNormalImpl()

template<typename T1 , typename T2 , typename T3 >
Tensor2<double, 3, 3> CylinderRigidBody::getDiffNormalImpl ( Tensor1< T1, 3 > &  t_coords,
Tensor1< T2, 3 > &  t_disp,
Tensor1< T3, 3 > &  t_normal 
)

Definition at line 503 of file RigidBodies.hpp.

505  {
506 
507  Tensor2<double, 3, 3> diff_normal;
508  Tensor1<double, 3> norm_diff;
509  const double clean_gap = gAp;
510  const double eps = 1e-8;
511  for (int ii = 0; ii != 3; ++ii) {
512  Tensor1<double, 3> pert_disp{t_disp(0), t_disp(1), t_disp(2)};
513  pert_disp(ii) += eps;
514  auto pert_normal = getNormal(t_coords, pert_disp);
515  norm_diff(i) = (pert_normal(i) - t_normal(i)) / eps;
516  dGap(ii) = (gAp - clean_gap) / eps;
517 
518  for (int jj = 0; jj != 3; ++jj) {
519  diff_normal(jj, ii) = norm_diff(jj);
520  }
521  }
522  return diff_normal;
523  };
static const double eps
Tensor1< double, 3 > getNormal(Tensor1< TPack3, 3 > &t_coords, Tensor1< TPack1, 3 > &t_disp)
Index< 'i', 3 > i
Definition: RigidBodies.hpp:31

◆ getGap()

double CylinderRigidBody::getGap ( Tensor1< TPack3, 3 > &  t_coords)
virtual

Implements RigidBodyData.

Definition at line 386 of file RigidBodies.hpp.

386  {
387  return getGapImpl(t_coords);
388  }
double getGapImpl(Tensor1< T1, 3 > &t_coords)

◆ getGapImpl()

template<typename T1 >
double CylinderRigidBody::getGapImpl ( Tensor1< T1, 3 > &  t_coords)

Definition at line 525 of file RigidBodies.hpp.

525  {
526  // gAp = tNormal(i) * (pointCoords(i) - tNormal(i) * rAdius);
527  return gAp;
528  };

◆ getNormal() [1/3]

Tensor1<double, 3> CylinderRigidBody::getNormal ( Tensor1< double, 3 > &  t_coords,
Tensor1< double, 3 > &  t_disp 
)
virtual

Implements RigidBodyData.

Definition at line 375 of file RigidBodies.hpp.

376  {
377  return getNormalImpl(t_coords, t_disp);
378  }
Tensor1< double, 3 > & getNormalImpl(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp)

◆ getNormal() [2/3]

Tensor1<double, 3> CylinderRigidBody::getNormal ( Tensor1< TPack3, 3 > &  t_coords,
Tensor1< double, 3 > &  t_disp 
)
virtual

Implements RigidBodyData.

Definition at line 371 of file RigidBodies.hpp.

372  {
373  return getNormalImpl(t_coords, t_disp);
374  }

◆ getNormal() [3/3]

Tensor1<double, 3> CylinderRigidBody::getNormal ( Tensor1< TPack3, 3 > &  t_coords,
Tensor1< TPack1, 3 > &  t_disp 
)
virtual

Implements RigidBodyData.

Definition at line 367 of file RigidBodies.hpp.

368  {
369  return getNormalImpl(t_coords, t_disp);
370  }

◆ getNormalImpl()

template<typename T1 , typename T2 >
Tensor1<double, 3>& CylinderRigidBody::getNormalImpl ( Tensor1< T1, 3 > &  t_coords,
Tensor1< T2, 3 > &  t_disp 
)

Definition at line 396 of file RigidBodies.hpp.

397  {
398  auto t_d = getPointCoords(t_coords, t_disp);
399  pointCoords(i) = rotationMat(i, j) * t_d(j);
400  Tensor1<double, 3> c_center(0., 0., 0.);
401  Tensor1<double, 3> q_minus_c(0., 0., 0.);
402  Tensor1<double, 3> K(0., 0., 0.);
403  auto Q = pointCoords;
404  auto P_Q = pointCoords;
405  auto Q_K = pointCoords;
406 
407  // https://www.geometrictools.com/Documentation/DistanceToCircle3.pdf
408  const double half_height = hEight / 2.;
409  const double rad2 = rAdius * rAdius;
410  tNormal(i) = 0.;
411 
412  if (pointCoords(2) > half_height) {
413  tNormal(2) = 1;
414  if (pow(pointCoords(1), 2) + pow(pointCoords(0), 2) < rad2) {
415  gAp = tNormal(i) * pointCoords(i) - half_height;
416  } else {
417  c_center(2) = half_height;
418  Q(2) = c_center(2); // Q
419  q_minus_c(i) = Q(i) - c_center(i);
420  q_minus_c.normalize();
421  K(i) = c_center(i) + rAdius * q_minus_c(i);
422  P_Q(i) = pointCoords(i) - Q(i);
423  Q_K(i) = Q(i) - K(i);
424  gAp = sqrt(P_Q.l2() * P_Q.l2() + Q_K.l2() * Q_K.l2());
425  }
426  } else if (pointCoords(2) < -half_height) {
427  tNormal(2) = -1;
428  if (pow(pointCoords(1), 2) + pow(pointCoords(0), 2) < rad2) {
429  gAp = tNormal(i) * pointCoords(i) - half_height;
430  } else {
431  c_center(2) = -half_height;
432  Q(2) = c_center(2); // Q
433  q_minus_c(i) = Q(i) - c_center(i);
434  q_minus_c.normalize();
435  K(i) = c_center(i) + rAdius * q_minus_c(i);
436  P_Q(i) = pointCoords(i) - Q(i);
437  Q_K(i) = Q(i) - K(i);
438  gAp = sqrt(P_Q.l2() * P_Q.l2() + Q_K.l2() * Q_K.l2());
439  }
440  } else {
441  pointCoords(2) = 0;
442  tNormal(i) = pointCoords(i);
443  tNormal.normalize();
444  gAp = tNormal(i) * (pointCoords(i) - tNormal(i) * rAdius);
445  }
446 
447  // auto length = [&](auto v) { return sqrt(v(0) * v(0) + v(1) * v(1)); };
448 
449  // auto sd_cylinder = [&](Tensor1<double, 3> p) {
450  // typedef Tensor1<double, 2> vec2;
451  // Index<'i', 2> i;
452  // vec2 d(length(p) - rAdius, abs(p(2)) - half_height);
453  // vec2 mx_d(max(d(0), 0.0), max(d(1), 0.0));
454  // return min(max(d(0), d(1)), 0.0) + length(mx_d);
455  // };
456 
457  // const double eps = 1e-6;
458  // auto cyl_distance = sd_cylinder(pointCoords);
459  // gAp = cyl_distance;
460 
461  // auto get_d = [&](int a) {
462  // auto pt = pointCoords;
463  // pt(a) += eps;
464  // return (sd_cylinder(pt) - cyl_distance) / eps;
465  // };
466 
467  // for (int dd = 0; dd != 3; ++dd)
468  // tNormal(dd) = get_d(dd);
469 
470  // for debugging
471  // auto my_rand_mn = [](double M, double N) {
472  // return M + (rand() / (RAND_MAX / (N - M)));
473  // };
474  // auto t_off = getBodyOffset();
475  // std::ofstream afile("cylinder_pts.csv", std::ios::ate);
476  // if (afile.is_open()) {
477  // afile << "X,Y,Z\n";
478  // Tensor1<double, 3> p(0, 0, 0);
479  // for (int n = 0; n != 1e6; n++) {
480  // double x = my_rand_mn(-300, 300);
481  // double y = my_rand_mn(-300, 300);
482  // double z = my_rand_mn(-300, 300);
483  // Tensor1<double, 3> coords(x, y, z);
484  // if (sd_cylinder(coords) < 0) {
485 
486  // afile << x + t_off(0) << "," << y + t_off(1) << "," << z + t_off(2)
487  // << "\n";
488  // }
489  // }
490  // afile.close();
491  // }
492  // cout << "kill process " << endl;
493 
494  tNormal.normalize();
495 
496  auto normal_copy = tNormal;
497  tNormal(i) = rotationMat(j, i) * normal_copy(j);
498 
499  return tNormal;
500  };
static Index< 'K', 3 > K
Tensor1< double, 3 > pointCoords
Definition: RigidBodies.hpp:40
Tensor1< double, 3 > tNormal
Definition: RigidBodies.hpp:38
Tensor1< double, 3 > & getPointCoords(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp)
Tensor2< double, 3, 3 > rotationMat
Definition: RigidBodies.hpp:44
Index< 'j', 3 > j
Definition: RigidBodies.hpp:32

◆ getRollerDataForTag()

MoFEMErrorCode CylinderRigidBody::getRollerDataForTag ( )
virtual

Implements RigidBodyData.

Definition at line 359 of file RigidBodies.hpp.

359  {
361  bodyType = CYLINDER;
362  dataForTags[0] = rAdius;
363  dataForTags[2] = hEight;
365  }
@ CYLINDER
Definition: RigidBodies.hpp:20
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
array< double, 9 > dataForTags
Definition: RigidBodies.hpp:50

Member Data Documentation

◆ hEight

double CylinderRigidBody::hEight

Definition at line 355 of file RigidBodies.hpp.

◆ rAdius

double CylinderRigidBody::rAdius

Definition at line 354 of file RigidBodies.hpp.


The documentation for this struct was generated from the following file: