Public Member Functions | Public Attributes | Private Attributes | List of all members
ScaleCoupledBeam Class Reference
+ Inheritance diagram for ScaleCoupledBeam:

Public Member Functions

 ScaleCoupledBeam ()
 
void setupOomph ()
 
void setupMercury ()
 
void actionsBeforeSolve () override
 Write header of output file. More...
 
void actionsAfterSolve () override
 Write header of output file. More...
 
void actionsBeforeOomphTimeStep () override
 Each time step, compute deflection, elastic, kinetic and gravitational energy, and write to output file. More...
 
double getBeamDeflection () const
 Computes beam deflection. More...
 
 ScaleCoupledBeam ()
 
void setupOomph ()
 
void setupMercury ()
 
void actionsBeforeSolve () override
 Write header of output file. More...
 
void actionsAfterSolve () override
 Write header of output file. More...
 
void actionsBeforeOomphTimeStep () override
 Each time step, compute deflection, elastic, kinetic and gravitational energy, and write to output file. More...
 
double getBeamDeflection () const
 Computes beam deflection. More...
 
- Public Member Functions inherited from ScaleCoupling< M, O >
const std::vector< CoupledElement > & getCoupledElements ()
 get coupled element More...
 
const std::vector< CoupledParticle > & getCoupledParticles ()
 get coupled particle More...
 
void setPenalty (double penalty)
 Sets penalty parameter. More...
 
void setCouplingWeight (const std::function< double(double x, double y, double z)> &couplingWeight)
 Sets weight function (determines which nodes/el_pts are in coupling zone) More...
 
void solveScaleCoupling ()
 Computes nStep, the ratio of FEM and DEM time steps, then calls solveSurfaceCoupling(nStep) More...
 
- Public Member Functions inherited from BaseCoupling< M, O >
 BaseCoupling ()=default
 
void setName (std::string name)
 
std::string getName () const
 
void removeOldFiles () const
 
void writeEneTimeStep (std::ostream &os) const override
 
void writeEneHeader (std::ostream &os) const override
 
void solveOomph ()
 
void solveMercury (unsigned long nt)
 
void setCGWidth (const double &width)
 
double getCGWidth ()
 
bool useCGMapping ()
 
CGFunctions::LucyXYZ getCGFunction ()
 

Public Attributes

std::ofstream out
 output file stream More...
 
std::ofstream outFile
 output file stream More...
 

Private Attributes

double length = 20.0
 
double distance = 0.4
 
double bulkDensity = 1309
 
double elasticModulus = 1e8
 
double overlapLength = 10.0*distance
 
double penalty = 8e9
 
double velocity = 1e-1
 

Detailed Description

Define a coupled problem

Constructor & Destructor Documentation

◆ ScaleCoupledBeam() [1/2]

ScaleCoupledBeam::ScaleCoupledBeam ( )
inline
47  {
48  //set name
49  setName("ScaleCoupledBeam");
50  //remove existing output files
51  removeOldFiles();
52 
53  // setup steps
54  setupOomph();
55  setupMercury(); //sets timeMax
56 
57  // set weight function
58  auto couplingWeight = [this] (double x,double y,double z) {
59  if (x<0.5*(length-overlapLength)) {
60  // DEM region
61  return 1.0;
62  } else if (x>0.5*(length+overlapLength)) {
63  // FEM region
64  return 0.0;
65  } else {
66  // coupled region
67  //return (0.5*(length+overlapLength)-x)/overlapLength;
68  double unit = (0.5*(length+overlapLength)-x)/overlapLength;
69  return 0.5+0.4*std::sin(constants::pi*(unit-0.5));
70  }
71  };
72  setCouplingWeight(couplingWeight);
73  setPenalty(penalty);
74 
75  // setup time
76  double waveSpeed = sqrt(elasticModulus/bulkDensity);
77  double propagationTime = length/waveSpeed;
78  setTimeMax(1.25*propagationTime);
79  setOomphTimeStep(getTimeStep());
80  //setOomphTimeStep(propagationTime/100.0);
81  setSaveCount(getTimeMax()/getTimeStep()/100.0);
82  logger(INFO,"TimeMax: %, nTimeSteps %", getTimeMax(), getTimeMax()/getTimeStep());
83 
84  // solve
85  writeToVTK();
86  solveScaleCoupling();
87  saveSolidMesh();
88  }
logger
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
Log::INFO
@ INFO
BaseCoupling::removeOldFiles
void removeOldFiles() const
Definition: BaseCoupling.h:68
BaseCoupling::setName
void setName(std::string name)
Definition: BaseCoupling.h:59
ScaleCoupledBeam::elasticModulus
double elasticModulus
Definition: ScaleCoupledBeam.cpp:40
ScaleCoupledBeam::bulkDensity
double bulkDensity
Definition: ScaleCoupledBeam.cpp:39
ScaleCoupledBeam::overlapLength
double overlapLength
Definition: ScaleCoupledBeam.cpp:41
ScaleCoupledBeam::penalty
double penalty
Definition: ScaleCoupledBeam.cpp:42
ScaleCoupledBeam::setupOomph
void setupOomph()
Definition: ScaleCoupledBeam.cpp:90
ScaleCoupledBeam::setupMercury
void setupMercury()
Definition: ScaleCoupledBeam.cpp:102
ScaleCoupledBeam::length
double length
Definition: ScaleCoupledBeam.cpp:37
ScaleCoupling::setCouplingWeight
void setCouplingWeight(const std::function< double(double x, double y, double z)> &couplingWeight)
Sets weight function (determines which nodes/el_pts are in coupling zone)
Definition: ScaleCoupling.h:120
ScaleCoupling::solveScaleCoupling
void solveScaleCoupling()
Computes nStep, the ratio of FEM and DEM time steps, then calls solveSurfaceCoupling(nStep)
Definition: ScaleCoupling.h:125
ScaleCoupling::setPenalty
void setPenalty(double penalty)
Sets penalty parameter.
Definition: ScaleCoupling.h:115
constants::pi
const Mdouble pi
Definition: ExtendedMath.h:45
mathsFunc::sin
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44

References bulkDensity, elasticModulus, INFO, length, logger, overlapLength, penalty, constants::pi, BaseCoupling< M, O >::removeOldFiles(), ScaleCoupling< M, O >::setCouplingWeight(), BaseCoupling< M, O >::setName(), ScaleCoupling< M, O >::setPenalty(), setupMercury(), setupOomph(), mathsFunc::sin(), and ScaleCoupling< M, O >::solveScaleCoupling().

◆ ScaleCoupledBeam() [2/2]

ScaleCoupledBeam::ScaleCoupledBeam ( )
inline
47  {
48  //set name
49  setName("ScaleCoupledBeamUnitTest");
50  //remove existing output files
51  removeOldFiles();
52 
53  // setup steps
54  setupOomph();
55  setupMercury(); //sets timeMax
56 
57  // set weight function
58  auto couplingWeight = [this] (double x,double y,double z) {
59  if (x<0.5*(length-overlapLength)) {
60  // DEM region
61  return 1.0;
62  } else if (x>0.5*(length+overlapLength)) {
63  // FEM region
64  return 0.0;
65  } else {
66  // coupled region
67  //return (0.5*(length+overlapLength)-x)/overlapLength;
68  double unit = (0.5*(length+overlapLength)-x)/overlapLength;
69  return 0.5+0.4*std::sin(constants::pi*(unit-0.5));
70  }
71  };
72  setCouplingWeight(couplingWeight);
73  setPenalty(penalty);
74 
75  // setup time
76  setTimeMax(2.0*getTimeStep()); // different to ScaleCoupledBeam
77  setOomphTimeStep(getTimeStep());
78  setSaveCount(getTimeMax()/getTimeStep()/100.0);
79  logger(INFO,"TimeMax: %, nTimeSteps %", getTimeMax(), getTimeMax()/getTimeStep());
80 
81  // solve
82  writeToVTK();
83  solveScaleCoupling();
84  saveSolidMesh();
85  }

References INFO, length, logger, overlapLength, penalty, constants::pi, BaseCoupling< M, O >::removeOldFiles(), ScaleCoupling< M, O >::setCouplingWeight(), BaseCoupling< M, O >::setName(), ScaleCoupling< M, O >::setPenalty(), setupMercury(), setupOomph(), mathsFunc::sin(), and ScaleCoupling< M, O >::solveScaleCoupling().

Member Function Documentation

◆ actionsAfterSolve() [1/2]

void ScaleCoupledBeam::actionsAfterSolve ( )
inlineoverride

Write header of output file.

150  {
151  out.close();
152  }
ScaleCoupledBeam::out
std::ofstream out
output file stream
Definition: ScaleCoupledBeam.cpp:178

References out.

◆ actionsAfterSolve() [2/2]

void ScaleCoupledBeam::actionsAfterSolve ( )
inlineoverride

Write header of output file.

142  {
143  outFile.close();
144  // checking errors
145  Vec3D couplingForce = getCoupledParticles().back().couplingForce;
146  helpers::check(couplingForce.X,-15917,0.1,"Coupling force");
147  }
ScaleCoupledBeam::outFile
std::ofstream outFile
output file stream
Definition: ScaleCoupledBeamUnitTest.cpp:173
ScaleCoupling::getCoupledParticles
const std::vector< CoupledParticle > & getCoupledParticles()
get coupled particle
Definition: ScaleCoupling.h:110
Vec3D
Definition: Vector.h:51
Vec3D::X
Mdouble X
the vector components
Definition: Vector.h:66
helpers::check
void check(double real, double ideal, double error, std::string errorMessage)
Definition: TestHelpers.cc:37

References helpers::check(), ScaleCoupling< M, O >::getCoupledParticles(), outFile, and Vec3D::X.

◆ actionsBeforeOomphTimeStep() [1/2]

void ScaleCoupledBeam::actionsBeforeOomphTimeStep ( )
inlineoverride

Each time step, compute deflection, elastic, kinetic and gravitational energy, and write to output file.

155  {
156  double mass, elasticEnergy, kineticEnergy;
157  Vector<double> com(3), linearMomentum(3), angularMomentum(3);
158  getMassMomentumEnergy(mass, com, linearMomentum, angularMomentum, elasticEnergy, kineticEnergy);
159  static double comZ0 = com[2];
160  double gravEnergy = 9.8*mass*(com[2]-comZ0);
161  out << getOomphTime() << ' ' << getBeamDeflection() << ' ' << elasticEnergy << ' ' << kineticEnergy << ' ' << gravEnergy << std::endl;
162  std::cout << getOomphTime() << ' ' << getBeamDeflection() << ' ' << elasticEnergy << ' ' << kineticEnergy << ' ' << gravEnergy << '\n';
163  }
ScaleCoupledBeam::getBeamDeflection
double getBeamDeflection() const
Computes beam deflection.
Definition: ScaleCoupledBeam.cpp:166

References getBeamDeflection(), and out.

◆ actionsBeforeOomphTimeStep() [2/2]

void ScaleCoupledBeam::actionsBeforeOomphTimeStep ( )
inlineoverride

Each time step, compute deflection, elastic, kinetic and gravitational energy, and write to output file.

150  {
151  double mass, elasticEnergy, kineticEnergy;
152  Vector<double> com(3), linearMomentum(3), angularMomentum(3);
153  getMassMomentumEnergy(mass, com, linearMomentum, angularMomentum, elasticEnergy, kineticEnergy);
154  static double comZ0 = com[2];
155  double gravEnergy = 9.8*mass*(com[2]-comZ0);
156  outFile << getOomphTime() << ' ' << getBeamDeflection() << ' ' << elasticEnergy << ' ' << kineticEnergy << ' ' << gravEnergy << std::endl;
157  std::cout << getOomphTime() << ' ' << getBeamDeflection() << ' ' << elasticEnergy << ' ' << kineticEnergy << ' ' << gravEnergy << '\n';
158  }

References getBeamDeflection(), and outFile.

◆ actionsBeforeSolve() [1/2]

void ScaleCoupledBeam::actionsBeforeSolve ( )
inlineoverride

Write header of output file.

142  {
143  helpers::writeToFile(getName()+".gnu", "set key autotitle columnheader\n"
144  "p 'SolidBeamUnsteady.out' u 1:3 w l, '' u 1:4 w l, '' u 1:5 w l, '' u 1:($3+$4+$5) t 'totalEnergy' w l");
145  out.open(getName()+".out");
146  out << "time deflection elasticEnergy kineticEnergy gravEnergy\n";
147  }
BaseCoupling::getName
std::string getName() const
Definition: BaseCoupling.h:64
helpers::writeToFile
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: FileIOHelpers.cc:58

References BaseCoupling< M, O >::getName(), out, and helpers::writeToFile().

◆ actionsBeforeSolve() [2/2]

void ScaleCoupledBeam::actionsBeforeSolve ( )
inlineoverride

Write header of output file.

134  {
135  helpers::writeToFile(getName()+".gnu", "set key autotitle columnheader\n"
136  "p 'SolidBeamUnsteady.out' u 1:3 w l, '' u 1:4 w l, '' u 1:5 w l, '' u 1:($3+$4+$5) t 'totalEnergy' w l");
137  outFile.open(getName()+".out");
138  outFile << "time deflection elasticEnergy kineticEnergy gravEnergy\n";
139  }

References BaseCoupling< M, O >::getName(), outFile, and helpers::writeToFile().

◆ getBeamDeflection() [1/2]

double ScaleCoupledBeam::getBeamDeflection ( ) const
inline

Computes beam deflection.

166  {
167  std::array<double, 3> min, max;
168  getDomainSize(min, max);
169 
170  Vector<double> xi(3);
171  xi[0] = max[0];
172  xi[1] = 0.5 * (max[1] + min[1]);
173  xi[2] = 0.5 * (max[2] + min[2]);
174  return getDeflection(xi, 2);
175  }

Referenced by actionsBeforeOomphTimeStep().

◆ getBeamDeflection() [2/2]

double ScaleCoupledBeam::getBeamDeflection ( ) const
inline

Computes beam deflection.

161  {
162  std::array<double, 3> min, max;
163  getDomainSize(min, max);
164 
165  Vector<double> xi(3);
166  xi[0] = max[0];
167  xi[1] = 0.5 * (max[1] + min[1]);
168  xi[2] = 0.5 * (max[2] + min[2]);
169  return getDeflection(xi, 2);
170  }

◆ setupMercury() [1/2]

void ScaleCoupledBeam::setupMercury ( )
inline
103  {
104  // positive overlap means adhesive forces
105  double overlap = 0.0*distance; //overlap between particles
106  double radius = 0.5*(distance+overlap); // particle radius
107  double particleDensity = bulkDensity*mathsFunc::cubic(distance) / (constants::pi/6.0*mathsFunc::cubic(2.0*radius));
108  double stiffness = elasticModulus*distance; // particle stiffness
109 
110  setDomain(Vec3D(-radius,-radius,0),Vec3D(radius,radius,0.5*(length+overlapLength)));
111  setParticlesWriteVTK(true);
112 
113  auto species = speciesHandler.copyAndAddObject(LinearViscoelasticReversibleAdhesiveSpecies());
114  species->setDensity(particleDensity);
115  species->setStiffness(stiffness);
116  species->setAdhesionStiffness(stiffness);
117  species->setAdhesionForceMax(stiffness*overlap);
118 
119  SphericalParticle p(species);
120  p.setRadius(radius);
121  //p.setVelocity(Vec3D(1,0,0));
122  double dx = distance;
123  auto n = (unsigned)(0.5*(length+overlapLength)/dx);
124  for (int i = 0; i<n; ++i) {
125  for (int j = 0; j<2; ++j) {
126  for (int k = 0; k<2; ++k) {
127  p.setPosition(Vec3D(1+2*i, j==0?-1:1, k==0?-1:1)*0.5*distance);
128  particleHandler.copyAndAddObject(p);
129  if (i==0) {
130  particleHandler.getLastObject()->setVelocity(Vec3D(velocity,0,0));
131  //particleHandler.getLastObject()->fixParticle();
132  //particleHandler.getLastObject()->setPrescribedVelocity(velocityAtBoundary);
133  }
134  }
135  }
136  }
137 
138  setTimeStep(0.5*species->computeTimeStep(p.getMass()));
139  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:32
LinearViscoelasticReversibleAdhesiveSpecies
Species< LinearViscoelasticNormalSpecies, EmptyFrictionSpecies, ReversibleAdhesiveSpecies > LinearViscoelasticReversibleAdhesiveSpecies
Definition: LinearViscoelasticReversibleAdhesiveSpecies.h:35
ScaleCoupledBeam::distance
double distance
Definition: ScaleCoupledBeam.cpp:38
ScaleCoupledBeam::velocity
double velocity
Definition: ScaleCoupledBeam.cpp:43
SphericalParticle
A spherical particle is the most simple particle used in MercuryDPM.
Definition: SphericalParticle.h:37
constants::i
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
mathsFunc::cubic
T cubic(const T val)
calculates the cube of a number
Definition: ExtendedMath.h:115

References bulkDensity, mathsFunc::cubic(), distance, elasticModulus, BaseParticle::getMass(), constants::i, length, n, overlapLength, constants::pi, BaseInteractable::setPosition(), BaseParticle::setRadius(), and velocity.

Referenced by ScaleCoupledBeam().

◆ setupMercury() [2/2]

void ScaleCoupledBeam::setupMercury ( )
inline
100  {
101  // positive overlap means adhesive forces
102  double overlap = 0.0*distance; //overlap between particles
103  double radius = 0.5*(distance+overlap); // particle radius
104  double particleDensity = bulkDensity*mathsFunc::cubic(distance) / (constants::pi/6.0*mathsFunc::cubic(2.0*radius));
105  double stiffness = elasticModulus*distance; // particle stiffness
106 
107  setDomain(Vec3D(-radius,-radius,0),Vec3D(radius,radius,0.5*(length+overlapLength)));
108  setParticlesWriteVTK(true);
109 
110  auto species = speciesHandler.copyAndAddObject(LinearViscoelasticReversibleAdhesiveSpecies());
111  species->setDensity(particleDensity);
112  species->setStiffness(stiffness);
113  species->setAdhesionStiffness(stiffness);
114  species->setAdhesionForceMax(stiffness*overlap);
115 
116  SphericalParticle p(species);
117  p.setRadius(radius);
118  p.setVelocity(Vec3D(1,0,0)); // different to ScaleCoupledBeam
119  double dx = distance;
120  auto n = (unsigned)(0.5*(length+overlapLength)/dx);
121  for (int i = 0; i<n; ++i) {
122  for (int j = 0; j<2; ++j) {
123  for (int k = 0; k<2; ++k) {
124  p.setPosition(Vec3D(1+2*i, j==0?-1:1, k==0?-1:1)*0.5*distance);
125  particleHandler.copyAndAddObject(p); // different to ScaleCoupledBeam
126  }
127  }
128  }
129 
130  setTimeStep(0.5*species->computeTimeStep(p.getMass()));
131  }

References bulkDensity, mathsFunc::cubic(), distance, elasticModulus, BaseParticle::getMass(), constants::i, length, n, overlapLength, constants::pi, BaseInteractable::setPosition(), BaseParticle::setRadius(), and BaseInteractable::setVelocity().

◆ setupOomph() [1/2]

void ScaleCoupledBeam::setupOomph ( )
inline
90  {
91  setElasticModulus(elasticModulus);
92  setDensity(bulkDensity);
93  setSolidCubicMesh(round(0.5*(length+overlapLength)/distance), 2, 2,
94  0.5*(length-overlapLength), length, -distance, distance, -distance, distance);
95  //pinBoundaries({Beam::Boundary::X_MIN});
96  setNewtonSolverTolerance(1e-2);
97  prepareForSolve();
98  linear_solver_pt()->disable_doc_time();
99  //disable_info_in_newton_solve();
100  }
helpers::round
Mdouble round(Mdouble value, unsigned int precision)
rounds a floating point number with a given precision
Definition: MathHelpers.cc:28

References bulkDensity, distance, elasticModulus, length, overlapLength, and helpers::round().

Referenced by ScaleCoupledBeam().

◆ setupOomph() [2/2]

void ScaleCoupledBeam::setupOomph ( )
inline
87  {
88  setElasticModulus(elasticModulus);
89  setDensity(bulkDensity);
90  setSolidCubicMesh(round(0.5*(length+overlapLength)/distance), 2, 2,
91  0.5*(length-overlapLength), length, -distance, distance, -distance, distance);
92  //pinBoundaries({Beam::Boundary::X_MIN});
93  setNewtonSolverTolerance(1e-2);
94  prepareForSolve();
95  linear_solver_pt()->disable_doc_time();
96  //disable_info_in_newton_solve();
97  }

References bulkDensity, distance, elasticModulus, length, overlapLength, and helpers::round().

Member Data Documentation

◆ bulkDensity

double ScaleCoupledBeam::bulkDensity = 1309
private

Referenced by ScaleCoupledBeam(), setupMercury(), and setupOomph().

◆ distance

double ScaleCoupledBeam::distance = 0.4
private

Referenced by setupMercury(), and setupOomph().

◆ elasticModulus

double ScaleCoupledBeam::elasticModulus = 1e8
private

Referenced by ScaleCoupledBeam(), setupMercury(), and setupOomph().

◆ length

double ScaleCoupledBeam::length = 20.0
private

Referenced by ScaleCoupledBeam(), setupMercury(), and setupOomph().

◆ out

std::ofstream ScaleCoupledBeam::out

output file stream

Referenced by actionsAfterSolve(), actionsBeforeOomphTimeStep(), and actionsBeforeSolve().

◆ outFile

std::ofstream ScaleCoupledBeam::outFile

output file stream

Referenced by actionsAfterSolve(), actionsBeforeOomphTimeStep(), and actionsBeforeSolve().

◆ overlapLength

double ScaleCoupledBeam::overlapLength = 10.0*distance
private

Referenced by ScaleCoupledBeam(), setupMercury(), and setupOomph().

◆ penalty

double ScaleCoupledBeam::penalty = 8e9
private

Referenced by ScaleCoupledBeam().

◆ velocity

double ScaleCoupledBeam::velocity = 1e-1
private

Referenced by setupMercury().

The documentation for this class was generated from the following files: