SCoupling< M, O > Class Template Reference

#include <SCoupling.h>

+ Inheritance diagram for SCoupling< M, O >:

Classes

struct  SCoupledElement
 

Public Types

typedef O::ELEMENT ELEMENT
 

Public Member Functions

 SCoupling ()=default
 
void solveSurfaceCoupling ()
 
void solveSurfaceCoupling (unsigned nStep)
 
void solveSurfaceCouplingFixedSolid ()
 
TriangleWallcreateTriangleWall (std::array< Vec3D, 3 > vertex)
 
void updateTriangleWall (TriangleWall *&wall, std::array< Vec3D, 3 > vertex)
 
void computeOneTimeStepForSCoupling (const unsigned &nStepsMercury)
 
void createDPMWallsFromFiniteElems ()
 
void updateDPMWallsFromFiniteElems ()
 
void updateTractionOnFiniteElems ()
 
bool computeSCouplingForcesFromTriangles (ELEMENT *const elem_pt, const unsigned &nTriangles, unsigned &wallID, Vector< Vector< double > > &nodalCouplingForces)
 
void getElementBoundingBox (ELEMENT *&elem_pt, Vec3D &min, Vec3D &max)
 
void getSCoupledElements ()
 
void coupleBoundary (unsigned b)
 
void coupleBoundaries (std::vector< unsigned > b)
 
void disableLogSurfaceCoupling ()
 
- 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 ()
 

Private Attributes

Vector< SCoupledElementsCoupledElements_
 List of surface-coupled elements. More...
 
std::vector< TriangleWall * > triangleWalls_
 List of triangle walls used for the surface coupling. More...
 
std::vector< unsigned > coupledBoundaries_
 
bool logSurfaceCoupling = true
 

Member Typedef Documentation

◆ ELEMENT

template<class M , class O >
typedef O::ELEMENT SCoupling< M, O >::ELEMENT

Constructor & Destructor Documentation

◆ SCoupling()

template<class M , class O >
SCoupling< M, O >::SCoupling ( )
default

Member Function Documentation

◆ computeOneTimeStepForSCoupling()

template<class M , class O >
void SCoupling< M, O >::computeOneTimeStepForSCoupling ( const unsigned &  nStepsMercury)
inline

Solve surface-coupled problem for one oomph time step (n mercury time steps)

238  {
239  auto t0 = std::chrono::system_clock::now();
241  auto t1 = std::chrono::system_clock::now();
242  BaseCoupling<M,O>::solveMercury(nStepsMercury);
243  auto t2 = std::chrono::system_clock::now();
245  auto t3 = std::chrono::system_clock::now();
247  auto t4 = std::chrono::system_clock::now();
248  if (logSurfaceCoupling) logger(INFO, "time % Elapsed time: FEM->DEM %, DEM %, DEM->FEM %, FEM %", M::getTime(),
249  (t1-t0).count(), (t2-t1).count(), (t3-t2).count(), (t4-t3).count());
250  }
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ INFO
void solveOomph()
Definition: BaseCoupling.h:141
void solveMercury(unsigned long nt)
Definition: BaseCoupling.h:150
void updateTractionOnFiniteElems()
Definition: SCoupling.h:345
void updateDPMWallsFromFiniteElems()
Definition: SCoupling.h:294
bool logSurfaceCoupling
Definition: SCoupling.h:708

References INFO, logger, BaseCoupling< M, O >::solveMercury(), and BaseCoupling< M, O >::solveOomph().

◆ computeSCouplingForcesFromTriangles()

template<class M , class O >
bool SCoupling< M, O >::computeSCouplingForcesFromTriangles ( ELEMENT *const  elem_pt,
const unsigned &  nTriangles,
unsigned &  wallID,
Vector< Vector< double > > &  nodalCouplingForces 
)
inline

Computes nodal scoupling forces from triangles.

Parameters
[in]elem_ptpointer to the element for which forces are computed
[in]nTrianglesnumber of triangles (equals number of nodes of the traction element)
[in,out]wallID
[out]nodalCouplingForces[l]= sum_i f * psi(l,cp) for each shape function l, interaction i, evaluated at contact point cp
Returns
410  {
411  // whether the element is coupled to particles
412  bool elemIsCoupled = false;
413 
414  // get number of nodes in the bulk element
415  const unsigned nnode = elem_pt->nnode();
416  // get dimension of the problem
417  const unsigned dim = elem_pt->dim();
418  // initialize the coupling force vector
419  nodalCouplingForces.resize(nnode, Vector<double>(dim, 0.0));
420 
421  // get number of TriangleWalls created from an oomph face element
422  unsigned n = 0;
423  // loop over TriangleWalls
424  while (n++ < nTriangles)
425  {
426  // get pointer to the wall
427  TriangleWall* w = triangleWalls_[wallID++];
428 
429  // skip the wall if not interactions with it
430  if (w->getInteractions().size() == 0) continue;
431 
432  // Set up memory for the shape/test functions
433  Shape psi(nnode);
434 
435  // loop over interactions with the wall
436  for (auto inter : w->getInteractions())
437  {
438  if (!inter->getForce().isZero())
439  {
440  // scale contact point and forces from DEM to FEM units
441  Vec3D xc = inter->getContactPoint();
442  Vec3D fc = inter->getForce();
443 
444  Vector<double> x(3, 0.0), f(3, 0.0);
445  x[0] = xc.getX();
446  x[1] = xc.getY();
447  x[2] = xc.getZ();
448  f[0] = fc.getX();
449  f[1] = fc.getY();
450  f[2] = fc.getZ();
451 
452  // get the local coordinate s if xc is located in the finite element
453  Vector<double> s(3, 0.0);
454  GeomObject* geom_obj_pt = 0;
455  elem_pt->locate_zeta(x, geom_obj_pt, s);
456 
457  // Get shape/test fcts
458  elem_pt->shape(s, psi);
459  // Loop over the test functions
460  for (unsigned l = 0; l < nnode; l++)
461  {
462  //Loop over the force components
463  for (unsigned i = 0; i < dim; i++)
464  {
465  // add contribution to the nodal coupling force
466  nodalCouplingForces[l][i] += f[i] * psi(l);
467  }
468  }
469  // set the flag to true
470  elemIsCoupled = true;
471  }
472  }
473  }
474  return elemIsCoupled;
475  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:32
const std::vector< BaseInteraction * > & getInteractions() const
Returns a list of interactions which belong to this interactable.
Definition: BaseInteractable.h:277
std::vector< TriangleWall * > triangleWalls_
List of triangle walls used for the surface coupling.
Definition: SCoupling.h:700
A TriangleWall is convex polygon defined as an intersection of InfiniteWall's.
Definition: TriangleWall.h:57
Definition: Vector.h:51
Mdouble getZ() const
Definition: Vector.h:408
Mdouble getX() const
Definition: Vector.h:402
Mdouble getY() const
Definition: Vector.h:405
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51

References BaseInteractable::getInteractions(), Vec3D::getX(), Vec3D::getY(), Vec3D::getZ(), constants::i, and n.

◆ coupleBoundaries()

template<class M , class O >
void SCoupling< M, O >::coupleBoundaries ( std::vector< unsigned >  b)
inline
686  {
687  coupledBoundaries_ = std::move(b);
688  }
std::vector< unsigned > coupledBoundaries_
Definition: SCoupling.h:706

◆ coupleBoundary()

template<class M , class O >
void SCoupling< M, O >::coupleBoundary ( unsigned  b)
inline

◆ createDPMWallsFromFiniteElems()

template<class M , class O >
void SCoupling< M, O >::createDPMWallsFromFiniteElems ( )
inline
253  {
254  // loop over bulk elements at boundary
255  for (auto sCoupledElement : sCoupledElements_)
256  {
257  // loop over nodes in element at boundary
258  Vector<Vector<double*> > position_pt = sCoupledElement.surface_node_position_pt;
259  // reordering vertices of oomph face element (default = {0,1,3,2})
260  swap(position_pt[2], position_pt[3]);
261 
262  // get global coordinate at the center
263  Vec3D center;
264  Vector<double> x(3, 0.0);
265  sCoupledElement.bulk_elem_pt->interpolated_x(sCoupledElement.center_local, x);
266  center.setX(x[0]);
267  center.setY(x[1]);
268  center.setZ(x[2]);
269 
270  // create TriangleWalls from oomph face element
271  const unsigned nTriangles = position_pt.size();
272  unsigned n = 0;
273  while (n < nTriangles)
274  {
275  // get vertices of TriangleWall (multiply vertex position with the length scale of the O)
276  std::array<Vec3D, 3> vertex;
277  // one vertex at the center
278  vertex[0] = center;
279  // two vertices from the O<element,TIMESTEPPER>
280  vertex[1] = Vec3D(*position_pt[0][0],*position_pt[0][1],*position_pt[0][2]);
281  vertex[2] = Vec3D(*position_pt[1][0],*position_pt[1][1],*position_pt[1][2]);
282 
283  // create triangle facet
284  TriangleWall* w = createTriangleWall(vertex);
285  triangleWalls_.push_back(w);
286 
287  // rotate forward by one element
288  rotate(position_pt.begin(), position_pt.begin() + 1, position_pt.end());
289  n++;
290  }
291  }
292  }
TriangleWall * createTriangleWall(std::array< Vec3D, 3 > vertex)
Definition: SCoupling.h:192
Vector< SCoupledElement > sCoupledElements_
List of surface-coupled elements.
Definition: SCoupling.h:697
void setY(Mdouble y)
Definition: Vector.h:396
void setZ(Mdouble z)
Definition: Vector.h:399
void setX(Mdouble x)
Definition: Vector.h:393

References n, Vec3D::setX(), Vec3D::setY(), and Vec3D::setZ().

◆ createTriangleWall()

template<class M , class O >
TriangleWall* SCoupling< M, O >::createTriangleWall ( std::array< Vec3D, 3 >  vertex)
inline

Create a triangle wall from a set of vertices

Parameters
vertex
Todo:
Why we need to set GroupID. I don't remember anymore.
Returns
193  {
194  TriangleWall wall;
195  auto species = M::speciesHandler.getObject(0);
196  wall.setSpecies(species);
197  wall.setGroupId(100);
198  wall.setVertices(vertex[0], vertex[1], vertex[2]);
199  auto w = M::wallHandler.copyAndAddObject(wall);
200  return w;
201  }
void setGroupId(unsigned groupId)
Definition: BaseObject.h:131
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:169
void setVertices(Vec3D A, Vec3D B, Vec3D C)
Sets member variables such that the wall represents a triangle with vertices A, B,...
Definition: TriangleWall.cc:165

References BaseObject::setGroupId(), BaseWall::setSpecies(), and TriangleWall::setVertices().

◆ disableLogSurfaceCoupling()

template<class M , class O >
void SCoupling< M, O >::disableLogSurfaceCoupling ( )
inline
690  {
691  logSurfaceCoupling = false;
692  }

◆ getElementBoundingBox()

template<class M , class O >
void SCoupling< M, O >::getElementBoundingBox ( ELEMENT *&  elem_pt,
Vec3D min,
Vec3D max 
)
inline

needed in computeSCouplingForcesFromCG to decide which particles affect a given element

568  {
569  // three arrays that contain the x, y and z coordinates of the bulk element
570  Vector<double> listOfCoordX;
571  Vector<double> listOfCoordY;
572  Vector<double> listOfCoordZ;
573 
574  // get the x, y and z coordinates of the bulk element
575  const unsigned nnode = elem_pt->nnode();
576  for (unsigned n = 0; n < nnode; n++)
577  {
578  listOfCoordX.push_back(elem_pt->node_pt(n)->x(0));
579  listOfCoordY.push_back(elem_pt->node_pt(n)->x(1));
580  listOfCoordZ.push_back(elem_pt->node_pt(n)->x(2));
581  }
582 
583  // get the bounding box of the bulk element
584  min.X = *min_element(listOfCoordX.begin(), listOfCoordX.end());
585  min.Y = *min_element(listOfCoordY.begin(), listOfCoordY.end());
586  min.Z = *min_element(listOfCoordZ.begin(), listOfCoordZ.end());
587  max.X = *max_element(listOfCoordX.begin(), listOfCoordX.end());
588  max.Y = *max_element(listOfCoordY.begin(), listOfCoordY.end());
589  max.Z = *max_element(listOfCoordZ.begin(), listOfCoordZ.end());
590 
591  // extend the bounding box if construct mapping with coarse graining
592  logger.assert_always(M::particleHandler.getLargestParticle(), "No particles detected");
593  //todo does it make sense with -2R here?
594  min -= Vec3D(1.0, 1.0, 1.0) * ( BaseCoupling<M,O>::getCGWidth() - 2 * M::particleHandler.getLargestParticle()->getRadius() );
595  max += Vec3D(1.0, 1.0, 1.0) * ( BaseCoupling<M,O>::getCGWidth() - 2 * M::particleHandler.getLargestParticle()->getRadius() );
596  }
double getCGWidth()
Definition: BaseCoupling.h:176
Mdouble Y
Definition: Vector.h:66
Mdouble Z
Definition: Vector.h:66
Mdouble X
the vector components
Definition: Vector.h:66

References BaseCoupling< M, O >::getCGWidth(), logger, n, Vec3D::X, Vec3D::Y, and Vec3D::Z.

◆ getSCoupledElements()

template<class M , class O >
void SCoupling< M, O >::getSCoupledElements ( )
inline

Get bulk elements along boundaries (SCoupling). The results are stored in sCoupledElements_.

603  {
604  // first clear bulk elements (for refineable problem)
605  sCoupledElements_.clear();
606 
607  // loop over all boundaries
608  for (unsigned b : coupledBoundaries_)
609  {
610  // we only need to couple the elements on the upper boundary
611  logger(INFO,"Coupling boundary %", b);
612 
613  // number of bulk elements adjacent to boundary b
614  unsigned n_element = this->solid_mesh_pt()->nboundary_element(b);
615 
616  // loop over the bulk elements adjacent to boundary b
617  for (unsigned e = 0; e < n_element; e++)
618  {
619  // initialize the struct SCoupledElement
620  SCoupledElement sCoupledElement;
621 
622  // get pointer to the bulk element that is adjacent to boundary b
623  sCoupledElement.bulk_elem_pt = dynamic_cast<ELEMENT*>(this->solid_mesh_pt()->boundary_element_pt(b, e));
624 
625  // get the index of the face of element e along boundary b
626  sCoupledElement.face_index = this->solid_mesh_pt()->face_index_at_boundary(b, e);
627 
628  // create temporary traction element
629  SolidTractionElement<ELEMENT> trac_elem(sCoupledElement.bulk_elem_pt, sCoupledElement.face_index);
630 
631  // number of nodes on traction element
632  unsigned n_node = trac_elem.nnode();
633 
634  // allocate space to store the local coordinates of the traction element's vertices
635  sCoupledElement.surface_node_position_pt.reserve(n_node);
636 
637  // store the local coordinates of the traction element's center
638  sCoupledElement.center_local.resize(3, 0.0);
639 
640  // assign addresses of nodal coordinates to the struct SCoupledElement
641  for (unsigned n = 0; n < n_node; n++)
642  {
643  // store pointer to solid nodes
644  sCoupledElement.node_pt.push_back(dynamic_cast<CoupledSolidNode*>(trac_elem.node_pt(n)));
645 
646  Vector<double> s(2, 0.0);
647  trac_elem.local_coordinate_of_node(n, s);
648 
649  Vector<double> s_bulk(3, 0.0);
650  trac_elem.get_local_coordinate_in_bulk(s, s_bulk);
651 
652  // pass the addresses of x coordinates at boundary to x_pt
653  Vector<double*> x_pt;
654  for (unsigned i = 0; i < 3; i++)
655  {
656  x_pt.push_back(&trac_elem.node_pt(n)->x(i));
657  sCoupledElement.center_local[i] += s_bulk[i];
658  }
659  /*
660  // debug that the addresses of x coordinates are passed to surface_node_position_pt
661  cout << "Address of x(1) of node on traction element\t " << &trac_elem.node_pt(n)->x(0) << endl;
662  int nodeIDOfBulkElement = sCoupledElement.bulk_elem_pt->get_node_number(trac_elem.node_pt(n));
663  cout << "Address of x(1) of the same node in bulk element "
664  << &sCoupledElement.bulk_elem_pt->node_pt(nodeIDOfBulkElement)->x(0) << endl;
665  cout << "Address of x(1) saved in Oomph interface object\t " << x_pt[0] << endl;
666  */
667  // save the addresses of Lagrange coordinates in each element at boundary
668  sCoupledElement.surface_node_position_pt.push_back(x_pt);
669  }
670 
671  // local coordinate of the center of the face element
672  for (unsigned i = 0; i < 3; i++) sCoupledElement.center_local[i] /= n_node;
673 
674  // save bulk elements and surface nodal position pointers
675  sCoupledElements_.push_back(sCoupledElement);
676  }
677  }
678  }
Definition: CoupledSolidNodes.h:16
Definition: SCoupledElement.h:39

References SCoupling< M, O >::SCoupledElement::bulk_elem_pt, SCoupling< M, O >::SCoupledElement::center_local, SCoupling< M, O >::SCoupledElement::face_index, constants::i, INFO, logger, n, SCoupling< M, O >::SCoupledElement::node_pt, and SCoupling< M, O >::SCoupledElement::surface_node_position_pt.

◆ solveSurfaceCoupling() [1/2]

template<class M , class O >
void SCoupling< M, O >::solveSurfaceCoupling ( )
inline
73  {
74  // compute nStep
75  unsigned nStep = O::getOomphTimeStep()/M::getTimeStep();
76  if (nStep==0) {
77  // if oomph has a smaller time step than Mercury
78  logger(INFO, "Set nStep %, change mercuryTimeStep from % to %",
79  nStep, M::getTimeStep(), O::getOomphTimeStep());
80  nStep = 1;
81  M::setTimeStep(O::getOomphTimeStep());
82  } else {
83  logger(INFO, "Set nStep %, change oomphTimeStep from % to %",
84  nStep, O::getOomphTimeStep(), nStep * M::getTimeStep());
85  O::setOomphTimeStep(nStep * M::getTimeStep());
86  }
87 
88  // call solve routine
89  solveSurfaceCoupling(nStep);
90  }
void solveSurfaceCoupling()
Definition: SCoupling.h:72

References INFO, and logger.

Referenced by CoupledBeam::CoupledBeam(), CoupledProblem::CoupledProblem(), and main().

◆ solveSurfaceCoupling() [2/2]

template<class M , class O >
void SCoupling< M, O >::solveSurfaceCoupling ( unsigned  nStep)
inline
94  {
95  // check whether time steps are set
96  logger.assert_always(O::getOomphTimeStep()>0,"Oomph time step not initialised");
97  logger.assert_always(M::getTimeStep()>0,"Mercury time step not initialised");
98  // check whether isCoupled is set
99  logger.assert_always(!coupledBoundaries_.empty(), "isCoupled needs to be set, e.g. via setIsCoupled([](unsigned b) { return b == Boundary::Y_MIN; })");
100 
101  // first part of the Mercury solve routine, containing the actions before time-stepping
102  M::initialiseSolve();
103 
104  // read Mercury time step
105  double mercury_dt = M::getTimeStep();
106  logger(INFO, "Mercury time step: %", mercury_dt);
107 
108  // set oomph time step
109  double oomph_dt = nStep * mercury_dt;
110  logger(INFO, "Oomph time step: %", oomph_dt);
111 
112  // set oomph initial conditions
113  O::assign_initial_values_impulsive(oomph_dt);
114 
115  // read oomph-mesh
116  logger(INFO, "Set up oomph mesh: % elements (% with traction)",
117  O::solid_mesh_pt()?O::solid_mesh_pt()->nelement():0, O::traction_mesh_pt()?O::traction_mesh_pt()->nelement():0);
118 
119  // get list of bulk elements along the surface-coupled boundaries
121 
122  // create DPM triangle walls from bulk finite elements
124 
125  // this is the main loop advancing over time
126  unsigned nDone = 0; //< last written file number
127  while (M::getTime() < M::getTimeMax())
128  {
129  this->actionsBeforeOomphTimeStep();
130  // solve the coupled problem for one time step
132  // write outputs of the oomphProb; this is slaved to the vtk output of Mercury, i.e. an oomph-lib output get written everytime a Mercury vtk file gets written
133  if (M::getParticlesWriteVTK() && M::getVtkWriter()->getFileCounter() > nDone)
134  {
135  O::writeToVTK();
136  nDone = M::getVtkWriter()->getFileCounter();
137  }
138  }
139 
140  // close output files of mercuryProb
141  M::finaliseSolve();
142  }
void getSCoupledElements()
Definition: SCoupling.h:602
void computeOneTimeStepForSCoupling(const unsigned &nStepsMercury)
Definition: SCoupling.h:237
void createDPMWallsFromFiniteElems()
Definition: SCoupling.h:252

References INFO, and logger.

◆ solveSurfaceCouplingFixedSolid()

template<class M , class O >
void SCoupling< M, O >::solveSurfaceCouplingFixedSolid ( )
inline
146  {
147  // first part of the Mercury solve routine, containing the actions before time-stepping
148  O::initialiseSolve();
149  logger.assert_always(!coupledBoundaries_.empty(), "isCoupled needs to be set, e.g. via setIsCoupled([](unsigned b) { return b == Boundary::Y_MIN; })");
150 
151  // read Mercury time step
152  double mercury_dt = M::getTimeStep();
153  logger(INFO, "Mercury time step: %", mercury_dt);
154 
155  // set oomph time step
156  logger(INFO, "Solid position fixed");
157 
158  // set oomph initial conditions
159  O::set_initial_conditions(0);
160 
161  // read oomph-mesh
162  logger(INFO, "Set up oomph mesh: % elements (% with traction)",
163  O::solid_mesh_pt()?O::solid_mesh_pt()->nelement():0, O::traction_mesh_pt()?O::traction_mesh_pt()->nelement():0);
164 
165  // get list of bulk elements along the surface-coupled boundaries
167 
168  // create DPM triangle walls from bulk finite elements
170 
171  // this is the main loop advancing over time
172  unsigned nDone = 0; //< last written file number
173  while (M::getTime() < M::getTimeMax())
174  {
175  O::actionsBeforeOomphTimeStep();
176  M::computeOneTimeStep();
177  //if (getParticlesWriteVTK() && getVtkWriter()->getFileCounter() > nDone) {
178  // writeToVTK();
179  // nDone = getVtkWriter()->getFileCounter();
180  //}
181  }
182  // close output files of mercuryProb
183  M::finaliseSolve();
184  }

References INFO, and logger.

◆ updateDPMWallsFromFiniteElems()

template<class M , class O >
void SCoupling< M, O >::updateDPMWallsFromFiniteElems ( )
inline
295  {
296  // loop over bulk elements at boundary
297  unsigned wallID = 0;
298  for (auto sCoupledElement : sCoupledElements_)
299  {
300  // get members of a SCoupledElement
301  Vector<Vector<double*> > position_pt = sCoupledElement.surface_node_position_pt;
302  Vector<CoupledSolidNode*> node_pt = sCoupledElement.node_pt;
303 
304  // reordering vertices of oomph face element (default = {0,1,3,2})
305  // \todo can we do this swap in getSCoupledElements, then we don't have to create a local copy position_pt all the time?
306  swap(position_pt[2], position_pt[3]);
307 
308  // get global coordinate at the center
309  Vector<double> x(3, 0.0);
310  sCoupledElement.bulk_elem_pt->interpolated_x(sCoupledElement.center_local, x);
311  Vec3D center {x[0], x[1], x[2]};
312 
313  // get number of TriangleWalls per oomph face element
314  const unsigned nTriangles = position_pt.size();
315  unsigned n = 0;
316  while (n < nTriangles)
317  {
318  // get vertices of TriangleWall (multiply vertex position with the length scale of the O)
319  std::array<Vec3D, 3> vertex;
320  // one vertex at the center
321  vertex[0] = center;
322  // two vertices from the O<element,TIMESTEPPER>
323  vertex[1] = Vec3D(*position_pt[0][0],
324  *position_pt[0][1],
325  *position_pt[0][2]);;
326  vertex[2] = Vec3D(*position_pt[1][0],
327  *position_pt[1][1],
328  *position_pt[1][2]);
329 
330  // update vertices of triangle facet (multiply vertex position with the length scale of the O<element,TIMESTEPPER>)
331  updateTriangleWall(triangleWalls_[wallID], vertex);
332 
333  // rotate forward by one element
334  rotate(position_pt.begin(), position_pt.begin() + 1, position_pt.end());
335  rotate(node_pt.begin(), node_pt.begin() + 1, node_pt.end());
336  n++;
337  wallID++;
338  }
339  }
340  }
void updateTriangleWall(TriangleWall *&wall, std::array< Vec3D, 3 > vertex)
Definition: SCoupling.h:208

References n.

◆ updateTractionOnFiniteElems()

template<class M , class O >
void SCoupling< M, O >::updateTractionOnFiniteElems ( )
inline

sets nodal_coupling_residual in each scoupled element

346  {
347  // if construct mapping with FEM basis functions
349  {
350  // tracks the id of the triangle walls (get incremented in computeSCouplingForcesFromTriangles)
351  unsigned wallID = 0;
352 
353  // loop over scoupled elements
354  for (auto sCoupledElement : sCoupledElements_)
355  {
356  // set up memory for nodal coupling force
357  Vector<Vector<double> > nodalCouplingForces;
358  // returns whether the element is coupled to particles
359  bool elemIsCoupled = computeSCouplingForcesFromTriangles(sCoupledElement.bulk_elem_pt,
360  sCoupledElement.surface_node_position_pt.size(),
361  wallID,
362  nodalCouplingForces);
363 
364  // assign nodal coupling force to the element to be used by element::fill_in_contribution_to_residuals(...)
365  sCoupledElement.bulk_elem_pt->set_nodal_coupling_residual(elemIsCoupled, nodalCouplingForces);
366  }
367  }
368  else
369  {
370 // // how many bulk elements in total
371 // unsigned n_element = O::solid_mesh_pt()->nelement();
372 //
373 // // loop over the bulk elements
374 // for (unsigned e = 0; e < n_element; e++)
375 // {
376 // // get pointer to the bulk element
377 // ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(O::solid_mesh_pt()->element_pt(e));
378 //
379 // // set up memory for nodal coupling force
380 // Vector<Vector<double> > nodalCouplingForces;
381 // // whether the element is coupled to particles
382 // bool elemIsCoupled = computeSCouplingForcesFromCG(elem_pt, nodalCouplingForces);
383 //
384 // // assign nodal coupling force to the element to be used by element::fill_in_contribution_to_residuals(...)
385 // elem_pt->set_nodal_coupling_residual(elemIsCoupled, nodalCouplingForces);
386 // }
387 // // reset the sum of the evaluated CG values
388 // for (auto inter : M::interactionHandler)
389 // {
390 // if (inter->getI()->getName() == "TriangleWall")
391 // {
394 // inter->resetTotalCGEval();
395 // }
396 // }
397  }
398  }
Definition: BaseCoupling.h:48
bool computeSCouplingForcesFromTriangles(ELEMENT *const elem_pt, const unsigned &nTriangles, unsigned &wallID, Vector< Vector< double > > &nodalCouplingForces)
Definition: SCoupling.h:408

◆ updateTriangleWall()

template<class M , class O >
void SCoupling< M, O >::updateTriangleWall ( TriangleWall *&  wall,
std::array< Vec3D, 3 >  vertex 
)
inline

Update the position of a triangle wall (used in updateDPMWallsFromFiniteElems)

Todo:

why is it using setPrescribedPosition

we need to set velocity

209  {
210  double time0 = M::getTime();
211  double dTime = O::getOomphTimeStep();
212  std::array<Vec3D,3> vertex0 = wall->getVertices();
213  std::array<Vec3D,3> dVertex = {
214  vertex[0] - vertex0[0],
215  vertex[1] - vertex0[1],
216  vertex[2] - vertex0[2]};
217  wall->setPrescribedPosition( [time0, dTime, vertex0, dVertex, wall]( double time ) {
218  double f = ( time - time0 ) / dTime;
219  std::array<Vec3D, 3> vertex = {
220  vertex0[0] + dVertex[0] * f,
221  vertex0[1] + dVertex[1] * f,
222  vertex0[2] + dVertex[2] * f };
223  wall->setVertices( vertex[0], vertex[1], vertex[2] );
224  //logger(INFO,"p %",vertex[0]);
225  return wall->getPosition();
226  } );
227  Vec3D velocity = (dVertex[0]+dVertex[1]+dVertex[2])/3./dTime;
228  wall->setPrescribedVelocity([velocity] (double time) {
229  //logger(INFO,"v %",velocity);
230  return velocity;
231  });
232  }
void setPrescribedPosition(const std::function< Vec3D(double)> &prescribedPosition)
Allows the position of an infinite mass interactable to be prescribed.
Definition: BaseInteractable.cc:413
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Definition: BaseInteractable.h:218
void setPrescribedVelocity(const std::function< Vec3D(double)> &prescribedVelocity)
Allows the velocity of an infinite mass interactable to be prescribed.
Definition: BaseInteractable.cc:444
std::array< Vec3D, 3 > getVertices() const
Definition: TriangleWall.h:105

References BaseInteractable::getPosition(), TriangleWall::getVertices(), BaseInteractable::setPrescribedPosition(), BaseInteractable::setPrescribedVelocity(), and TriangleWall::setVertices().

Member Data Documentation

◆ coupledBoundaries_

template<class M , class O >
std::vector<unsigned> SCoupling< M, O >::coupledBoundaries_
private

Function to determine whether boundary is coupled. Needs to be set before solveSurfaceCoupling is called.

◆ logSurfaceCoupling

template<class M , class O >
bool SCoupling< M, O >::logSurfaceCoupling = true
private

◆ sCoupledElements_

template<class M , class O >
Vector<SCoupledElement> SCoupling< M, O >::sCoupledElements_
private

List of surface-coupled elements.

◆ triangleWalls_

template<class M , class O >
std::vector<TriangleWall*> SCoupling< M, O >::triangleWalls_
private

List of triangle walls used for the surface coupling.


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