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 "triangle.h"
#include "node.h"
#include "tetra.h"
#include "settings.h"
 
namespace Mesh
    {
    using Edge = std::pair<int, int>;
 
class mesh
    {
public:
    mesh(Settings &mySets ,
         bool simplified = false )
        : paramTetra(mySets.paramTetra), volumeRegions(mySets.paramTetra.size())
        {
        readMesh(mySets);
        if (simplified)
            { return; }
        if (!controlTriangles())
            { exit(1); }
 
        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);
 
        totalMagVol = 0;
        // Compute the per-region volumes and the total volume.
        std::for_each(tet.begin(), tet.end(),
                [this](const Tetra::Tet &te)
                    {
                    double vol_tet = te.calc_vol();
                    paramTetra[te.idxPrm].volume += vol_tet;
                    if(isMagnetic(te))
                        totalMagVol += vol_tet;
                    });
        vol = std::transform_reduce(paramTetra.begin(), paramTetra.end(), 0.0,
                std::plus<>(), [](const Tetra::prm &region){ return region.volume; });
 
        // Build the list of tetrahedrons for each region.
        std::for_each(tet.begin(), tet.end(), [this](const Tetra::Tet &te)
            {
            volumeRegions[te.idxPrm].push_back(te.idx);
            });
 
        // 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](const Tetra::Tet &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
        magNode.resize(node.size());
        std::fill(magNode.begin(),magNode.end(),false);
        std::for_each(tet.begin(),tet.end(),[this](Tetra::Tet &te)
                {
                if(isMagnetic(te))
                    {
                    magTet.push_back(te.idx);
                    for (int i=0;i<4;i++)
                        { magNode[te.ind[i]] = true; }
                    }
                });
 
        for(unsigned int i=0;i<tri.size();i++)
            {
            if (isMagnetic(tri[i]) && !mySets.paramTriangle[tri[i].idxPrm].suppress_charges)
                { magTri.push_back(i); }
            }
 
        if(mySets.getFieldType() == R4toR3)
            { setExtSpaceField(mySets); }
        }
 
    mesh(const mesh &) = delete;
 
    mesh &operator=(const mesh &) = delete;
 
    inline int getNbNodes(void) const { return node.size(); }
 
    inline int getNbTris(void) const { return tri.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, const Eigen::Vector3d &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(const int i, const double vp, const 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;
 
    double totalMagVol;
 
    std::vector<Triangle::Tri> tri;
 
    std::vector<Tetra::Tet> tet;
 
    std::vector<Tetra::prm> &paramTetra;
 
    double readSol(bool VERBOSE ,
                   const std::string& fileName );
 
    inline void init_distrib(const Settings &mySets )
        {
        for (int nodeIdx = 0; nodeIdx < int(node.size()); ++nodeIdx)
            {
            Nodes::Node &n = node[nodeIdx];
            n.d[Nodes::NEXT].phi = 0.;
            n.d[Nodes::NEXT].phiv = 0.;
 
            // A non-magnetic node's magnetization is a vector of NAN.
            if (!magNode[nodeIdx])
                {
                n.d[Nodes::CURRENT].u = Eigen::Vector3d(NAN, NAN, NAN);
                n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
                continue;
                }
 
            // If the initial magnetization depends only on the node position, we do not need to
            // build the list of region names.
            if (mySets.getMagType() == POSITION_ONLY)
                {
                n.d[Nodes::CURRENT].u = mySets.getMagnetization(n.p);
                n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
                continue;
                }
 
            // Get the list of region indices this node belongs to.
            std::set<int> nodeRegions;
            // skip region 0 (__default__) which should be empty
            for (size_t regIdx = 1; regIdx < volumeRegions.size(); ++regIdx)
                {
                const std::vector<int> &regionTetras = volumeRegions[regIdx];
                bool node_in_region = std::any_of(EXEC_POL,
                    regionTetras.begin(), regionTetras.end(),
                    [this, nodeIdx](const int tetIdx)
                    {
                    const Tetra::Tet &tetrahedron = tet[tetIdx];
                    for (int i = 0; i < Tetra::N; ++i) // node of tetrahedron tetIdx
                        {
                        if (tetrahedron.ind[i] == nodeIdx)
                            { return true; }
                        }
                    return false;
                    });
                if (node_in_region)
                    { nodeRegions.insert(regIdx); }
                }
 
            // Get the list of region names this node belongs to.
            std::vector<std::string> region_names;
            region_names.resize(nodeRegions.size());
            std::transform(nodeRegions.begin(), nodeRegions.end(), region_names.begin(),
                [this](const int regIdx){ return paramTetra[regIdx].regName; });
 
            n.d[Nodes::CURRENT].u = mySets.getMagnetization(n.p, region_names);
            n.d[Nodes::NEXT].u = n.d[Nodes::CURRENT].u;
            }
        }
 
    double thiele(const int region ) const;
 
    double avg(const 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,
                [](const double a, const double b){ return std::min(a, b); },
                [this](const Edge &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 ,
             const std::string &metadata ,
             bool withSpinAcc ,
             std::vector<Eigen::Vector3d> &s ) const;
 
    inline void set(const int i ,
                    const std::function<void(Nodes::Node &, const double)>& what_to_set ,
                    const double val )
        { what_to_set(node[i], val); }
 
    inline int getNodeIndex(const int i) const { return node_index[i]; }
 
    inline bool isMagnetic(const Tetra::Tet &t) const { return (paramTetra[t.idxPrm].Ms > 0); }
 
    inline bool isMagnetic(const Triangle::Tri &f)
        { return (magNode[f.ind[0]] && magNode[f.ind[1]] && magNode[f.ind[2]]); }
 
    bool controlTriangles();
 
    std::vector<Edge> edges;
 
    std::vector<bool> magNode;
 
    std::vector<int> magTet;
 
    std::vector<int> magTri;
 
    std::vector< Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> > extSpaceField;
 
private:
    std::vector<Nodes::Node> node;
 
    std::vector<int> node_index;
 
    std::vector<std::vector<int>> volumeRegions;
 
    void checkMeshFile(const Settings &mySets );
    
    void readNodes(const Settings &mySets );
 
    void readTetraedrons(const Settings &mySets );
 
    void readTriangles(const Settings &mySets );
 
    void readMesh(const Settings &mySets );
 
    double doOnNodes(const double init_val ,
                     const Nodes::index coord ,
                     const std::function<bool(double, double)>& whatToDo ) const;
 
    inline double minNodes(const Nodes::index coord ) const
        {
        return doOnNodes(__DBL_MAX__, coord, [](const double a, const double b)
                        { return a < b; });
        }
 
    inline double maxNodes(const Nodes::index coord ) const
        {
        return doOnNodes(__DBL_MIN__, coord, [](const double a, const double b)
                        { return a > b; });
        }
 
    void sortNodes(Nodes::index long_axis );
 
    bool findErrInTriangle(const int nbSurfTri, const int nbTetraFaces,
                         const std::pair<int,int> pairIdVolRegs, const int idSurfReg);
 
    double surface(std::vector<int> &triIndices) const;
 
    void setExtSpaceField(Settings &s );
    }; // end class mesh
    }  // end namespace Mesh
#endif