dev_fdaPDE/Mesh_imp.h at master · simonepanzeri/dev_fdaPDE · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
//
// Created by simonepanzeri on 25/11/2021.
//

#ifndef DEV_FDAPDE_MESH_IMP_H
#define DEV_FDAPDE_MESH_IMP_H

template <UInt ORDER, UInt mydim, UInt ndim>
MeshHandler<ORDER, mydim, ndim>::MeshHandler(const RNumericMatrix& points, const RIntegerMatrix& sides,
                                             const RIntegerMatrix& elements, const RIntegerMatrix& neighbors,
                                             UInt search) :
    points_(points), sides_(sides), elements_(elements), neighbors_(neighbors), search_(search) {
    //if(search == 2)
        //tree_ptr_ = make_unique<const ADTree<meshElement>> (Rmesh);
}

template <UInt ORDER, UInt mydim, UInt ndim>
Point<ndim> MeshHandler<ORDER, mydim, ndim>::getPoint(const UInt id) const {
    return Point<ndim>(id, points_);
}

template <UInt ORDER, UInt mydim, UInt ndim>
typename MeshHandler<ORDER, mydim, ndim>::meshElement MeshHandler<ORDER, mydim, ndim>::getElement(const UInt id) const {
    typename meshElement::elementPoints elPoints;
    for (int j = 0; j < how_many_nodes(ORDER, mydim); ++j)
        elPoints[j] = getPoint(elements_(id,j));
    return meshElement(id, elPoints);
}

template <UInt ORDER, UInt mydim, UInt ndim>
typename MeshHandler<ORDER, mydim, ndim>::meshElement MeshHandler<ORDER, mydim, ndim>::getNeighbors(const UInt id_element, const UInt number) const {
    int id_neighbor{neighbors_(id_element, number)};
    //return empty element if "neighbor" not present (out of boundary!)
    return (id_neighbor == -1) ? meshElement() : getElement(id_neighbor);
}

template <UInt ORDER, UInt mydim, UInt ndim>
template <bool isManifold>
typename std::enable_if<!isManifold, typename MeshHandler<ORDER, mydim, ndim>::meshElement>::type
MeshHandler<ORDER, mydim, ndim>::findLocation(const Point<ndim>& point) const {
    if(search_ == 2)
        return findLocationTree(point);
    else if(search_ == 3)
        return findLocationWalking(point, getElement(0));
    else
        return findLocationNaive(point);
}

template <UInt ORDER, UInt mydim, UInt ndim>
template <bool isManifold>
typename std::enable_if<isManifold, typename MeshHandler<ORDER, mydim, ndim>::meshElement>::type
MeshHandler<ORDER, mydim, ndim>::findLocation(const Point<ndim>& point) const {
    if(search_ == 2)
        return findLocationTree(point);
    else
        return findLocationNaive(point);
}

template <UInt ORDER, UInt mydim, UInt ndim>
typename MeshHandler<ORDER, mydim, ndim>::meshElement MeshHandler<ORDER, mydim, ndim>::findLocationNaive(const Point<ndim>& point) const {
    for(UInt id = 0; id < elements_.nrows(); ++id){
        meshElement current_element{getElement(id)};
        if(current_element.isPointInside(point))
            return current_element;
    }
    return meshElement(); //default element with NVAL ID
}

//Visibility walk algorithm which uses barycentric coordinate [Sundareswara et al]
//Starting triangles usually n^(1/3) points
template <UInt ORDER, UInt mydim, UInt ndim>
typename MeshHandler<ORDER, mydim, ndim>::meshElement MeshHandler<ORDER, mydim, ndim>::findLocationWalking(const Point<ndim>& point, const meshElement& starting_element) const {
    static_assert(mydim == ndim, "ERROR! WALKING SEARCH CANNOT BE USED ON MANIFOLD MESHES! See mesh_imp.h");
    meshElement current_element{starting_element};
    //Test for found Element, or out of border
    while(current_element.hasValidId() && !current_element.isPointInside(point))
        current_element = getNeighbors(current_element.getId(), current_element.getPointDirection(point));

    return current_element;
}

template <UInt ORDER, UInt mydim, UInt ndim>
typename MeshHandler<ORDER, mydim, ndim>::meshElement MeshHandler<ORDER, mydim, ndim>::findLocationTree(const Point<ndim>& point) const {
    std::set<int> found;
    std::vector<Real> region;
    region.reserve(2*ndim);

    for (UInt i = 0; i < ndim; ++i){
        region.push_back(point[i]);
    }
    for (UInt i = 0; i < ndim; ++i){
        region.push_back(point[i]);
    }

    if(!tree_ptr_->search(region, found)) {
        return meshElement();
    }

    for (const auto &i : found) {
        const UInt index = tree_ptr_->pointId(i);
        meshElement tmp = getElement(index);
        if(tmp.isPointInside(point)) {
            return tmp;
        }
    }
    return meshElement();
}

template <UInt ORDER, UInt mydim, UInt ndim>
void MeshHandler<ORDER, mydim, ndim>::printPoints(std::ostream& os) const
{
    os << "# Nodes: " << points_.nrows() << std::endl;
    for(UInt i = 0; i < points_.nrows(); ++i)
        os << getPoint(i);
}

template <UInt ORDER, UInt mydim, UInt ndim>
void MeshHandler<ORDER, mydim, ndim>::printElements(std::ostream& os) const
{
    os << "# Triangles: " << elements_.nrows() << std::endl;
    for (UInt i = 0; i < elements_.nrows(); ++i )
        os << getElement(i);
}

template <UInt ORDER, UInt mydim, UInt ndim>
void MeshHandler<ORDER, mydim, ndim>::printNeighbors(std::ostream& os) const
{
    os << "# Neighbors list: " << elements_.nrows() << std::endl;
    for (UInt i = 0; i < elements_.nrows(); ++i){
        for( UInt j = 0; j < mydim+1; ++j)
            os << neighbors(i,j) << " ";
        os << std::endl;
    }
}

template <UInt ORDER, UInt mydim, UInt ndim>
void MeshHandler<ORDER, mydim, ndim>::printTree(std::ostream & os) const
{
    os << "# Tree characteristic: " << std::endl;
    if (tree_ptr_)
        os << *tree_ptr_ << std::endl;
    else
        os << "No tree!" << std::endl;
}

#endif //DEV_FDAPDE_MESH_IMP_H