SolidProblem< ELEMENT_TYPE > Class Template Reference

#include <SolidProblem.h>

+ Inheritance diagram for SolidProblem< ELEMENT_TYPE >:

Classes

class  SolidCubicMesh
 

Public Types

enum  Boundary : unsigned {
  Z_MIN = 0 , Y_MIN = 1 , X_MAX = 2 , Y_MAX = 3 ,
  X_MIN = 4 , Z_MAX = 5
}
 

Public Member Functions

 SolidProblem ()
 Constructor: set default constitutive law and time stepper. More...
 
void setName (const std::string &name)
 set function for name_ More...
 
std::string getName () const
 get function for name_ More...
 
void setElasticModulus (double elasticModulus)
 set function for elasticModulus_ More...
 
double getElasticModulus () const
 get function for elasticModulus_ More...
 
void setOomphGravity (double gravity)
 set function for elasticModulus_ More...
 
double getOomphGravity () const
 get function for gravity_ More...
 
void setPoissonRatio (double poissonRatio)
 set function for poissonRatio_ More...
 
double getPoissonRatio () const
 get function for poissonRatio_ More...
 
void setDensity (double density)
 set function for density_ More...
 
double getDensity () const
 get function for density_ More...
 
void setBodyForceAsGravity (double gravity=9.8)
 set function for body_force_pt More...
 
void setIsPinned (std::function< bool(SolidNode *, unsigned)> isPinned)
 set function for isPinned_ More...
 
void pinBoundary (unsigned b)
 
void pinBoundaries (std::vector< unsigned > b)
 
std::enable_if< std::is_base_of< RefineableQDPVDElement< 3, 2 >, ELEMENT >::value, void > setDissipation (double dissipation)
 Sets the dissipation coefficient for all elements. More...
 
void setNewtonSolverTolerance (double Newton_solver_tolerance)
 set function for Newton_solver_tolerance More...
 
void setMaxNewtonIterations (unsigned Max_newton_iterations)
 set function for Max_newton_iterations More...
 
void setOomphTimeStep (double dt)
 set function for time step More...
 
double getOomphTimeStep () const
 get function for time step More...
 
double getOomphTime () const
 get function for current time More...
 
SolidMesh *& solid_mesh_pt ()
 Get function for the solid mesh pointer. More...
 
SolidMesh *& traction_mesh_pt ()
 Get function for the traction mesh pointer. More...
 
SolidMesh *const & solid_mesh_pt () const
 Get function for the solid mesh pointer. More...
 
SolidMesh *const & traction_mesh_pt () const
 Get function for the traction mesh pointer. More...
 
void getDomainSize (std::array< double, 3 > &min, std::array< double, 3 > &max) const
 Computes the domain size (min/max of the nodal positions in x/y/z) More...
 
void prepareForSolve ()
 
void countPinned ()
 returns statistics about pinned nodes to the console More...
 
void solveSteady ()
 
virtual void actionsBeforeSolve ()
 
virtual void actionsAfterSolve ()
 
virtual void actionsBeforeOomphTimeStep ()
 
void solveUnsteady (double timeMax, double dt, unsigned saveCount=10)
 
void removeOldFiles () const
 
void get_x (const Vector< double > &xi, Vector< double > &x) const
 
double getDeflection (Vector< double > xi, unsigned d) const
 
void setSolidCubicMesh (const unsigned &nx, const unsigned &ny, const unsigned &nz, const double &xMin, const double &xMax, const double &yMin, const double &yMax, const double &zMin, const double &zMax)
 
void saveSolidMesh ()
 
void loadSolidMesh (std::string name)
 
void writeToVTK ()
 
void getMassMomentumEnergy (double &mass, Vector< double > &com, Vector< double > &linearMomentum, Vector< double > &angularMomentum, double &elasticEnergy, double &kineticEnergy)
 See PVDEquationsBase<DIM>::get_energy. More...
 

Protected Types

typedef ELEMENT_TYPE ELEMENT
 

Protected Attributes

std::string name_
 
double elasticModulus_ = 0
 Elastic modulus (set via setSolid) More...
 
double poissonRatio_ = 0
 Poisson's ratio (set via setSolid) More...
 
double density_ = 0
 Density. More...
 
double gravity_ = 0
 Density. More...
 
void(* body_force_fct )(const double &time, const Vector< double > &xi, Vector< double > &b) = nullptr
 Pointer to the body force function. More...
 
ConstitutiveLaw * constitutive_law_pt = nullptr
 Pointer to constitutive law (should be set in constructor) More...
 
SolidMeshSolid_mesh_pt = nullptr
 Pointer to solid mesh. More...
 
SolidMeshTraction_mesh_pt = nullptr
 Pointer to mesh of traction elements. More...
 
std::function< bool(SolidNode *, unsigned)> isPinned_
 Function to determine which nodal positions are pinned. More...
 

Detailed Description

template<class ELEMENT_TYPE>
class SolidProblem< ELEMENT_TYPE >

Base class for (surface-coupled) problems with a solid body.

Assumptions:

  • The element type is derived from QPVDElement, i.e. a solid.
  • The mesh type is derived from SolidMesh.
  • The constitutive law is derived from Hookean, with parameters elasticModulus, poissonRatio, density.

The solid properties (elasticModulus, poissonRatio, density, constitutive_law_pt, body_force_fct, etc.) are stored in member variables, not in globals as typical in oomph-lib, and accessed by setters and getters.

Additional functionality:

  • setIsPinned allows one to define the pinned nodes using a function.

Member Typedef Documentation

◆ ELEMENT

template<class ELEMENT_TYPE >
typedef ELEMENT_TYPE SolidProblem< ELEMENT_TYPE >::ELEMENT
protected

Member Enumeration Documentation

◆ Boundary

template<class ELEMENT_TYPE >
enum SolidProblem::Boundary : unsigned
Enumerator
Z_MIN 
Y_MIN 
X_MAX 
Y_MAX 
X_MIN 
Z_MAX 
520  {
521  Z_MIN = 0, Y_MIN = 1, X_MAX = 2, Y_MAX = 3, X_MIN = 4, Z_MAX = 5
522  };
@ Z_MIN
Definition: SolidProblem.h:521
@ X_MIN
Definition: SolidProblem.h:521
@ Z_MAX
Definition: SolidProblem.h:521
@ Y_MAX
Definition: SolidProblem.h:521
@ Y_MIN
Definition: SolidProblem.h:521
@ X_MAX
Definition: SolidProblem.h:521

Constructor & Destructor Documentation

◆ SolidProblem()

template<class ELEMENT_TYPE >
SolidProblem< ELEMENT_TYPE >::SolidProblem ( )
inline

Constructor: set default constitutive law and time stepper.

119  {
120  logger(INFO, "Set default constitutive law (GeneralisedHookean) and time stepper (Newmark<2>)");
121 
122  // Set Newton solver tolerance and maximum number of Newton iterations
123  Max_newton_iterations = 20;
124  Newton_solver_tolerance = 1e-10;
125  Max_residuals = constants::inf;
126 
127  // Set constitutive law
128  constitutive_law_pt = new GeneralisedHookean(&poissonRatio_, &elasticModulus_);
129 
130  // Allocate the timestepper
131  add_time_stepper_pt(new Newmark<2>);
132 
133  //build_mesh();
134 
135  // setup boundary conditions
136  //pin_and_assign_eqn_numbers();
137  };
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ INFO
double elasticModulus_
Elastic modulus (set via setSolid)
Definition: SolidProblem.h:89
double poissonRatio_
Poisson's ratio (set via setSolid)
Definition: SolidProblem.h:92
ConstitutiveLaw * constitutive_law_pt
Pointer to constitutive law (should be set in constructor)
Definition: SolidProblem.h:104
const Mdouble inf
Definition: GeneralDefine.h:44

References constants::inf, INFO, and logger.

Member Function Documentation

◆ actionsAfterSolve()

template<class ELEMENT_TYPE >
virtual void SolidProblem< ELEMENT_TYPE >::actionsAfterSolve ( )
inlinevirtual

Reimplemented in Beam.

424 {}

◆ actionsBeforeOomphTimeStep()

template<class ELEMENT_TYPE >
virtual void SolidProblem< ELEMENT_TYPE >::actionsBeforeOomphTimeStep ( )
inlinevirtual

Reimplemented in Beam, SolidBag, and Beam.

426 {}

◆ actionsBeforeSolve()

template<class ELEMENT_TYPE >
virtual void SolidProblem< ELEMENT_TYPE >::actionsBeforeSolve ( )
inlinevirtual

Reimplemented in Beam, and Beam.

422 {}

◆ countPinned()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::countPinned ( )
inline

returns statistics about pinned nodes to the console

392  {
393  // count pinned
394  std::array<unsigned,3> countPinned {0,0,0};
395  unsigned countAll = 0;
396  for (unsigned n = 0; n < solid_mesh_pt()->nnode(); n++)
397  {
398  for (unsigned i = 0; i < 3; i++)
399  {
400  countPinned[i] += solid_mesh_pt()->node_pt(n)->position_is_pinned(i);
401  countAll++;
402  }
403  }
404  unsigned countPinnedAll = countPinned[0] + countPinned[1] +countPinned[2];
405  logger(INFO, "Pinned % of % positions (% free): % in x, % in y, % in z", countPinnedAll, countAll, countAll - countPinnedAll, countPinned[0], countPinned[1], countPinned[2]);
406  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:32
SolidMesh *& solid_mesh_pt()
Get function for the solid mesh pointer.
Definition: SolidProblem.h:284
void countPinned()
returns statistics about pinned nodes to the console
Definition: SolidProblem.h:391
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51

References constants::i, INFO, logger, and n.

◆ get_x()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::get_x ( const Vector< double > &  xi,
Vector< double > &  x 
) const
inline

Returns the x value at a certain xi value

484  {
485  if (OOMPH_MPI_PROCESSOR_NUM>1) {
486  logger(INFO, "get_x does not work with MPI");
487  for (int i = 0; i < 3; ++i) {
488  x[i] = xi[i];
489  }
490  } else {
491  Vector<double> s(3);
492  GeomObject *geom_obj_pt = nullptr;
493  const unsigned long nelement = solid_mesh_pt()->nelement();
494  for (unsigned long i = 0; i < nelement; i++) {
495  auto el_pt = dynamic_cast<ELEMENT *>(solid_mesh_pt()->element_pt(i));
496  el_pt->locate_zeta(xi, geom_obj_pt, s);
497  if (geom_obj_pt) {
498  //logger(INFO,"Point % % % is in element % at % % %",
499  // xi[0],xi[1],xi[2],i,s[0], s[1], s[2]);
500  el_pt->interpolated_x(s, x); //deformed coordinate
501  return;
502  }
503  }
504  logger(ERROR, "x(xi) could not be found");
505  }
506  }
@ ERROR
#define OOMPH_MPI_PROCESSOR_NUM
Definition: SolidProblem.h:56

References ERROR, constants::i, INFO, logger, and OOMPH_MPI_PROCESSOR_NUM.

◆ getDeflection()

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::getDeflection ( Vector< double xi,
unsigned  d 
) const
inline

Returns the difference between the x and xi value in dimension d at a certain xi value

512  {
513  Vector<double> x(3);
514  get_x(xi, x);
515  return x[d] - xi[d];
516  }
void get_x(const Vector< double > &xi, Vector< double > &x) const
Definition: SolidProblem.h:483

◆ getDensity()

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::getDensity ( ) const
inline

get function for density_

200  {
201  return density_;
202  }
double density_
Density.
Definition: SolidProblem.h:95

Referenced by main().

◆ getDomainSize()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::getDomainSize ( std::array< double, 3 > &  min,
std::array< double, 3 > &  max 
) const
inline

Computes the domain size (min/max of the nodal positions in x/y/z)

297  {
298  min[0] = min[1] = min[2] = constants::inf;
299  max[0] = max[1] = max[2] = -constants::inf;
300  logger.assert_always(solid_mesh_pt(),"mesh not found");
301  for (unsigned i = 0; i < solid_mesh_pt()->nnode(); i++)
302  {
303  const auto n = solid_mesh_pt()->node_pt(i);
304  for (int j = 0; j < 3; ++j)
305  {
306  min[j] = std::min(min[j], n->xi(j));
307  max[j] = std::max(max[j], n->xi(j));
308  }
309  }
310  }

References constants::i, constants::inf, logger, and n.

Referenced by CoupledProblem::setupMercury().

◆ getElasticModulus()

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::getElasticModulus ( ) const
inline

get function for elasticModulus_

161  {
162  return elasticModulus_;
163  }

Referenced by main().

◆ getMassMomentumEnergy()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::getMassMomentumEnergy ( double mass,
Vector< double > &  com,
Vector< double > &  linearMomentum,
Vector< double > &  angularMomentum,
double elasticEnergy,
double kineticEnergy 
)
inline

See PVDEquationsBase<DIM>::get_energy.

868  {
869  // Initialise mass
870  mass = 0;
871  // Initialise center of mass
872  com.initialise(0);
873  // Initialise momentum
874  linearMomentum.initialise(0);
875  angularMomentum.initialise(0);
876  // Initialise energy
877  elasticEnergy = 0;
878  kineticEnergy = 0;
879 
880  // For each element
881  const unsigned long nelement = this->solid_mesh_pt()->nelement();
882  for (unsigned long e = 0; e < nelement; e++) {
883 
884  // Get the pointer to the element
885  ELEMENT* e_pt = dynamic_cast<ELEMENT*>(Solid_mesh_pt->element_pt(e));
886 
887  const unsigned DIM = e_pt->dim();
888 
889  //Find out how many integration points there are
890  unsigned n_intpt = e_pt->integral_pt()->nweight();
891 
892  //Set the Vector to hold local coordinates
893  Vector<double> s(DIM);
894 
895  //Find out how many nodes there are
896  const unsigned n_node = e_pt->nnode();
897 
898  //Find out how many positional dofs there are
899  const unsigned n_position_type = e_pt->nnodal_position_type();
900 
901  //Set up memory for the shape functions
902  Shape psi(n_node, n_position_type);
903  DShape dpsidxi(n_node, n_position_type, DIM);
904 
905  // Timescale ratio (non-dim density)
906  double lambda_sq = e_pt->lambda_sq();
907 
908  //Loop over the integration points
909  for (unsigned ipt = 0; ipt < n_intpt; ipt++) {
910  //Assign local coordinate s
911  for (unsigned i = 0; i < DIM; i++) { s[i] = e_pt->integral_pt()->knot(ipt, i); }
912 
913  //Get the integral weight
914  double w = e_pt->integral_pt()->weight(ipt);
915 
916  //Evaluate the shape function and its derivatives, and get Jacobian
917  double J = e_pt->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
918 
919  //Get mass and the coupling weight at the integration point
920  double coupling_w = 0;
921  //for (unsigned l = 0; l < n_node; l++)
922  //{
923  // double nodal_coupling_w = dynamic_cast<SolidNode*>(e_pt->node_pt(l))->get_coupling_weight();
924  // double psi_ = psi(l);
925  // coupling_w += psi_ * nodal_coupling_w;
926  //}
927 
928  //Get the coordinates of the integration point
929  Vector<double> x(DIM, 0.0);
930  e_pt->interpolated_x(s,x);
931 
932  //Storage for Lagrangian coordinates and velocity (initialised to zero)
933  Vector<double> interpolated_xi(DIM, 0.0);
934  Vector<double> veloc(DIM, 0.0);
935 
936  //Calculate lagrangian coordinates
937  for (unsigned l = 0; l < n_node; l++) {
938  //Loop over positional dofs
939  for (unsigned k = 0; k < n_position_type; k++) {
940  //Loop over displacement components (deformed position)
941  for (unsigned i = 0; i < DIM; i++) {
942  //Calculate the Lagrangian coordinates
943  interpolated_xi[i] += e_pt->lagrangian_position_gen(l, k, i) * psi(l, k);
944 
945  //Calculate the velocity components (if unsteady solve)
946  if (e_pt->is_inertia_enabled()) {
947  veloc[i] += e_pt->dnodal_position_gen_dt(l, k, i) * psi(l, k);
948  }
949  }
950  }
951  }
952 
953  //Get isotropic growth factor
954  double gamma = 1.0;
955  e_pt->get_isotropic_growth(ipt, s, interpolated_xi, gamma);
956 
957  //Premultiply the undeformed volume ratio (from the isotropic
958  // growth), the integral weights, the coupling weights, and the Jacobian
959  double W = gamma * w * (1.0 - coupling_w) * J;
960 
961  DenseMatrix<double> sigma(DIM, DIM);
962  DenseMatrix<double> strain(DIM, DIM);
963 
964  //Now calculate the stress tensor from the constitutive law
965  e_pt->get_stress(s, sigma);
966 
967  // Add pre-stress
968  for (unsigned i = 0; i < DIM; i++) {
969  for (unsigned j = 0; j < DIM; j++) {
970  sigma(i, j) += e_pt->prestress(i, j, interpolated_xi);
971  }
972  }
973 
974  //get the strain
975  e_pt->get_strain(s, strain);
976 
977  // Initialise
978  double localElasticEnergy = 0;
979  double velocitySquared = 0;
980 
981  // Compute integrals
982  for (unsigned i = 0; i < DIM; i++) {
983  for (unsigned j = 0; j < DIM; j++) {
984  localElasticEnergy += sigma(i, j) * strain(i, j);
985  }
986  velocitySquared += veloc[i] * veloc[i];
987  }
988 
989  // Mass
990  mass += lambda_sq * W;
991  // Linear momentum and angular momentum
992  Vector<double> cross_product(DIM, 0);
993  VectorHelpers::cross(interpolated_xi, veloc, cross_product);
994  for (unsigned i = 0; i < DIM; i++) {
995  com[i] += lambda_sq * W * x[i];
996  linearMomentum[i] += lambda_sq * veloc[i] * W;
997  angularMomentum[i] += lambda_sq * cross_product[i] * W;
998  }
999  // Potential energy
1000  elasticEnergy += 0.5 * localElasticEnergy * W;
1001  // Kinetic energy
1002  kineticEnergy += lambda_sq * 0.5 * velocitySquared * W;
1003  }
1004  }
1005  for (unsigned i = 0; i < com.size(); i++) {
1006  com[i] /= mass;
1007  }
1008  }
SolidMesh * Solid_mesh_pt
Pointer to solid mesh.
Definition: SolidProblem.h:107
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:137

References mathsFunc::gamma(), and constants::i.

◆ getName()

template<class ELEMENT_TYPE >
std::string SolidProblem< ELEMENT_TYPE >::getName ( ) const
inline

get function for name_

148  {
149  return name_;
150  }
std::string name_
Definition: SolidProblem.h:86

◆ getOomphGravity()

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::getOomphGravity ( ) const
inline

get function for gravity_

174  {
175  return gravity_;
176  }
double gravity_
Density.
Definition: SolidProblem.h:98

◆ getOomphTime()

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::getOomphTime ( ) const
inline

get function for current time

279  {
280  return time_pt()->time();
281  }

◆ getOomphTimeStep()

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::getOomphTimeStep ( ) const
inline

get function for time step

274  {
275  return time_pt()->dt();
276  }

◆ getPoissonRatio()

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::getPoissonRatio ( ) const
inline

get function for poissonRatio_

187  {
188  return poissonRatio_;
189  }

◆ loadSolidMesh()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::loadSolidMesh ( std::string  name)
inline
582  {
583  //if (not std::is_base_of<QElement<3, 2>, ELEMENT>::value) {
584  // logger(ERROR, "This kind of mesh requires a QElement.");
585  //}
586 
587  std::ifstream mesh(name);
588  logger.assert_always(mesh.good(),"Mesh file % could not be opened",name);
589  unsigned nx, ny, nz;
590  mesh >> nx >> ny >> nz;
591  logger(INFO, "Loaded %x%x% cubic mesh from %", nx, ny, nz, name);
592  solid_mesh_pt() = new SolidCubicMesh(nx, ny, nz, 0, 1, 0, 1, 0, 1, time_stepper_pt());
593 
594  double xi, x;
595  bool pin;
596  for (int i = 0; i < solid_mesh_pt()->nnode(); ++i)
597  {
598  SolidNode* n = solid_mesh_pt()->node_pt(i);
599  for (int j = 0; j < 3; ++j)
600  {
601  mesh >> xi >> x >> pin;
602  n->xi(j) = xi;
603  n->x(j) = x;
604  if (pin) n->pin_position(j);
605  }
606  }
607  }
Simple cubic mesh upgraded to become a solid mesh.
Definition: glued_mesh_stuff.h:208
std::string name
Definition: MercuryProb.h:48

References constants::i, INFO, logger, n, and units::name.

◆ pinBoundaries()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::pinBoundaries ( std::vector< unsigned >  b)
inline
235  {
236  for (const auto a: b) {
237  logger(INFO, "Pinning nodes on boundary %", a);
238  }
239  isPinned_ = [b](SolidNode *n, unsigned d) {
240  for (const auto a : b) {
241  if (n->is_on_boundary(a)) return true;
242  }
243  return false;
244  };
245  }
std::function< bool(SolidNode *, unsigned)> isPinned_
Function to determine which nodal positions are pinned.
Definition: SolidProblem.h:113

References INFO, logger, and n.

Referenced by main().

◆ pinBoundary()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::pinBoundary ( unsigned  b)
inline
227  {
228  logger(INFO, "Pinning nodes on boundary %", b);
229  isPinned_ = [b](SolidNode *n, unsigned d) {
230  return n->is_on_boundary(b);
231  };
232  }

References INFO, logger, and n.

Referenced by main().

◆ prepareForSolve()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::prepareForSolve ( )
inline
  • Asserts all parameters are set
  • Assign pointers for all elements
  • Builds global mesh
  • Pins nodes
  • Pins solid pressure
  • Assigns equation numbers

This function should be called before solve(), after creating the mesh, defining body force and constitutive law.

323  {
324  // check certain values are set
325  logger.assert_always(!name_.empty(), "Set name via setName(..)");
326  logger.assert_always(elasticModulus_>0, "Set elastic modulus via setElasticModulus(..)");
327  logger.assert_always(density_>0, "Set density via setDensity(..)");
328  logger.assert_always(solid_mesh_pt(), "Set solid mesh via e.g. setSolidCubicMesh(..)");
329 
330  // Assign constitutive_law_pt and body_force_fct_pt of each element
331  logger(INFO, "Assign constitutive_law, body_force, density to all elements");
332  for (unsigned i = 0; i < solid_mesh_pt()->nelement(); i++)
333  {
334  //Cast to a solid element
335  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(solid_mesh_pt()->element_pt(i));
336  // Set the constitutive law
337  el_pt->constitutive_law_pt() = constitutive_law_pt;
338  // Set body force
339  el_pt->body_force_fct_pt() = body_force_fct;
340  // Set density
341  el_pt->lambda_sq_pt() = &density_;
342  }
343 
344  // Construct the traction element mesh
345  // Traction_mesh_pt = new SolidMesh;
346  // create_traction_elements();
347 
348  //Build combined "global" mesh
349  add_sub_mesh(solid_mesh_pt());
350  if (traction_mesh_pt()) {
351  add_sub_mesh(traction_mesh_pt());
352  logger(INFO, "Built global mesh from solid mesh and traction mesh");
353  } else {
354  logger(INFO, "Built global mesh from solid mesh");
355  }
356  build_global_mesh();
357 
358  //Pin nodes
359  if (isPinned_)
360  {
361  for ( unsigned n = 0; n < solid_mesh_pt()->nnode(); n++ )
362  {
363  SolidNode* n_pt = solid_mesh_pt()->node_pt( n );
364  //Pin all nodes
365  for ( unsigned i = 0; i < 3; i++ )
366  {
367  if ( isPinned_( n_pt, i ) )
368  {
369  n_pt->pin_position( i );
370  }
371  else
372  {
373  n_pt->unpin_position( i );
374  }
375  }
376  }
377  }
378  countPinned();
379 
380  // Pin the redundant solid pressures (if any)
381  PVDEquationsBase<3>::pin_redundant_nodal_solid_pressures(
382  solid_mesh_pt()->element_pt());
383  logger(INFO, "Pinned redundant nodal solid pressures");
384 
385  // Attach the boundary conditions to the mesh
386  unsigned n_eq = assign_eqn_numbers();
387  logger(INFO, "Assigned % equation numbers", n_eq);
388  }
SolidMesh *& traction_mesh_pt()
Get function for the traction mesh pointer.
Definition: SolidProblem.h:287
void(* body_force_fct)(const double &time, const Vector< double > &xi, Vector< double > &b)
Pointer to the body force function.
Definition: SolidProblem.h:101

References constants::i, INFO, logger, and n.

Referenced by main().

◆ removeOldFiles()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::removeOldFiles ( ) const
inline

Removes old output files (vtu) with the same name as this problem.

471  {
472  for (int i = 0; true; ++i)
473  {
474  std::string fileName = name_ + "FEM_" + std::to_string(i) + ".vtu";
475  if (remove(fileName.c_str())) break;
476  std::cout << "Deleted " << fileName << '\n';
477  }
478  }

References constants::i.

◆ saveSolidMesh()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::saveSolidMesh ( )
inline
558  {
559  auto solid_cubic_mesh_pt = dynamic_cast<SolidCubicMesh*>(solid_mesh_pt());
560  if (not solid_cubic_mesh_pt) {
561  logger(INFO,"Solid mesh can only be stored if based on cubic constructor");
562  return;
563  }
564 
565  std::ofstream mesh(name_ + ".mesh");
566  mesh << solid_cubic_mesh_pt->nx() << ' ';
567  mesh << solid_cubic_mesh_pt->ny() << ' ';
568  mesh << solid_cubic_mesh_pt->nz() << '\n';
569  for (int i = 0; i < solid_mesh_pt()->nnode(); ++i)
570  {
571  SolidNode* n = solid_mesh_pt()->node_pt(i);
572  for (int j = 0; j < 3; ++j)
573  {
574  mesh << n->xi(j) << ' ' << n->x(j) << ' ' << n->position_is_pinned(j) << ' ';
575  }
576  mesh << '\n';
577  }
578  logger(INFO, "Saved %x%x% mesh to %.mesh",solid_cubic_mesh_pt->nx(), solid_cubic_mesh_pt->ny(), solid_cubic_mesh_pt->nz(),name_);
579  }

References constants::i, INFO, logger, and n.

Referenced by main().

◆ setBodyForceAsGravity()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setBodyForceAsGravity ( double  gravity = 9.8)
inline

set function for body_force_pt

206  {
207  logger(INFO, "Setting oomph-gravity in z-direction");
208  gravity_ = gravity;
209  // define a static body force
210  static double& Density = density_;
211  static double& Gravity = gravity_;
212  body_force_fct = [](const double& time, const Vector<double>& xi, Vector<double>& b) {
213  b[0] = 0.0;
214  b[1] = 0.0;
215  b[2] = -Gravity * Density;
216  };
217  }

References INFO, and logger.

Referenced by main().

◆ setDensity()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setDensity ( double  density)
inline

set function for density_

193  {
194  density_ = density;
195  logger(INFO, "Density: %", density_);
196  }

References INFO, and logger.

Referenced by main().

◆ setDissipation()

template<class ELEMENT_TYPE >
std::enable_if<std::is_base_of<RefineableQDPVDElement<3,2>, ELEMENT>::value, void> SolidProblem< ELEMENT_TYPE >::setDissipation ( double  dissipation)
inline

Sets the dissipation coefficient for all elements.

249  {
250  for (int i = 0; i < solid_mesh_pt()->nelement(); ++i)
251  {
252  dynamic_cast<RefineableQDPVDElement<3,2>*>(solid_mesh_pt()->element_pt(i))->setDissipation(dissipation);
253  }
254  }
std::enable_if< std::is_base_of< RefineableQDPVDElement< 3, 2 >, ELEMENT >::value, void > setDissipation(double dissipation)
Sets the dissipation coefficient for all elements.
Definition: SolidProblem.h:248
Definition: RefineableQDPVDElement.h:19

References constants::i.

◆ setElasticModulus()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setElasticModulus ( double  elasticModulus)
inline

set function for elasticModulus_

154  {
155  elasticModulus_ = elasticModulus;
156  logger(INFO, "Elastic Modulus: %", elasticModulus_);
157  }

References INFO, and logger.

Referenced by main().

◆ setIsPinned()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setIsPinned ( std::function< bool(SolidNode *, unsigned)>  isPinned)
inline

set function for isPinned_

221  {
222  isPinned_ = std::move(isPinned);
223  logger(INFO, "Setting which positions on which nodes are pinned");
224  }

References INFO, and logger.

◆ setMaxNewtonIterations()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setMaxNewtonIterations ( unsigned  Max_newton_iterations)
inline

set function for Max_newton_iterations

264  {
265  this->Max_newton_iterations = Max_newton_iterations;
266  }

◆ setName()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setName ( const std::string &  name)
inline

set function for name_

141  {
142  name_ = name;
143  logger(INFO, "Name: %", name_);
144  }

References INFO, logger, and units::name.

Referenced by main().

◆ setNewtonSolverTolerance()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setNewtonSolverTolerance ( double  Newton_solver_tolerance)
inline

set function for Newton_solver_tolerance

258  {
259  this->Newton_solver_tolerance = Newton_solver_tolerance;
260  }

Referenced by main().

◆ setOomphGravity()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setOomphGravity ( double  gravity)
inline

set function for elasticModulus_

167  {
168  gravity_ = gravity;
169  logger(INFO, "Elastic Modulus: %", gravity_);
170  }

References INFO, and logger.

◆ setOomphTimeStep()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setOomphTimeStep ( double  dt)
inline

set function for time step

269  {
270  time_pt()->initialise_dt(dt);
271  }

Referenced by main().

◆ setPoissonRatio()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setPoissonRatio ( double  poissonRatio)
inline

set function for poissonRatio_

180  {
181  poissonRatio_ = poissonRatio;
182  logger(INFO, "Poisson Ratio: %", poissonRatio_);
183  }

References INFO, and logger.

◆ setSolidCubicMesh()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::setSolidCubicMesh ( const unsigned &  nx,
const unsigned &  ny,
const unsigned &  nz,
const double xMin,
const double xMax,
const double yMin,
const double yMax,
const double zMin,
const double zMax 
)
inline
549  {
550  // Create mesh
551  solid_mesh_pt() = new SolidCubicMesh(nx, ny, nz, xMin, xMax, yMin, yMax, zMin, zMax, time_stepper_pt());
552 
553  logger(INFO, "Created %x%x% cubic mesh on domain [%,%]x[%,%]x[%,%]",
554  nx, ny, nz, xMin, xMax, yMin, yMax, zMin, zMax);
555  }

References INFO, and logger.

Referenced by main().

◆ solid_mesh_pt() [1/2]

template<class ELEMENT_TYPE >
SolidMesh*& SolidProblem< ELEMENT_TYPE >::solid_mesh_pt ( )
inline

Get function for the solid mesh pointer.

284 { return Solid_mesh_pt; }

◆ solid_mesh_pt() [2/2]

template<class ELEMENT_TYPE >
SolidMesh* const& SolidProblem< ELEMENT_TYPE >::solid_mesh_pt ( ) const
inline

Get function for the solid mesh pointer.

290 { return Solid_mesh_pt; }

◆ solveSteady()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::solveSteady ( )
inline

Solves a steady problem.

412  {
413  logger.assert_always(mesh_pt(), "Mesh pointer not set; did you call prepareForSolve?");
414  logger(INFO, "Solve steady-state problem");
416  newton_solve();
418  writeToVTK();
419  saveSolidMesh();
420  }
void saveSolidMesh()
Definition: SolidProblem.h:557
virtual void actionsBeforeSolve()
Definition: SolidProblem.h:422
void writeToVTK()
Definition: SolidProblem.h:610
virtual void actionsAfterSolve()
Definition: SolidProblem.h:424

References INFO, and logger.

Referenced by main().

◆ solveUnsteady()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::solveUnsteady ( double  timeMax,
double  dt,
unsigned  saveCount = 10 
)
inline

Solves an unsteady problem.

432  {
433  logger.assert_always(mesh_pt(), "Mesh pointer not set; did you call prepareForSolve?");
434 
435  std::cout << "Solving oomph with dt=" << dt << " until timeMax=" << timeMax << std::endl;
436 
437  // Setup initial conditions. Default is a non-impulsive start
438  assign_initial_values_impulsive(dt);
440 
441  // This is the main loop over advancing time
442  double& time = time_stepper_pt()->time();
443  unsigned count = 0;
444  const unsigned countMax = std::ceil(timeMax/dt);
445  while (time < timeMax)
446  {
447  logger(INFO, "Time %s of %s (% of %)", time, timeMax, count, countMax);
449  // solve the oomphProb for one time step (this also increments time)
450  adaptive_unsteady_newton_solve(dt,2e-7);
451  // increase count
452  count++;
453  // write outputs of the oomphProb
454  if (count % saveCount == 0 or time + dt > timeMax) {
455  writeToVTK();
456  }
457  // abort if problem
458  if (Max_res.back()==0) {
459  logger(WARN,"Maximum residual is 0; aborting time loop");
460  break;
461  }
462  }
463  saveSolidMesh();
465  }
@ WARN
virtual void actionsBeforeOomphTimeStep()
Definition: SolidProblem.h:426

References INFO, logger, and WARN.

Referenced by main().

◆ traction_mesh_pt() [1/2]

template<class ELEMENT_TYPE >
SolidMesh*& SolidProblem< ELEMENT_TYPE >::traction_mesh_pt ( )
inline

Get function for the traction mesh pointer.

287 { return Traction_mesh_pt; }
SolidMesh * Traction_mesh_pt
Pointer to mesh of traction elements.
Definition: SolidProblem.h:110

◆ traction_mesh_pt() [2/2]

template<class ELEMENT_TYPE >
SolidMesh* const& SolidProblem< ELEMENT_TYPE >::traction_mesh_pt ( ) const
inline

Get function for the traction mesh pointer.

293 { return Traction_mesh_pt; }

◆ writeToVTK()

template<class ELEMENT_TYPE >
void SolidProblem< ELEMENT_TYPE >::writeToVTK ( )
inline
611  {
612  //set local coordinates list
613  std::vector<std::vector<double>> sList0;
614  // order of nodes
615  std::vector<unsigned> nList;
616  // vtkFormat for given shape
617  // https://kitware.github.io/vtk-examples/site/VTKFileFormats/
618  unsigned vtkFormat;
619 
620  // define differently for tetrahedral and hexahedral (quad) elements
621  if (std::is_base_of<SolidTElement<3, 2>, ELEMENT>::value)
622  {
623  vtkFormat = 10; // TETRA
624 
625  sList0 = {
626  {1, 0, 0},
627  {0, 1, 0},
628  {0, 0, 1},
629  {0, 0, 0}
630  };
631 
632  nList = {0, 1, 2, 3};
633  }
634  else if (std::is_base_of<SolidQElement<3, 2>, ELEMENT>::value)
635  {
636  vtkFormat = 12; // HEXAHEDRON
637 
638  sList0 = {
639  {-1, -1, -1},
640  {-1, -1, +1},
641  {-1, +1, +1},
642  {-1, +1, -1},
643  {+1, -1, -1},
644  {+1, -1, +1},
645  {+1, +1, +1},
646  {+1, +1, -1}
647  };
648  // order of nodes
649  nList = {0, 4, 6, 2, 1, 5, 7, 3};
650  } else {
651  logger(ERROR,"Element type unknown");
652  }
653 
654  // convert to Vector
655  std::vector<Vector<double>> sList;
656  for (auto s0 : sList0)
657  {
658  Vector<double> s(3);
659  s[0] = s0[0];
660  s[1] = s0[1];
661  s[2] = s0[2];
662  sList.push_back(s);
663  }
664 
665  // number of cells
666  unsigned nel = solid_mesh_pt()->nelement();
667 
668  // set up vtu Points
669  struct Point
670  {
671  Vector<double> coordinate;
672  struct Data
673  {
674  std::string name;
675  std::vector<double> value;
676  };
677  std::vector<Data> data;
678  };
679  std::vector<Point> points;
680  points.reserve(nel * sList.size());
681 
682  // set up vtu Cells
683  struct Cell
684  {
685  std::vector<unsigned long> connectivity;
686  unsigned long offset;
687  unsigned type;
688  };
689  std::vector<Cell> cells;
690  points.reserve(nel);
691 
692  //for all elements
693  for (unsigned e = 0; e < nel; e++)
694  {
695  // Get pointer to element
696  auto el_pt = dynamic_cast<ELEMENT*>(
697  solid_mesh_pt()->element_pt(e));
698 
699  std::vector<unsigned long> connectivity;
700  unsigned ix = 0; // node index
701  for (const auto& s : sList)
702  {
703  // pointer to node info
704  auto n = nList[ix];
705  auto n_pt = dynamic_cast<SolidNode*>(el_pt->node_pt(n));
706  // Get Eulerian coordinates
707  Vector<double> x(3);
708  el_pt->interpolated_x(s, x);
709  // get velocity
710  Vector<double> dxdt(3);
711  el_pt->interpolated_dxdt(s, 1, dxdt);
712  // get displacement
713  Vector<double> xi(3);
714  el_pt->interpolated_xi(s, xi);
715  // getBodyForce
716  Vector<double> body_force(3);
717  auto bodyForceFct = el_pt->body_force_fct_pt();
718  if (bodyForceFct) bodyForceFct(time(), xi, body_force);
719  // get stress/pressure
720  DenseMatrix<double> sigma(3,3);
721  el_pt->get_stress(s, sigma);
722  double pressure = (sigma(0,0)+sigma(1,1)+sigma(2,2))/3;
723  // get strain
724  DenseMatrix<double> strain(3,3);
725  el_pt->get_strain(s, strain);
726  // get velocity (fails)
727  std::vector<double> dudt
728  {el_pt->dnodal_position_dt(n, 0),
729  el_pt->dnodal_position_dt(n, 1),
730  el_pt->dnodal_position_dt(n, 2)};
731  // get coupling force
732  std::vector<double> coupledForce(3,0.0);
733  double couplingWeight = 0.0;
734  auto c_el_pt = dynamic_cast<ScaleCoupledElement<RefineableQDPVDElement<3, 2>>*>(el_pt);
735  if (c_el_pt != nullptr) {
736  double nnode = c_el_pt->nnode();
737  double dim = c_el_pt->dim();
738  Shape psi(nnode);
739  c_el_pt->shape(s, psi);
740  for (int n=0; n<nnode; ++n) {
741  for (int d=0; d<dim; ++d) {
742  coupledForce[d] += c_el_pt->get_coupling_residual(n, d) * psi(n);
743  }
744  couplingWeight += c_el_pt->get_coupling_weight(n) * psi(n);
745  }
746  }
747  // get boundary (works)
748  std::set<unsigned>* boundaries_pt;
749  n_pt->get_boundaries_pt(boundaries_pt);
750  double b = boundaries_pt ? *boundaries_pt->begin() : -1;
751  std::vector<double> pin {(double) n_pt->position_is_pinned(0),
752  (double) n_pt->position_is_pinned(1),
753  (double) n_pt->position_is_pinned(2)};
754  points.push_back(
755  {x, {
756  {"Velocity", {dxdt[0], dxdt[1], dxdt[2]}},
757  {"Displacement", {x[0] - xi[0], x[1] - xi[1], x[2] - xi[2]}},
758  {"BodyForce", {body_force[0], body_force[1], body_force[2]}},
759  {"Pin", pin},
760  {"Velocity2", dudt},
761  {"CoupledForce",coupledForce},
762  {"CouplingWeight", {couplingWeight}},
763 // {"Pressure", {pressure}},
764 // {"Stress", {sigma(0,0), sigma(0,1), sigma(0,2),
765 // sigma(1,0), sigma(1,1), sigma(1,2),
766 // sigma(2,0), sigma(2,1), sigma(2,2)}},
767 // {"Strain", {strain(0,0), strain(0,1), strain(0,2),
768 // strain(1,0), strain(1,1), strain(1,2),
769 // strain(2,0), strain(2,1), strain(2,2)}},
770  //\todo undo this comment, but it causes a problem in the output
771  {"Boundary", {b}}
772  }});
773 
774  // add to connectivity
775  connectivity.push_back(points.size() - 1);
776  ix++;
777  }
778  cells.push_back({connectivity, points.size(), vtkFormat});
779  }
780 
781  //open vtk file
782  static unsigned count = 0;
783 #ifdef OOMPH_HAS_MPI
784  std::string vtkFileName = name_ + "FEM_" + std::to_string(OOMPH_MPI_PROCESSOR_ID) + "_" + std::to_string(count++) + ".vtu";
785 #else
786  std::string vtkFileName = name_ + "FEM_" + std::to_string(count++) + ".vtu";
787 #endif
788 
789  std::ofstream vtk(vtkFileName);
790 
791  //write vtk file
792  vtk << "<?xml version=\"1.0\"?>\n"
793  "<!-- time " << getOomphTime() << "-->\n"
794  "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n"
795  "<UnstructuredGrid>\n"
796  "<Piece NumberOfPoints=\""
797  << points.size()
798  << "\" NumberOfCells=\""
799  << cells.size()
800  << "\">\n"
801  "\n"
802  "<Points>\n"
803  "<DataArray type=\"Float32\" Name=\"Position\" NumberOfComponents=\"3\" format=\"ascii\">\n";
804  for (auto point : points)
805  {
806  vtk << point.coordinate[0] << " " << point.coordinate[1] << " " << point.coordinate[2] << "\n";
807  }
808  vtk << "</DataArray>\n"
809  "</Points>\n\n";
810  if (not points.empty())
811  {
812  vtk << "<PointData Vectors=\"vector\">\n";
813  for (int i = 0; i < points.front().data.size(); ++i)
814  {
815  auto data = points.front().data[i];
816  vtk << R"(<DataArray type="Float32" Name=")"
817  << points.front().data[i].name
818  << R"(" NumberOfComponents=")"
819  << points.front().data[i].value.size()
820  << R"(" format="ascii">)" << "\n";
821  for (const Point& point : points)
822  {
823  for (const auto& value : point.data[i].value)
824  {
825  vtk << value << " ";
826  }
827  }
828  vtk << "\n"
829  "</DataArray>\n";
830  }
831  vtk << "</PointData>\n"
832  "\n";
833  }
834  vtk << "<Cells>\n"
835  "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
836  for (const Cell& cell : cells)
837  {
838  for (auto point : cell.connectivity)
839  {
840  vtk << point << " ";
841  }
842  }
843  vtk << "</DataArray>\n"
844  "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n";
845  for (const Cell& cell : cells)
846  {
847  vtk << cell.offset << " ";
848  }
849  vtk << "</DataArray>\n"
850  "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n";
851  for (const Cell& cell : cells)
852  {
853  vtk << cell.type << " ";
854  }
855  vtk << "\n"
856  "</DataArray>\n"
857  "</Cells>\n"
858  "\n"
859  "</Piece>\n"
860  "</UnstructuredGrid>\n"
861  "</VTKFile>\n";
862  vtk.close();
863  logger(INFO, "Written %", vtkFileName);
864  }
#define OOMPH_MPI_PROCESSOR_ID
Definition: SolidProblem.h:62
double getOomphTime() const
get function for current time
Definition: SolidProblem.h:279
Wrapper class for solid elements to be coupled with discrete solid particles on the surfaces.
Definition: ScaleCoupledElement.h:37

References ERROR, constants::i, INFO, logger, n, units::name, and OOMPH_MPI_PROCESSOR_ID.

Member Data Documentation

◆ body_force_fct

template<class ELEMENT_TYPE >
void(* SolidProblem< ELEMENT_TYPE >::body_force_fct) (const double &time, const Vector< double > &xi, Vector< double > &b) = nullptr
protected

Pointer to the body force function.

◆ constitutive_law_pt

template<class ELEMENT_TYPE >
ConstitutiveLaw* SolidProblem< ELEMENT_TYPE >::constitutive_law_pt = nullptr
protected

Pointer to constitutive law (should be set in constructor)

◆ density_

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::density_ = 0
protected

Density.

◆ elasticModulus_

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::elasticModulus_ = 0
protected

Elastic modulus (set via setSolid)

◆ gravity_

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::gravity_ = 0
protected

Density.

◆ isPinned_

template<class ELEMENT_TYPE >
std::function<bool(SolidNode*, unsigned)> SolidProblem< ELEMENT_TYPE >::isPinned_
protected

Function to determine which nodal positions are pinned.

◆ name_

template<class ELEMENT_TYPE >
std::string SolidProblem< ELEMENT_TYPE >::name_
protected

◆ poissonRatio_

template<class ELEMENT_TYPE >
double SolidProblem< ELEMENT_TYPE >::poissonRatio_ = 0
protected

Poisson's ratio (set via setSolid)

◆ Solid_mesh_pt

template<class ELEMENT_TYPE >
SolidMesh* SolidProblem< ELEMENT_TYPE >::Solid_mesh_pt = nullptr
protected

Pointer to solid mesh.

◆ Traction_mesh_pt

template<class ELEMENT_TYPE >
SolidMesh* SolidProblem< ELEMENT_TYPE >::Traction_mesh_pt = nullptr
protected

Pointer to mesh of traction elements.


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