mesh.h

Go to the documentation of this file.

#ifndef mesh_h
#define mesh_h
 
#include <algorithm>
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#include <execution>
#pragma GCC diagnostic pop
 
#include "facette.h"
#include "node.h"
#include "tetra.h"
#include "feellgoodSettings.h"
 
namespace Mesh
    {
 
class mesh
    {
public:
    inline mesh(Settings const &mySets )
        {
        readMesh(mySets);
        indexReorder(mySets.paramTetra);
 
        if (mySets.verbose)
            { std::cout << "  reindexed\n"; }
 
        double xmin = minNodes(Nodes::IDX_X);
        double xmax = maxNodes(Nodes::IDX_X);
 
        double ymin = minNodes(Nodes::IDX_Y);
        double ymax = maxNodes(Nodes::IDX_Y);
 
        double zmin = minNodes(Nodes::IDX_Z);
        double zmax = maxNodes(Nodes::IDX_Z);
 
        l = Eigen::Vector3d(xmax - xmin, ymax - ymin, zmax - zmin);
        c = Eigen::Vector3d(0.5 * (xmax + xmin), 0.5 * (ymax + ymin), 0.5 * (zmax + zmin));
 
        // Find the longest axis of the sample.
        Nodes::index long_axis;
        if (l.x() > l.y())
            {
            if (l.x() > l.z())
                long_axis = Nodes::IDX_X;
            else
                long_axis = Nodes::IDX_Z;
            }
        else
            {
            if (l.y() > l.z())
                long_axis = Nodes::IDX_Y;
            else
                long_axis = Nodes::IDX_Z;
            }
        sortNodes(long_axis);
 
        vol = std::transform_reduce(EXEC_POL, tet.begin(), tet.end(), 0.0, std::plus{},
                                    [](Tetra::Tet const &te) { return te.calc_vol(); });
        }
 
    inline int getNbNodes(void) const { return node.size(); }
 
    inline int getNbFacs(void) const { return fac.size(); }
 
    inline int getNbTets(void) const { return tet.size(); }
 
    inline const Eigen::Vector3d getNode_p(const int i) const { return node[i].p; }
    
    inline const Eigen::Vector3d getNode_u(const int i) const { return node[i].get_u(Nodes::NEXT); }
 
    inline const Eigen::Vector3d getNode_v(const int i) const { return node[i].get_v(Nodes::NEXT); }
 
    inline double getProj_ep(const int i) const {return node[i].proj_ep();}
 
    inline double getProj_eq(const int i) const {return node[i].proj_eq();}
 
    inline void set_node_u0(const int i, Eigen::Vector3d const &val)
        { node[i].d[Nodes::CURRENT].u = val; }
 
    inline void set_node_zero_v(const int i) { node[i].d[Nodes::NEXT].v.setZero(); }
 
    void infos(void) const;
 
    inline void setBasis(const double r)
        {
        std::for_each(EXEC_POL, node.begin(), node.end(),
                      [&r](Nodes::Node &nod) { nod.setBasis(r); });
        }
 
    void updateNodes(Eigen::Ref<Eigen::VectorXd> X, const double dt);
 
    inline void evolution(void)
        {
        std::for_each(EXEC_POL, node.begin(), node.end(),
                      [](Nodes::Node &nod) { nod.evolution(); });
        }
 
    Eigen::Vector3d c;
 
    Eigen::Vector3d l;
 
    double vol;
 
    std::vector<Facette::Fac> fac;
 
    std::vector<Tetra::Tet> tet;
 
    double readSol(bool VERBOSE ,
                   const std::string fileName );
 
    inline void init_distrib(Settings const &mySets )
        {
        std::for_each(node.begin(),node.end(),[&mySets](Nodes::Node &n)
                                              {
                                              n.d[Nodes::CURRENT].u = mySets.getMagnetization(n.p);
                                              n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
                                              n.d[Nodes::NEXT].phi = 0.;
                                              n.d[Nodes::NEXT].phiv = 0.;
                                              } );
        }
 
    double avg(std::function<double(Nodes::Node, Nodes::index)> getter ,
               Nodes::index d ) const;
 
    void savesol(const int precision ,
                 const std::string fileName ,
                 std::string const &metadata ) const;
 
    bool savesol(const int precision ,
                 const std::string fileName ,
                 std::string const &metadata ,
                 std::vector<double> const &val ) const;
 
    inline void set(const int i ,
                    std::function<void(Nodes::Node &, const double)> what_to_set ,
                    const double val )
        { what_to_set(node[i], val); }
 
private:
    std::vector<Nodes::Node> node;
 
    std::vector<int> node_index;
 
    void checkMeshFile(Settings const &mySets );
    
    void readNodes(Settings const &mySets );
 
    void readTetraedrons(Settings const &mySets );
 
    void readTriangles(Settings const &mySets );
 
    void readMesh(Settings const &mySets );
 
    double doOnNodes(const double init_val ,
                     const Nodes::index coord ,
                     std::function<bool(double, double)> whatToDo ) const;
 
    inline double minNodes(const Nodes::index coord ) const
        {
        return doOnNodes(__DBL_MAX__, coord, [](double a, double b) { return a < b; });
        }
 
    inline double maxNodes(const Nodes::index coord ) const
        {
        return doOnNodes(__DBL_MIN__, coord, [](double a, double b) { return a > b; });
        }
 
    void indexReorder(std::vector<Tetra::prm> const &prmTetra);
 
    void sortNodes(Nodes::index long_axis );
 
    }; // end class mesh
 
    }  // end namespace Mesh
#endif