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

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

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

Public Member Functions

 ConeRigidBody (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 >
double getGapImpl (Tensor1< T1, 3 > &t_coords)
 
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 , typename T2 >
Tensor1< double, 3 > getdGapImpl (Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_normal)
 
virtual 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

shared_ptr< ArbitrarySurfaceDatacP
 
double sLope
 
double hEight
 
double rAdius
 
double smallRadius
 
double oFfset
 
double aNgle
 
- 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 956 of file RigidBodies.hpp.

Constructor & Destructor Documentation

◆ ConeRigidBody()

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

Definition at line 965 of file RigidBodies.hpp.

966  : RigidBodyData(c_coords, roller_disp, id) {
967  // , cP(R, r, 1)
968  sLope = 1;
969  hEight = 1;
970  rAdius = 1;
971  oFfset = 0;
972  // compute rotation matrix here
973  }
RigidBodyData()=delete

Member Function Documentation

◆ getBodyOptions()

virtual MoFEMErrorCode ConeRigidBody::getBodyOptions ( )
virtual

Implements RigidBodyData.

Definition at line 1174 of file RigidBodies.hpp.

1174  {
1176  PetscBool rflg;
1177  int nb_dirs = 3;
1178  double angle = 45;
1179  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
1180  string param = "-radius" + to_string(iD);
1181  CHKERR PetscOptionsScalar(param.c_str(), "set cone radius", "", rAdius,
1182  &rAdius, PETSC_NULL);
1183  param = "-height" + to_string(iD);
1184  CHKERR PetscOptionsScalar(param.c_str(), "set cone height", "", hEight,
1185  &hEight, PETSC_NULL);
1186  param = "-angle" + to_string(iD);
1187  CHKERR PetscOptionsScalar(param.c_str(), "set cone half angle", "", angle,
1188  &angle, PETSC_NULL);
1189  param = "-direction" + to_string(iD);
1190  CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(), &oRientation(0),
1191  &nb_dirs, &rflg);
1192  ierr = PetscOptionsEnd();
1193  CHKERRG(ierr);
1194  aNgle = angle;
1195  sLope = std::tan(aNgle * M_PI / 180);
1196  oFfset = (rAdius / sLope) - hEight;
1198 
1199  if (rflg)
1201 
1202  cP = make_shared<ArbitrarySurfaceData>(1, sLope);
1203 
1205  }
#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)
shared_ptr< ArbitrarySurfaceData > cP
double smallRadius
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> ConeRigidBody::getdGap ( Tensor1< TPack3, 3 > &  t_coords,
Tensor1< TPack1, 3 > &  t_normal 
)
virtual

Implements RigidBodyData.

Definition at line 1006 of file RigidBodies.hpp.

1007  {
1008  return getdGapImpl(t_coords, t_normal);
1009  }
Tensor1< double, 3 > getdGapImpl(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_normal)

◆ getdGapImpl()

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

Definition at line 1169 of file RigidBodies.hpp.

1170  {
1171  return dGap;
1172  };
Tensor1< double, 3 > dGap
Definition: RigidBodies.hpp:39

◆ getDiffNormal()

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

Implements RigidBodyData.

Definition at line 997 of file RigidBodies.hpp.

999  {
1000  return getDiffNormalImpl(t_coords, t_disp, t_normal);
1001  }
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> ConeRigidBody::getDiffNormalImpl ( Tensor1< T1, 3 > &  t_coords,
Tensor1< T2, 3 > &  t_disp,
Tensor1< T3, 3 > &  t_normal 
)

Definition at line 1147 of file RigidBodies.hpp.

1149  {
1150  Tensor2<double, 3, 3> diff_normal;
1151  Tensor1<double, 3> norm_diff;
1152  const double clean_gap = gAp;
1153  const double eps = 1e-8;
1154  for (int ii = 0; ii != 3; ++ii) {
1155  Tensor1<double, 3> pert_disp{t_disp(0), t_disp(1), t_disp(2)};
1156  pert_disp(ii) += eps;
1157  auto pert_normal = getNormal(t_coords, pert_disp);
1158  norm_diff(i) = (pert_normal(i) - t_normal(i)) / eps;
1159  dGap(ii) = (gAp - clean_gap) / eps;
1160 
1161  for (int jj = 0; jj != 3; ++jj) {
1162  diff_normal(jj, ii) = norm_diff(jj);
1163  }
1164  }
1165  return diff_normal;
1166  };
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 ConeRigidBody::getGap ( Tensor1< TPack3, 3 > &  t_coords)
virtual

Implements RigidBodyData.

Definition at line 1003 of file RigidBodies.hpp.

1003  {
1004  return getGapImpl(t_coords);
1005  }
double getGapImpl(Tensor1< T1, 3 > &t_coords)

◆ getGapImpl()

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

Definition at line 1142 of file RigidBodies.hpp.

1142  {
1143  return gAp;
1144  };

◆ getNormal() [1/3]

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

Implements RigidBodyData.

Definition at line 992 of file RigidBodies.hpp.

993  {
994  return getNormalImpl(t_coords, t_disp);
995  }
Tensor1< double, 3 > & getNormalImpl(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp)

◆ getNormal() [2/3]

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

Implements RigidBodyData.

Definition at line 988 of file RigidBodies.hpp.

989  {
990  return getNormalImpl(t_coords, t_disp);
991  }

◆ getNormal() [3/3]

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

Implements RigidBodyData.

Definition at line 984 of file RigidBodies.hpp.

985  {
986  return getNormalImpl(t_coords, t_disp);
987  }

◆ getNormalImpl()

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

Definition at line 1012 of file RigidBodies.hpp.

1013  {
1014 
1015  auto t_d = getPointCoords(t_coords, t_disp);
1016  pointCoords(i) = rotationMat(i, j) * t_d(j);
1017  t_d(i) = pointCoords(i);
1018 
1019  // Tensor1<double, 3> c_center(0., 0., 0.);
1020  // Tensor1<double, 3> q_minus_c(0., 0., 0.);
1021  // Tensor1<double, 3> K(0., 0., 0.);
1022  // auto Q = pointCoords;
1023  // auto P_Q = pointCoords;
1024  // auto Q_K = pointCoords;
1025 
1026  const double half_height = hEight / 2.;
1027  // const double rad2 = rAdius * rAdius;
1028  // const double small_rad2 = smallRadius * smallRadius;
1029 
1030  // tNormal(i) = 0.;
1031  // pointCoords(2) += (oFfset + half_height);
1032  // t_d(2) += oFfset + half_height;
1033  // if (t_d(2) > half_height) {
1034  // tNormal(2) = 1;
1035  // if (pow(t_d(1), 2) + pow(t_d(0), 2) < rad2) {
1036  // gAp = tNormal(i) * t_d(i) - half_height;
1037  // } else {
1038  // c_center(2) = half_height;
1039  // Q(2) = c_center(2); // Q
1040  // q_minus_c(i) = Q(i) - c_center(i);
1041  // q_minus_c.normalize();
1042  // K(i) = c_center(i) + rAdius * q_minus_c(i);
1043  // P_Q(i) = t_d(i) - Q(i);
1044  // Q_K(i) = Q(i) - K(i);
1045  // gAp = sqrt(P_Q.l2() * P_Q.l2() + Q_K.l2() * Q_K.l2());
1046  // }
1047  // } else if (t_d(2) < -half_height) {
1048  // tNormal(2) = -1;
1049  // if (pow(t_d(1), 2) + pow(t_d(0), 2) < small_rad2) {
1050  // gAp = tNormal(i) * t_d(i) - half_height;
1051  // } else {
1052  // c_center(2) = -half_height;
1053  // Q(2) = c_center(2); // Q
1054  // q_minus_c(i) = Q(i) - c_center(i);
1055  // q_minus_c.normalize();
1056  // K(i) = c_center(i) + smallRadius * q_minus_c(i);
1057  // P_Q(i) = t_d(i) - Q(i);
1058  // Q_K(i) = Q(i) - K(i);
1059  // gAp = sqrt(P_Q.l2() * P_Q.l2() + Q_K.l2() * Q_K.l2());
1060  // }
1061  // } else {
1062  // auto closest_pt_cone = cP->getClosestPointToCone(t_d);
1063  // closestPoint(i) = closest_pt_cone(i);
1064 
1065  // tNormal(j) = t_d(j) - closestPoint(j);
1066  // auto test_normal = tNormal.l2();
1067  // if (cP->isPointInsideCone(t_d))
1068  // tNormal(i) *= -1;
1069  // tNormal.normalize();
1070 
1071  // gAp = tNormal(i) * (t_d(i) - closestPoint(i));
1072  // }
1073 
1074  auto clamp = [](auto x, auto min, auto max) {
1075  if (x < min)
1076  x = min;
1077  else if (x > max)
1078  x = max;
1079  return x;
1080  };
1081 
1082  auto sd_capped_cone = [&](Tensor1<double, 3> p, double h, double r1,
1083  double r2) {
1084  typedef Tensor1<double, 2> vec2;
1085  Index<'i', 2> i;
1086  vec2 cb;
1087  vec2 q(sqrt(p(i) * p(i)), p(2));
1088  vec2 k1(r2, h);
1089  vec2 k2(r2 - r1, 2.0 * h);
1090  vec2 ca(q(0) - std::min(q(0), (q(1) < 0.0) ? r1 : r2), abs(q(1)) - h);
1091  vec2 k1_q(k1(0) - q(0), k1(1) - q(1));
1092  cb(i) = q(i) - k1(i) +
1093  k2(i) * clamp((k1_q(i) * k2(i)) / (k2(i) * k2(i)), 0.0, 1.0);
1094  double s = (cb(0) < 0.0 && ca(1) < 0.0) ? -1.0 : 1.0;
1095  return s * sqrt(std::min((ca(i) * ca(i)), (cb(i) * cb(i))));
1096  };
1097 
1098  const double eps = 1e-6;
1099  auto cone_distance =
1100  sd_capped_cone(pointCoords, half_height, rAdius, smallRadius);
1101 
1102  auto get_d = [&](int a) {
1103  auto pt = pointCoords;
1104  pt(a) += eps;
1105  return (sd_capped_cone(pt, half_height, rAdius, smallRadius) -
1106  cone_distance) /
1107  eps;
1108  };
1109 
1110  // for debugging
1111  // auto my_rand_mn = [](double M, double N) {
1112  // return M + (rand() / (RAND_MAX / (N - M)));
1113  // };
1114  // auto t_off = getBodyOffset();
1115  // std::ofstream afile("cone_pts.csv", std::ios::ate);
1116  // if (afile.is_open()) {
1117  // afile << "X,Y,Z\n";
1118  // Tensor1<double, 3> p(0, 0, 0);
1119  // for (int n = 0; n != 1e6; n++) {
1120  // double x = my_rand_mn(-3, 3);
1121  // double y = my_rand_mn(-3, 3);
1122  // double z = my_rand_mn(-3, 3);
1123  // Tensor1<double, 3> coords(x, y, z);
1124  // if (sd_capped_cone(coords, half_height, rAdius, smallRadius) < 0) {
1125 
1126  // afile << x + t_off(0) << "," << y + t_off(1) << "," << z + t_off(2)
1127  // << "\n";
1128  // }
1129  // }
1130  // afile.close();
1131  // }
1132 
1133  Tensor1<double, 3> my_normal(get_d(0), get_d(1), get_d(2));
1134  my_normal.normalize();
1135 
1136  gAp = cone_distance;
1137  tNormal(i) = rotationMat(j, i) * my_normal(j);
1138 
1139  return tNormal;
1140  };
static Index< 'p', 3 > p
constexpr double a
Definition: single.cpp:4
double h
convective heat coefficient
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 ConeRigidBody::getRollerDataForTag ( )
virtual

Implements RigidBodyData.

Definition at line 975 of file RigidBodies.hpp.

975  {
977  bodyType = CONE;
978  dataForTags[0] = rAdius;
979  dataForTags[2] = hEight;
980  dataForTags[3] = aNgle;
982  }
@ CONE
Definition: RigidBodies.hpp:21
#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

◆ aNgle

double ConeRigidBody::aNgle

Definition at line 964 of file RigidBodies.hpp.

◆ cP

shared_ptr<ArbitrarySurfaceData> ConeRigidBody::cP

Definition at line 957 of file RigidBodies.hpp.

◆ hEight

double ConeRigidBody::hEight

Definition at line 959 of file RigidBodies.hpp.

◆ oFfset

double ConeRigidBody::oFfset

Definition at line 962 of file RigidBodies.hpp.

◆ rAdius

double ConeRigidBody::rAdius

Definition at line 960 of file RigidBodies.hpp.

◆ sLope

double ConeRigidBody::sLope

Definition at line 958 of file RigidBodies.hpp.

◆ smallRadius

double ConeRigidBody::smallRadius

Definition at line 961 of file RigidBodies.hpp.


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