ScaleCoupling< M, O > Class Template Reference

#include <ScaleCoupling.h>

Inheritance diagram for ScaleCoupling< M, O >:

Classes
struct	CoupledElement
	Stores properties of a coupled element: pointer to the element and list of particles coupled to it. More...
struct	CoupledParticle
	For all particles, stores coupling properties: coupling force, pointer to coupled element and location in coupled element. More...

Public Member Functions
const std::vector< CoupledElement > &	getCoupledElements ()
	get coupled element More...
const std::vector< CoupledParticle > &	getCoupledParticles ()
	get coupled particle More...
void	setPenalty (double penalty)
	Sets penalty parameter. More...
void	setCouplingWeight (const std::function< double(double x, double y, double z)> &couplingWeight)
	Sets weight function (determines which nodes/el_pts are in coupling zone) More...
void	solveScaleCoupling ()
	Computes nStep, the ratio of FEM and DEM time steps, then calls solveSurfaceCoupling(nStep) More...
Public Member Functions inherited from BaseCoupling< M, O >
	BaseCoupling ()=default
void	setName (std::string name)
std::string	getName () const
void	removeOldFiles () const
void	writeEneTimeStep (std::ostream &os) const override
void	writeEneHeader (std::ostream &os) const override
void	solveOomph ()
void	solveMercury (unsigned long nt)
void	setCGWidth (const double &width)
double	getCGWidth ()
bool	useCGMapping ()
CGFunctions::LucyXYZ	getCGFunction ()

Private Types
typedef O::ELEMENT	ELEMENT
	type of el_pt in O (e.g. RefineableQDPVDElement<3, 2>) More...

Private Member Functions
void	solveScaleCoupling (unsigned nStep)
	Solve scale-coupled problem. More...
void	initialiseCoupledElements ()
	initialises the coupledElements vector More...
void	computeOneTimeStepScaleCoupling (unsigned nStep)
	What to do in each time step. More...
void	updateCoupledElements ()
	Updates the coupledElements and coupledParticles vectors: finds which particles interact with which element. More...
Vector< Vector< double > >	computeProjectionMatrix (const CoupledElement &coupledElement)
	Computes projectionMatrix. More...
void	computeCouplingForce ()
	Computes coupling force for each element and particle. More...
void	computeNodalCouplingForces (const CoupledElement &coupledElement, const Vector< Vector< double >> &projectionMatrix)
	Computes coupling force for each node. More...
void	computeExternalForces (BaseParticle *p) override
	Applies coupling force to MercuryDPM in each time step. More...

Private Attributes
std::function< double(double x, double y, double z)>	couplingWeight_
std::vector< CoupledElement >	coupledElements_
	Vector of all coupled elements. More...
std::vector< CoupledParticle >	coupledParticles_
	The i-th element of this vector describes the coupling properties of the i-th DPM particle. More...
double	penalty_ = constants::NaN
	Penalty parameter, i.e., proportionality constant between velocity difference and coupling force. More...

Member Typedef Documentation

ELEMENT

template<class M , class O >

typedef O::ELEMENT ScaleCoupling< M, O >::ELEMENT

private

type of el_pt in O (e.g. RefineableQDPVDElement<3, 2>)

Member Function Documentation

computeCouplingForce()

template<class M , class O >

void ScaleCoupling< M, O >::computeCouplingForce ( )

inlineprivate

Computes coupling force for each element and particle.

    {
        // first loop over the coupled bulk el_pts
        for (const auto& coupledElement : coupledElements_)
        {
            // construct projection matrix, defined as projectionMatrix_{n,p} = N_{n,p}*V_p / sum_p(N_{n,p}*V_p) 
            // (shape function N, volume V, node n, particle p)
            Vector<Vector<double>> projectionMatrix = computeProjectionMatrix(coupledElement);
            // get pointer to the bulk el_pt
            ELEMENT* el_pt = coupledElement.el_pt;
            const unsigned nParticles = coupledElement.coupledParticles.size();
            // get number of nodes in the bulk el_pt
            const unsigned nnode = coupledElement.el_pt->nnode();
            // prepare vector of nodal velocities projected from particle centers
            const unsigned dim = el_pt->dim();
            // get coupling force at nodal positions, based on penalizing the difference in velocity
            computeNodalCouplingForces(coupledElement, projectionMatrix);
            // project coupling forces from nodes to particles (note projectionMatrix^{-1} = projectionMatrix^T)
            for (unsigned p=0; p < nParticles; p++) {
                BaseParticle* particle = coupledElement.coupledParticles[p];
                CoupledParticle& coupledParticle = coupledParticles_[particle->getIndex()];
                for (unsigned d=0; d < dim; d++) {
                    double forceComponent = 0;
                    for (unsigned n=0; n < nnode; n++) {
                        // note projectionMatrix^{-1} = projectionMatrix^T
                        forceComponent += projectionMatrix[n][p] * el_pt->get_coupling_residual(n, d);
                    }
                    coupledParticle.couplingForce.setComponent(d,forceComponent);
                }
                logger(VERBOSE,"Coupling force % on particle %", coupledParticle.couplingForce,p);
            }
        }
    }

References ScaleCoupling< M, O >::CoupledParticle::couplingForce, BaseObject::getIndex(), logger, n, Vec3D::setComponent(), and VERBOSE.

computeExternalForces()

template<class M , class O >

void ScaleCoupling< M, O >::computeExternalForces ( BaseParticle *  p )

inlineoverrideprivate

Applies coupling force to MercuryDPM in each time step.

    {
        if (!p->isFixed())
        {
            // Applying the force due to gravity (F = m.g)
            p->addForce(M::getGravity() * p->getMass());
            // add particle coupling force 1/w_i*f_i^c to the total force
            const CoupledParticle& cp = coupledParticles_[p->getIndex()];
            double couplingWeight = couplingWeight_(p->getPosition().X,p->getPosition().Y,p->getPosition().Z);
            p->addForce(cp.couplingForce / couplingWeight);
        }
    }

References BaseInteractable::addForce(), ScaleCoupling< M, O >::CoupledParticle::couplingForce, BaseObject::getIndex(), BaseParticle::getMass(), BaseInteractable::getPosition(), BaseParticle::isFixed(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

computeNodalCouplingForces()

template<class M , class O >

void ScaleCoupling< M, O >::computeNodalCouplingForces ( const CoupledElement &  coupledElement, const Vector< Vector< double >> &  projectionMatrix  )

inlineprivate

Computes coupling force for each node.

1) compute velocity difference at nodal positions

2) compute nodal coupling force, penalizing the difference in velocity

3) assign it to the coupled bulk element to be used by ELEMENT::fill_in_contribution_to_residuals(...)

    {
        // return if no particle in cell
        if (coupledElement.coupledParticles.empty()) return;
 
        // a few shortcut variables
        const ELEMENT* el_pt = coupledElement.el_pt;
        const unsigned nnode = el_pt->nnode();
        const unsigned dim = el_pt->dim();
        const unsigned nParticles = coupledElement.coupledParticles.size();
        
        Vector<Vector<double>> nodalVelocityDifference(nnode,Vector<double>(dim,0.0));
        // get projected nodal velocity field from particle velocities
        for (unsigned n=0; n < nnode; n++) {
            for (unsigned d=0; d < dim; d++) {
                double nodalVelocityDPM = 0;
                for (unsigned p=0; p < nParticles; p++) {
                    nodalVelocityDPM += projectionMatrix[n][p] * coupledElement.coupledParticles[p]->getVelocity().getComponent(d);
                }
                double nodalVelocityFEM = el_pt->dnodal_position_gen_dt(n,0,d);
                nodalVelocityDifference[n][d] = nodalVelocityDPM-nodalVelocityFEM;
            }
        }
 
        Vector<Vector<double> > nodalCouplingForces(nnode, Vector<double>(dim, 0.0));
        
        // Set up memory for the shape functions
        Shape psi(nnode);
        DShape dpsidxi(nnode,dim);
        const unsigned n_in = el_pt->integral_pt()->nweight();
        // Loop over the integration points
        for (unsigned in=0; in<n_in; in++)
        {
            // Get the integral weight
            double w = el_pt->integral_pt()->weight(in);
            // Call the derivatives of the shape functions (and get Jacobian)
            double J = el_pt->dshape_lagrangian_at_knot(in,psi,dpsidxi);
            // Loop over the test functions, nodes of the element
            for (unsigned n=0; n < nnode; n++)
            {
                // Loop over the nodes of the element again
                for (unsigned nn=0; nn < nnode; nn++)
                {
                    double volume = penalty_ * psi(nn) * psi(n) * w * J;
 
                    for (unsigned d=0; d < dim; d++)
                    {
                        //nodalVelocityDifference[n][d] = 1/penalty_/time_pt()->dt(0);
                        double displacement = nodalVelocityDifference[n][d] * O::time_pt()->dt(0);
                        nodalCouplingForces[n][d] -= volume * displacement;
                        //logger(INFO, "w % J % psi % % displacement % nodalCouplingForces % dt %", w, J, psi(nn), psi(n), displacement, nodalCouplingForces[n][d], time_pt()->dt(0));
                    }
                }
            }
        }
 
        coupledElement.el_pt->set_coupling_residual(nodalCouplingForces);
 
        logger(VERBOSE,"VelocityDifference % % %, coupling force % % % on element %, node 0, particles %",
               nodalVelocityDifference[0][0], nodalVelocityDifference[0][1], nodalVelocityDifference[0][2],
               nodalCouplingForces[0][0], nodalCouplingForces[0][1], nodalCouplingForces[0][2],
               el_pt, coupledElement.coupledParticles.size());
    }

References ScaleCoupling< M, O >::CoupledElement::coupledParticles, ScaleCoupling< M, O >::CoupledElement::el_pt, logger, n, and VERBOSE.

computeOneTimeStepScaleCoupling()

template<class M , class O >

void ScaleCoupling< M, O >::computeOneTimeStepScaleCoupling ( unsigned  nStep )

inlineprivate

What to do in each time step.

                                                         {
        updateCoupledElements();
        computeCouplingForce();
        this->solveOomph();
        this->solveMercury(nStep);
    };

computeProjectionMatrix()

template<class M , class O >

Vector<Vector<double> > ScaleCoupling< M, O >::computeProjectionMatrix ( const CoupledElement &  coupledElement )

inlineprivate

Computes projectionMatrix.

Construct mapping with FEM basis functions

    {
        Vector<Vector<double>> projectionMatrix;
        // get constants
        const unsigned nParticles = coupledElement.coupledParticles.size();
        const unsigned nnode = coupledElement.el_pt->nnode();
        const unsigned dim = coupledElement.el_pt->dim();
 
        // prepare projection matrix
        projectionMatrix.resize(nnode,Vector<double>(nParticles,0.0));
 
        // storage for position s and shape psi
        Vector<double> pos(dim,0.0);
        Shape psi(nnode);
        // loop over shape functions at the nodes
        // loop over shape functions at the particles
        for (unsigned p = 0; p < nParticles; p++) {
            BaseParticle* particle = coupledElement.coupledParticles[p];
            // get the shape function evaluated at the center of particle p
            coupledElement.el_pt->shape(coupledParticles_[particle->getIndex()].s, psi);
            for (unsigned n=0; n < nnode; n++) {
                projectionMatrix[n][p] = psi(n) * particle->getVolume();
            }
        }
        // normalize each row of the projection matrix to sum to one
        for (unsigned n=0; n < nnode; n++) {
            double sum = 0;
            for (unsigned p = 0; p<nParticles; p++) {
                sum += projectionMatrix[n][p];
            }
            if (sum != 0) {
                for (unsigned p = 0; p < nParticles; p++) {
                    projectionMatrix[n][p] /= sum;
                }
            }
        }
        return projectionMatrix;
    }

References ScaleCoupling< M, O >::CoupledElement::coupledParticles, ScaleCoupling< M, O >::CoupledElement::el_pt, BaseObject::getIndex(), BaseParticle::getVolume(), and n.

getCoupledElements()

template<class M , class O >

const std::vector<CoupledElement>& ScaleCoupling< M, O >::getCoupledElements ( )

inline

get coupled element

Todo:

check if set

                                                          {
        return coupledElements_;
    }

getCoupledParticles()

template<class M , class O >

const std::vector<CoupledParticle>& ScaleCoupling< M, O >::getCoupledParticles ( )

inline

get coupled particle

                                                            {
        return coupledParticles_;
    }

Referenced by ScaleCoupledBeam::actionsAfterSolve().

initialiseCoupledElements()

template<class M , class O >

void ScaleCoupling< M, O >::initialiseCoupledElements ( )

inlineprivate

initialises the coupledElements vector

                                     {
        // loop through all el_pts, check if any nodes have a fractional coupling weight.
        // If yes, add the el_pt to the coupledElements_
        auto nelement = O::mesh_pt()->nelement();
        for (int e=0; e < nelement; ++e) {
            bool coupledElement = false;
            auto* el_pt = dynamic_cast<ELEMENT*>(O::mesh_pt()->finite_element_pt(e));
            auto nnode = el_pt->nnode();
            Vector<double> nodalCouplingWeights(nnode);
            for (int n=0; n < nnode; ++n) {
                auto* n_pt = dynamic_cast<SolidNode*>(el_pt->node_pt(n));
                nodalCouplingWeights[n] = couplingWeight_(n_pt->xi(0), n_pt->xi(1), n_pt->xi(2));
                if (nodalCouplingWeights[n]>0) { // check for nodalCouplingWeights[n]<1 ?
                    coupledElement = true;
                    logger(VERBOSE,"Coupled element at % % %, weight %",n_pt->xi(0), n_pt->xi(1), n_pt->xi(2), nodalCouplingWeights[n]);
                }
            }
            if (coupledElement) {
                coupledElements_.push_back({el_pt});
                el_pt->set_coupling_weight(nodalCouplingWeights);
            }
        }
        logger(INFO,"% of % el_pts are coupled", coupledElements_.size(), nelement);
    }

References INFO, logger, n, and VERBOSE.

setCouplingWeight()

template<class M , class O >

void ScaleCoupling< M, O >::setCouplingWeight ( const std::function< double(double x, double y, double z)> &  couplingWeight )

inline

Sets weight function (determines which nodes/el_pts are in coupling zone)

                                                                                                  {
        couplingWeight_ = couplingWeight;
    }

Referenced by ScaleCoupledBeam::ScaleCoupledBeam().

setPenalty()

template<class M , class O >

void ScaleCoupling< M, O >::setPenalty ( double  penalty )

inline

Sets penalty parameter.

                                    {
        penalty_ = penalty;
    }

Referenced by ScaleCoupledBeam::ScaleCoupledBeam().

solveScaleCoupling() [1/2]

template<class M , class O >

void ScaleCoupling< M, O >::solveScaleCoupling ( )

inline

Computes nStep, the ratio of FEM and DEM time steps, then calls solveSurfaceCoupling(nStep)

    {
        // compute nStep
        unsigned nStep = O::getOomphTimeStep()/M::getTimeStep();
        if (nStep==0) {
            // if oomph has a smaller time step than Mercury
            logger(INFO, "Set nStep %, change mercuryTimeStep from % to %",
                   nStep, M::getTimeStep(), O::getOomphTimeStep());
            nStep = 1;
            M::setTimeStep(O::getOomphTimeStep());
        } else {
            logger(INFO, "Set nStep %, change oomphTimeStep from % to %",
                   nStep, O::getOomphTimeStep(), nStep * M::getTimeStep());
            O::setOomphTimeStep(nStep * M::getTimeStep());
        }
        // call solve routine
        solveScaleCoupling(nStep);
    }

References INFO, and logger.

Referenced by ScaleCoupledBeam::ScaleCoupledBeam().

solveScaleCoupling() [2/2]

template<class M , class O >

void ScaleCoupling< M, O >::solveScaleCoupling ( unsigned  nStep )

inlineprivate

Solve scale-coupled problem.

    {
        // Check main parameters are set
        logger.assert_always(O::getOomphTimeStep()>0, "Oomph time step not initialised");
        logger.assert_always(M::getTimeStep()>0, "Mercury time step not initialised");
        logger.assert_always(std::isfinite(penalty_), "Penalty parameter not set");
        logger.assert_always((bool)couplingWeight_, "Coupling weight is not set");
 
        // Initialise the coupled element vector
        initialiseCoupledElements();
        
        // This is the main loop advancing over time
        unsigned nDone = 0; //< last written file number
        while (M::getTime() < M::getTimeMax())
        {
            O::actionsBeforeOomphTimeStep();
            // solve the coupled problem for one time step
            computeOneTimeStepScaleCoupling(nStep);
            // 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
            if (M::getParticlesWriteVTK() && M::getVtkWriter()->getFileCounter() > nDone)
            {
                O::writeToVTK();
                nDone = M::getVtkWriter()->getFileCounter();
            }
        }
 
        // Write final output files
        M::finaliseSolve();
    }

References logger.

updateCoupledElements()

template<class M , class O >

void ScaleCoupling< M, O >::updateCoupledElements ( )

inlineprivate

Updates the coupledElements and coupledParticles vectors: finds which particles interact with which element.

1) Find min and max of coupled region

2) update the coupledParticles vector, and the particles vector of each coupledElement

    {
        Vector<double> min {constants::inf, constants::inf, constants::inf};
        Vector<double> max {-constants::inf, -constants::inf, -constants::inf};
        // loop over all coupled elements
        for (auto& coupledElement : coupledElements_)
        {
            ELEMENT* el_pt = coupledElement.el_pt;
            const unsigned nnode = el_pt->nnode();
            const unsigned dim = el_pt->dim(); // /todo I need to use el_pt->ndim=3, not el_pt->ndim=0, why?
            // Find min and max
            for (int n=0; n<nnode; ++n) {
                auto* n_pt = dynamic_cast<SolidNode *>(el_pt->node_pt(n));
                for (int d = 0; d < dim; ++d) {
                    min[d] = std::min(min[d], n_pt->x(d));
                    max[d] = std::max(max[d], n_pt->x(d));
                }
            }
        }
        logger(VERBOSE,"Coupling zone: % < x < %, % < y < %, % < z < % ",min[0],max[0],min[1],max[1],min[2],max[2]);
 
        // resize the coupledParticles vector
        coupledParticles_.resize(M::particleHandler.getSize());
        // now set the coupledElement pointer for each coupledParticles, and the particles vector for each coupledElement
        for (int p=0; p<coupledParticles_.size(); ++p) {
            CoupledParticle& coupledParticle = coupledParticles_[p];
            Vec3D pos = M::particleHandler.getObject(p)->getPosition();
            // this is the pointer we need to set
            if (pos.X>=min[0] && pos.X<=max[0] && pos.Y>=min[1] && pos.Y<=max[1] && pos.Z>=min[2] && pos.Z<=max[2]) {
                logger(VERBOSE,"Particle % of % is in coupling zone",p, M::particleHandler.getSize());
                // if near the coupling zone
                Vector<double> x {pos.X, pos.Y, pos.Z};
                Vector<double>& s = coupledParticle.s;
                GeomObject* geom_obj_pt;
                // if el_pt is already set and still valid, keep
                if (coupledParticle.coupledElement_pt != nullptr) {
                    coupledParticle.coupledElement_pt->el_pt->locate_zeta(x, geom_obj_pt, s);
                    if (geom_obj_pt!=nullptr) {
                        logger(VERBOSE,"Particle % remains coupled with element %",p, coupledParticle.coupledElement_pt->el_pt);
                        continue;
                    }
                }
                // otherwise, search for the right element 
                for (auto& coupledElement : coupledElements_) {
                    coupledElement.el_pt->locate_zeta(x,geom_obj_pt,s);
                    if (geom_obj_pt!=nullptr) {
                        logger(VERBOSE,"Coupled particle % to element % at % % %",p, &coupledElement, s[0], s[1], s[2]);
                        //\todo If a containing element is found, we stop; note though, particles on the boundary could belong to multiple elements; also not clear to me whether force should be applied to global basis function or local
                        coupledParticle.setCoupledElement(&coupledElement, M::particleHandler.getObject(p));
                        break;
                    }
                }
            } else {
                // if away from the coupling, set coupledElement to null
                if (coupledParticle.coupledElement_pt) logger(VERBOSE, "Removed particle % from element %", p, coupledParticle.coupledElement_pt);
                coupledParticle.removeCoupledElement(M::particleHandler.getObject(p));
            }
        }
    }

References ScaleCoupling< M, O >::CoupledParticle::coupledElement_pt, ScaleCoupling< M, O >::CoupledElement::el_pt, constants::inf, logger, n, ScaleCoupling< M, O >::CoupledParticle::removeCoupledElement(), ScaleCoupling< M, O >::CoupledParticle::s, ScaleCoupling< M, O >::CoupledParticle::setCoupledElement(), VERBOSE, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Member Data Documentation

coupledElements_

template<class M , class O >

std::vector<CoupledElement> ScaleCoupling< M, O >::coupledElements_

private

Vector of all coupled elements.

coupledParticles_

template<class M , class O >

std::vector<CoupledParticle> ScaleCoupling< M, O >::coupledParticles_

private

The i-th element of this vector describes the coupling properties of the i-th DPM particle.

couplingWeight_

template<class M , class O >

std::function<double(double x,double y,double z)> ScaleCoupling< M, O >::couplingWeight_

private

coupling weight: 0: pure FEM, 1: pure DEM, (0,1): overlap region

penalty_

template<class M , class O >

double ScaleCoupling< M, O >::penalty_ = constants::NaN

private

Penalty parameter, i.e., proportionality constant between velocity difference and coupling force.

---

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

• ScaleCoupling.h