VCoupling.h
Go to the documentation of this file.
1 //Copyright (c) 2013-2023, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #ifndef VOLUME_COUPLING_H
27 #define VOLUME_COUPLING_H
28 
29 #include "BaseCoupling.h"
30 
32 {
33 public:
34  typedef VolumeCoupledElement<OomphProblem::ELEMENT> VELEMENT;
35 
40  {
41  // the pointer to a bulk element
43  // the vector that contains pointers to the coupled particles
44  Vector<BaseParticle*> listOfCoupledParticles;
45  Vector<BaseParticle*> listOfCoupledParticlesExt;
46  // the list of local coordinates of particle centers
47  Vector<Vector<double>> listOfParticleCentersLoc;
48  // the list of shape functions evaluated at particle centers
49  Vector<Vector<double>> listOfShapesAtParticleCenters;
50  };
51 
52  // default constructor
53  VolumeCoupling() = default;
54 
55  // constructor passing in the minimal coupling weight and flag for using the implicit or explicit scheme
56  VolumeCoupling(const double& wMin, const bool& flag)
57  {
58  setExplicit(flag);
59  weightMin_ = wMin;
60  }
61 
66  {
67  //Checks that the current particle is not "fixed"
68  //and hence infinitely massive!
69  if (!CI->isFixed())
70  {
71  // Applying the force due to bodyForce_fct (F = m.g)
72  CI->addForce(getGravity() * CI->getMass());
73  // add particle coupling force 1/w_i*f_i^c to the total force
74  CI->addForce(CI->getCouplingForce());
75  }
76  }
77 
82  void actionsAfterTimeStep() override
83  {
84  // Reset coupling weights and forces to default values
85  for (BaseParticle* p : particleHandler)
86  {
87  if (!p->isFixed())
88  {
89  // reset the coupling force to zero for the next time step
90  p->resetCouplingForce();
91  p->resetInvCouplingWeight();
92  p->resetCouplingMass();
93  }
94  }
95 
96  // check how much equilibrium is broken
97  const auto p = particleHandler.getLastObject();
98  Vec3D f = p->getMass() * getGravity();
99  Vec3D a = p->getForce();
100  if (a.getLengthSquared() > 1.0e-5 * f.getLengthSquared())
101  {
102  logger(INFO, "time %, body force: %, total force: %, position %",
103  getTime(), f, a, p->getPosition());
104  }
105 
106  // check if each particle-wall interaction is unique
107  for (auto p : particleHandler)
108  {
109  unsigned n = 0;
110  for (auto i : p->getInteractions())
111  {
112  if (i->getI()->getName() == "TriangleWall") n++;
113  }
114  if (n > 1)
115  {
116  logger(INFO, "particle % is interacting with % TriangleWalls", p->getId(), n);
117  for (auto i : p->getInteractions())
118  {
119  if (i->getI()->getName() == "TriangleWall")
120  {
121  logger(INFO, "particle % is interacting with wall %", p->getId(), i->getI()->getId());
122  }
123  }
124  }
125  }
126  }
127 
128  // solve volume coupled OomphMercuryProblem using CG for micro-macro mapping
129  void solveVolumeCoupling(const double& width)
130  {
131  // set the coarse-grainning width w.r.t. length scale
132  setCGWidth(width * particleHandler.getLargestParticle()->getRadius());
133  // solve volume coupled OomphMercuryProblem
135  }
136 
137  // solve volume coupled OomphMercuryProblem
139  {
140  initialiseSolve();
141 
142  // get the time step of mercuryProb
143  double mercury_dt = getTimeStep();
144  /*
145  // %TODO use the DPM time scale to non-dimensionalise the OomphProblem<SELEMENT,TIMESTEPPER>?
146  Global_Physical_Variables::timeScale = mercury_dt;
147  */
148 
149  /* // get number of mercury-steps per Oomph-step (critical time step = intrinsic_time or from_eigen)
150  double oomph_dt = 0.2*Global_Physical_Variables::intrinsic_time;
151  double oomph_dt = 0.2*OomphProblem::get_critical_timestep_from_eigen();
152  unsigned nstep = unsigned(floor(oomph_dt/mercury_dt));
153  */
154  unsigned nstep = 1;
155  double oomph_dt = nstep * mercury_dt;
156 
157  // setup initial conditions
158  OomphProblem::set_initial_conditions(oomph_dt);
159  // MercuryBeam::set_initial_conditions();
160 
161  // This is the main loop over advancing time
162  unsigned nDone = 0;
163  while (getTime() < getTimeMax())
164  {
165  // solve the coupled problem for one time step
167  // write outputs of the oomphProb
168  if (getVtkWriter()->getFileCounter() > nDone)
169  {
170  OomphProblem::doc_solution();
171  nDone = getVtkWriter()->getFileCounter();
172  }
173  }
174  // close output files of mercuryProb
175  finaliseSolve();
176  }
177 
178  // solve volume coupled OomphMercuryProblem for one time step
179  void computeOneTimeStepForVCoupling(const unsigned& nstep)
180  {
183  // setNodalWeightOnCubicMesh();
186  if (explicit_)
187  {
189  solve();
190  solveMercury(nstep);
191  }
192  else
193  {
196  solveOomph();
199  solveMercury(nstep - 1);
200  }
201  }
202 
207  {
208  logger(DEBUG, "starting solveHalfTimeStep()");
209 
210  logger(DEBUG, "about to call writeOutputFiles()");
211  writeOutputFiles(); //everything is written at the beginning of the time step!
212 
213  logger(DEBUG, "about to call hGridActionsBeforeIntegration()");
214  hGridActionsBeforeIntegration();
215 
216  //Computes the half-time step velocity and full time step position and updates the particles accordingly
217  logger(DEBUG, "about to call integrateBeforeForceComputation()");
218  integrateBeforeForceComputation();
219  //New positions require the MPI and parallel periodic boundaries to do things
220  logger(DEBUG, "about to call performGhostParticleUpdate()");
221  performGhostParticleUpdate();
222 
226 
227  logger(DEBUG, "about to call checkInteractionWithBoundaries()");
228  checkInteractionWithBoundaries(); // INSERTION boundaries handled
229 
230  logger(DEBUG, "about to call hGridActionsAfterIntegration()");
231  hGridActionsAfterIntegration();
232 
233  // Compute forces
235  // INSERTION/DELETION boundary flag change
236  for (BaseBoundary* b : boundaryHandler)
237  {
238  b->checkBoundaryBeforeTimeStep(this);
239  }
240 
241  logger(DEBUG, "about to call actionsBeforeTimeStep()");
242  actionsBeforeTimeStep();
243 
244  logger(DEBUG, "about to call checkAndDuplicatePeriodicParticles()");
245  checkAndDuplicatePeriodicParticles();
246 
247  logger(DEBUG, "about to call hGridActionsBeforeTimeStep()");
248  hGridActionsBeforeTimeStep();
249  }
250 
255  {
256  //Creates and updates interactions and computes forces based on these
257  logger(DEBUG, "about to call computeAllForces()");
258  computeAllForces();
259 
260  logger(DEBUG, "about to call removeDuplicatePeriodicParticles()");
261  removeDuplicatePeriodicParticles();
262 
263  //Computes new velocities and updates the particles accordingly
264  logger(DEBUG, "about to call integrateAfterForceComputation()");
265  integrateAfterForceComputation();
266 
267  logger(DEBUG, "about to call actionsAfterTimeStep()");
269 
270  //erase interactions that have not been used during the last time step
271  logger(DEBUG, "about to call interactionHandler.eraseOldInteractions(getNumberOfTimeSteps())");
272  interactionHandler.eraseOldInteractions(getNumberOfTimeSteps());
273  logger(DEBUG, "about to call interactionHandler.actionsAfterTimeStep()");
274  interactionHandler.actionsAfterTimeStep();
275  particleHandler.actionsAfterTimeStep();
276 
277  continueInTime();
278 
279  logger(DEBUG, "finished computeOneTimeStep()");
280  }
281 
283  {
284  // how many bulk elements in total
285  unsigned n_element = OomphProblem::solid_mesh_pt()->nelement();
286 
287  // loop over the bulk elements
288  for (unsigned e = 0; e < n_element; e++)
289  {
290  // get pointer to the bulk element
291  auto bulk_elem_pt = dynamic_cast<VELEMENT*>(OomphProblem::solid_mesh_pt()->element_pt(e));
292 
293  // get number of nodes in the bulk element
294  const unsigned nnode = bulk_elem_pt->nnode();
295 
296  // Set up memory for the shape/test functions
297  Shape psi(nnode);
298 
299  // three arrays that contain the x, y and z coordinates of the bulk element
300  Vector<double> listOfCoordX;
301  Vector<double> listOfCoordY;
302  Vector<double> listOfCoordZ;
303 
304  // get the x, y and z coordinates of the bulk element
305  for (unsigned n = 0; n < nnode; n++)
306  {
307  listOfCoordX.push_back(bulk_elem_pt->node_pt(n)->x(0) * Global_Physical_Variables::lenScale);
308  listOfCoordY.push_back(bulk_elem_pt->node_pt(n)->x(1) * Global_Physical_Variables::lenScale);
309  listOfCoordZ.push_back(bulk_elem_pt->node_pt(n)->x(2) * Global_Physical_Variables::lenScale);
310  }
311 
312  // get the bounding box of the bulk element
313  Vec3D min;
314  min.X = *min_element(listOfCoordX.begin(), listOfCoordX.end());
315  min.Y = *min_element(listOfCoordY.begin(), listOfCoordY.end());
316  min.Z = *min_element(listOfCoordZ.begin(), listOfCoordZ.end());
317  Vec3D max;
318  max.X = *max_element(listOfCoordX.begin(), listOfCoordX.end());
319  max.Y = *max_element(listOfCoordY.begin(), listOfCoordY.end());
320  max.Z = *max_element(listOfCoordZ.begin(), listOfCoordZ.end());
321 
322  // extend the bounding box if construct mapping with coarse graining
323  min -= Vec3D(1.0, 1.0, 1.0) * ( getCGWidth() - 2 * particleHandler.getLargestParticle()->getRadius() );
324  max += Vec3D(1.0, 1.0, 1.0) * ( getCGWidth() - 2 * particleHandler.getLargestParticle()->getRadius() );
325 
326  // get particles if there are in the bounding box
327  Vector<BaseParticle*> pList;
328  getParticlesInCell(min, max, pList);
329  if (pList.size() != 0)
330  {
331  DPMVCoupledElement DPM_coupled_elem;
332  // save the extended particle list if construct mapping with coarse graining functions
333  if (useCGMapping())
334  { DPM_coupled_elem.listOfCoupledParticlesExt = pList; }
335  // check if a particle resides within the bulk element
336  for (auto it = pList.begin(); it != pList.end(); it++)
337  {
338  Vector<double> s(3, 0.0), x(3, 0.0);
339  GeomObject* geom_obj_pt = 0;
340  // get global coordinates of particle center (DPM)
341  Vec3D xp = ( *it )->getPosition();
343  x[0] = xp.getX();
344  x[1] = xp.getY();
345  x[2] = xp.getZ();
346  // get local coordinates of particle center (FEM)
347  bulk_elem_pt->locate_zeta(x, geom_obj_pt, s);
348  // If true, the iterator is decremented after it is passed to erase() but before erase() is executed
349  if (geom_obj_pt == nullptr)
350  { pList.erase(it--); }
351  else
352  {
353  // save local coordinates of particle center
354  DPM_coupled_elem.listOfParticleCentersLoc.push_back(s);
355  // evaluate shape functions at particle center
356  bulk_elem_pt->shape(s, psi);
357  Vector<double> shapes;
358  for (unsigned l = 0; l < nnode; l++) shapes.push_back(psi(l));
359  // save shape functions at particle center
360  DPM_coupled_elem.listOfShapesAtParticleCenters.push_back(shapes);
361  }
362  }
363  // if there are particles in the bulk element (check with locate_zeta)
364  if (pList.size() != 0)
365  {
366  // save pointer to the bulk element and the pointer list to particles
367  DPM_coupled_elem.bulk_elem_pt = bulk_elem_pt;
368  DPM_coupled_elem.listOfCoupledParticles = pList;
369  listOfDPMVCoupledElements_.push_back(DPM_coupled_elem);
370  }
371  }
372  }
373  }
374 
379  {
380  // get a full list of all nodes that are coupled
381  Vector<CoupledSolidNode*> listOfSolidNodes;
383  {
384  // get pointers to all nodes of the bulk element
385  const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
386  for (unsigned n = 0; n < nnode; n++)
387  {
388  CoupledSolidNode* node_pt = dynamic_cast<CoupledSolidNode*>(coupled_elem.bulk_elem_pt->node_pt(n));
389  // if it is not a boundary node, add it to the list
390  if (!node_pt->is_on_boundary())
391  {
392  listOfSolidNodes.push_back(node_pt);
393  }
394  // set nodal weights to one if it is a boundary node
395  else
396  { node_pt->set_coupling_weight(0.5); }
397  }
398  }
399  // count number of occurrences of a node and set nodal weights
400  for (CoupledSolidNode* node_pt : listOfSolidNodes)
401  {
402  unsigned n = count(listOfSolidNodes.begin(), listOfSolidNodes.end(), node_pt);
403  if (n == pow(2, getSystemDimensions()))
404  { node_pt->set_coupling_weight(0.5); }
405  else
406  { node_pt->set_coupling_weight(0.5); }
407  // cout << "The solid node " << node_pt << " is shared by " << n << " coupled elements. ";
408  // cout << "Set weight to " << node_pt->get_coupling_weight() << endl;
409  }
410  }
411 
416  {
417  unsigned direction = 1;
418  // store the position[direction] of all coupled nodes in a vector
419  Vector<double> listOfPos;
421  {
422  // get pointers to all nodes of the bulk element
423  const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
424  for (unsigned n = 0; n < nnode; n++)
425  {
426  CoupledSolidNode* node_pt = dynamic_cast<CoupledSolidNode*>(coupled_elem.bulk_elem_pt->node_pt(n));
427  listOfPos.push_back(node_pt->x(direction));
428  }
429  }
430  // get the minimum and maximum coordinates[direction]
431  auto posMin = *min_element(listOfPos.begin(), listOfPos.end());
432  auto posMax = *max_element(listOfPos.begin(), listOfPos.end());
438  double weightMax = 1.0 - weightMin_;
440  {
441  // get pointers to all nodes of the bulk element
442  const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
443  for (unsigned n = 0; n < nnode; n++)
444  {
445  CoupledSolidNode* node_pt = dynamic_cast<CoupledSolidNode*>(coupled_elem.bulk_elem_pt->node_pt(n));
446  // assign nodal coupling weight using a linear function of distance
447  double weight = weightMin_
448  + ( weightMax - weightMin_ ) / ( posMax - posMin ) * ( node_pt->x(direction) - posMin );
449  // // assign nodal coupling weight using a cosine function of distance
450  // double amp = 0.5*(weightMax - weightMin_);
451  // double omega = M_PI/(posMax-posMin);
452  // double weight = weightMin_ + amp - amp*cos(omega*(node_pt->x(direction)-posMin));
453  node_pt->set_coupling_weight(weight);
454  // cout << "Set weight to " << node_pt->get_coupling_weight();
455  // cout << " at a distance of " << node_pt->x(direction)-posMin <<endl;
456  }
457  }
458  }
459 
461  Vector<Vector<double>>& couplingMatrix, Vector<Vector<double>>& nodalDisplDPM)
462  {
463  // get coupled particles in the bulk element
464  Vector<BaseParticle*> particles;
465  if (!useCGMapping())
466  { particles = coupled_elem.listOfCoupledParticles; }
467  else
468  { particles = coupled_elem.listOfCoupledParticlesExt; }
469  // get number of coupled particles in the bulk element
470  const unsigned nParticles = particles.size();
471  // get nodal coordinates of the bulk element
472  const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
473  // get dimension of the problem
474  const unsigned dim = coupled_elem.bulk_elem_pt->dim();
475  // prepare vector of particle displacement
476  Vector<Vector<double>> particleDispl(nParticles, Vector<double>(dim, 0));
477  // fill in the projection matrix of the coupled bulk element
478  getProjectionMatrix(coupled_elem, couplingMatrix);
479  // get particle displacements
480  for (unsigned m = 0; m < nParticles; m++)
481  {
482  Vec3D displ = particles[m]->getDisplacement() / Global_Physical_Variables::lenScale;
483  particleDispl[m][0] = displ.getX();
484  particleDispl[m][1] = displ.getY();
485  particleDispl[m][2] = displ.getZ();
486  }
487  // get projected nodal displacement field from particle displacement at the center
488  for (unsigned i = 0; i < nnode; i++)
489  {
490  for (unsigned k = 0; k < dim; k++)
491  {
492  for (unsigned j = 0; j < nParticles; j++)
493  {
494  nodalDisplDPM[i][k] += couplingMatrix[i][j] * particleDispl[j][k];
495  }
496  }
497  }
498  }
499 
500  void getProjectionMatrix(const DPMVCoupledElement& coupled_elem, Vector<Vector<double>>& couplingMatrix)
501  {
502  // get number of coupled particles in the bulk element
503  const unsigned nParticles = couplingMatrix[0].size();
504  // get nodal coordinates of the bulk element
505  const unsigned nnode = couplingMatrix.size();
506  const unsigned dim = coupled_elem.bulk_elem_pt->dim();
507  // if construct mapping with FEM basis functions
508  if (!useCGMapping())
509  {
510  // loop over shape functions at the nodes
511  for (unsigned l = 0; l < nnode; l++)
512  {
513  double sum = 0;
514  // loop over shape functions at the particles
515  for (unsigned m = 0; m < nParticles; m++)
516  {
517  // get the l-th shape function evaluated at the center of particle m
518  double shape = coupled_elem.listOfShapesAtParticleCenters[m][l];
519  // shape function weighted by particle volume N_{l,m}*V_m
520  shape *= coupled_elem.listOfCoupledParticles[m]->getVolume();
521  couplingMatrix[l][m] = shape;
522  sum += shape;
523  }
524  // normalize each row of the projection matrix to sum to one
525  for (unsigned m = 0; m < nParticles; m++)
526  {
527  // projection rule defined as N_{l,m}*V_m / sum_m N_{l,m}*V_m
528  if (sum != 0)
529  { couplingMatrix[l][m] /= sum; }
530  }
531  }
532  }
533  else
534  {
535  // set up memory for the shape functions
536  Shape psi(nnode);
537  // set up memory for the local and global coordinates
538  Vector<double> s(dim), x(dim);
539  // get the element pointer
540  auto el_pt = coupled_elem.bulk_elem_pt;
541  // get the number of integration points
542  const unsigned n_intpt = el_pt->integral_pt()->nweight();
543  // loop over the integration points
544  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
545  {
546  // get the integral weight
547  double w = el_pt->integral_pt()->weight(ipt);
548  // find the shape functions at the integration points r_i
549  el_pt->shape_at_knot(ipt, psi);
550  // get the value of the local coordinates at the integration point
551  for (unsigned i = 0; i < dim; i++)
552  { s[i] = el_pt->integral_pt()->knot(ipt, i); }
553  // get the value of the global coordinates at the integration point
554  el_pt->interpolated_x(s, x);
555  // loop over the nodes
556  for (unsigned l = 0; l < nnode; l++)
557  {
558  // loop over shape functions at the particles
559  for (unsigned m = 0; m < nParticles; m++)
560  {
561  // set CG coordinates
562  CGCoordinates::XYZ coord;
563  coord.setXYZ(Vec3D(x[0], x[1], x[2]));
564  // evaluate the value of CG function around particle m at CGcoords \phi(\vec r_i-r_m)
565  double phi = getCGFunction().evaluateCGFunction(
566  coupled_elem.listOfCoupledParticlesExt[m]->getPosition() /
568  // CG mapping defined as \tilde{N_{l,m}}_ipt = w_ipt * \phi(\vec r_i-r_m) * N_l(r_i)
569  double shape = w * phi * psi(l);
570  // CG mapping weighted by particle volume \tilde{N_{l,m}}*V_m
571  shape *= coupled_elem.listOfCoupledParticlesExt[m]->getVolume();
572  // sum over the integration points
573  couplingMatrix[l][m] += shape;
574  }
575  }
576  }
577  // normalize each row of the projection matrix to sum to one
578  for (unsigned l = 0; l < nnode; l++)
579  {
580  double sum = 0;
581  for (unsigned m = 0; m < nParticles; m++)
582  { sum += couplingMatrix[l][m]; }
583  for (unsigned m = 0; m < nParticles; m++)
584  {
585  // projection rule defined as \tilde{N_{l,m}}*V_m / sum_m \tilde{N_{l,m}}*V_m
586  if (sum != 0)
587  { couplingMatrix[l][m] /= sum; }
588  }
589  }
590  }
591  }
592 
594  {
595  // first loop over the coupled bulk elements
597  {
598  // get pointer to the bulk element
599  VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
600  // get particles and the number of coupled particles in the bulk element
601  Vector<BaseParticle*> particles;
602  if (!useCGMapping())
603  { particles = coupled_elem.listOfCoupledParticles; }
604  else
605  { particles = coupled_elem.listOfCoupledParticlesExt; }
606  const unsigned nParticles = particles.size();
607  // get number of nodes in the bulk element
608  const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
609  // prepare projection matrix
610  Vector<Vector<double>> couplingMatrix(nnode, Vector<double>(nParticles, 0.0));
611  // prepare vector of nodal displacements projected from particle centers
612  const unsigned dim = elem_pt->dim();
613  Vector<Vector<double>> nodalDisplDPM(nnode, Vector<double>(dim, 0.0));
614  // compute projected nodal displacements
615  getProjectionAndProjected(coupled_elem, couplingMatrix, nodalDisplDPM);
616  // compute nodal coupling forces by penalizing the difference in displacement
617  Vector<Vector<double>> nodalCouplingForces(nnode, Vector<double>(dim, 0.0));
618  computeNodalCouplingForces(elem_pt, nnode, dim, nodalDisplDPM, nodalCouplingForces);
619  // get coupling forces projected from coupled nodes to particles (note couplingMatrix^{-1} = couplingMatrix^T)
620  Vector<Vector<double>> particleCouplingForces(nParticles, Vector<double>(dim, 0.0));
621  for (unsigned i = 0; i < nParticles; i++)
622  {
623  for (unsigned k = 0; k < dim; k++)
624  {
625  for (unsigned j = 0; j < nnode; j++)
626  {
627  // note that particle and nodal coupling forces are in opposite directions
628  // the negative sign on the FEM side is taken care of in by element
629  particleCouplingForces[i][k] += couplingMatrix[j][i] * nodalCouplingForces[j][k];
630  }
631  }
632  // add particle coupling forces to be used by actionsAfterTimeStep(...)
633  auto p = particles[i];
634  Vec3D cForce = Vec3D(particleCouplingForces[i][0],
635  particleCouplingForces[i][1],
636  particleCouplingForces[i][2]);
637  // if construct mapping with coarse graining, add elements' contribution to particle coupling forces
638  if (!useCGMapping())
639  {
640  p->setCouplingForce(cForce * Global_Physical_Variables::forceScale());
641  }
642  // otherwise, set particle coupling force since the operation is locally within each element
643  else
644  {
645  p->addCouplingForce(cForce * Global_Physical_Variables::forceScale());
646  }
647  }
648  // assign nodal coupling force to the element to be used by element::fill_in_contribution_to_residuals(...)
649  elem_pt->set_nodal_coupling_residual(nodalCouplingForces);
650  }
651  }
652 
654  {
655  // first loop over the coupled bulk elements
657  {
658  // get pointer to the bulk element
659  VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
660 
661  // get particles and the number of coupled particles in the bulk element
662  Vector<BaseParticle*> particles;
663  if (!useCGMapping())
664  { particles = coupled_elem.listOfCoupledParticles; }
665  else
666  { particles = coupled_elem.listOfCoupledParticlesExt; }
667  const unsigned nParticles = particles.size();
668 
669  // get number of nodes in the bulk element
670  const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
671  const unsigned dim = elem_pt->dim();
672 
673  // prepare projection matrix \Pi_{i,\alpha}
674  Vector<Vector<double>> couplingMatrix(nnode, Vector<double>(nParticles, 0.0));
675  // fill in the projection matrix of the coupled bulk element
676  getProjectionMatrix(coupled_elem, couplingMatrix);
677 
678  // prepare the coupled stiffness matrix of the bulk element, penalty*V_{ij}
679  Vector<Vector<double>> bulkStiffness(nnode, Vector<double>(nnode, 0.0));
680  // fill in the volume coupled stiffness matrix
681  getBulkStiffness(elem_pt, bulkStiffness);
682  // assign penalty*\Pi_{\alpha,i}*V_{i,j} to the VCoupled element
683  Vector<Vector<double>> couplingStiffness(nParticles, Vector<double>(nnode, 0.0));
684  for (unsigned m = 0; m < nParticles; m++)
685  {
686  for (unsigned j = 0; j < nnode; j++)
687  {
688  for (unsigned i = 0; i < nnode; i++)
689  {
690  couplingStiffness[m][j] += couplingMatrix[i][m] * bulkStiffness[i][j];
691  }
692  }
693  }
694  elem_pt->setCouplingStiffness(couplingStiffness);
695 
696  // prepare the coupling residual projected from the particles to the node
697  Vector<Vector<double>> nodalResidualDPM(nnode, Vector<double>(dim, 0.0));
698 
699  // get displacement u(t+dt) = v(t+0.5dt)*dt after integrateBeforeForceComputation(...)
700  for (unsigned m = 0; m < nParticles; m++)
701  {
702  Vec3D cForce;
703  for (unsigned l = 0; l < nnode; l++)
704  {
705  for (unsigned ll = 0; ll < nnode; ll++)
706  {
707  for (unsigned mm = 0; mm < nParticles; mm++)
708  {
709  // change finite element coupling stiffness into the DEM scale
710  cForce += couplingMatrix[l][m] * bulkStiffness[l][ll] *
711  Global_Physical_Variables::stiffScale() * couplingMatrix[ll][mm] *
712  particles[mm]->getDisplacement();
713  }
714  }
715  }
716  // add the first particle coupling force: -penalty*\PI_{\alpha i}*V_{ij}*\PI_{j \beta}*u_\beta
717  particles[m]->addCouplingForce(-cForce);
718  }
719  // compute the residual projected from the particles to the nodes
720  for (unsigned i = 0; i < nnode; i++)
721  {
722  for (unsigned j = 0; j < nnode; j++)
723  {
724  for (unsigned m = 0; m < nParticles; m++)
725  {
726  Vec3D displ = particles[m]->getDisplacement() / Global_Physical_Variables::lenScale;
727  double particleStiffness = bulkStiffness[i][j] * couplingMatrix[j][m];
728  nodalResidualDPM[i][0] -= particleStiffness * displ.getX();
729  nodalResidualDPM[i][1] -= particleStiffness * displ.getY();
730  nodalResidualDPM[i][2] -= particleStiffness * displ.getZ();
731  }
732  }
733  }
734  // assign residual to the coupled bulk element to be used by element::fill_in_contribution_to_residuals(...)
735  elem_pt->set_nodal_coupling_residual(nodalResidualDPM);
736 
737  // assign Jacobian to the coupled bulk element to be used by element::fill_in_contribution_to_residuals(...)
738  elem_pt->set_nodal_coupling_jacobian(bulkStiffness);
739  }
740  }
741 
743  {
744  // first loop over the coupled bulk elements
746  {
747  // get pointer to the bulk element
748  VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
749 
750  // get particles and the number of coupled particles in the bulk element
751  Vector<BaseParticle*> particles;
752  if (!useCGMapping())
753  { particles = coupled_elem.listOfCoupledParticles; }
754  else
755  { particles = coupled_elem.listOfCoupledParticlesExt; }
756  const unsigned nParticles = particles.size();
757 
758  // get number of nodes in the bulk element
759  const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
760  const unsigned dim = elem_pt->dim();
761 
762  // loop over the particles
763  for (unsigned m = 0; m < nParticles; m++)
764  {
765  // get "residual" from the bulk element
766  Vec3D cForce;
767  for (unsigned l = 0; l < nnode; l++)
768  {
769  double coupling = elem_pt->getCouplingStiffness(m, l);
770  // get the nodal displacement
771  Vec3D force = coupling * Vec3D(elem_pt->nodal_position(0, l, 0) - elem_pt->nodal_position(1, l, 0),
772  elem_pt->nodal_position(0, l, 1) - elem_pt->nodal_position(1, l, 1),
773  elem_pt->nodal_position(0, l, 2) - elem_pt->nodal_position(1, l, 2));
774  // add contributions to the coupling force
775  cForce += force;
776  }
777  // add the second particle coupling force: penalty*\PI_{\alpha i}*V_{ij}*u_j
778  particles[m]->addCouplingForce(cForce * Global_Physical_Variables::forceScale());
779  }
780  }
781  }
782 
783  void computeNodalCouplingForces(VELEMENT*& elem_pt, const unsigned& nnode, const unsigned& dim,
784  const Vector<Vector<double>>& nodalDisplDPM,
785  Vector<Vector<double>>& nodalCouplingForces)
786  {
787  // get the difference between nodal displacement from FEM and projected ones from DPM
788  Vector<Vector<double>> nodalDisplDiff(nnode, Vector<double>(dim, 0));
789  for (unsigned l = 0; l < nnode; l++)
790  {
791  for (unsigned i = 0; i < dim; i++)
792  {
793  double displFEM = elem_pt->nodal_position(0, l, i) - elem_pt->nodal_position(1, l, i);
794  double displDPM = nodalDisplDPM[l][i];
795  nodalDisplDiff[l][i] = displDPM - displFEM;
796  }
797  }
798  // Set up memory for the shape functions
799  Shape psi(nnode);
800  DShape dpsidxi(nnode, dim);
801  const unsigned n_intpt = elem_pt->integral_pt()->nweight();
802  // Loop over the integration points
803  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
804  {
805  // Get the integral weight
806  double w = elem_pt->integral_pt()->weight(ipt);
807  // Call the derivatives of the shape functions (and get Jacobian)
808  double J = elem_pt->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
809  // Loop over the test functions, nodes of the element
810  for (unsigned l = 0; l < nnode; l++)
811  {
812  // Loop over the nodes of the element again
813  for (unsigned ll = 0; ll < nnode; ll++)
814  {
815  double vol = Global_Physical_Variables::penalty * psi(ll) * psi(l) * w * J;
816  /*
817  for (unsigned i=0; i<dim; i++)
818  {
819  vol += Global_Physical_Variables::penalty*2*Global_Physical_Variables::L/3.0 * dpsidxi(ll,i)*dpsidxi(l,i)*w*J;
820  }
821  */
822  for (unsigned i = 0; i < dim; i++)
823  {
824  nodalCouplingForces[l][i] += -vol * nodalDisplDiff[l][i];
825  }
826  }
827  }
828  }
829  }
830 
831 
832  void getBulkStiffness(VELEMENT*& elem_pt, Vector<Vector<double>>& bulkStiffness)
833  {
834  // get number of nodes in the bulk element
835  const unsigned nnode = elem_pt->nnode();
836  // get dimension of the problem
837  const unsigned dim = elem_pt->dim();
838  // Set up memory for the shape functions
839  Shape psi(nnode);
840  DShape dpsidxi(nnode, dim);
841  const unsigned n_intpt = elem_pt->integral_pt()->nweight();
842  // Loop over the integration points
843  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
844  {
845  // Get the integral weight
846  double w = elem_pt->integral_pt()->weight(ipt);
847  // Call the derivatives of the shape functions (and get Jacobian)
848  double J = elem_pt->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
849  // Loop over the test functions, nodes of the element
850  for (unsigned l = 0; l < nnode; l++)
851  {
852  // Loop over the nodes of the element again
853  for (unsigned ll = 0; ll < nnode; ll++)
854  {
855  bulkStiffness[l][ll] += Global_Physical_Variables::penalty * psi(l) * psi(ll) * w * J;
856  }
857  }
858  }
859  }
860 
861  void getInvCoupledMassPerTime(const DPMVCoupledElement& elem, const Vector<Vector<double>>& couplingMatrix,
862  const Vector<Vector<double>>& bulkStiffness, Vector<double>& invCMassPerTime)
863  {
864  // get number of nodes in the bulk element
865  const unsigned nnode = bulkStiffness.size();
866  // get number of coupled particles in the bulk element
867  const unsigned nParticles = invCMassPerTime.size();
868  // get coupled particles in the bulk element
869  Vector<BaseParticle*> particles;
870  if (!useCGMapping())
871  { particles = elem.listOfCoupledParticles; }
872  else
873  { particles = elem.listOfCoupledParticlesExt; }
874  // loop over the particles
875  for (unsigned m = 0; m < nParticles; m++)
876  {
877  // get particle mass
878  double particleMass = particles[m]->getMass() / particles[m]->getInvCouplingWeight();
879  // initialize additional coupling mass from the bulk element
880  double cMassPerTimeSquared = 0.0;
881  // loop over the nodes
882  for (unsigned l = 0; l < nnode; l++)
883  {
884  // loop over the nodes again
885  for (unsigned ll = 0; ll < nnode; ll++)
886  {
887  // add contributions to the coupling stiffness or mass over time squared
888  cMassPerTimeSquared += couplingMatrix[l][m] * bulkStiffness[l][ll] * couplingMatrix[ll][m];
889  }
890  }
891  cMassPerTimeSquared *= 0.5;
892  // add coupling mass to the particle (in the DPM scale)
893  double cMass = cMassPerTimeSquared * Global_Physical_Variables::stiffScale() * pow(getTimeStep(), 2);
894  particles[m]->setCouplingMass(cMass);
895  // set coupling mass/time only along the diagonal (lumped version)
896  double scale = Global_Physical_Variables::massScale() / Global_Physical_Variables::timeScale;
897  double oomph_dt = time_pt()->dt(0);
898  invCMassPerTime[m] = 1.0 / ( particleMass / getTimeStep() / scale + cMassPerTimeSquared * oomph_dt );
899  }
900  }
901 
902  /*
903  * used in computeOneTimeStepForVCoupling to compute weight on particles
904  */
906  {
907  // first loop over the coupled bulk elements
909  {
910  // get pointer to the bulk element
911  VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
912  // Find out how many nodes there are
913  const unsigned n_node = elem_pt->nnode();
914  // Set up memory for the shape/test functions
915  Shape psi(n_node);
916  // loop over the particles in the bulk element
917  for (unsigned i = 0; i < coupled_elem.listOfCoupledParticles.size(); i++)
918  {
919  auto shapes = coupled_elem.listOfShapesAtParticleCenters[i];
920  // Compute weight at particle center
921  double weightAtCenter = 0;
922  for (unsigned l = 0; l < n_node; l++)
923  {
924  CoupledSolidNode* node_pt = dynamic_cast<CoupledSolidNode*>(elem_pt->node_pt(l));
925  weightAtCenter += shapes[l] * node_pt->get_coupling_weight();
926  }
927  double invWeightAtCenter = 1.0 / weightAtCenter;
928 
929  // get point to particle i in the coupled bulk element and set the coupling weight
930  auto p = coupled_elem.listOfCoupledParticles[i];
931  p->setInvCouplingWeight(invWeightAtCenter);
932 
933  // skip particle i if no interactions
934  if (p->getInteractions().size() == 0) continue;
935 
936  if (explicit_)
937  {
938  // loop over interactions with particle i (FIXME accessing coupled interactions is n squared...)
939  for (auto inter : p->getInteractions())
940  {
941  // skip if the coupling weight is already computed
942  if (inter->isCoupled())
943  { continue; }
944  else
945  {
946  // check if a contact point resides within the bulk element
947  Vector<double> s(3, 0.0), x(3, 0.0);
948  GeomObject* geom_obj_pt = 0;
949  // get global coordinates of contact point (DPM)
950  Vec3D xc = inter->getContactPoint();
952  x[0] = xc.getX();
953  x[1] = xc.getY();
954  x[2] = xc.getZ();
955  // get local coordinates of contact point (FEM)
956  elem_pt->locate_zeta(x, geom_obj_pt, s);
957  // if true compute weight at contact point
958  if (geom_obj_pt != nullptr)
959  {
960  // Get shape/test functions
961  elem_pt->shape(s, psi);
962  double weightAtContactPoint = 0;
963  for (unsigned l = 0; l < n_node; l++)
964  {
965  CoupledSolidNode* node_pt = dynamic_cast<CoupledSolidNode*>(elem_pt->node_pt(l));
966  weightAtContactPoint += node_pt->get_coupling_weight() * psi(l);
967  }
968  inter->setCouplingWeight(weightAtContactPoint);
969  inter->setCoupled();
970  }
971  }
972  }
973  }
974  }
975  }
976  }
977 
978  void setExplicit(const bool& flag)
979  { explicit_ = flag; }
980 
981 private:
982  // vector of bulk elements coupled with discrete particles
983  Vector<DPMVCoupledElement> listOfDPMVCoupledElements_;
984  // minimal weighting factor
985  double weightMin_ = 0.0;
986  // flag to use implicit or explicit coupling scheme (VolumeCoupling)
987  bool explicit_ = false;
988 };
989 
990 #endif //VOLUME_COUPLING_H
const unsigned n
Definition: CG3DPackingUnitTest.cpp:32
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ INFO
@ DEBUG
Definition: BaseBoundary.h:49
Definition: BaseCoupling.h:48
bool useCGMapping()
Definition: BaseCoupling.h:181
void solveOomph()
Definition: BaseCoupling.h:141
double getCGWidth()
Definition: BaseCoupling.h:176
CGFunctions::LucyXYZ getCGFunction()
Definition: BaseCoupling.h:184
void setCGWidth(const double &width)
Definition: BaseCoupling.h:158
void solveMercury(unsigned long nt)
Definition: BaseCoupling.h:150
void addForce(const Vec3D &addForce)
Adds an amount to the force on this BaseInteractable.
Definition: BaseInteractable.cc:116
Definition: BaseParticle.h:54
bool isFixed() const override
Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Ma...
Definition: BaseParticle.h:93
Mdouble getMass() const
Returns the particle's mass.
Definition: BaseParticle.h:322
Defines the position of the CGPoint, in the non-averaged directions, i.e. all directions on which spa...
Definition: XYZ.h:75
void setXYZ(Vec3D p)
Returns the position of the current CGPoint, in the non-averaged directions.
Definition: XYZ.cc:61
Mdouble evaluateCGFunction(const Vec3D &position, const Coordinates &r)
Evaluates the coarse-graining function.
Definition: Vector.h:51
Mdouble Y
Definition: Vector.h:66
Mdouble Z
Definition: Vector.h:66
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:332
Mdouble X
the vector components
Definition: Vector.h:66
Mdouble getZ() const
Definition: Vector.h:408
Mdouble getX() const
Definition: Vector.h:402
Mdouble getY() const
Definition: Vector.h:405
Definition: VCoupling.h:32
void computeOneTimeStepForVCoupling(const unsigned &nstep)
Definition: VCoupling.h:179
void setNodalWeightByDistance()
Definition: VCoupling.h:415
void computeCouplingForce()
Definition: VCoupling.h:593
VolumeCoupling(const double &wMin, const bool &flag)
Definition: VCoupling.h:56
void actionsAfterTimeStep() override
Definition: VCoupling.h:82
void computeCouplingOnDPM()
Definition: VCoupling.h:742
void solveFirstHalfTimeStep()
Definition: VCoupling.h:206
void computeNodalCouplingForces(VELEMENT *&elem_pt, const unsigned &nnode, const unsigned &dim, const Vector< Vector< double >> &nodalDisplDPM, Vector< Vector< double >> &nodalCouplingForces)
Definition: VCoupling.h:783
void checkParticlesInFiniteElems()
Definition: VCoupling.h:282
void getProjectionMatrix(const DPMVCoupledElement &coupled_elem, Vector< Vector< double >> &couplingMatrix)
Definition: VCoupling.h:500
bool explicit_
Definition: VCoupling.h:987
void solveVolumeCoupling()
Definition: VCoupling.h:138
void computeCouplingOnFEM()
Definition: VCoupling.h:653
VolumeCoupledElement< OomphProblem::ELEMENT > VELEMENT
Definition: VCoupling.h:34
VolumeCoupling()=default
void setExplicit(const bool &flag)
Definition: VCoupling.h:978
void solveVolumeCoupling(const double &width)
Definition: VCoupling.h:129
void computeExternalForces(BaseParticle *CI) override
Definition: VCoupling.h:65
void getBulkStiffness(VELEMENT *&elem_pt, Vector< Vector< double >> &bulkStiffness)
Definition: VCoupling.h:832
double weightMin_
Definition: VCoupling.h:985
void getInvCoupledMassPerTime(const DPMVCoupledElement &elem, const Vector< Vector< double >> &couplingMatrix, const Vector< Vector< double >> &bulkStiffness, Vector< double > &invCMassPerTime)
Definition: VCoupling.h:861
void solveSecondHalfTimeStep()
Definition: VCoupling.h:254
void getProjectionAndProjected(const DPMVCoupledElement &coupled_elem, Vector< Vector< double >> &couplingMatrix, Vector< Vector< double >> &nodalDisplDPM)
Definition: VCoupling.h:460
Vector< DPMVCoupledElement > listOfDPMVCoupledElements_
Definition: VCoupling.h:983
void setNodalWeightOnCubicMesh()
Definition: VCoupling.h:378
void computeWeightOnParticles()
Definition: VCoupling.h:905
double lenScale
Length scale.
Definition: VCoupledElement.h:31
double timeScale
Time scale.
Definition: VCoupledElement.h:34
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
Definition: VCoupling.h:40
Vector< BaseParticle * > listOfCoupledParticlesExt
Definition: VCoupling.h:45
Vector< Vector< double > > listOfParticleCentersLoc
Definition: VCoupling.h:47
VELEMENT * bulk_elem_pt
Definition: VCoupling.h:42
Vector< BaseParticle * > listOfCoupledParticles
Definition: VCoupling.h:44
Vector< Vector< double > > listOfShapesAtParticleCenters
Definition: VCoupling.h:49