BaseCoupling.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef BASE_COUPLING_H
#define BASE_COUPLING_H
 
#include "Math/Vector.h"
#include "Walls/TriangleWall.h"
#include "Species/LinearViscoelasticSpecies.h"
#include "CG/Functions/Gauss.h"
#include "CG/Functions/Lucy.h"
#include "CG/Functions/Heaviside.h"
#include "Interactions/BaseInteraction.h"
#include "Logger.h"
 
template<class M, class O>
class BaseCoupling : public M, public O
{
//protected:
//    // A reference to the Mercury problem
//    Mercury3D& m = this;
//    // A reference to the Oomph problem
//    SolidProblem<typename OomphProblem::ELEMENT_TYPE>& o = this;
 
public:
    
    BaseCoupling() = default;
    
    void setName(std::string name) {
        M::setName(name);
        O::setName(name);
    }
 
    std::string getName() const {
        return M::getName();
    }
 
    void removeOldFiles() const {
        M::removeOldFiles();
        O::removeOldFiles();
    }
    
    void writeEneTimeStep(std::ostream& os) const override
    {
        if (M::eneFile.getCounter() == 1 || M::eneFile.getFileType() == FileType::MULTIPLE_FILES ||
                M::eneFile.getFileType() == FileType::MULTIPLE_FILES_PADDED)
        {
            writeEneHeader(os);
        }
 
        const Mdouble m = M::particleHandler.getMass();
        const Vec3D com = M::particleHandler.getMassTimesPosition();
        //Ensure the numbers fit into a constant width column: for this we need the precision given by the operating system,
        //plus a few extra characters for characters like a minus and scientific notation.
        const static long int width = os.precision() + 6;
        os << std::setw(width) << M::getTime()
//           << " " << std::setw(width) << getCoupledMass()
           << " " << std::setw(width) << M::particleHandler.getMomentum().getX()
           << " " << std::setw(width) << M::particleHandler.getMomentum().getY()
           << " " << std::setw(width) << M::particleHandler.getMomentum().getZ()
           << " " << std::setw(width) << M::particleHandler.getAngularMomentum().getX()
           << " " << std::setw(width) << M::particleHandler.getAngularMomentum().getY()
           << " " << std::setw(width) << M::particleHandler.getAngularMomentum().getZ()
           << " " << std::setw(width) << -Vec3D::dot(M::getGravity(), com)
           << " " << std::setw(width) << M::particleHandler.getKineticEnergy()
           << " " << std::setw(width) << M::particleHandler.getRotationalEnergy()
           //<< " " << std::setw(width) << M::getElasticEnergyCoupled()
           // we need to write x, y and z coordinates separately, otherwise the width of the columns is incorrect
           << " " << std::setw(width)
           << ( m == 0 ? constants::NaN : com.X / m ) //set to nan because 0/0 implementation in gcc and clang differs
           << " " << std::setw(width) << ( m == 0 ? constants::NaN : com.Y / m )
           << " " << std::setw(width) << ( m == 0 ? constants::NaN : com.Z / m )
           << std::endl;
    }
    
    void writeEneHeader(std::ostream& os) const override
    {
        //only write if we don't restart
        if (M::getAppend()) {
            return;
        }
        
        long width = os.precision() + 6;
        os << std::setw(width)
           << "time " << std::setw(width)
           << "mass " << std::setw(width)
           << "momentum_X " << std::setw(width)
           << "momentum_Y " << std::setw(width)
           << "momentum_Z " << std::setw(width)
           << "angMomentum_X " << std::setw(width)
           << "angMomentum_Y " << std::setw(width)
           << "angMomentum_Z " << std::setw(width)
           << "gravitEnergy " << std::setw(width) //gravitational potential energy
           << "traKineticEnergy " << std::setw(width) //translational kinetic energy
           << "rotKineticEnergy " << std::setw(width) //rotational kE
           << "elasticEnergy " << std::setw(width)
           << "centerOfMassX " << std::setw(width)
           << "centerOfMassY " << std::setw(width)
           << "centerOfMassZ\n";
    }
    
    void solveOomph(int max_adapt = 0)
    {
        O::actionsBeforeOomphTimeStep();
        // the coupled codes seem to not work if newton_solve is used
        //this->newton_solve(this->time_pt()->dt());
        if(max_adapt <= 0)
        {
            this->unsteady_newton_solve(this->time_pt()->dt());
        }
        else
        {
            this->unsteady_newton_solve(this->time_pt()->dt(), max_adapt, false);
        }
        checkResidual();
    }
 
    // check whether the maximum residual is zero. If so, the computation has failed and the code will exit.
    void checkResidual(){
        oomph::DoubleVector residual;
        O::get_residuals(residual);
        double residualNorm = residual.max();
        if (residualNorm == 0) {
            logger(ERROR, "Maximum residuals %; the computation has failed and the code will exit", residualNorm);
        }
    }
 
    void solveMercury(unsigned long nt)
    {
        for (int n = 0; n < nt; ++n)
        {
            M::computeOneTimeStep();
        }
    }
    
    inline void setCGWidth(const double& width)
    {
        if (width == 0)
        {
            CGMapping_ = false;
        }
        else
        {
            CGMapping_ = true;
            // note CG function needs to be nondimensionalized
            // fixme CG width is not set with respect to particle radius. Should we do that?
            cgFunction_.setWidth(width);
        }
    }
 
    inline double getCGWidth()
    {
        return cgFunction_.getWidth();
    }
    
    inline bool useCGMapping()
    { return CGMapping_; }
    
    inline CGFunctions::LucyXYZ getCGFunction()
    { return cgFunction_; }
 
private:
    // flag to construct mapping with FEM basis functions (default) or DPM coarse graining
    bool CGMapping_ = false;
    // coarse-graining functions and coordinates
    CGFunctions::LucyXYZ cgFunction_;
};
 
#endif //BASE_COUPLING_H