26 #ifndef VOLUME_COUPLING_H
27 #define VOLUME_COUPLING_H
34 typedef VolumeCoupledElement<OomphProblem::ELEMENT>
VELEMENT;
74 CI->
addForce(CI->getCouplingForce());
90 p->resetCouplingForce();
91 p->resetInvCouplingWeight();
92 p->resetCouplingMass();
97 const auto p = particleHandler.getLastObject();
98 Vec3D f = p->getMass() * getGravity();
99 Vec3D a = p->getForce();
102 logger(
INFO,
"time %, body force: %, total force: %, position %",
103 getTime(), f, a, p->getPosition());
107 for (
auto p : particleHandler)
110 for (
auto i : p->getInteractions())
112 if (
i->getI()->getName() ==
"TriangleWall")
n++;
116 logger(
INFO,
"particle % is interacting with % TriangleWalls", p->getId(),
n);
117 for (
auto i : p->getInteractions())
119 if (
i->getI()->getName() ==
"TriangleWall")
121 logger(
INFO,
"particle % is interacting with wall %", p->getId(),
i->getI()->getId());
132 setCGWidth(width * particleHandler.getLargestParticle()->getRadius());
143 double mercury_dt = getTimeStep();
155 double oomph_dt = nstep * mercury_dt;
158 OomphProblem::set_initial_conditions(oomph_dt);
163 while (getTime() < getTimeMax())
168 if (getVtkWriter()->getFileCounter() > nDone)
170 OomphProblem::doc_solution();
171 nDone = getVtkWriter()->getFileCounter();
213 logger(
DEBUG,
"about to call hGridActionsBeforeIntegration()");
214 hGridActionsBeforeIntegration();
217 logger(
DEBUG,
"about to call integrateBeforeForceComputation()");
218 integrateBeforeForceComputation();
220 logger(
DEBUG,
"about to call performGhostParticleUpdate()");
221 performGhostParticleUpdate();
227 logger(
DEBUG,
"about to call checkInteractionWithBoundaries()");
228 checkInteractionWithBoundaries();
230 logger(
DEBUG,
"about to call hGridActionsAfterIntegration()");
231 hGridActionsAfterIntegration();
238 b->checkBoundaryBeforeTimeStep(
this);
241 logger(
DEBUG,
"about to call actionsBeforeTimeStep()");
242 actionsBeforeTimeStep();
244 logger(
DEBUG,
"about to call checkAndDuplicatePeriodicParticles()");
245 checkAndDuplicatePeriodicParticles();
247 logger(
DEBUG,
"about to call hGridActionsBeforeTimeStep()");
248 hGridActionsBeforeTimeStep();
260 logger(
DEBUG,
"about to call removeDuplicatePeriodicParticles()");
261 removeDuplicatePeriodicParticles();
264 logger(
DEBUG,
"about to call integrateAfterForceComputation()");
265 integrateAfterForceComputation();
267 logger(
DEBUG,
"about to call actionsAfterTimeStep()");
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();
285 unsigned n_element = OomphProblem::solid_mesh_pt()->nelement();
288 for (
unsigned e = 0; e < n_element; e++)
291 auto bulk_elem_pt =
dynamic_cast<VELEMENT*
>(OomphProblem::solid_mesh_pt()->element_pt(e));
294 const unsigned nnode = bulk_elem_pt->nnode();
300 Vector<double> listOfCoordX;
301 Vector<double> listOfCoordY;
302 Vector<double> listOfCoordZ;
305 for (
unsigned n = 0;
n < nnode;
n++)
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());
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());
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() );
327 Vector<BaseParticle*> pList;
328 getParticlesInCell(min, max, pList);
329 if (pList.size() != 0)
336 for (
auto it = pList.begin(); it != pList.end(); it++)
338 Vector<double> s(3, 0.0), x(3, 0.0);
339 GeomObject* geom_obj_pt = 0;
341 Vec3D xp = ( *it )->getPosition();
347 bulk_elem_pt->locate_zeta(x, geom_obj_pt, s);
349 if (geom_obj_pt ==
nullptr)
350 { pList.erase(it--); }
356 bulk_elem_pt->shape(s, psi);
357 Vector<double> shapes;
358 for (
unsigned l = 0; l < nnode; l++) shapes.push_back(psi(l));
364 if (pList.size() != 0)
381 Vector<CoupledSolidNode*> listOfSolidNodes;
385 const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
386 for (
unsigned n = 0;
n < nnode;
n++)
388 CoupledSolidNode* node_pt =
dynamic_cast<CoupledSolidNode*
>(coupled_elem.bulk_elem_pt->node_pt(
n));
390 if (!node_pt->is_on_boundary())
392 listOfSolidNodes.push_back(node_pt);
396 { node_pt->set_coupling_weight(0.5); }
400 for (CoupledSolidNode* node_pt : listOfSolidNodes)
402 unsigned n = count(listOfSolidNodes.begin(), listOfSolidNodes.end(), node_pt);
403 if (
n == pow(2, getSystemDimensions()))
404 { node_pt->set_coupling_weight(0.5); }
406 { node_pt->set_coupling_weight(0.5); }
417 unsigned direction = 1;
419 Vector<double> listOfPos;
423 const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
424 for (
unsigned n = 0;
n < nnode;
n++)
426 CoupledSolidNode* node_pt =
dynamic_cast<CoupledSolidNode*
>(coupled_elem.bulk_elem_pt->node_pt(
n));
427 listOfPos.push_back(node_pt->x(direction));
431 auto posMin = *min_element(listOfPos.begin(), listOfPos.end());
432 auto posMax = *max_element(listOfPos.begin(), listOfPos.end());
442 const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
443 for (
unsigned n = 0;
n < nnode;
n++)
445 CoupledSolidNode* node_pt =
dynamic_cast<CoupledSolidNode*
>(coupled_elem.bulk_elem_pt->node_pt(
n));
448 + ( weightMax -
weightMin_ ) / ( posMax - posMin ) * ( node_pt->x(direction) - posMin );
453 node_pt->set_coupling_weight(weight);
461 Vector<Vector<double>>& couplingMatrix, Vector<Vector<double>>& nodalDisplDPM)
464 Vector<BaseParticle*> particles;
470 const unsigned nParticles = particles.size();
472 const unsigned nnode = coupled_elem.
bulk_elem_pt->nnode();
476 Vector<Vector<double>> particleDispl(nParticles, Vector<double>(dim, 0));
480 for (
unsigned m = 0; m < nParticles; m++)
483 particleDispl[m][0] = displ.
getX();
484 particleDispl[m][1] = displ.
getY();
485 particleDispl[m][2] = displ.
getZ();
488 for (
unsigned i = 0;
i < nnode;
i++)
490 for (
unsigned k = 0; k < dim; k++)
492 for (
unsigned j = 0; j < nParticles; j++)
494 nodalDisplDPM[
i][k] += couplingMatrix[
i][j] * particleDispl[j][k];
503 const unsigned nParticles = couplingMatrix[0].size();
505 const unsigned nnode = couplingMatrix.size();
511 for (
unsigned l = 0; l < nnode; l++)
515 for (
unsigned m = 0; m < nParticles; m++)
521 couplingMatrix[l][m] = shape;
525 for (
unsigned m = 0; m < nParticles; m++)
529 { couplingMatrix[l][m] /= sum; }
538 Vector<double> s(dim), x(dim);
542 const unsigned n_intpt = el_pt->integral_pt()->nweight();
544 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
547 double w = el_pt->integral_pt()->weight(ipt);
549 el_pt->shape_at_knot(ipt, psi);
551 for (
unsigned i = 0;
i < dim;
i++)
552 { s[
i] = el_pt->integral_pt()->knot(ipt,
i); }
554 el_pt->interpolated_x(s, x);
556 for (
unsigned l = 0; l < nnode; l++)
559 for (
unsigned m = 0; m < nParticles; m++)
569 double shape = w * phi * psi(l);
573 couplingMatrix[l][m] += shape;
578 for (
unsigned l = 0; l < nnode; l++)
581 for (
unsigned m = 0; m < nParticles; m++)
582 { sum += couplingMatrix[l][m]; }
583 for (
unsigned m = 0; m < nParticles; m++)
587 { couplingMatrix[l][m] /= sum; }
599 VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
601 Vector<BaseParticle*> particles;
603 { particles = coupled_elem.listOfCoupledParticles; }
605 { particles = coupled_elem.listOfCoupledParticlesExt; }
606 const unsigned nParticles = particles.size();
608 const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
610 Vector<Vector<double>> couplingMatrix(nnode, Vector<double>(nParticles, 0.0));
612 const unsigned dim = elem_pt->dim();
613 Vector<Vector<double>> nodalDisplDPM(nnode, Vector<double>(dim, 0.0));
617 Vector<Vector<double>> nodalCouplingForces(nnode, Vector<double>(dim, 0.0));
620 Vector<Vector<double>> particleCouplingForces(nParticles, Vector<double>(dim, 0.0));
621 for (
unsigned i = 0;
i < nParticles;
i++)
623 for (
unsigned k = 0; k < dim; k++)
625 for (
unsigned j = 0; j < nnode; j++)
629 particleCouplingForces[
i][k] += couplingMatrix[j][
i] * nodalCouplingForces[j][k];
633 auto p = particles[
i];
635 particleCouplingForces[
i][1],
636 particleCouplingForces[
i][2]);
640 p->setCouplingForce(cForce * Global_Physical_Variables::forceScale());
645 p->addCouplingForce(cForce * Global_Physical_Variables::forceScale());
649 elem_pt->set_nodal_coupling_residual(nodalCouplingForces);
659 VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
662 Vector<BaseParticle*> particles;
664 { particles = coupled_elem.listOfCoupledParticles; }
666 { particles = coupled_elem.listOfCoupledParticlesExt; }
667 const unsigned nParticles = particles.size();
670 const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
671 const unsigned dim = elem_pt->dim();
674 Vector<Vector<double>> couplingMatrix(nnode, Vector<double>(nParticles, 0.0));
679 Vector<Vector<double>> bulkStiffness(nnode, Vector<double>(nnode, 0.0));
683 Vector<Vector<double>> couplingStiffness(nParticles, Vector<double>(nnode, 0.0));
684 for (
unsigned m = 0; m < nParticles; m++)
686 for (
unsigned j = 0; j < nnode; j++)
688 for (
unsigned i = 0;
i < nnode;
i++)
690 couplingStiffness[m][j] += couplingMatrix[
i][m] * bulkStiffness[
i][j];
694 elem_pt->setCouplingStiffness(couplingStiffness);
697 Vector<Vector<double>> nodalResidualDPM(nnode, Vector<double>(dim, 0.0));
700 for (
unsigned m = 0; m < nParticles; m++)
703 for (
unsigned l = 0; l < nnode; l++)
705 for (
unsigned ll = 0; ll < nnode; ll++)
707 for (
unsigned mm = 0; mm < nParticles; mm++)
710 cForce += couplingMatrix[l][m] * bulkStiffness[l][ll] *
711 Global_Physical_Variables::stiffScale() * couplingMatrix[ll][mm] *
712 particles[mm]->getDisplacement();
717 particles[m]->addCouplingForce(-cForce);
720 for (
unsigned i = 0;
i < nnode;
i++)
722 for (
unsigned j = 0; j < nnode; j++)
724 for (
unsigned m = 0; m < nParticles; m++)
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();
735 elem_pt->set_nodal_coupling_residual(nodalResidualDPM);
738 elem_pt->set_nodal_coupling_jacobian(bulkStiffness);
748 VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
751 Vector<BaseParticle*> particles;
753 { particles = coupled_elem.listOfCoupledParticles; }
755 { particles = coupled_elem.listOfCoupledParticlesExt; }
756 const unsigned nParticles = particles.size();
759 const unsigned nnode = coupled_elem.bulk_elem_pt->nnode();
760 const unsigned dim = elem_pt->dim();
763 for (
unsigned m = 0; m < nParticles; m++)
767 for (
unsigned l = 0; l < nnode; l++)
769 double coupling = elem_pt->getCouplingStiffness(m, l);
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));
778 particles[m]->addCouplingForce(cForce * Global_Physical_Variables::forceScale());
784 const Vector<Vector<double>>& nodalDisplDPM,
785 Vector<Vector<double>>& nodalCouplingForces)
788 Vector<Vector<double>> nodalDisplDiff(nnode, Vector<double>(dim, 0));
789 for (
unsigned l = 0; l < nnode; l++)
791 for (
unsigned i = 0;
i < dim;
i++)
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;
800 DShape dpsidxi(nnode, dim);
801 const unsigned n_intpt = elem_pt->integral_pt()->nweight();
803 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
806 double w = elem_pt->integral_pt()->weight(ipt);
808 double J = elem_pt->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
810 for (
unsigned l = 0; l < nnode; l++)
813 for (
unsigned ll = 0; ll < nnode; ll++)
815 double vol = Global_Physical_Variables::penalty * psi(ll) * psi(l) * w * J;
822 for (
unsigned i = 0;
i < dim;
i++)
824 nodalCouplingForces[l][
i] += -vol * nodalDisplDiff[l][
i];
835 const unsigned nnode = elem_pt->nnode();
837 const unsigned dim = elem_pt->dim();
840 DShape dpsidxi(nnode, dim);
841 const unsigned n_intpt = elem_pt->integral_pt()->nweight();
843 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
846 double w = elem_pt->integral_pt()->weight(ipt);
848 double J = elem_pt->dshape_lagrangian_at_knot(ipt, psi, dpsidxi);
850 for (
unsigned l = 0; l < nnode; l++)
853 for (
unsigned ll = 0; ll < nnode; ll++)
855 bulkStiffness[l][ll] += Global_Physical_Variables::penalty * psi(l) * psi(ll) * w * J;
862 const Vector<Vector<double>>& bulkStiffness, Vector<double>& invCMassPerTime)
865 const unsigned nnode = bulkStiffness.size();
867 const unsigned nParticles = invCMassPerTime.size();
869 Vector<BaseParticle*> particles;
875 for (
unsigned m = 0; m < nParticles; m++)
878 double particleMass = particles[m]->getMass() / particles[m]->getInvCouplingWeight();
880 double cMassPerTimeSquared = 0.0;
882 for (
unsigned l = 0; l < nnode; l++)
885 for (
unsigned ll = 0; ll < nnode; ll++)
888 cMassPerTimeSquared += couplingMatrix[l][m] * bulkStiffness[l][ll] * couplingMatrix[ll][m];
891 cMassPerTimeSquared *= 0.5;
893 double cMass = cMassPerTimeSquared * Global_Physical_Variables::stiffScale() * pow(getTimeStep(), 2);
894 particles[m]->setCouplingMass(cMass);
897 double oomph_dt = time_pt()->dt(0);
898 invCMassPerTime[m] = 1.0 / ( particleMass / getTimeStep() / scale + cMassPerTimeSquared * oomph_dt );
911 VELEMENT* elem_pt = coupled_elem.bulk_elem_pt;
913 const unsigned n_node = elem_pt->nnode();
917 for (
unsigned i = 0;
i < coupled_elem.listOfCoupledParticles.size();
i++)
919 auto shapes = coupled_elem.listOfShapesAtParticleCenters[
i];
921 double weightAtCenter = 0;
922 for (
unsigned l = 0; l < n_node; l++)
924 CoupledSolidNode* node_pt =
dynamic_cast<CoupledSolidNode*
>(elem_pt->node_pt(l));
925 weightAtCenter += shapes[l] * node_pt->get_coupling_weight();
927 double invWeightAtCenter = 1.0 / weightAtCenter;
930 auto p = coupled_elem.listOfCoupledParticles[
i];
931 p->setInvCouplingWeight(invWeightAtCenter);
934 if (p->getInteractions().size() == 0)
continue;
939 for (
auto inter : p->getInteractions())
942 if (inter->isCoupled())
947 Vector<double> s(3, 0.0), x(3, 0.0);
948 GeomObject* geom_obj_pt = 0;
950 Vec3D xc = inter->getContactPoint();
956 elem_pt->locate_zeta(x, geom_obj_pt, s);
958 if (geom_obj_pt !=
nullptr)
961 elem_pt->shape(s, psi);
962 double weightAtContactPoint = 0;
963 for (
unsigned l = 0; l < n_node; l++)
965 CoupledSolidNode* node_pt =
dynamic_cast<CoupledSolidNode*
>(elem_pt->node_pt(l));
966 weightAtContactPoint += node_pt->get_coupling_weight() * psi(l);
968 inter->setCouplingWeight(weightAtContactPoint);
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.
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.
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
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