OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
MeshGenerator.cpp
Go to the documentation of this file.
1#include "Structure/MeshGenerator.h"
2#include "Physics/Physics.h"
3#include "Utilities/Util.h"
4#include "AbsBeamline/Bend2D.h"
5#include "AbsBeamline/RBend3D.h"
6#include "AbsBeamline/Multipole.h"
7#include "AbstractObjects/OpalData.h"
8
9#include <iostream>
10#include <fstream>
11
12extern Inform *gmsg;
13
14MeshGenerator::MeshGenerator():
15 elements_m()
16{ }
17
18void MeshGenerator::add(const ElementBase &element) {
19 double start = 0.0;
20
21 MeshData mesh;
22 if (element.getType() == ElementType::SBEND ||
23 element.getType() == ElementType::RBEND) {
24
25 const Bend2D* dipole = static_cast<const Bend2D*>(&element);
26 mesh = dipole->getSurfaceMesh();
27 mesh.type_m = DIPOLE;
28 } else if (element.getType() == ElementType::RBEND3D) {
29 const RBend3D* dipole = static_cast<const RBend3D*>(&element);
30 mesh = dipole->getSurfaceMesh();
31 mesh.type_m = DIPOLE;
32 } else if (element.getType() == ElementType::DRIFT) {
33 return;
34 } else {
35 double end, length;
36 element.getElementDimensions(start, end);
37 length = end - start;
38 auto apert = element.getAperture();
39
40 switch (apert.first) {
41 case ApertureType::RECTANGULAR:
42 case ApertureType::CONIC_RECTANGULAR:
43 mesh = getBox(length, apert.second[0], apert.second[1], apert.second[2]);
44 break;
45 case ApertureType::ELLIPTICAL:
46 case ApertureType::CONIC_ELLIPTICAL:
47 mesh = getCylinder(length, apert.second[0], apert.second[1], apert.second[2]);
48 break;
49 default:
50 return;
51 }
52
53 switch(element.getType()) {
54 case ElementType::MULTIPOLE:
55 switch(static_cast<const Multipole*>(&element)->getMaxNormalComponentIndex()) {
56 case 1:
57 mesh.type_m = DIPOLE;
58 break;
59 case 2:
60 mesh.type_m = QUADRUPOLE;
61 break;
62 case 3:
63 mesh.type_m = SEXTUPOLE;
64 break;
65 case 4:
66 mesh.type_m = OCTUPOLE;
67 break;
68 default:
69 break;
70 }
71 break;
72 case ElementType::RFCAVITY:
73 case ElementType::VARIABLERFCAVITY:
74 case ElementType::TRAVELINGWAVE:
75 mesh.type_m = RFCAVITY;
76 break;
77 case ElementType::SOLENOID:
78 mesh.type_m = SOLENOID;
79 break;
80 default:
81 mesh.type_m = OTHER;
82 }
83 }
84
85 CoordinateSystemTrafo trafo = element.getCSTrafoGlobal2Local().inverted();
86 Vector_t z = trafo.rotateTo(Vector_t(0, 0, 1));
87 for (unsigned int i = 0; i < mesh.vertices_m.size(); ++ i) {
88 mesh.vertices_m[i] = trafo.transformTo(mesh.vertices_m[i]) + start * z;
89 }
90
91 for (unsigned int i = 0; i < mesh.decorations_m.size(); ++ i) {
92 mesh.decorations_m[i].first = trafo.transformTo(mesh.decorations_m[i].first) + start * z;
93 mesh.decorations_m[i].second = trafo.transformTo(mesh.decorations_m[i].second) + start * z;
94 }
95
96 elements_m.push_back(mesh);
97}
98
99#include <boost/iostreams/filtering_streambuf.hpp>
100#include <boost/iostreams/copy.hpp>
101#include <boost/iostreams/filter/zlib.hpp>
102
103void MeshGenerator::write(const std::string &fname) {
104 std::string filename = Util::combineFilePath({
105 OpalData::getInstance()->getAuxiliaryOutputDirectory(),
106 fname + "_ElementPositions.py"
107 });
108 std::ofstream out(filename);
109 const char *buffer;
110 const std::string indent(" ");
111
112 out << std::fixed << std::setprecision(6);
113
114 out << "import os, sys, argparse, math, base64, zlib, struct\n\n";
115 out << "if sys.version_info < (3,0):\n";
116 out << " range = xrange\n\n";
117
118 std::stringstream vertices_ascii;
119 std::ostringstream vertices_compressed;
120 out << "vertices = []\n";
121 out << "numVertices = [";
122 for (auto &element: elements_m) {
123 auto &vertices = element.vertices_m;
124 out << vertices.size() << ", ";
125 for (unsigned int i = 0; i < vertices.size(); ++ i) {
126 buffer = reinterpret_cast<const char*>(&(vertices[i](0)));
127 vertices_ascii.write(buffer, 3 * sizeof(double));
128 }
129 }
130 out.seekp(-2, std::ios_base::end);
131 out << "]\n";
132
133 {
134 boost::iostreams::filtering_streambuf<boost::iostreams::input> in;
135 in.push(boost::iostreams::zlib_compressor());
136 in.push(vertices_ascii);
137 boost::iostreams::copy(in, vertices_compressed);
138 }
139
140 std::string vertices_base64 = Util::base64_encode(vertices_compressed.str());
141 out << "vertices_base64 = \"" << vertices_base64 << "\"\n";
142 out << "triangles = [";
143 for (auto &element: elements_m) {
144 out << "[ ";
145 auto &triangles = element.triangles_m;
146 for (unsigned int i = 0; i < triangles.size(); ++ i) {
147 out << triangles[i](0) << ", " << triangles[i](1) << ", " << triangles[i](2) << ", ";
148 }
149 out.seekp(-2, std::ios_base::end);
150 out << "], ";
151 }
152 out.seekp(-2, std::ios_base::end);
153 out << "]\n";
154
155
156 out << "decoration = [";
157 for (auto &element: elements_m) {
158 out << "[ ";
159 for (auto &decoration: element.decorations_m) {
160 out << (decoration.first)(0) << ", " << (decoration.first)(1) << ", " << (decoration.first)(2) << ", "
161 << (decoration.second)(0) << ", " << (decoration.second)(1) << ", " << (decoration.second)(2) << ", ";
162 }
163 if (!element.decorations_m.empty())
164 out.seekp(-2, std::ios_base::end);
165 out << "], ";
166 }
167 if (!elements_m.empty())
168 out.seekp(-2, std::ios_base::end);
169 out << "]\n\n";
170
171 out << "color = [";
172 for (auto &element: elements_m) {
173 out << element.type_m << ", ";
174 }
175 out.seekp(-2, std::ios_base::end);
176 out << "]\n\n";
177
178 std::stringstream index_ascii;
179 std::ostringstream index_compressed;
180 index_ascii << "<!DOCTYPE html>\n"
181 << "<html>\n"
182 << " <head>\n"
183 << " <meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n"
184 << "\n"
185 << " <title>Babylon.js sample code</title>\n"
186 << " <!-- Babylon.js -->\n"
187 << " <script src=\"http://cdn.babylonjs.com/hand.minified-1.2.js\"></script>\n"
188 << " <script src=\"http://cdn.babylonjs.com/cannon.js\"></script>\n"
189 << " <script src=\"http://cdn.babylonjs.com/oimo.js\"></script>\n"
190 << " <script src=\"http://cdn.babylonjs.com/babylon.js\"></script>\n"
191 << " <style>\n"
192 << " html, body {\n"
193 << " overflow: hidden;\n"
194 << " width: 100%;\n"
195 << " height: 100%;\n"
196 << " margin: 0;\n"
197 << " padding: 0;\n"
198 << " }\n"
199 << "\n"
200 << " #renderCanvas {\n"
201 << " width: 100%;\n"
202 << " height: 100%;\n"
203 << " touch-action: none;\n"
204 << " }\n"
205 << " </style>\n"
206 << " </head>\n"
207 << "<body>\n"
208 << " <canvas id=\"renderCanvas\"></canvas>\n"
209 << " <script>\n"
210 << " var canvas = document.getElementById(\"renderCanvas\");\n"
211 << " var engine = new BABYLON.Engine(canvas, true);\n"
212 << "\n"
213 << " var createScene = function () {\n"
214 << " var scene = new BABYLON.Scene(engine);\n"
215 << "\n"
216 << " //Adding a light\n"
217 << " var light = new BABYLON.PointLight(\"Omni\", new BABYLON.Vector3(0, -100, 0), scene);\n"
218 << " //light.diffuse = new BABYLON.Color3(0.8, 0.8, 0.8);\n"
219 << " light.specular = new BABYLON.Color3(0.1, 0.1, 0.1);\n"
220 << " //light.groundColor = new BABYLON.Color3(0, 0, 0);\n"
221 << "\n"
222 << " //Adding an Arc Rotate Camera\n"
223 << " var camera = new BABYLON.ArcRotateCamera(\"Camera\", 0.0, Math.PI, 50, BABYLON.Vector3.Zero(), scene);\n"
224 << " camera.attachControl(canvas, false);\n"
225 << " camera.wheelPrecision = 20;\n"
226 << "\n"
227 << " var mymesh = ##DATA##;\n"
228 << " // The first parameter can be used to specify which mesh to import. Here we import all meshes\n"
229 << " BABYLON.SceneLoader.ImportMesh(\"\", \"\", mymesh, scene, function (newMeshes) {\n"
230 << " // Set the target of the camera to the first imported mesh\n"
231 << " camera.target = newMeshes[0];\n"
232 << "\n"
233 << " var material1 = new BABYLON.StandardMaterial(\"texture1\", scene);\n"
234 << " //material1.wireframe = true;\n"
235 << " newMeshes[0].material = material1;\n"
236 << " newMeshes[0].material.emissiveColor = new BABYLON.Color3(0.2, 0.2, 0.2);\n"
237 << " });\n"
238 << "\n"
239 << " return scene;\n"
240 << " }\n"
241 << "\n"
242 << "\n"
243 << " var scene = createScene();\n"
244 << "\n"
245 << " engine.runRenderLoop(function () {\n"
246 << " scene.render();\n"
247 << " });\n"
248 << "\n"
249 << " // Resize\n"
250 << " window.addEventListener(\"resize\", function () {\n"
251 << " engine.resize();\n"
252 << " });\n"
253 << " </script>\n"
254 << "</body>\n"
255 << "</html>";
256 {
257 boost::iostreams::filtering_streambuf<boost::iostreams::input> in;
258 in.push(boost::iostreams::zlib_compressor());
259 in.push(index_ascii);
260 boost::iostreams::copy(in, index_compressed);
261 }
262
263 out << "index_base64 = '" << Util::base64_encode(index_compressed.str()) << "'\n\n";
264
265 out << "def decodeVertices():\n";
266 out << indent << "vertices_binary = zlib.decompress(base64.b64decode(vertices_base64))\n";
267 out << indent << "k = 0\n";
268 out << indent << "for i in range(len(numVertices)):\n";
269 out << indent << indent << "current = []\n";
270 out << indent << indent << "for j in range(numVertices[i] * 3):\n";
271 out << indent << indent << indent << "current.append(float(struct.unpack('=d', vertices_binary[k:k+8])[0]))\n";
272 out << indent << indent << indent << "k += 8\n";
273 out << indent << indent << "vertices.append(current)\n\n";
274
275 out << "def dot(a, b):\n";
276 out << indent << "return sum(a[i]*b[i] for i in range(len(a)))\n\n";
277
278 out << "def normalize(a):\n";
279 out << indent << "length = math.sqrt(dot(a, a))\n";
280 out << indent << "for i in range(3):\n";
281 out << indent << indent << "a[i]/=length\n";
282 out << indent << "return a\n\n";
283
284 out << "def distance(a, b):\n";
285 out << indent << "c = [0] * 2\n";
286 out << indent << "for i in range(len(a)):\n";
287 out << indent << indent << "c[i] = a[i] - b[i]\n";
288 out << indent << "return math.sqrt(dot(c, c))\n\n";
289
290 out << "def cross(a, b):\n";
291 out << indent << "c = [0, 0, 0]\n";
292 out << indent << "c[0] = a[1]*b[2] - a[2]*b[1]\n";
293 out << indent << "c[1] = a[2]*b[0] - a[0]*b[2]\n";
294 out << indent << "c[2] = a[0]*b[1] - a[1]*b[0]\n";
295 out << indent << "return c\n\n";
296
297 out << "class Quaternion:\n";
298 out << indent << "def __init__(self, values):\n";
299 out << indent << indent << "self.R = values\n\n";
300
301 out << indent << "def __str__(self):\n";
302 out << indent << indent << "return str(self.R)\n\n";
303
304 out << indent << "def __mul__(self, other):\n";
305 out << indent << indent << "imagThis = self.imag()\n";
306 out << indent << indent << "imagOther = other.imag()\n";
307 out << indent << indent << "cro = cross(imagThis, imagOther)\n\n";
308
309 out << indent << indent << "ret = ([self.R[0] * other.R[0] - dot(imagThis, imagOther)] +\n";
310 out << indent << indent << indent << " [self.R[0] * imagOther[0] + other.R[0] * imagThis[0] + cro[0],\n";
311 out << indent << indent << indent << indent << "self.R[0] * imagOther[1] + other.R[0] * imagThis[1] + cro[1],\n";
312 out << indent << indent << indent << indent << "self.R[0] * imagOther[2] + other.R[0] * imagThis[2] + cro[2]])\n\n";
313
314 out << indent << indent << "return Quaternion(ret)\n\n";
315
316 out << indent << "def imag(self):\n";
317 out << indent << indent << "return self.R[1:]\n\n";
318
319 out << indent << "def conjugate(self):\n";
320 out << indent << indent << "ret = [0] * 4\n";
321 out << indent << indent << "ret[0] = self.R[0]\n";
322 out << indent << indent << "ret[1] = -self.R[1]\n";
323 out << indent << indent << "ret[2] = -self.R[2]\n";
324 out << indent << indent << "ret[3] = -self.R[3]\n\n";
325
326 out << indent << indent << "return Quaternion(ret)\n\n";
327
328 out << indent << "def inverse(self):\n";
329 out << indent << indent << "return self.conjugate()\n\n";
330
331 out << indent << "def rotate(self, vec3D):\n";
332 out << indent << indent << "quat = Quaternion([0.0] + vec3D)\n";
333 out << indent << indent << "conj = self.conjugate()\n\n";
334
335 out << indent << indent << "ret = self * (quat * conj)\n";
336 out << indent << indent << "return ret.R[1:]\n\n";
337
338 out << "def getQuaternion(u, ref):\n";
339 out << indent << "u = normalize(u)\n";
340 out << indent << "ref = normalize(ref)\n";
341 out << indent << "axis = cross(u, ref)\n";
342 out << indent << "normAxis = math.sqrt(dot(axis, axis))\n\n";
343
344 out << indent << "if normAxis < 1e-12:\n";
345 out << indent << indent << "if math.fabs(dot(u, ref) - 1.0) < 1e-12:\n";
346 out << indent << indent << indent << "return Quaternion([1.0, 0.0, 0.0, 0.0])\n\n";
347
348 out << indent << indent << "return Quaternion([0.0, 1.0, 0.0, 0.0])\n\n";
349
350 out << indent << "axis = normalize(axis)\n";
351 out << indent << "cosAngle = math.sqrt(0.5 * (1 + dot(u, ref)))\n";
352 out << indent << "sinAngle = math.sqrt(1 - cosAngle * cosAngle)\n\n";
353
354 out << indent << "return Quaternion([cosAngle, sinAngle * axis[0], sinAngle * axis[1], sinAngle * axis[2]])\n\n";
355
356 out << "def exportVTK():\n";
357 out << indent << "vertices_str = \"\"\n";
358 out << indent << "triangles_str = \"\"\n";
359 out << indent << "color_str = \"\"\n";
360 out << indent << "cellTypes_str = \"\"\n";
361 out << indent << "startIdx = 0\n";
362 out << indent << "vertexCounter = 0\n";
363 out << indent << "cellCounter = 0\n";
364 out << indent << "lookup_table = []\n";
365 out << indent << "lookup_table.append([0.5, 0.5, 0.5, 0.5])\n";
366 out << indent << "lookup_table.append([1.0, 0.847, 0.0, 1.0])\n";
367 out << indent << "lookup_table.append([1.0, 0.0, 0.0, 1.0])\n";
368 out << indent << "lookup_table.append([0.537, 0.745, 0.525, 1.0])\n";
369 out << indent << "lookup_table.append([0.5, 0.5, 0.0, 1.0])\n";
370 out << indent << "lookup_table.append([1.0, 0.541, 0.0, 1.0])\n";
371 out << indent << "lookup_table.append([0.0, 0.0, 1.0, 1.0])\n\n";
372
373 out << indent << "decodeVertices()\n\n";
374
375 out << indent << "for i in range(len(vertices)):\n";
376 out << indent << indent << "for j in range(0, len(vertices[i]), 3):\n";
377 out << indent << indent << indent << "vertices_str += (\"%f %f %f\\n\" %(vertices[i][j], vertices[i][j+1], vertices[i][j+2]))\n";
378 out << indent << indent << indent << "vertexCounter += 1\n\n";
379
380 out << indent << indent << "for j in range(0, len(triangles[i]), 3):\n";
381 out << indent << indent << indent << "triangles_str += (\"3 %d %d %d\\n\" % (triangles[i][j] + startIdx, triangles[i][j+1] + startIdx, triangles[i][j+2] + startIdx))\n";
382 out << indent << indent << indent << "cellTypes_str += \"5\\n\"\n";
383 out << indent << indent << indent << "tmp_color = lookup_table[color[i]]\n";
384 out << indent << indent << indent << "color_str += (\"%f %f %f %f\\n\" % (tmp_color[0], tmp_color[1], tmp_color[2], tmp_color[3]))\n";
385 out << indent << indent << indent << "cellCounter += 1\n";
386
387 out << indent << indent << "startIdx = vertexCounter\n\n";
388
389 out << indent << "fh = open('" << fname << "_ElementPositions.vtk','w')\n";
390 out << indent << "fh.write(\"# vtk DataFile Version 2.0\\n\")\n";
391 out << indent << "fh.write(\"test\\nASCII\\n\\n\")\n";
392 out << indent << "fh.write(\"DATASET UNSTRUCTURED_GRID\\n\")\n";
393 out << indent << "fh.write(\"POINTS \" + str(vertexCounter) + \" float\\n\")\n";
394 out << indent << "fh.write(vertices_str)\n";
395 out << indent << "fh.write(\"CELLS \" + str(cellCounter) + \" \" + str(cellCounter * 4) + \"\\n\")\n";
396 out << indent << "fh.write(triangles_str)\n";
397 out << indent << "fh.write(\"CELL_TYPES \" + str(cellCounter) + \"\\n\")\n";
398 out << indent << "fh.write(cellTypes_str)\n";
399 out << indent << "fh.write(\"CELL_DATA \" + str(cellCounter) + \"\\n\")\n";
400 out << indent << "fh.write(\"COLOR_SCALARS type 4\\n\")\n";
401 out << indent << "fh.write(color_str + \"\\n\")\n";
402 out << indent << "fh.close()\n\n";
403
404 out << "def getNormal(tri_vertices):\n";
405 out << indent << "vec1 = [tri_vertices[1][0] - tri_vertices[0][0],\n";
406 out << indent << " tri_vertices[1][1] - tri_vertices[0][1],\n";
407 out << indent << " tri_vertices[1][2] - tri_vertices[0][2]]\n";
408 out << indent << "vec2 = [tri_vertices[2][0] - tri_vertices[0][0],\n";
409 out << indent << " tri_vertices[2][1] - tri_vertices[0][1],\n";
410 out << indent << " tri_vertices[2][2] - tri_vertices[0][2]]\n";
411 out << indent << "return normalize(cross(vec1,vec2))\n\n";
412
413 out << "def exportWeb(bgcolor):\n";
414 // out << indent << "if not os.path.exists('scenes'):\n";
415 // out << indent << indent << "os.makedirs('scenes')\n";
416 // out << indent << "fh = open('scenes/" << fname << "_ElementPositions.babylon','w')\n";
417 // out << indent << "indent = \"\";\n\n";
418
419 out << indent << "lookup_table = []\n";
420 out << indent << "lookup_table.append([0.5, 0.5, 0.5])\n";
421 out << indent << "lookup_table.append([1.0, 0.847, 0.0])\n";
422 out << indent << "lookup_table.append([1.0, 0.0, 0.0])\n";
423 out << indent << "lookup_table.append([0.537, 0.745, 0.525])\n";
424 out << indent << "lookup_table.append([0.5, 0.5, 0.0])\n";
425 out << indent << "lookup_table.append([1.0, 0.541, 0.0])\n";
426 out << indent << "lookup_table.append([0.0, 0.0, 1.0])\n\n";
427
428 out << indent << "decodeVertices()\n\n";
429
430 out << indent << "mesh = \"'data:\"\n";
431 out << indent << "mesh += \"{\"\n";
432 out << indent << "mesh += \"\\\"autoClear\\\":true,\"\n";
433 out << indent << "mesh += \"\\\"clearColor\\\":[0.0000,0.0000,0.0000],\"\n";
434 out << indent << "mesh += \"\\\"ambientColor\\\":[0.0000,0.0000,0.0000],\"\n";
435 out << indent << "mesh += \"\\\"gravity\\\":[0.0000,-9.8100,0.0000],\"\n";
436 out << indent << "mesh += \"\\\"cameras\\\":[\"\n";
437 out << indent << "mesh += \"{\"\n";
438 out << indent << "mesh += \"\\\"name\\\":\\\"Camera\\\",\"\n";
439 out << indent << "mesh += \"\\\"id\\\":\\\"Camera\\\",\"\n";
440 out << indent << "mesh += \"\\\"position\\\":[21.7936,2.2312,-85.7292],\"\n";
441 out << indent << "mesh += \"\\\"rotation\\\":[0.0432,-0.1766,-0.0668],\"\n";
442 out << indent << "mesh += \"\\\"fov\\\":0.8578,\"\n";
443 out << indent << "mesh += \"\\\"minZ\\\":10.0000,\"\n";
444 out << indent << "mesh += \"\\\"maxZ\\\":10000.0000,\"\n";
445 out << indent << "mesh += \"\\\"speed\\\":1.0000,\"\n";
446 out << indent << "mesh += \"\\\"inertia\\\":0.9000,\"\n";
447 out << indent << "mesh += \"\\\"checkCollisions\\\":false,\"\n";
448 out << indent << "mesh += \"\\\"applyGravity\\\":false,\"\n";
449 out << indent << "mesh += \"\\\"ellipsoid\\\":[0.2000,0.9000,0.2000]\"\n";
450 out << indent << "mesh += \"}],\"\n";
451 out << indent << "mesh += \"\\\"activeCamera\\\":\\\"Camera\\\",\"\n";
452 out << indent << "mesh += \"\\\"lights\\\":[\"\n";
453 out << indent << "mesh += \"{\"\n";
454 out << indent << "mesh += \"\\\"name\\\":\\\"Lamp\\\",\"\n";
455 out << indent << "mesh += \"\\\"id\\\":\\\"Lamp\\\",\"\n";
456 out << indent << "mesh += \"\\\"type\\\":0.0000,\"\n";
457 out << indent << "mesh += \"\\\"position\\\":[4.0762,34.9321,-63.5788],\"\n";
458 out << indent << "mesh += \"\\\"intensity\\\":1.0000,\"\n";
459 out << indent << "mesh += \"\\\"diffuse\\\":[1.0000,1.0000,1.0000],\"\n";
460 out << indent << "mesh += \"\\\"specular\\\":[1.0000,1.0000,1.0000]\"\n";
461 out << indent << "mesh += \"}],\"\n";
462 out << indent << "mesh += \"\\\"materials\\\":[],\"\n";
463 out << indent << "mesh += \"\\\"meshes\\\":[\"\n\n";
464
465 out << indent << "for i in range(len(triangles)):\n";
466 out << indent << indent << "vertex_list = []\n";
467 out << indent << indent << "indices_list = []\n";
468 out << indent << indent << "normals_list = []\n";
469 out << indent << indent << "color_list = []\n";
470 out << indent << indent << "for j in range(0, len(triangles[i]), 3):\n";
471 out << indent << indent << indent << "tri_vertices = []\n";
472 out << indent << indent << indent << "idcs = triangles[i][j:j + 3]\n";
473 out << indent << indent << indent << "tri_vertices.append(vertices[i][3 * idcs[0]:3 * (idcs[0] + 1)])\n";
474 out << indent << indent << indent << "tri_vertices.append(vertices[i][3 * idcs[1]:3 * (idcs[1] + 1)])\n";
475 out << indent << indent << indent << "tri_vertices.append(vertices[i][3 * idcs[2]:3 * (idcs[2] + 1)])\n";
476 out << indent << indent << indent << "indices_list.append(','.join(str(n) for n in range(len(vertex_list),len(vertex_list) + 3)))\n";
477 out << indent << indent << indent << "# left hand order!\n";
478 out << indent << indent << indent << "vertex_list.append(','.join(\"%.5f\" % (round(n,5) + 0.0) for n in tri_vertices[0]))\n";
479 out << indent << indent << indent << "vertex_list.append(','.join(\"%.5f\" % (round(n,5) + 0.0) for n in tri_vertices[2]))\n";
480 out << indent << indent << indent << "vertex_list.append(','.join(\"%.5f\" % (round(n,5) + 0.0) for n in tri_vertices[1]))\n";
481 out << indent << indent << indent << "normal = getNormal(tri_vertices)\n";
482 out << indent << indent << indent << "normals_list.append(','.join(\"%.5f\" % (round(n,5) + 0.0) for n in normal * 3))\n";
483 out << indent << indent << indent << "color_list.append(','.join([str(n) for n in lookup_table[color[i]]] * 3))\n";
484 out << indent << indent << "mesh += \"{\"\n";
485 out << indent << indent << "mesh += \"\\\"name\\\":\\\"element_\" + str(i) + \"\\\",\"\n";
486 out << indent << indent << "mesh += \"\\\"id\\\":\\\"element_\" + str(i) + \"\\\",\"\n";
487 out << indent << indent << "mesh += \"\\\"position\\\":[0.0,0.0,0.0],\"\n";
488 out << indent << indent << "mesh += \"\\\"rotation\\\":[0.0,0.0,0.0],\"\n";
489 out << indent << indent << "mesh += \"\\\"scaling\\\":[1.0,1.0,1.0],\"\n";
490 out << indent << indent << "mesh += \"\\\"isVisible\\\":true,\"\n";
491 out << indent << indent << "mesh += \"\\\"isEnabled\\\":true,\"\n";
492 out << indent << indent << "mesh += \"\\\"useFlatShading\\\":false,\"\n";
493 out << indent << indent << "mesh += \"\\\"checkCollisions\\\":false,\"\n";
494 out << indent << indent << "mesh += \"\\\"billboardMode\\\":0,\"\n";
495 out << indent << indent << "mesh += \"\\\"receiveShadows\\\":false,\"\n";
496 out << indent << indent << "mesh += \"\\\"positions\\\":[\" + ','.join(vertex_list) + \"],\"\n";
497 out << indent << indent << "mesh += \"\\\"normals\\\":[\" + ','.join(normals_list) + \"],\"\n";
498 out << indent << indent << "mesh += \"\\\"indices\\\":[\" + ','.join(indices_list) + \"],\"\n";
499 out << indent << indent << "mesh += \"\\\"colors\\\":[\" + ','.join(color_list) + \"],\"\n";
500 out << indent << indent << "mesh += \"\\\"subMeshes\\\":[\"\n";
501 out << indent << indent << "mesh += \"{\"\n";
502 out << indent << indent << "mesh += \"\\\"materialIndex\\\":0,\"\n";
503 out << indent << indent << "mesh += \" \\\"verticesStart\\\":0,\"\n";
504 out << indent << indent << "mesh += \" \\\"verticesCount\\\":\" + str(len(triangles[i])) + \",\"\n";
505 out << indent << indent << "mesh += \" \\\"indexStart\\\":0,\"\n";
506 out << indent << indent << "mesh += \" \\\"indexCount\\\":\" + str(len(triangles[i])) + \"\"\n";
507 out << indent << indent << "mesh += \"}]\"\n";
508 out << indent << indent << "mesh += \"}\"\n";
509 out << indent << indent << "mesh += \",\"\n\n";
510
511 out << indent << indent << "del normals_list[:]\n";
512 out << indent << indent << "del vertex_list[:]\n";
513 out << indent << indent << "del color_list[:]\n";
514 out << indent << indent << "del indices_list[:]\n\n";
515
516 out << indent << "mesh = mesh[:-1] + \"]\"\n";
517 out << indent << "mesh += \"}'\"\n";
518
519 out << indent << "index_compressed = base64.b64decode(index_base64)\n";
520 out << indent << "index = str(zlib.decompress(index_compressed))\n";
521 out << indent << "if (len(bgcolor) == 3):\n";
522 out << indent << indent << "mesh += \";\\n \"\n";
523 out << indent << indent << "mesh += \"scene.clearColor = new BABYLON.Color3(%f, %f, %f)\" % (bgcolor[0], bgcolor[1], bgcolor[2])\n\n";
524 out << indent << "index = index.replace('##DATA##', mesh)\n";
525 out << indent << "fh = open('" << fname << "_ElementPositions.html','w')\n";
526 out << indent << "fh.write(index)\n";
527 out << indent << "fh.close()\n\n";
528
529 out << "def computeMinAngle(idx, curAngle, positions, connections, check):\n";
530 out << indent << "matrix = [[-math.cos(curAngle), -math.sin(curAngle)], [math.sin(curAngle), -math.cos(curAngle)]]\n\n";
531
532 out << indent << "minAngle = 2 * math.pi\n";
533 out << indent << "nextIdx = -1\n\n";
534
535 out << indent << "for j in connections[idx]:\n";
536 out << indent << indent << "direction = [positions[j][0] - positions[idx][0],\n";
537 out << indent << indent << indent << " positions[j][1] - positions[idx][1]]\n";
538 out << indent << indent << "direction = [dot(matrix[0],direction), dot(matrix[1],direction)]\n\n";
539
540 out << indent << indent << "if math.fabs(dot([1.0, 0.0], direction) / distance(positions[j], positions[idx]) - 1.0) < 1e-6: continue\n\n";
541
542 out << indent << indent << "angle = math.fmod(math.atan2(direction[1],direction[0]) + 2 * math.pi, 2 * math.pi)\n\n";
543
544 out << indent << indent << "if angle < minAngle:\n";
545 out << indent << indent << indent << "nextIdx = j\n";
546 out << indent << indent << indent << "minAngle = angle\n";
547 out << indent << indent << "if angle == minAngle and check:\n";
548 out << indent << indent << indent << "dire = [positions[j][0] - positions[idx][0],\n";
549 out << indent << indent << indent << " positions[j][1] - positions[idx][1]]\n";
550 out << indent << indent << indent << "minA0 = math.atan2(dire[1], dire[0])\n";
551 out << indent << indent << indent << "minA1 = computeMinAngle(nextIdx, minA0, positions, connections, False)\n";
552 out << indent << indent << indent << "minA2 = computeMinAngle(j, minA0, positions, connections, False)\n";
553 out << indent << indent << indent << "if minA2[1] < minA2[1]:\n";
554 out << indent << indent << indent << indent << "nextIdx = j\n\n";
555
556 out << indent << "if nextIdx == -1:\n";
557 out << indent << indent << "nextIdx = connections[idx][0]\n\n";
558
559 out << indent << "return (nextIdx, minAngle)\n\n";
560
561 out << "def squashVertices(positionDict, connections):\n";
562 out << indent << "removedItems = []\n";
563 out << indent << "indices = [int(k) for k in positionDict.keys()]\n";
564 out << indent << "idxChanges = indices\n";
565 out << indent << "for i in indices:\n";
566 out << indent << indent << "if i in removedItems:\n";
567 out << indent << indent << indent << "continue\n";
568 out << indent << indent << "for j in indices:\n";
569 out << indent << indent << indent << "if j in removedItems or j == i:\n";
570 out << indent << indent << indent << indent << "continue\n\n";
571
572 out << indent << indent << indent << "if distance(positionDict[i], positionDict[j]) < 1e-6:\n";
573 out << indent << indent << indent << indent << "connections[j] = list(set(connections[j]))\n";
574 out << indent << indent << indent << indent << "if i in connections[j]:\n";
575 out << indent << indent << indent << indent << indent << "connections[j].remove(i)\n";
576 out << indent << indent << indent << indent << "if j in connections[i]:\n";
577 out << indent << indent << indent << indent << indent << "connections[i].remove(j)\n\n";
578
579 out << indent << indent << indent << indent << "connections[i].extend(connections[j])\n";
580 out << indent << indent << indent << indent << "connections[i] = list(set(connections[i]))\n";
581 out << indent << indent << indent << indent << "connections[i].sort()\n\n";
582
583 out << indent << indent << indent << indent << "for k in connections.keys():\n";
584 out << indent << indent << indent << indent << indent << "if k == i: continue\n";
585 out << indent << indent << indent << indent << indent << "if j in connections[k]:\n";
586 out << indent << indent << indent << indent << indent << indent << "idx = connections[k].index(j)\n";
587 out << indent << indent << indent << indent << indent << indent << "connections[k][idx] = i\n\n";
588
589 out << indent << indent << indent << indent << "idxChanges[j] = i\n";
590 out << indent << indent << indent << indent << "removedItems.append(j)\n\n";
591
592 out << indent << "for i in removedItems:\n";
593 out << indent << indent << "del positionDict[i]\n";
594 out << indent << indent << "del connections[i]\n\n";
595
596 out << indent << "for i in connections.keys():\n";
597 out << indent << indent << "connections[i] = list(set(connections[i]))\n";
598 out << indent << indent << "connections[i].sort()\n\n";
599
600 out << indent << "return idxChanges\n\n";
601
602 out << "def computeLineEquations(positions, connections):\n";
603 out << indent << "cons = set()\n";
604 out << indent << "for i in connections.keys():\n";
605 out << indent << indent << "for j in connections[i]:\n";
606 out << indent << indent << indent << "cons.add((min(i, j), max(i, j)))\n\n";
607
608 out << indent << "lineEquations = {}\n";
609 out << indent << "for item in cons:\n";
610 out << indent << indent << "a = (positions[item[1]][1] - positions[item[0]][1])\n";
611 out << indent << indent << "b = -(positions[item[1]][0] - positions[item[0]][0])\n\n";
612
613 out << indent << indent << "xm = 0.5 * (positions[item[0]][0] + positions[item[1]][0])\n";
614 out << indent << indent << "ym = 0.5 * (positions[item[0]][1] + positions[item[1]][1])\n";
615 out << indent << indent << "c = -(a * xm + b * ym)\n\n";
616
617 out << indent << indent << "lineEquations[item] = (a, b, c)\n\n";
618
619 out << indent << "return lineEquations\n\n";
620
621 out << "def checkPossibleSegmentIntersection(segment, positions, lineEquations, minAngle, lastIntersection):\n";
622 out << indent << "item1 = (min(segment), max(segment))\n\n";
623
624 out << indent << "(a1, b1, c1) = (0,0,0)\n";
625 out << indent << "A = [0]*2\n";
626 out << indent << "B = A\n\n";
627
628 out << indent << "if segment[0] == None:\n";
629 out << indent << indent << "A = lastIntersection\n";
630 out << indent << indent << "B = positions[segment[1]]\n\n";
631
632 out << indent << indent << "a1 = B[1] - A[1]\n";
633 out << indent << indent << "b1 = -(B[0] - A[0])\n";
634 out << indent << indent << "xm = 0.5 * (A[0] + B[0])\n";
635 out << indent << indent << "ym = 0.5 * (A[1] + B[1])\n";
636 out << indent << indent << "c1 = -(a1 * xm + b1 * ym)\n\n";
637
638 out << indent << "else:\n";
639 out << indent << indent << "A = positions[segment[0]]\n";
640 out << indent << indent << "B = positions[segment[1]]\n\n";
641
642 out << indent << indent << "(a1, b1, c1) = lineEquations[item1]\n\n";
643
644 out << indent << indent << "if segment[1] < segment[0]:\n";
645 out << indent << indent << indent << "(a1, b1, c1) = (-a1, -b1, -c1)\n\n";
646
647 out << indent << "curAngle = math.atan2(a1, -b1)\n";
648 out << indent << "matrix = [[-math.cos(curAngle), -math.sin(curAngle)], [math.sin(curAngle), -math.cos(curAngle)]]\n\n";
649
650 out << indent << "origMinAngle = minAngle\n\n";
651
652 out << indent << "segment1 = [B[0] - A[0], B[1] - A[1], 0.0]\n\n";
653
654 out << indent << "intersectingSegments = []\n";
655 out << indent << "distanceAB = distance(A, B)\n";
656 out << indent << "for item2 in lineEquations.keys():\n";
657 out << indent << indent << "item = item2\n";
658 out << indent << indent << "C = positions[item[0]]\n";
659 out << indent << indent << "D = positions[item[1]]\n\n";
660
661 out << indent << indent << "if segment[0] == None:\n";
662 out << indent << indent << indent << "if (segment[1] == item[0] or\n";
663 out << indent << indent << indent << indent << "segment[1] == item[1]): continue\n";
664 out << indent << indent << "else:\n";
665 out << indent << indent << indent << "if (item1[0] == item[0] or\n";
666 out << indent << indent << indent << indent << "item1[1] == item[0] or\n";
667 out << indent << indent << indent << indent << "item1[0] == item[1] or\n";
668 out << indent << indent << indent << indent << "item1[1] == item[1]): continue\n\n";
669
670 out << indent << indent << "nodes = set([item1[0],item1[1],item[0],item[1]])\n";
671 out << indent << indent << "if len(nodes) < 4: continue # share same vertex\n\n";
672
673 out << indent << indent << "(a2, b2, c2) = lineEquations[item]\n\n";
674
675 out << indent << indent << "segment2 = [C[0] - A[0], C[1] - A[1], 0.0]\n";
676 out << indent << indent << "segment3 = [D[0] - A[0], D[1] - A[1], 0.0]\n\n";
677
678 out << indent << indent << "# check that C and D aren't on the same side of AB\n";
679 out << indent << indent << "if cross(segment1, segment2)[2] * cross(segment1, segment3)[2] > 0.0: continue\n\n";
680
681 out << indent << indent << "if cross(segment1, segment2)[2] < 0.0 or cross(segment1, segment3)[2] > 0.0:\n";
682 out << indent << indent << indent << "(C, D, a2, b2, c2) = (D, C, -a2, -b2, -c2)\n";
683 out << indent << indent << indent << "item = (item[1], item[0])\n\n";
684
685 out << indent << indent << "denominator = a1 * b2 - b1 * a2\n";
686 out << indent << indent << "if math.fabs(denominator) < 1e-9: continue\n\n";
687
688 out << indent << indent << "px = (b1 * c2 - b2 * c1) / denominator\n";
689 out << indent << indent << "py = (a2 * c1 - a1 * c2) / denominator\n\n";
690
691 out << indent << indent << "distanceCD = distance(C, D)\n\n";
692
693 out << indent << indent << "distanceAP = distance(A, [px, py])\n";
694 out << indent << indent << "distanceBP = distance(B, [px, py])\n";
695 out << indent << indent << "distanceCP = distance(C, [px, py])\n";
696 out << indent << indent << "distanceDP = distance(D, [px, py])\n\n";
697
698 out << indent << indent << "# check intersection is between AB and CD\n";
699 out << indent << indent << "check1 = (distanceAP - 1e-6 < distanceAB and distanceBP - 1e-6 < distanceAB)\n";
700 out << indent << indent << "check2 = (distanceCP - 1e-6 < distanceCD and distanceDP - 1e-6 < distanceCD)\n";
701 out << indent << indent << "if not check1 or not check2: continue\n\n";
702
703 out << indent << indent << "if math.fabs(dot(segment1, [D[0] - C[0], D[1] - C[1], 0.0]) / (distanceAB * distanceCD) + 1.0) < 1e-9: continue\n\n";
704
705 out << indent << indent << "if ((distanceAP < 1e-6) or\n";
706 out << indent << indent << indent << "(distanceBP < 1e-6 and distanceCP < 1e-6) or\n";
707 out << indent << indent << indent << "(distanceDP < 1e-6)):\n";
708 out << indent << indent << indent << "continue\n\n";
709
710 out << indent << indent << "direction = [D[0] - C[0], D[1] - C[1]]\n";
711 out << indent << indent << "direction = [dot(matrix[0], direction), dot(matrix[1], direction)]\n";
712 out << indent << indent << "angle = math.fmod(math.atan2(direction[1], direction[0]) + 2 * math.pi, 2 * math.pi)\n\n";
713
714 out << indent << indent << "newSegment = ([px, py], item[1], distanceAP, angle)\n\n";
715
716 out << indent << indent << "if distanceCP < 1e-6 and angle > origMinAngle: continue\n";
717 out << indent << indent << "if distanceBP > 1e-6 and angle > math.pi: continue\n\n";
718
719 out << indent << indent << "if len(intersectingSegments) == 0:\n";
720 out << indent << indent << indent << "intersectingSegments.append(newSegment)\n";
721 out << indent << indent << indent << "minAngle = angle\n";
722 out << indent << indent << "else:\n";
723 out << indent << indent << indent << "if intersectingSegments[0][2] - 1e-9 > distanceAP:\n";
724 out << indent << indent << indent << indent << "intersectingSegments[0] = newSegment\n";
725 out << indent << indent << indent << indent << "minAngle = angle\n";
726 out << indent << indent << indent << "elif intersectingSegments[0][2] + 1e-9 > distanceAP and angle < minAngle:\n";
727 out << indent << indent << indent << indent << "intersectingSegments[0] = newSegment\n";
728 out << indent << indent << indent << indent << "minAngle = angle\n\n";
729
730 out << indent << "return intersectingSegments\n\n";
731
732 out << "def projectToPlane(normal):\n";
733 out << indent << "fh = open('" << fname << "_ElementPositions.gpl','w')\n";
734 // out << indent << "fg = open('test2.dat', 'w')\n";
735 out << indent << "normal = normalize(normal)\n";
736 out << indent << "ori = getQuaternion(normal, [0, 0, 1])\n\n";
737
738 out << indent << "left2Right = [0, 0, 1]\n";
739 out << indent << "if math.fabs(math.fabs(dot(normal, left2Right)) - 1) < 1e-9:\n";
740 out << indent << indent << "left2Right = [1, 0, 0]\n";
741 out << indent << "rotL2R = ori.rotate(left2Right)\n";
742 out << indent << "deviation = math.atan2(rotL2R[1], rotL2R[0])\n";
743 out << indent << "rotBack = Quaternion([math.cos(0.5 * deviation), 0, 0, -math.sin(0.5 * deviation)])\n";
744 out << indent << "ori = rotBack * ori\n\n";
745
746 out << indent << "decodeVertices()\n\n";
747
748 out << indent << "for i in range(len(vertices)):\n";
749 out << indent << indent << "positions = {}\n";
750 out << indent << indent << "connections = {}\n";
751 out << indent << indent << "for j in range(0, len(vertices[i]), 3):\n";
752 out << indent << indent << indent << "nextPos3D = ori.rotate(vertices[i][j:j+3])\n";
753 out << indent << indent << indent << "nextPos2D = nextPos3D[0:2]\n";
754 out << indent << indent << indent << "positions[j/3] = nextPos2D\n\n";
755
756 out << indent << indent << "if len(positions) == 0:\n";
757 out << indent << indent << indent << "continue\n\n";
758
759 out << indent << indent << "idx = 0\n";
760 out << indent << indent << "maxX = positions[0][0]\n";
761 out << indent << indent << "for j in positions.keys():\n";
762 out << indent << indent << indent << "if positions[j][0] > maxX:\n";
763 out << indent << indent << indent << indent << "maxX = positions[j][0]\n";
764 out << indent << indent << indent << indent << "idx = j\n";
765 out << indent << indent << indent << "if positions[j][0] == maxX and positions[j][1] > positions[idx][1]:\n";
766 out << indent << indent << indent << indent << "idx = j\n\n";
767
768 out << indent << indent << "for j in range(0, len(triangles[i]), 3):\n";
769 out << indent << indent << indent << "for k in range(0, 3):\n";
770 out << indent << indent << indent << indent << "vertIdx = triangles[i][j + k]\n";
771 out << indent << indent << indent << indent << "if not vertIdx in connections:\n";
772 out << indent << indent << indent << indent << indent << "connections[vertIdx] = []\n";
773 out << indent << indent << indent << indent << "for l in range(1, 3):\n";
774 out << indent << indent << indent << indent << indent << "connections[vertIdx].append(triangles[i][j + ((k + l) % 3)])\n\n";
775
776 out << indent << indent << "numConnections = 0\n";
777 out << indent << indent << "for j in connections.keys():\n";
778 out << indent << indent << indent << "connections[j] = list(set(connections[j]))\n";
779 out << indent << indent << indent << "numConnections += len(connections[j])\n\n";
780 out << indent << indent << "numConnections /= 2\n";
781
782 out << indent << indent << "idChanges = squashVertices(positions, connections)\n\n";
783
784 out << indent << indent << "lineEquations = computeLineEquations(positions, connections)\n\n";
785
786 // out << indent << indent << "for j in positions.keys():\n";
787 // out << indent << indent << indent << "for k in connections[j]:\n";
788 // out << indent << indent << indent << indent << "if j < k:\n";
789 // out << indent << indent << indent << indent << indent << "fg.write(\"%.8f %.8f \\n%.8f %.8f\\n# %d %d\\n\\n\" %(positions[j][0], positions[j][1], positions[k][0], positions[k][1], j, k))\n";
790 // out << indent << indent << indent << "fg.write(\"\\n\")\n\n";
791
792 out << indent << indent << "idx = idChanges[int(idx)]\n";
793 out << indent << indent << "fh.write(\"%.6f %.6f\\n\" % (positions[idx][0], positions[idx][1]))\n\n";
794
795 out << indent << indent << "curAngle = math.pi\n";
796 out << indent << indent << "origin = idx\n";
797 out << indent << indent << "count = 0\n";
798 out << indent << indent << "passed = []\n";
799 out << indent << indent << "while (count == 0 or distance(positions[origin], positions[idx]) > 1e-9) and not count > numConnections:\n";
800 out << indent << indent << indent << "nextGen = computeMinAngle(idx, curAngle, positions, connections, False)\n";
801 out << indent << indent << indent << "nextIdx = nextGen[0]\n";
802 out << indent << indent << indent << "direction = [positions[nextIdx][0] - positions[idx][0],\n";
803 out << indent << indent << indent << indent << " positions[nextIdx][1] - positions[idx][1]]\n";
804 out << indent << indent << indent << "curAngle = math.atan2(direction[1], direction[0])\n\n";
805
806 out << indent << indent << indent << "intersections = checkPossibleSegmentIntersection((idx, nextIdx), positions, lineEquations, nextGen[1], [])\n";
807 out << indent << indent << indent << "if len(intersections) > 0:\n";
808 out << indent << indent << indent << indent << "while len(intersections) > 0 and not count > numConnections:\n";
809 out << indent << indent << indent << indent << indent << "fh.write(\"%.6f %.6f\\n\" %(intersections[0][0][0], intersections[0][0][1]))\n";
810 out << indent << indent << indent << indent << indent << "count += 1\n";
811 out << indent << "\n";
812 out << indent << indent << indent << indent << indent << "idx = intersections[0][1]\n";
813 out << indent << indent << indent << indent << indent << "direction = [positions[idx][0] - intersections[0][0][0],\n";
814 out << indent << indent << indent << indent << indent << " positions[idx][1] - intersections[0][0][1]]\n";
815 out << indent << indent << indent << indent << indent << "curAngle = math.atan2(direction[1], direction[0])\n";
816 out << indent << indent << indent << indent << indent << "nextGen = computeMinAngle(idx, curAngle, positions, connections, False)\n";
817 out << indent << "\n";
818 out << indent << indent << indent << indent << indent << "intersections = checkPossibleSegmentIntersection((None, idx), positions, lineEquations, nextGen[1], intersections[0][0])\n";
819 out << indent << indent << indent << "else:\n";
820 out << indent << indent << indent << indent << "idx = nextIdx\n";
821 out << indent << "\n";
822 out << indent << indent << indent << "fh.write(\"%.6f %.6f\\n\" % (positions[idx][0], positions[idx][1]))\n";
823 out << indent << "\n";
824 out << indent << indent << indent << "if idx in passed:\n";
825 out << indent << indent << indent << indent << "direction1 = [positions[idx][0] - positions[passed[-1]][0],\n";
826 out << indent << indent << indent << indent << indent << " positions[idx][1] - positions[passed[-1]][1]]\n";
827 out << indent << indent << indent << indent << "direction2 = [positions[origin][0] - positions[passed[-1]][0],\n";
828 out << indent << indent << indent << indent << indent << " positions[origin][1] - positions[passed[-1]][1]]\n";
829 out << indent << indent << indent << indent << "dist1 = distance(positions[idx], positions[passed[-1]])\n";
830 out << indent << indent << indent << indent << "dist2 = distance(positions[origin], positions[passed[-1]])\n";
831 out << indent << indent << indent << indent << "if dist1 * dist2 > 0.0 and math.fabs(math.fabs(dot(direction1, direction2) / (dist1 * dist2)) - 1.0) > 1e-9:\n";
832 out << indent << indent << indent << indent << indent << "sys.stderr.write(\"error: projection cycling on element id: %d, vertex id: %d\\n\" %(i, idx))\n";
833 out << indent << indent << indent << indent << "break\n\n";
834
835 out << indent << indent << indent << "passed.append(idx)\n";
836 out << indent << indent << indent << "count += 1\n";
837 out << indent << indent << "fh.write(\"\\n\")\n\n";
838
839 out << indent << indent << "if count > numConnections:\n";
840 out << indent << indent << indent << "sys.stderr.write(\"error: projection cycling on element id: %d\\n\" % i)\n\n";
841
842 out << indent << indent << "for j in range(0, len(decoration[i]), 6):\n";
843 out << indent << indent << indent << "for k in range(j, j + 6, 3):\n";
844 out << indent << indent << indent << indent << "nextPos3D = ori.rotate(decoration[i][k:k+3])\n";
845 out << indent << indent << indent << indent << "fh.write(\"%.6f %.6f\\n\" % (nextPos3D[0], nextPos3D[1]))\n";
846 out << indent << indent << indent << "fh.write(\"\\n\")\n\n";
847
848 out << indent << "fh.close()\n\n";
849
850 out << "if __name__ == \"__main__\":\n";
851 out << indent << "parser = argparse.ArgumentParser()\n";
852 out << indent << "parser.add_argument('--export-vtk', action='store_true')\n";
853 out << indent << "parser.add_argument('--export-web', action='store_true')\n";
854 out << indent << "parser.add_argument('--background', nargs=3, type=float)\n";
855 out << indent << "parser.add_argument('--project-to-plane', action='store_true')\n";
856 out << indent << "parser.add_argument('--normal', nargs=3, type=float)\n";
857 out << indent << "args = parser.parse_args()\n\n";
858
859 out << indent << "if (args.export_vtk):\n";
860 out << indent << indent << "exportVTK()\n";
861 out << indent << indent << "sys.exit()\n\n";
862
863 out << indent << "if (args.export_web):\n";
864 out << indent << indent << "bgcolor = []\n";
865 out << indent << indent << "if (args.background):\n";
866 out << indent << indent << indent << "validBackground = True\n";
867 out << indent << indent << indent << "for comp in bgcolor:\n";
868 out << indent << indent << indent << indent << "if comp < 0.0 or comp > 1.0:\n";
869 out << indent << indent << indent << indent << indent << "validBackground = False\n";
870 out << indent << indent << indent << indent << indent << "break\n";
871 out << indent << indent << indent << "if (validBackground):\n";
872 out << indent << indent << indent << indent << "bgcolor = args.background\n";
873 out << indent << indent << "exportWeb(bgcolor)\n";
874 out << indent << indent << "sys.exit()\n\n";
875
876 out << indent << "if (args.project_to_plane):\n";
877 out << indent << indent << "normal = [0.0, 1.0, 0.0]\n";
878 out << indent << indent << "if (args.normal):\n";
879 out << indent << indent << indent << "normal = args.normal\n";
880 out << indent << indent << "projectToPlane(normal)\n";
881 out << indent << indent << "sys.exit()\n\n";
882
883 out << indent << "parser.print_help()\n";
884}
885
886MeshData MeshGenerator::getCylinder(double length,
887 double minor,
888 double major,
889 double formFactor,
890 const unsigned int numSegments) {
891 double angle = 0;
892 double dAngle = Physics::two_pi / numSegments;
893
894 MeshData mesh;
895 mesh.vertices_m.push_back(Vector_t(0.0));
896 for (unsigned int i = 0; i < numSegments; ++ i, angle += dAngle) {
897 Vector_t node(major * cos(angle), minor * sin(angle), 0);
898 mesh.vertices_m.push_back(node);
899
900 unsigned int next = (i + 1) % numSegments;
901 Vektor<unsigned int, 3> baseTriangle(0u, next + 1, i + 1);
902 mesh.triangles_m.push_back(baseTriangle);
903
904 Vektor<unsigned int, 3> sideTriangle1(i + 1, next + 1, i + numSegments + 2);
905 mesh.triangles_m.push_back(sideTriangle1);
906
907 Vektor<unsigned int, 3> sideTriangle2(next + 1, next + numSegments + 2, i + numSegments + 2);
908 mesh.triangles_m.push_back(sideTriangle2);
909 }
910
911 mesh.vertices_m.push_back(Vector_t(0.0, 0.0, length));
912 for (unsigned int i = 0; i < numSegments; ++ i, angle += dAngle) {
913 Vector_t node(formFactor * major * cos(angle), formFactor * minor * sin(angle), length);
914 mesh.vertices_m.push_back(node);
915
916 unsigned int next = (i + 1) % numSegments;
917 Vektor<unsigned int, 3> topTriangle(numSegments + 1, i + numSegments + 2, next + numSegments + 2);
918 mesh.triangles_m.push_back(topTriangle);
919 }
920
921 mesh.decorations_m.push_back(std::make_pair(Vector_t(0.0), Vector_t(0, 0, length)));
922
923 return mesh;
924}
925
926MeshData MeshGenerator::getBox(double length,
927 double width,
928 double height,
929 double formFactor) {
930
931 MeshData mesh;
932 mesh.vertices_m.push_back(Vector_t(width, height, 0.0));
933 mesh.vertices_m.push_back(Vector_t(-width, height, 0.0));
934 mesh.vertices_m.push_back(Vector_t(-width, -height, 0.0));
935 mesh.vertices_m.push_back(Vector_t(width, -height, 0.0));
936
937 mesh.vertices_m.push_back(Vector_t(formFactor * width, formFactor * height, length));
938 mesh.vertices_m.push_back(Vector_t(-formFactor * width, formFactor * height, length));
939 mesh.vertices_m.push_back(Vector_t(-formFactor * width, -formFactor * height, length));
940 mesh.vertices_m.push_back(Vector_t(formFactor * width, -formFactor * height, length));
941
942 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(0, 2, 1));
943 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(0, 3, 2));
944
945 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(0, 1, 4));
946 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(4, 1, 5));
947 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(1, 2, 5));
948 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(5, 2, 6));
949 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(2, 3, 6));
950 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(6, 3, 7));
951 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(3, 0, 7));
952 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(7, 0, 4));
953
954 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(4, 5, 6));
955 mesh.triangles_m.push_back(Vektor<unsigned int, 3>(4, 6, 7));
956
957 mesh.decorations_m.push_back(std::make_pair(Vector_t(0.0), Vector_t(0, 0, length)));
958
959 return mesh;
960}
end
PartBunchBase< T, Dim >::ConstIterator end(PartBunchBase< T, Dim > const &bunch)
Definition: PartBunchBase.h:736
gmsg
Inform * gmsg
Definition: Main.cpp:61
cos
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
sin
Tps< T > sin(const Tps< T > &x)
Sine.
Definition: TpsMath.h:111
Physics::two_pi
constexpr double two_pi
The value of.
Definition: Physics.h:33
Util::combineFilePath
std::string combineFilePath(std::initializer_list< std::string > ilist)
Definition: Util.cpp:196
Util::base64_encode
std::string base64_encode(const std::string &string_to_encode)
Definition: Util.cpp:364
OpalData::getInstance
static OpalData * getInstance()
Definition: OpalData.cpp:196
OpalData::getAuxiliaryOutputDirectory
std::string getAuxiliaryOutputDirectory() const
get the name of the the additional data directory
Definition: OpalData.cpp:661
Bend2D
Definition: Bend2D.h:51
Bend2D::getSurfaceMesh
MeshData getSurfaceMesh() const
Definition: Bend2D.cpp:1387
ElementBase::getAperture
std::pair< ApertureType, std::vector< double > > getAperture() const
Definition: ElementBase.h:524
ElementBase::getCSTrafoGlobal2Local
CoordinateSystemTrafo getCSTrafoGlobal2Local() const
Definition: ElementBase.h:497
ElementBase::getType
virtual ElementType getType() const =0
Get element type std::string.
ElementBase::getElementDimensions
virtual void getElementDimensions(double &begin, double &end) const
Definition: ElementBase.h:173
Multipole
Interface for general multipole.
Definition: Multipole.h:47
RBend3D
Interface for solenoids.
Definition: RBend3D.h:39
RBend3D::getSurfaceMesh
MeshData getSurfaceMesh() const
Definition: RBend3D.cpp:236
CoordinateSystemTrafo::transformTo
Vector_t transformTo(const Vector_t &r) const
Definition: CoordinateSystemTrafo.h:78
CoordinateSystemTrafo::rotateTo
Vector_t rotateTo(const Vector_t &r) const
Definition: CoordinateSystemTrafo.h:88
CoordinateSystemTrafo::inverted
CoordinateSystemTrafo inverted() const
Definition: CoordinateSystemTrafo.h:70
MeshData::decorations_m
std::vector< std::pair< Vector_t, Vector_t > > decorations_m
Definition: MeshGenerator.h:26
MeshData::triangles_m
std::vector< Vektor< unsigned int, 3 > > triangles_m
Definition: MeshGenerator.h:25
MeshData::vertices_m
std::vector< Vector_t > vertices_m
Definition: MeshGenerator.h:24
MeshGenerator::add
void add(const ElementBase &element)
Definition: MeshGenerator.cpp:18
MeshGenerator::getBox
static MeshData getBox(double length, double width, double height, double formFactor)
Definition: MeshGenerator.cpp:926
MeshGenerator::elements_m
std::vector< MeshData > elements_m
Definition: MeshGenerator.h:60
MeshGenerator::write
void write(const std::string &fname)
Definition: MeshGenerator.cpp:103
MeshGenerator::getCylinder
static MeshData getCylinder(double length, double minor, double major, double formFactor, const unsigned int numSegments=36)
Definition: MeshGenerator.cpp:886
Inform
Definition: Inform.h:42
Vector_t
Vektor< double, 3 > Vector_t
Definition: Vektor.h:6
Generated on Thu Oct 20 2022 17:40:11 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.3