FiniteElementMesh.h 18 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

15
class FiniteElementMeshInterface {
16
public:
17
    FiniteElementMeshInterface() {};
JasonPries's avatar
JasonPries committed
18

19
    FiniteElementMeshInterface(std::vector<XY> nodes, std::vector<std::shared_ptr<DiscreteRegion<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : Nodes(nodes), Regions(r), Boundaries(b) {
20 21
        sort_boundaries();
        sort_regions();
22
    };
23

24
    virtual ~FiniteElementMeshInterface() {};
25

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

28
    auto &boundaries() { return Boundaries; };
JasonPries's avatar
JasonPries committed
29

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

32
    std::vector<std::shared_ptr<DiscreteBoundary<2>>> boundary(std::shared_ptr<Curve const> const &c) const {
33 34
        // Boundary may be discontinuous (e.g. multiple overlapping curves). Therefore, multiple DiscreteBoundaries may be returned (in general, upper != lower)

35 36
        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()); };
37

38 39
        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);
40

41 42
        if (lower != Boundaries.end()) {
            return std::vector<std::shared_ptr<DiscreteBoundary<2>>>(lower, upper);
43
        } else {
44
            return std::vector<std::shared_ptr<DiscreteBoundary<2>>>();
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 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 205 206 207 208 209 210 211 212 213 214 215
    // TODO: virtual std::shared_ptr<DiscreteRegion<2>> region(std::shared_ptr<Region const> const &r) const = 0;

    virtual std::vector<std::shared_ptr<DiscreteBoundary<2>>> make_discontinuous(std::shared_ptr<Curve const> const &c) = 0;

    auto const &boundaries() const { return Boundaries; };

    auto const &nodes() const { return Nodes; };

    auto const &node(size_t i) const { return Nodes[i]; };

    auto &regions() { return Regions; };

    auto &region(size_t i) { return Regions[i]; };

    auto const &regions() const { return Regions; };

    auto const &region(size_t i) const { return Regions[i]; };

    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>>();
        }
    }

    size_t nodes_size() const { return Nodes.size(); };

    virtual size_t elements_size() const = 0;

    virtual void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const = 0;

    virtual void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, std::string path, std::string file_name) const = 0;

    virtual DiagonalMatrixGroup determinant() const = 0;

    virtual SparseMatrixGroup basis() const = 0;

    virtual DerivativeMatrixGroup derivative() const = 0;

    virtual std::vector<std::vector<XY>> quadrature_points() const = 0;

protected:
    std::vector<XY> Nodes;

    std::vector<std::shared_ptr<DiscreteRegion<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
    std::vector<std::shared_ptr<DiscreteBoundary<2>>> Boundaries;

    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);
    }
};

template<size_t ElementOrder, size_t QuadratureOrder>
class FiniteElementMesh : public FiniteElementMeshInterface {
public:
    FiniteElementMesh() {};

    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<ElementOrder>> tris, std::vector<std::shared_ptr<DiscreteRegion<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : FiniteElementMeshInterface{nodes,r,b}, Triangles(tris) {};

    FiniteElementMesh(Mesh &m) {
        // std::vector<XY> nodes
        // std::vector<Triangle<Order>> tris
        // std::vector<std::shared_ptr<Region<2>>> r
        // std::vector<std::shared_ptr<Boundary<2>>> b

        Nodes.reserve(m.size_points()); // TODO: This will be larger when Order > 1
        for (size_t i = 0; i!= m.size_points(); ++i) {
            Nodes.emplace_back(m.point(i));
        }

        /* TODO: Changes to rotational motion implementation
            Initially, perform a straightforward conversion as currently implemented
            Define a new Boundary type called (called something like SlidingBoundary)
                SlidingBoundary holds a reference to the two boundaries and ensures the correct mapping
            In FiniteElementMesh, add a method called make_sliding(...) or something similar
                This method will implement a procedure which duplicates the nodes on the boundary and renumbers the nodes in the appropriate elements
            Then, the SlidingInterface boundary condition holds a pointer to the SlidingBoundary instead of two Boundary objects

            Perhaps a similar approach can be taken for implementing periodic boundary conditions (e.g. MappedBoundary holds two (or several arrays of) boundaries and enforces the proper invariant)
                e.g. make_mapped(...)
            Then PeriodicBoundaryCondition holds a pointer to the MappedBoundary instead of two Boundary objects

            TODO: To make these easier to implement from the user side, Boundary should hold a shared_ptr to the defining curve (from the Mesh object)
                Then selection can be implemented by passing the shared_ptr that the user has (similar to how the Mesh operations work)
         */

        Boundaries.reserve(m.size_boundary_constraints());
        for (size_t i = 0; i!= m.size_boundary_constraints(); ++i) {
            Boundaries.push_back(std::make_shared<DiscreteBoundary<2>>(m.curve(i), m.boundary_nodes(i)));
        }

        Regions.reserve(m.size_contours());
        for (size_t i = 0; i != m.size_contours(); ++i) {
            Regions.push_back(std::make_shared<DiscreteRegion<2>>(m.contour(i))); // TODO: Assign material properties here
        }

        Triangles.reserve(m.size_triangles());
        for (size_t i = 0; i != m.size_triangles(); ++i) {
            auto &dt = m.triangle(i);
            Regions[dt.contour()]->push_back(i);
            Triangles.emplace_back(i, m.nodes_from_triangle(dt));
        }

        sort_boundaries();
        sort_regions();
    };

    std::vector<Triangle<ElementOrder>> const &triangles() const { return Triangles; };

    Triangle<ElementOrder> const &triangle(size_t i) const { return Triangles[i]; };

    size_t elements_size() const override { return Triangles.size(); };

    DiagonalMatrixGroup determinant() const override {
        DiagonalMatrixGroup mat(TriangleQuadrature<QuadratureOrder>::size, Triangles.size());

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

        return mat;
    }

    SparseMatrixGroup basis() const override {
        SparseMatrixGroup mat(TriangleQuadrature<QuadratureOrder>::size, Nodes.size(), Triangles.size(), Triangle<ElementOrder>::NumNodes);

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

        return mat;
    }

    DerivativeMatrixGroup derivative() const override {
        DerivativeMatrixGroup df(TriangleQuadrature<QuadratureOrder>::size, Nodes.size(), Triangles.size(), Triangle<ElementOrder>::NumNodes);

        for (size_t i = 0; i != Triangles.size(); ++i) {
            Triangles[i]. template derivative<QuadratureOrder>(df, Nodes);
        }

        return df;
    }

    std::vector<std::vector<XY>> quadrature_points() const override {
        std::vector<std::vector<XY>> qp;
        qp.reserve(Triangles.size());

        for (size_t i = 0; i != Triangles.size(); ++i) {
            qp.push_back(Triangles[i]. template quadrature_points<QuadratureOrder>(Nodes));
        }

        return qp;
    }

    // TODO: std::shared_ptr<DiscreteRegion<2>> region(std::shared_ptr<Region const> const &r) const override {...}

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

219 220 221 222 223
        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
224 225

        // Copy nodes
226 227 228 229 230 231
        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]]);
232 233
        }

234 235 236 237 238 239 240 241
        // 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)
242 243 244

        // 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
245 246 247 248 249
        // 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();
250

251
            for (Triangle<ElementOrder> &t : Triangles) {
252
                for (size_t k = 0; k != 3; ++k) {
253
                    if (t.node(k) == boundary_nodes[i]) {
254 255 256 257
                        size_t k0 = t.node(0);
                        size_t k1 = t.node(1);
                        size_t k2 = t.node(2);

258 259
                        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;
260 261 262

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

267
                    if (t.node(k) == boundary_nodes[ii]) {
268 269 270 271
                        size_t k0 = t.node(0);
                        size_t k1 = t.node(1);
                        size_t k2 = t.node(2);

272 273
                        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;
274 275 276

                        double_t cross = dx0 * dy1 - dy0 * dx1;
                        if (cross < 0) {
277
                            t.node(k, discontinuous_nodes[ii]);
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 310 311 312 313 314 315 316 317 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 357 358 359 360 361 362
        // 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;
363 364
    }

365
    void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const override {
366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
        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
384 385
        for (Triangle<ElementOrder> const &t : Triangles) {
            for (size_t i = 0; i != Triangle<ElementOrder>::NumNodes - 1; ++i) {
386 387
                fs << t.node(i) << ',';
            }
388
            fs << t.node(Triangle<ElementOrder>::NumNodes - 1) << std::endl;
389 390 391 392 393 394 395 396
        }

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

        // Write values
397
        for (size_t i = 0; i != v.size(); ++i) {
398 399 400 401 402 403 404
            fs << v(i) << std::endl;
        }


        fs.close();
    }

405
    void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, std::string path, std::string file_name) const override {
406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
        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
424 425
        for (Triangle<ElementOrder> const &t : Triangles) {
            for (size_t i = 0; i != Triangle<ElementOrder>::NumNodes - 1; ++i) {
426 427
                fs << t.node(i) << ',';
            }
428
            fs << t.node(Triangle<ElementOrder>::NumNodes - 1) << std::endl;
429 430 431 432 433 434 435 436 437 438 439 440 441 442 443
        }

        // 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();
    }

444
protected:
445
    std::vector<Triangle<ElementOrder>> Triangles;
446 447 448
};

#endif //OERSTED_FINITEELEMENTMESH_H