1 #include "../config.h"
2 #include <stdio.h>
3 #ifdef WIN32
4 #include <direct.h>
5 #else
6 #include <unistd.h>
7 #endif
8
9 #include <iostream>
10 #include <math.h>
11 #include "Mesh.h"
12 #include "GeometryNode.h"
13 #include "GeometryEdge.h"
14 #include "Connect.h"
15 #include "BoundaryLayer.h"
16 #include "VoronoiVertex.h"
17 #include "VoronoiSegment.h"
18 #include "SSMFVoronoiSegment.h"
19 #include "SSSFVoronoiSegment.h"
20 #include "Border.h"
21 #include "Body.h"
22 #include "QuadLayer.h"
23 #include "TriangleNELayer.h"
24
25 #include "Node.h"
26 #include "Element.h"
27
28 #include <vector>
29
Mesh()30 Mesh::Mesh()
31 {
32 }
33
34 void Mesh::
convertEdgeFormat(MeshParser & parser)35 convertEdgeFormat( MeshParser& parser )
36 {
37 int i, len;
38
39 std::map<int, std::list<EdgeToken*> > nodeEdges;
40 std::map<int, std::list<LayerToken*> > nodeLayers;
41 std::map<int, int> edgeBody;
42
43 std::vector<Layer *> layers;
44
45 len = parser.nodes.size();
46 for( i = 0; i < len; ++i )
47 {
48 NodeToken *nt = parser.nodes[i];
49 GeometryNode *nd = new GeometryNode( nt->tag, nt->x, nt->y);
50 nd->boundaryTag = nt->boundaryTag;
51 // nd->setDelta( nt->delta );
52 geometryNodes[ nd->tag ] = nd;
53 }
54
55 len = parser.edges.size();
56 for( i = 0; i < len; ++i )
57 {
58 EdgeToken *et = parser.edges[i];
59
60 GeometryEdge *ed = new GeometryEdge( et->tag, et->segCount );
61 ed->setBoundaryTag(et->boundaryTag);
62
63 int sz = et->nodes.size();
64 for(int j = 0; j < sz; ++j)
65 {
66 nodeEdges[et->nodes[j]].push_back(et);
67
68 GeometryNodeMap::iterator it = geometryNodes.find(et->nodes[j]);
69 if (it == geometryNodes.end())
70 blm_error("Nonexisting node in edge", et->tag);
71
72 ed->addNode( it->second );
73 }
74
75 geometryEdges[ ed->tag ] = ed;
76 }
77
78 len = parser.bodies.size();
79 for( i = 0; i < len; ++i )
80 {
81 BodyToken *bt = parser.bodies[i];
82
83 Body *body = new Body( bt->tag, bt->parabolic );
84 int layerLen = bt->layers.size();
85 for( int m = 0; m < layerLen; ++m )
86 {
87 LayerToken *ly = bt->layers[m];
88
89 if (ly->delta <= 0.0)
90 ly->delta = bt->delta;
91
92 BGMesh *bgmesh = NULL;
93 Layer *layer;
94
95 // Here comes the selection depending on type
96 if( ly->type == QUAD_GRID )
97 layer = new QuadLayer( bt->tag );
98 else if( ly->type == TRIANGLE_NE_GRID )
99 layer = new TriangleNELayer( bt->tag );
100 else if( ly->type == TRIANGLE_NW_GRID )
101 layer = new TriangleNWLayer( bt->tag );
102 else if( ly->type == TRIANGLE_UJ_NE_GRID )
103 layer = new TriangleUJNELayer( bt->tag );
104 else if( ly->type == TRIANGLE_UJ_NW_GRID )
105 layer = new TriangleUJNWLayer( bt->tag );
106 else if( ly->type == TRIANGLE_FB_NE_GRID )
107 layer = new TriangleFBNELayer( bt->tag );
108 else if( ly->type == TRIANGLE_FB_NW_GRID )
109 layer = new TriangleFBNWLayer( bt->tag );
110 else
111 {
112 if( ly->type == CONNECT )
113 {
114 bgmesh = new BGTriangleMesh;
115 layer = new Connect( bt->tag, bgmesh );
116 }
117 else if ( ly->type == BOUNDARY_MESH )
118 {
119 bgmesh = new BGTriangleMesh;
120 layer = new BoundaryLayer( bt->tag, bgmesh );
121 }
122 else
123 {
124 if ( ly->bg->type == "Grid" )
125 bgmesh = new BGGridMesh(2.0);
126 else if ( ly->bg->type == "DenseGrid" )
127 bgmesh = new BGGridMesh(1.0);
128 else if ( ly->bg->type == "SparseGrid" )
129 bgmesh = new BGGridMesh(3.0);
130 else if ( ly->bg->type == "Delaunay" || ly->bg->type == "Explicit" )
131 bgmesh = new BGTriangleMesh;
132 else
133 blm_error("Unknown background mesh type:", ly->bg->type.c_str());
134
135 if( ly->bg->type == "Explicit" )
136 {
137 std::vector< GeometryNode * > bgnodes;
138 std::vector< GeometryEdge * > dummy;
139
140 for(int ind = 0; ind < ly->bg->nodes.size(); ++ind)
141 {
142 NodeToken nt = ly->bg->nodes[ind];
143 GeometryNode *nd = new GeometryNode( 0, nt.x, nt.y);
144 nd->setDelta( nt.delta * parser.scale );
145 bgnodes.push_back( nd );
146 }
147
148 bgmesh->initialize(bgnodes, dummy);
149 }
150
151 if( ly->type == VORONOI_VERTEX )
152 layer = new VoronoiVertex( bt->tag, bgmesh );
153 else if( ly->type == VORONOI_SEGMENT )
154 layer = new VoronoiSegment( bt->tag, bgmesh );
155 else if( ly->type == SSMF_VORONOI_SEGMENT )
156 layer = new SSMFVoronoiSegment( bt->tag, bgmesh );
157 else if( ly->type == SSSF_VORONOI_SEGMENT )
158 layer = new SSSFVoronoiSegment( bt->tag, bgmesh );
159 else
160 blm_error("Unknown layer type code!!!!");
161 }
162
163 // Insert the fixed nodes to the layer
164 for (int i = 0; i < ly->fixedNodes.size(); i++)
165 {
166 NodeToken *nt = ly->fixedNodes[i];
167 GeometryNode *nd = new GeometryNode( nt->tag, nt->x, nt->y);
168
169 if (nt->delta <= 0.0)
170 nt->delta = ly->delta > 0.0 ? ly->delta : parser.delta;
171
172 std::cout << "Node " << nt->tag << ": " << nt->delta * parser.scale << std::endl;
173 nd->setDelta( nt->delta * parser.scale );
174
175 layer->addFixedNode( nd );
176 }
177 }
178
179 int seedDirection = 1;
180
181 Border *bor = new Border( ly->loops.size() );
182 int loopLen = ly->loops.size();
183 for( int j = 0; j < loopLen; ++j )
184 {
185 int edgeLen = ly->loops[j]->edges.size();
186
187 for( int k = 0; k < edgeLen; k++ )
188 {
189 int edgeTag;
190 if (ly->loops[j]->direction > 0)
191 edgeTag = ly->loops[j]->edges[k];
192 else
193 edgeTag = -ly->loops[j]->edges[edgeLen-1-k];
194
195 int dir = edgeTag > 0 ? 1 : -1;
196 edgeTag = abs( edgeTag );
197
198 GeometryEdgeMap::iterator it = geometryEdges.find(edgeTag);
199 if (it == geometryEdges.end())
200 blm_error("Nonexisting edge in layer", ly->tag);
201 GeometryEdge *ed = it->second;
202
203 if (ly->gridh > 0 && ly->gridv > 0)
204 ed->setSegments(k & 1 ? ly->gridv : ly->gridh);
205
206 std::map<int,int>::iterator mi;
207 if ((mi = edgeBody.find(edgeTag)) == edgeBody.end())
208 edgeBody[edgeTag] = bt->tag;
209 else if (mi->second == bt->tag)
210 ed->makeVirtual();
211 else
212 ed->makeInner();
213
214 // Add background mesh to edge
215 if (bgmesh != NULL)
216 ed->addBGMesh(bgmesh);
217
218 // Push layer to node layer list, last node from next edge
219 std::vector<GeometryNode*> nds;
220 ed->exportGeometryNodes(nds, dir);
221 for (int n = 0; n < nds.size() - 1; n++)
222 nodeLayers[nds[n]->tag].push_back(ly);
223
224 bor->addLoopEdge( j, dir, ed );
225
226 if (ly->seed && ly->seed->edge == edgeTag)
227 seedDirection = dir;
228 }
229 }
230
231 // Set the seed
232 if (ly->type == SSMF_VORONOI_SEGMENT || ly->type == SSSF_VORONOI_SEGMENT)
233 {
234 SSVoronoiSegment *vs = static_cast<SSVoronoiSegment *>( layer );
235 if( ly->seed->type == "Explicit" )
236 {
237 vs->makeExplicitSeed( ly->seed->nodes[0],
238 ly->seed->nodes[1], ly->seed->nodes[2] );
239 }
240 else if( ly->seed->type == "Implicit" )
241 {
242 vs->makeImplicitSeed( geometryEdges[ ly->seed->edge ], seedDirection );
243 }
244 else
245 blm_error("Unknown seed type:", ly->seed->type.c_str());
246 }
247
248 layer->setBounds( bor );
249 body->addLayer( layer );
250 layers.push_back(layer);
251 }
252
253 bodies[ body->tag ] = body;
254 }
255
256 len = parser.nodes.size();
257 for (i = 0; i < len; i++)
258 {
259 NodeToken *nt = parser.nodes[i];
260
261 if (nt->delta <= 0.0)
262 {
263 double mean = 1.0;
264 int n = 0;
265
266 std::list<EdgeToken*> &edges = nodeEdges[nt->tag];
267 std::list<EdgeToken*>::iterator it;
268 for (it = edges.begin(); it != edges.end(); it++)
269 {
270 if ((*it)->delta > 0.0)
271 {
272 mean *= (*it)->delta;
273 n++;
274 }
275 }
276
277 if (n == 0)
278 {
279 std::list<LayerToken*> &layers = nodeLayers[nt->tag];
280 std::list<LayerToken*>::iterator it;
281 for (it = layers.begin(); it != layers.end(); it++)
282 {
283 if ((*it)->delta > 0.0)
284 {
285 mean *= (*it)->delta;
286 n++;
287 }
288 }
289 }
290
291 if (n > 0)
292 nt->delta = pow(mean, 1.0 / n);
293 else
294 nt->delta = parser.delta;
295 }
296
297 std::cout << "Node " << nt->tag << ": " << nt->delta * parser.scale << std::endl;
298 geometryNodes[nt->tag]->setDelta(nt->delta * parser.scale);
299 }
300
301 for (std::vector<Layer*>::iterator it = layers.begin(); it != layers.end(); it++)
302 {
303 (*it)->initialize();
304 }
305
306 std::cout << "Conversion completed." << std::endl;
307 }
308
309 void Mesh::
discretize()310 discretize()
311 {
312 GeometryNodeMapIt nit;
313 for( nit = geometryNodes.begin(); nit != geometryNodes.end(); ++nit )
314 {
315 MeshNode *mnd = new MeshNode( *((*nit).second) );
316 fixedNodes[ mnd->tag ] = mnd;
317 mnd->boundarynode = true;
318 }
319
320 GeometryEdgeMapIt eit;
321 for( eit = geometryEdges.begin(); eit != geometryEdges.end(); ++eit )
322 {
323 GeometryEdge *ed = (*eit).second;
324 if (ed->isConstant()) ed->discretize( fixedNodes );
325 }
326
327 for( eit = geometryEdges.begin(); eit != geometryEdges.end(); ++eit )
328 {
329 GeometryEdge *ed = (*eit).second;
330 if (!ed->isConstant()) ed->discretize( fixedNodes );
331 }
332
333 BodyMapIt bit;
334 for( bit = bodies.begin(); bit != bodies.end(); ++bit )
335 {
336 Body *bd = (*bit).second;
337 bd->discretize( fixedNodes, meshNodes, elements );
338 std::cout << "Body " << bd->tag << " completed!" << std::endl;
339 }
340
341 for( eit = geometryEdges.begin(); eit != geometryEdges.end(); ++eit )
342 {
343 GeometryEdge *ed = (*eit).second;
344 ed->elements(boundaryElements, 0);
345 }
346
347 createMiddleNodes();
348
349 std::cout << "Nodes: " << meshNodes.size() << std::endl;
350 std::cout << "Elements: " << elements.size() << std::endl;
351 }
352
353 void Mesh::
createMiddleNodes()354 createMiddleNodes()
355 {
356 std::map<std::pair<int, int>, Node *> edgeNodeMap;
357
358 std::list<Element *>::iterator i;
359 for (i = elements.begin(); i != elements.end(); i++)
360 {
361 Element *e = *i;
362 if (bodies[e->partOfBody()]->isParabolic())
363 {
364 int size = e->size();
365 std::vector<Node *> newNodes;
366 for (int i = 0; i < size; i++)
367 {
368 Node *node;
369
370 Node *node0 = e->nodeAt(i), *node1 = e->nodeAt((i+1)%size);
371 std::pair<int, int> p(std::min(node0->tag, node1->tag), std::max(node0->tag, node1->tag));
372 std::map<std::pair<int, int>, Node *>::iterator ni;
373 if ((ni = edgeNodeMap.find(p)) == edgeNodeMap.end())
374 {
375 double x = (node0->x + node1->x) / 2;
376 double y = (node0->y + node1->y) / 2;
377 node = new MeshNode(x, y);
378 edgeNodeMap[p] = node;
379 meshNodes[node->tag] = node;
380 }
381 else
382 {
383 node = ni->second;
384 }
385
386 newNodes.push_back(node);
387 }
388
389 #if 0
390 if (e->elementType() > 400)
391 {
392 double x = (e->nodeAt(0)->x + e->nodeAt(1)->x + e->nodeAt(2)->x + e->nodeAt(3)->x) / 4.0;
393 double y = (e->nodeAt(0)->y + e->nodeAt(1)->y + e->nodeAt(2)->y + e->nodeAt(3)->y) / 4.0;
394 Node *node = new MeshNode(x, y);
395 meshNodes[node->tag] = node;
396 newNodes.push_back(node);
397 }
398 #endif
399
400 e->upgrade(newNodes);
401 }
402 }
403
404 std::vector<BoundaryElement *>::iterator bi;
405 for (bi = boundaryElements.begin(); bi != boundaryElements.end(); bi++)
406 {
407 BoundaryElement *e = *bi;
408 int t0 = e->from()->tag, t1 = e->to()->tag;
409 std::pair<int, int> p(std::min(t0, t1), std::max(t0, t1));
410 std::map<std::pair<int, int>, Node *>::iterator ni;
411 if ((ni = edgeNodeMap.find(p)) != edgeNodeMap.end())
412 e->addMiddleNode(ni->second);
413 }
414 }
415
output(std::ofstream & o)416 void Mesh::output(std::ofstream& o)
417 {
418 NodeMapIt noit;
419 std::list< Element * >::iterator elit;
420 GeometryEdgeMapIt git;
421
422 int index = 1;
423
424 o << meshNodes.size() << ' ' << elements.size() << std::endl;
425 for( noit = meshNodes.begin(); noit != meshNodes.end(); ++noit)
426 {
427 Node *nd = (*noit).second;
428 nd->tag = index++;
429 o << *nd;
430 }
431
432 for( elit = elements.begin(); elit != elements.end(); ++elit)
433 {
434 o << **elit;
435 }
436 }
437
outputHeader(std::ofstream & o)438 void Mesh::outputHeader(std::ofstream& o)
439 {
440 int ngnds = 0;
441 GeometryNodeMap::iterator it;
442 for (it = geometryNodes.begin(); it != geometryNodes.end(); it++)
443 if (it->second->boundaryTag > 0)
444 ngnds++;
445
446 o << meshNodes.size() << ' ' << elements.size() << ' ' << boundaryElements.size() + ngnds << '\n';
447
448 std::map<int, int> n;
449
450 std::list< Element* >::const_iterator e;
451 for (e = elements.begin(); e != elements.end(); e++)
452 n[(*e)->elementType()]++;
453
454 std::vector< BoundaryElement* >::const_iterator be;
455 for (be = boundaryElements.begin(); be != boundaryElements.end(); be++)
456 n[(*be)->middle() != NULL ? 203 : 202]++;
457
458 n[101] += ngnds;
459
460 o << n.size() << '\n';
461
462 std::map<int, int>::iterator i;
463 for (i = n.begin(); i != n.end(); i++)
464 o << i->first << ' ' << i->second << '\n';
465 }
466
outputNodes(std::ofstream & o)467 void Mesh::outputNodes(std::ofstream &o)
468 {
469 int tag = 1;
470 NodeMapIt noit;
471
472 // o.setf(ios::scientific);
473 o.precision(16);
474
475 for( noit = meshNodes.begin(); noit != meshNodes.end(); ++noit)
476 {
477 Node *nd = (*noit).second;
478 nd->tag = tag++;
479 o << *nd;
480 }
481 }
482
outputElements(std::ofstream & o)483 void Mesh::outputElements(std::ofstream &o)
484 {
485 std::list< Element * >::iterator elit;
486
487 for( elit = elements.begin(); elit != elements.end(); ++elit)
488 {
489 o << **elit;
490 }
491 }
492
outputBoundary(std::ofstream & o)493 void Mesh::outputBoundary(std::ofstream &o)
494 {
495 std::vector< BoundaryElement* >::iterator bel;
496 for (bel = boundaryElements.begin(); bel != boundaryElements.end(); bel++)
497 {
498 o << **bel;
499 }
500
501 int id = boundaryElements.size() + 1;
502 GeometryNodeMapIt n;
503 for (n = geometryNodes.begin(); n != geometryNodes.end(); n++)
504 {
505 if (n->second->boundaryTag > 0)
506 o << id++ << ' ' << n->second->boundaryTag << " -1 -1 101 " << n->first << '\n';
507 }
508 }
509