dev_fdaPDE/Domain_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
//
// Created by simonepanzeri on 25/11/2021.
//

#ifndef DEV_FDAPDE_DOMAIN_IMP_H
#define DEV_FDAPDE_DOMAIN_IMP_H

template<class T>
Real Domain<T>::tolerance_ = 1.e-3;

template<class T>
Real Domain<T>::mindiff_ = std::numeric_limits<Real>::min();

template<class T>
Domain<T>::Domain(){ //default domain: whole domain scale from [0,1]
    origin_.resize(T::dt());
    std::vector<Real> tmp(T::dt(),1);
    scalingfactors_ = tmp;
}

template<class T>
Domain<T>::Domain(std::vector<std::vector<Real>> const& coord) {
    int ndimp = T::dp();
    origin_.resize(T::dt());
    scalingfactors_.resize(T::dt());

    /* Find geometric limits.
     *
     * If loops are put outside for loops in order to improve performance.
     */
    if(ndimp == int(coord.size())) {
        // T is equal to Point<3>, Element<mydim, ndim>, Box<NDIMP>
        if(T::dp() == T::dt()) {
            // T is equal to Point<3>
            for(int i = 0; i < ndimp; ++i) {
                origin_[i] = *(std::min_element(coord[i].begin(), coord[i].end()));
                scalingfactors_[i] = *(std::max_element(coord[i].begin(), coord[i].end()));

                // Add the tolerance.
                Real delta = scalingfactors_[i] - origin_[i];
                origin_[i] -= delta*gettolerance();
                scalingfactors_[i] += delta*gettolerance();

                delta = scalingfactors_[i] - origin_[i];
                scalingfactors_[i] = 1./std::max(delta, getmindiff());
            }
        } else {
            // T is equal to Element<mydim, ndim>, Box<NDIMP>
            for(int i = 0; i < ndimp; ++i) {
                origin_[i] = *(std::min_element(coord[i].begin(), coord[i].end()));
                scalingfactors_[i] = *(std::max_element(coord[i].begin(), coord[i].end()));

                // Add the tolerance.
                double delta = scalingfactors_[i] - origin_[i];
                origin_[i] -= delta*gettolerance();
                scalingfactors_[i] += delta*gettolerance();

                delta = scalingfactors_[i] - origin_[i];
                scalingfactors_[i] = 1./std::max(delta, getmindiff());

                /* Repeat the limits because tree dimension is in fact 2 * physical space dimension
                 * because the tree contains triangle or tetrahedron bounding boxes.
                 */
                origin_[i + ndimp] = origin_[i];
                scalingfactors_[i + ndimp] = scalingfactors_[i];
            }
        }
    } else {
        //Special case of Point<2> because Point has default myDim=3
        if(T::dp() == T::dt()) {
            int realndimp = int(coord.size()); // in case of 2-dimen Point (default dimension is 3 and real dimension is 2)
            //need to resize the vectors of origin_, scalingfactors_
            origin_.resize(realndimp);
            scalingfactors_.resize(realndimp);

            for(int i = 0; i < realndimp; ++i) {
                origin_[i] = *(std::min_element(coord[i].begin(), coord[i].end()));
                scalingfactors_[i] = *(std::max_element(coord[i].begin(), coord[i].end()));

                // Add the tolerance.
                Real delta = scalingfactors_[i] - origin_[i];
                origin_[i] -= delta*gettolerance();
                scalingfactors_[i] += delta*gettolerance();

                delta = scalingfactors_[i] - origin_[i];
                scalingfactors_[i] = 1./std::max(delta, getmindiff());
            }
        }
    }
}

template<class T>
std::ostream& operator<<(std::ostream& ostr, Domain<T> const& d) {
    ostr << "Domain" << std::endl;
    ostr << "------" << std::endl;

    int dimp = d.origin_.size(); //not recommended to use T::dp() because of exception case of Point<2>
    for(int i = 0; i < dimp; ++i)
        ostr << "x" << i+1 << "_min = " << d.origin_[i] << std::endl;
    ostr << std::endl;

    for(int i = 0; i < dimp; ++i)
        ostr << "x" << i+1 << "_max = " << d.origin_[i]+1./d.scalingfactors_[i] << std::endl;
    ostr << "------" << std::endl;

    return ostr;
}

#endif //DEV_FDAPDE_DOMAIN_IMP_H