mesh.h

Go to the documentation of this file.

#ifndef mesh_h
#define mesh_h
 
#include <cmath>
#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 &mySets )
        : paramTetra(mySets.paramTetra)
        {
        readMesh(mySets);
        indexReorder();
 
        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);
 
        // Compute the per-region volumes and the total volume.
        std::for_each(tet.begin(), tet.end(),
                [this](Tetra::Tet const &te)
                    {
                    paramTetra[te.idxPrm].volume += te.calc_vol();
                    });
        vol = std::transform_reduce(paramTetra.begin(), paramTetra.end(), 0.0,
                std::plus<>(), [](const Tetra::prm &region){ return region.volume; });
 
        // Build the list of all the mesh edges.
        edges.reserve(tet.size() * 6);  // 6 (non unique) edges per tetrahedron
        std::for_each(tet.begin(), tet.end(),
                [this](Tetra::Tet const &te)
                    {
                    for (int i = 0; i < 3; ++i)
                        {
                        for (int j = i + 1; j < 4; ++j)
                            { edges.push_back(std::minmax(te.ind[i], te.ind[j])); }
                        }
                    });
        std::sort(EXEC_POL, edges.begin(), edges.end());
        auto last = std::unique(EXEC_POL, edges.begin(), edges.end());
        edges.erase(last, edges.end());
        edges.shrink_to_fit();  // save memory, as this could be quite big
        }
 
    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); });
        }
 
    inline void updateNode(int i, double vp, double vq, const double dt)
        { node[i].make_evol(vp*gamma0, vq*gamma0, 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;
 
    std::vector<Tetra::prm> &paramTetra;
 
    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 ,
               int region = -1 ) const;
 
    double max_angle() const
        {
        double min_dot_product = std::transform_reduce(EXEC_POL, edges.begin(), edges.end(), 1.0,
                [](double a, double b){ return std::min(a, b); },
                [this](const std::pair<int, int> edge)
                    {
                    Eigen::Vector3d m1 = getNode_u(edge.first);
                    Eigen::Vector3d m2 = getNode_u(edge.second);
                    return m1.dot(m2);
                    });
        return std::acos(min_dot_product);
        }
 
    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<std::pair<int, int>> edges;
 
    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();
 
    void sortNodes(Nodes::index long_axis );
 
    }; // end class mesh
 
    }  // end namespace Mesh
#endif