PLaSK library

Loading...
Searching...
No Matches
generator_triangular.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include "generator_triangular.hpp"
15#include "triangular2d.hpp"
16
17#include <triangle.h>
18 using namespace triangle;
19
20 namespace plask {
21
22 struct TrifreeCaller {
23 void operator()(void* ptr) const { trifree(ptr); }
24};
25
26 struct VecFuzzyCompare {
27 bool operator()(const typename Primitive<2>::DVec& a, const typename Primitive<2>::DVec& b) const {
28 return Primitive<2>::vecFuzzyCompare(a, b);
29 }
30};
31
32 shared_ptr<MeshD<2>> TriangleGenerator::generate(const shared_ptr<GeometryObjectD<2>>& geometry) {
33 typedef typename GeometryObjectD<2>::LineSegment LineSegment;
34 typedef typename Primitive<2>::DVec DVec;
35 typedef std::map<typename Primitive<2>::DVec, int, VecFuzzyCompare> PointMap;
36
37 std::set<LineSegment> lineSegments = geometry->getLineSegments();
38
39 // Make sets of points and segments removing duplicates
40 Box2D bb = geometry->getBoundingBox();
41
42 PointMap pointmap;
43 std::set<std::pair<int, int>> segmentsset;
44 if (full) {
45 pointmap[DVec(bb.left(), bb.bottom())] = 0;
46 pointmap[DVec(bb.right(), bb.bottom())] = 1;
47 pointmap[DVec(bb.right(), bb.top())] = 2;
48 pointmap[DVec(bb.left(), bb.top())] = 3;
49
50 if (pointmap.size() < 4)
51 throw BadMesh("TriangularGenerator", "geometry object is too small to create triangular mesh");
52
53 segmentsset.insert(std::make_pair(0, 1));
54 segmentsset.insert(std::make_pair(1, 2));
55 segmentsset.insert(std::make_pair(2, 3));
56 segmentsset.insert(std::make_pair(0, 3));
57 }
58
59 size_t n = pointmap.size();
60 for (const LineSegment& segment : lineSegments) {
61 int iseg[2];
62 for (int i = 0; i < 2; ++i) {
63 bool ok;
64 PointMap::iterator it;
65 PLASK_NO_CONVERSION_WARNING_BEGIN
66 std::tie(it, ok) = pointmap.insert(std::make_pair(segment[i], n));
67 PLASK_NO_WARNING_END
68 if (ok) ++n;
69 iseg[i] = it->second;
70 }
71 if (iseg[0] != iseg[1]) {
72 if (iseg[0] > iseg[1]) std::swap(iseg[0], iseg[1]);
73 segmentsset.insert(std::make_pair(iseg[0], iseg[1]));
74 }
75 }
76
77 triangulateio in = {}, out = {}; // are fields are nulled, so we will only fill fields we need
78
79 in.numberofpoints = pointmap.size();
80 std::unique_ptr<double[]> in_points(new double[2 * in.numberofpoints]);
81 in.pointlist = in_points.get();
82 for (auto pi : pointmap) {
83 in_points[2 * pi.second] = pi.first.c0;
84 in_points[2 * pi.second + 1] = pi.first.c1;
85 // std::cerr << format("{}: {:5.3f}, {:5.3f}\n", pi.second, pi.first.c0, pi.first.c1);
86 }
87
88 in.numberofsegments = segmentsset.size();
89 std::unique_ptr<int[]> in_segments(new int[2 * in.numberofsegments]);
90 in.segmentlist = in_segments.get();
91 n = 0;
92 for (auto s : segmentsset) {
93 in_segments[n] = s.first;
94 in_segments[n + 1] = s.second;
95 n += 2;
96 // std::cerr << s.first << "-" << s.second << "\n";
97 }
98
99 triangulate(const_cast<char*>(getSwitches().c_str()), &in, &out, nullptr);
100
101 // just for case we free memory which could be allocated by triangulate but we do not need (some of this can be
102 // nullptr):
103 trifree(out.pointattributelist);
104 trifree(out.pointmarkerlist);
105 trifree(out.triangleattributelist);
106 trifree(out.trianglearealist);
107 trifree(out.neighborlist);
108 trifree(out.segmentlist);
109 trifree(out.segmentmarkerlist);
110 trifree(out.holelist);
111 trifree(out.regionlist);
112 trifree(out.edgelist);
113 trifree(out.edgemarkerlist);
114 trifree(out.normlist);
115
116 // this will free rest of memory allocated by triangulate (even if an exception will be thrown):
117 std::unique_ptr<REAL[], TrifreeCaller> out_points(out.pointlist);
118 std::unique_ptr<int[], TrifreeCaller> out_triangles(out.trianglelist);
119
120 shared_ptr<TriangularMesh2D> result = make_shared<TriangularMesh2D>();
121 result->nodes.reserve(out.numberofpoints);
122 for (std::size_t i = 0; i < std::size_t(out.numberofpoints) * 2; i += 2)
123 result->nodes.emplace_back(out.pointlist[i], out.pointlist[i + 1]);
124 result->elementNodes.reserve(out.numberoftriangles);
125 for (std::size_t i = 0; i < std::size_t(out.numberoftriangles) * 3; i += 3)
126 result->elementNodes.push_back({std::size_t(out.trianglelist[i]), std::size_t(out.trianglelist[i + 1]),
127 std::size_t(out.trianglelist[i + 2])});
128 // it is also possible by triangle, to use 6 numbers per triangle, but we do not support such things at the moment
129 return result;
130}
131
132std::string TriangleGenerator::getSwitches() const {
133 std::ostringstream result;
134
135 // p - specifies vertices, segments, holes, regional attributes, and regional area constraints
136 // z - points (and other items) are numbered from zero
137 // Q - quiet
138 // B - suppresses boundary markers in the output
139 // P - suppresses the output .poly file
140 result << "pzQqBP";
141
142 // D - Conforming Delaunay triangulation
143 if (delaunay) result << 'D';
144
145 // a - imposes a maximum triangle area
146 if (maxTriangleArea) result << 'a' << std::fixed << *maxTriangleArea;
147
148 // q - adds vertices to the mesh to ensure that all angles are between given and 140 degrees.
149 if (minTriangleAngle) {
150 result << 'q';
151 if (!isnan(*minTriangleAngle)) result << std::fixed << *minTriangleAngle;
152 }
153
154 // TODO more configuration
155
156 return result.str();
157}
158
159shared_ptr<MeshGenerator> readTriangleGenerator(XMLReader& reader, const Manager&) {
160 shared_ptr<TriangleGenerator> result = make_shared<TriangleGenerator>();
161 if (reader.requireTagOrEnd("options")) {
162 result->maxTriangleArea = reader.getAttribute<double>("maxarea");
163 result->minTriangleAngle = reader.getAttribute<double>("minangle");
164 result->delaunay = reader.getAttribute<bool>("delaunay", false);
165 result->full = reader.getAttribute<bool>("full", false);
166 reader.requireTagEnd(); // end of options
167 reader.requireTagEnd(); // end of generator
168 }
169 return result;
170}
171
172 static RegisterMeshGeneratorReader trianglegenerator_reader("triangular2d.triangle", readTriangleGenerator);
173
174} // namespace plask