FiniteElementMesh.h 8.26 KB
Newer Older
1 2 3
#ifndef OERSTED_FINITEELEMENTMESH_H
#define OERSTED_FINITEELEMENTMESH_H

JasonPries's avatar
JasonPries committed
4
#include <memory>
5 6
#include <vector>

7 8
#include "Mesh.hpp"

JasonPries's avatar
JasonPries committed
9
#include "DiscreteBoundary.h"
10
#include "MatrixGroup.h"
11 12
#include "Node.h"
#include "Triangle.h"
JasonPries's avatar
JasonPries committed
13
#include "Region.h"
14

JasonPries's avatar
JasonPries committed
15
template<size_t Dimension, size_t Order>
16 17
class FiniteElementMesh {
};
18

JasonPries's avatar
JasonPries committed
19
template<size_t Order>
20
class FiniteElementMesh<2, Order> {
21
public:
JasonPries's avatar
JasonPries committed
22 23
    FiniteElementMesh() {};

JasonPries's avatar
JasonPries committed
24 25
    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<std::shared_ptr<Region<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {
        auto comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, std::shared_ptr<DiscreteBoundary<2>> const &y) { return (size_t) (x->curve().get()) < (size_t) (y->curve().get()); };
26 27
        std::sort(Boundaries.begin(), Boundaries.end(), comp);
    };
28 29

    FiniteElementMesh(Mesh &m);
30

JasonPries's avatar
JasonPries committed
31 32
    size_t size_nodes() const { return Nodes.size(); };

JasonPries's avatar
JasonPries committed
33 34
    size_t size_elements() const { return Triangles.size(); };

JasonPries's avatar
JasonPries committed
35
    std::vector<Triangle<Order>> const &triangles() const { return Triangles; };
JasonPries's avatar
JasonPries committed
36

JasonPries's avatar
JasonPries committed
37
    Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
JasonPries's avatar
JasonPries committed
38

39
    auto &boundaries() { return Boundaries; };
JasonPries's avatar
JasonPries committed
40

41
    auto &boundary(size_t i) { return Boundaries[i]; };
JasonPries's avatar
JasonPries committed
42

JasonPries's avatar
JasonPries committed
43 44
    std::shared_ptr<DiscreteBoundary<2>> boundary(std::shared_ptr<Curve const> const &c) const {
        auto comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, size_t const &y) { return (size_t) (x->curve().get()) < y; };
45 46 47 48 49 50

        auto result = std::lower_bound(Boundaries.begin(), Boundaries.end(), (size_t) c.get(), comp);

        if (result != Boundaries.end()) {
            return *result;
        } else {
JasonPries's avatar
JasonPries committed
51
            std::shared_ptr<DiscreteBoundary<2>> null;
52 53 54 55
            return null;
        }
    };

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
    std::shared_ptr<DiscontinuousBoundary<2>> make_discontinuous(std::shared_ptr<Curve const> const &c) {
        // TODO: Should this be implemented in the refineable mesh instead? It might be simpler.
        std::shared_ptr<DiscreteBoundary<2>> b = boundary(c);

        if (!b) {
            std::shared_ptr<DiscontinuousBoundary<2>> null;
            return null;
        }

        std::shared_ptr<DiscontinuousBoundary<2>> db = std::dynamic_pointer_cast<DiscontinuousBoundary<2>>(b);
        if (db) {
            return db;
        }

        // Copy nodes
        std::vector<size_t> bnodes = b->nodes();
        std::vector<size_t> dnodes;
        dnodes.reserve(bnodes.size());
        for (size_t i = 0; i != bnodes.size(); ++i) {
            dnodes.push_back(Nodes.size());
            Nodes.emplace_back(Nodes[bnodes[i]]);
        }

        // Construct new DiscontinuousBoundary and replace
        db = std::make_shared<DiscontinuousBoundary<2>>(bnodes, dnodes, std::vector<double_t>(dnodes.size(), 1.0), b->is_closed());
        std::replace(Boundaries.begin(), Boundaries.end(), b, std::dynamic_pointer_cast<DiscreteBoundary<2>>(db));

        // Renumber nodes in elements when the contains nodes with the same orientation as the boundary
        // TODO: This a brute renumbering algorithm with O(t*n) complexity where t = Triangles.size() and n = bnodes.size()
        // TODO: Sort + renumber should result in O(n*log(t)) complexity, but all present assumptions about this class assume the the ordering of the elements never changes
        for (size_t i = 0; i != (dnodes.size() - 1 + db->is_closed()); ++i) {
            size_t ii = (i + 1) % dnodes.size();
            double_t dx0 = Nodes[bnodes[ii]].x() - Nodes[bnodes[i]].x();
            double_t dy0 = Nodes[bnodes[ii]].y() - Nodes[bnodes[i]].y();

            for (Triangle<Order> &t : Triangles) {
                for (size_t k = 0; k != 3; ++k) {
                    if (t.node(k) == bnodes[i]) {
                        size_t k0 = t.node(0);
                        size_t k1 = t.node(1);
                        size_t k2 = t.node(2);

                        double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[bnodes[i]].x() + Nodes[bnodes[ii]].x()) / 2.0;
                        double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0  - (Nodes[bnodes[i]].y() + Nodes[bnodes[ii]].y()) / 2.0;

                        double_t cross = dx0 * dy1 - dy0 * dx1;
                        if (cross < 0) {
                            t.node(k, dnodes[i]);
                        }
                    }

                    if (t.node(k) == bnodes[ii]) {
                        size_t k0 = t.node(0);
                        size_t k1 = t.node(1);
                        size_t k2 = t.node(2);

                        double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[bnodes[i]].x() + Nodes[bnodes[ii]].x()) / 2.0;
                        double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0  - (Nodes[bnodes[i]].y() + Nodes[bnodes[ii]].y()) / 2.0;

                        double_t cross = dx0 * dy1 - dy0 * dx1;
                        if (cross < 0) {
                            t.node(k, dnodes[ii]);
                        }
                    }
                }
            }
        }

        return db;
    }

127
    auto const &boundaries() const { return Boundaries; };
JasonPries's avatar
JasonPries committed
128

129
    auto const &nodes() const { return Nodes; };
JasonPries's avatar
JasonPries committed
130

131
    auto const &node(size_t i) const { return Nodes[i]; };
JasonPries's avatar
JasonPries committed
132

133
    auto &regions() { return Regions; };
JasonPries's avatar
JasonPries committed
134

135
    auto &region(size_t i) { return Regions[i]; };
JasonPries's avatar
JasonPries committed
136

137 138 139
    auto const &regions() const { return Regions; };

    auto const &region(size_t i) const { return Regions[i]; };
JasonPries's avatar
JasonPries committed
140

141
    template<size_t Q>
142 143
    DiagonalMatrixGroup<TriangleQuadrature<Q>::size> determinant() const {
        DiagonalMatrixGroup<TriangleQuadrature<Q>::size> mat(Triangles.size());
144 145 146 147 148 149 150 151

        for (size_t i = 0; i != Triangles.size(); ++i) {
            Triangles[i].determinant<Q>(mat, Nodes);
        }

        return mat;
    }

JasonPries's avatar
JasonPries committed
152
    template<size_t Q>
153 154
    SparseMatrixGroup<TriangleQuadrature<Q>::size> basis() const {
        SparseMatrixGroup<TriangleQuadrature<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
JasonPries's avatar
JasonPries committed
155 156

        for (size_t i = 0; i != Triangles.size(); ++i) {
157
            Triangles[i].basis<Q>(mat, Nodes);
JasonPries's avatar
JasonPries committed
158 159 160
        }

        return mat;
161
    }
JasonPries's avatar
JasonPries committed
162 163

    template<size_t Q>
164 165
    DerivativeMatrixGroup<TriangleQuadrature<Q>::size> derivative() const {
        DerivativeMatrixGroup<TriangleQuadrature<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
JasonPries's avatar
JasonPries committed
166

167
        for (size_t i = 0; i != Triangles.size(); ++i) {
168
            Triangles[i].derivative<Q>(df, Nodes);
JasonPries's avatar
JasonPries committed
169
        }
170

JasonPries's avatar
JasonPries committed
171 172
        return df;
    }
173

174 175 176 177 178
    template<size_t Q>
    std::vector<std::vector<XY>> quadrature_points() const {
        std::vector<std::vector<XY>> qp;
        qp.reserve(Triangles.size());

179
        for (size_t i = 0; i != Triangles.size(); ++i) {
180 181 182 183 184 185
            qp.push_back(Triangles[i].quadrature_points<Q>(Nodes));
        }

        return qp;
    }

186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217
    void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const {
        if (!boost::filesystem::exists(path)) {
            boost::filesystem::create_directories(path);
        }

        std::fstream fs;
        fs.open(path + file_name + ".oesc", std::fstream::out);

        // Write header
        size_t loc = 3;
        fs << loc << ',' << Triangles.size() << std::endl;
        loc += Triangles.size();

        fs << loc << ',' << Nodes.size() << std::endl;
        loc += Nodes.size();

        fs << loc << ',' << v.size() << std::endl;

        // Write triangles
        for (Triangle<Order> const &t : Triangles) {
            for (size_t i = 0; i != Triangle<Order>::NumNodes - 1; ++i) {
                fs << t.node(i) << ',';
            }
            fs << t.node(Triangle<Order>::NumNodes - 1) << std::endl;
        }

        // Write nodes
        for (XY const &n : Nodes) {
            fs << n.x() << ',' << n.y() << std::endl;
        }

        // Write values
218
        for (size_t i = 0; i != v.size(); ++i) {
219 220 221 222 223 224 225
            fs << v(i) << std::endl;
        }


        fs.close();
    }

226 227
protected:
    std::vector<XY> Nodes;
JasonPries's avatar
JasonPries committed
228 229
    std::vector<Triangle<Order>> Triangles;

JasonPries's avatar
JasonPries committed
230
    std::vector<std::shared_ptr<Region<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
JasonPries's avatar
JasonPries committed
231
    std::vector<std::shared_ptr<DiscreteBoundary<2>>> Boundaries;
232 233 234
};

#endif //OERSTED_FINITEELEMENTMESH_H