FiniteElementMesh.h 14 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"
13
#include "DiscreteRegion.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() {};

24 25 26
    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<std::shared_ptr<DiscreteRegion<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {
        sort_boundaries();
        sort_regions();
27
    };
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

43 44 45
    std::vector<std::shared_ptr<DiscreteBoundary<2>>> boundary(std::shared_ptr<Curve const> const &c) const {
        auto lower_comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, size_t const &y) { return (size_t) (x->curve().get()) < y; };
        auto upper_comp = [](size_t const &y, std::shared_ptr<DiscreteBoundary<2>> const &x) { return y < (size_t) (x->curve().get()); };
46

47 48
        auto lower = std::lower_bound(Boundaries.begin(), Boundaries.end(), (size_t) c.get(), lower_comp);
        auto upper = std::upper_bound(Boundaries.begin(), Boundaries.end(), (size_t) c.get(), upper_comp);
49

50 51
        if (lower != Boundaries.end()) {
            return std::vector<std::shared_ptr<DiscreteBoundary<2>>>(lower, upper);
52
        } else {
53
            return std::vector<std::shared_ptr<DiscreteBoundary<2>>>();
54 55 56
        }
    };

57
    std::vector<std::shared_ptr<DiscreteBoundary<2>>> make_discontinuous(std::shared_ptr<Curve const> const &c) {
58
        // TODO: Should this be implemented in the refineable mesh instead? It might be simpler.
59
        std::vector<std::shared_ptr<DiscreteBoundary<2>>> b_vec = boundary(c);
60

61 62 63 64 65
        if (b_vec.size() == 0 || b_vec.size() == 2) { // Boundary curve not found || boundary pair found
            return b_vec;
        } else if (b_vec.size() > 2) { // should not happen
            throw std::runtime_error("FiniteElementMesh<2>::make_discontinuous expects to find at most 2 DiscreteBoundary<2> objects");
        } // else b.size() == 1 and we need to construct an new boundary
66 67

        // Copy nodes
68 69 70 71 72 73
        std::vector<size_t> boundary_nodes = b_vec[0]->nodes();
        std::vector<size_t> discontinuous_nodes;
        discontinuous_nodes.reserve(boundary_nodes.size());
        for (size_t i = 0; i != boundary_nodes.size(); ++i) {
            discontinuous_nodes.push_back(Nodes.size());
            Nodes.emplace_back(Nodes[boundary_nodes[i]]);
74 75
        }

76 77 78 79 80 81 82 83
        // Construct new DiscontinuousBoundary
        bool is_closed = b_vec[0]->is_closed();
        std::shared_ptr<DiscreteBoundary<2>> db = std::make_shared<DiscreteBoundary<2>>(b_vec[0]->curve(), discontinuous_nodes, is_closed);
        Boundaries.push_back(db);
        sort_boundaries();

        // Renumber nodes in elements when the element is on the "right" side of the plane defined by the discontinuous boundary edge
        // (cross product of edge vector with edge midpoint to incenter vector is negative)
84 85 86

        // 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
87 88 89 90 91
        // TODO: However, triangles do store an ID which could reverse the sorting
        for (size_t i = 0; i != (discontinuous_nodes.size() - 1 + is_closed); ++i) {
            size_t ii = (i + 1) % discontinuous_nodes.size();
            double_t dx0 = Nodes[boundary_nodes[ii]].x() - Nodes[boundary_nodes[i]].x();
            double_t dy0 = Nodes[boundary_nodes[ii]].y() - Nodes[boundary_nodes[i]].y();
92 93 94

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

100 101
                        double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[boundary_nodes[i]].x() + Nodes[boundary_nodes[ii]].x()) / 2.0;
                        double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0  - (Nodes[boundary_nodes[i]].y() + Nodes[boundary_nodes[ii]].y()) / 2.0;
102 103 104

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

109
                    if (t.node(k) == boundary_nodes[ii]) {
110 111 112 113
                        size_t k0 = t.node(0);
                        size_t k1 = t.node(1);
                        size_t k2 = t.node(2);

114 115
                        double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[boundary_nodes[i]].x() + Nodes[boundary_nodes[ii]].x()) / 2.0;
                        double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0  - (Nodes[boundary_nodes[i]].y() + Nodes[boundary_nodes[ii]].y()) / 2.0;
116 117 118

                        double_t cross = dx0 * dy1 - dy0 * dx1;
                        if (cross < 0) {
119
                            t.node(k, discontinuous_nodes[ii]);
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 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
        // Renumber other boundary nodes when the other boundary is on the right-hand side of the plane defined by the discontinuous boundary edges
        b_vec.push_back(db);
        for (auto &b : Boundaries) {
            if (b == b_vec[0] || b == b_vec[1]) {
                continue;
            }

            size_t k0, k1, k2;

            if (b->nodes().front() == boundary_nodes.front()) {
                k0 = b->node(0);
                k1 = b->node(1);
                k2 = boundary_nodes[1];

                double_t dx10 = Nodes[k1].x() - Nodes[k0].x();
                double_t dy10 = Nodes[k1].y() - Nodes[k0].y();

                double_t dx20 = Nodes[k2].x() - Nodes[k0].x();
                double_t dy20 = Nodes[k2].y() - Nodes[k0].y();

                double_t cross = dx10 * dy20 - dy10 * dx20;
                if (cross > 0) {
                    b->node(0, discontinuous_nodes.front());
                }
            }

            if (b->nodes().front() == boundary_nodes.back()) {
                k0 = b->node(0);
                k1 = b->node(1);
                k2 = boundary_nodes[boundary_nodes.size() - 2];

                double_t dx10 = Nodes[k1].x() - Nodes[k0].x();
                double_t dy10 = Nodes[k1].y() - Nodes[k0].y();

                double_t dx20 = Nodes[k2].x() - Nodes[k0].x();
                double_t dy20 = Nodes[k2].y() - Nodes[k0].y();

                double_t cross = dx10 * dy20 - dy10 * dx20;
                if (cross < 0) {
                    b->node(0, discontinuous_nodes.back());
                }
            }

            if (b->nodes().back() == boundary_nodes.front()) {
                k0 = b->node(b->nodes().size() - 1);
                k1 = b->node(b->nodes().size() - 2);
                k2 = boundary_nodes[1];

                double_t dx10 = Nodes[k1].x() - Nodes[k0].x();
                double_t dy10 = Nodes[k1].y() - Nodes[k0].y();

                double_t dx20 = Nodes[k2].x() - Nodes[k0].x();
                double_t dy20 = Nodes[k2].y() - Nodes[k0].y();

                double_t cross = dx10 * dy20 - dy10 * dx20;
                if (cross > 0) {
                    b->node(b->nodes().size() - 1, discontinuous_nodes.front());
                }
            }

            if (b->nodes().back() == boundary_nodes.back()) {
                k0 = b->node(b->nodes().size() - 1);
                k1 = b->node(b->nodes().size() - 2);
                k2 = boundary_nodes[boundary_nodes.size() - 2];

                double_t dx10 = Nodes[k1].x() - Nodes[k0].x();
                double_t dy10 = Nodes[k1].y() - Nodes[k0].y();

                double_t dx20 = Nodes[k2].x() - Nodes[k0].x();
                double_t dy20 = Nodes[k2].y() - Nodes[k0].y();

                double_t cross = dx10 * dy20 - dy10 * dx20;
                if (cross < 0) {
                    b->node(b->nodes().size() - 1, discontinuous_nodes.back());
                }
            }
        }

        return b_vec;
205 206
    }

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

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

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

213
    auto &regions() { return Regions; };
JasonPries's avatar
JasonPries committed
214

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

217 218 219
    auto const &regions() const { return Regions; };

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

221 222 223 224 225 226 227 228 229 230 231 232
    std::shared_ptr<DiscreteRegion<2>> region(std::shared_ptr<Contour const> const &c) const {
        auto comp = [](std::shared_ptr<DiscreteRegion<2>> const &x, size_t const &y) { return (size_t) (x->region().get()) < y; };

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

        if (iter != Regions.end()) {
            return *iter;
        } else {
            return std::make_shared<DiscreteRegion<2>>();
        }
    }

233
    template<size_t Q>
234 235
    DiagonalMatrixGroup<TriangleQuadrature<Q>::size> determinant() const {
        DiagonalMatrixGroup<TriangleQuadrature<Q>::size> mat(Triangles.size());
236 237 238 239 240 241 242 243

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

        return mat;
    }

JasonPries's avatar
JasonPries committed
244
    template<size_t Q>
245 246
    SparseMatrixGroup<TriangleQuadrature<Q>::size> basis() const {
        SparseMatrixGroup<TriangleQuadrature<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
JasonPries's avatar
JasonPries committed
247 248

        for (size_t i = 0; i != Triangles.size(); ++i) {
249
            Triangles[i].basis<Q>(mat, Nodes);
JasonPries's avatar
JasonPries committed
250 251 252
        }

        return mat;
253
    }
JasonPries's avatar
JasonPries committed
254 255

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

259
        for (size_t i = 0; i != Triangles.size(); ++i) {
260
            Triangles[i].derivative<Q>(df, Nodes);
JasonPries's avatar
JasonPries committed
261
        }
262

JasonPries's avatar
JasonPries committed
263 264
        return df;
    }
265

266 267 268 269 270
    template<size_t Q>
    std::vector<std::vector<XY>> quadrature_points() const {
        std::vector<std::vector<XY>> qp;
        qp.reserve(Triangles.size());

271
        for (size_t i = 0; i != Triangles.size(); ++i) {
272 273 274 275 276 277
            qp.push_back(Triangles[i].quadrature_points<Q>(Nodes));
        }

        return qp;
    }

278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309
    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
310
        for (size_t i = 0; i != v.size(); ++i) {
311 312 313 314 315 316 317
            fs << v(i) << std::endl;
        }


        fs.close();
    }

318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
    void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, 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 + ".oeve", 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 << ',' << Fx.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
        for (size_t i = 0; i != Fx.size(); ++i) {
            fs << Fx(i) << ',' << Fy(i) <<  std::endl;
        }

        fs.close();
    }

357 358
protected:
    std::vector<XY> Nodes;
JasonPries's avatar
JasonPries committed
359 360
    std::vector<Triangle<Order>> Triangles;

361
    std::vector<std::shared_ptr<DiscreteRegion<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
JasonPries's avatar
JasonPries committed
362
    std::vector<std::shared_ptr<DiscreteBoundary<2>>> Boundaries;
363 364 365 366 367 368 369 370 371 372

    void sort_boundaries() {
        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()); };
        std::sort(Boundaries.begin(), Boundaries.end(), comp);
    }

    void sort_regions() {
        auto comp = [](std::shared_ptr<DiscreteRegion<2>> const &x, std::shared_ptr<DiscreteRegion<2>> const &y) { return (size_t) (x->region().get()) < (size_t) (y->region().get()); };
        std::sort(Regions.begin(), Regions.end(), comp);
    }
373 374 375
};

#endif //OERSTED_FINITEELEMENTMESH_H