OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
trianglemesherinterface.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
36 #include "floatarray.h"
37 #include "intarray.h"
38 #include "mathfem.h"
39 
40 #include <set>
41 #include <cstdio>
42 
43 #ifdef __TRIANGLE_MODULE
44 extern "C" {
45  #include <triangle.h>
46 
47  // Function prototypes from triangle.c
48  void parsecommandline(int argc, char **argv, struct behavior *b);
49  void transfernodes(struct mesh *m, struct behavior *b, REAL *pointlist, REAL *pointattriblist, int *pointmarkerlist, int numberofpoints, int numberofpointattribs);
50  void formskeleton(struct mesh *m, struct behavior *b, int *segmentlist, int *segmentmarkerlist, int numberofsegments);
51  void highorder(struct mesh *m, struct behavior *b);
52  void numbernodes(struct mesh *m, struct behavior *b);
53  void writeelements(struct mesh *m, struct behavior *b, int **trianglelist, REAL **triangleattriblist);
54  void writepoly(struct mesh *m, struct behavior *b, int **segmentlist, int **segmentmarkerlist);
55  void writeedges(struct mesh *m, struct behavior *b, int **edgelist, int **edgemarkerlist);
56  void statistics(struct mesh *m, struct behavior *b);
57  void writenodes(struct mesh *m, struct behavior *b, REAL **pointlist, REAL **pointattriblist, int **pointmarkerlist);
58  void checkmesh(struct mesh *m, struct behavior *b);
59  void checkdelaunay(struct mesh *m, struct behavior *b);
60  void writeneighbors(struct mesh *m, struct behavior *b, int **neighborlist);
61  void writevoronoi(struct mesh *m, struct behavior *b, REAL **vpointlist, REAL **vpointattriblist, int **vpointmarkerlist, int **vedgelist, int **vedgemarkerlist, REAL **vnormlist);
62  int reconstruct(struct mesh *m, struct behavior *b, int *trianglelist, REAL *triangleattriblist, REAL *trianglearealist, int elements, int corners, int attribs, int *segmentlist, int *segmentmarkerlist, int numberofsegments);
63  subseg *subsegtraverse(struct mesh *m);
64  void regionplague(struct mesh *m, struct behavior *b, REAL attribute, REAL area);
65  void plague(struct mesh *m, struct behavior *b);
66  void poolinit(struct memorypool *pool, int bytecount, int itemcount, int firstitemcount, int alignment);
67  void pooldeinit(struct memorypool *pool);
68 
69  // Macros from triangle.c
70  #define setelemattribute(otri, attnum, value) \
71  ( ( REAL * ) ( otri ).tri ) [ m->elemattribindex + ( attnum ) ] = value
72  #define elemattribute(otri, attnum) \
73  ( ( REAL * ) ( otri ).tri ) [ m->elemattribindex + ( attnum ) ]
74 
75  #define stpivot(osub, otri) \
76  ptr = ( triangle ) ( osub ).ss [ 6 + ( osub ).ssorient ]; \
77  decode(ptr, otri)
78  #define ssymself(osub) \
79  ( osub ).ssorient = 1 - ( osub ).ssorient
80  #define mark(osub) ( * ( int * ) ( ( osub ).ss + 8 ) )
81  #define decode(ptr, otri) \
82  ( otri ).orient = ( int ) ( ( long ) ptr & ( long ) 3l ); \
83  ( otri ).tri = ( triangle * ) ( ( long ) ( ptr ) ^ ( long ) ( otri ).orient )
84  #define infect(otri) \
85  ( otri ).tri [ 6 ] = ( triangle ) ( ( long ) ( otri ).tri [ 6 ] | ( long ) 2l )
86  #define org(otri, vertexptr) \
87  vertexptr = ( vertex ) ( otri ).tri [ plus1mod3 [ ( otri ).orient ] + 3 ]
88  #define dest(otri, vertexptr) \
89  vertexptr = ( vertex ) ( otri ).tri [ minus1mod3 [ ( otri ).orient ] + 3 ]
90  #define apex(otri, vertexptr) \
91  vertexptr = ( vertex ) ( otri ).tri [ ( otri ).orient + 3 ]
92  #define sorg(osub, vertexptr) \
93  vertexptr = ( vertex ) ( osub ).ss [ 2 + ( osub ).ssorient ]
94  #define sdest(osub, vertexptr) \
95  vertexptr = ( vertex ) ( osub ).ss [ 3 - ( osub ).ssorient ]
96 
97 
98  #define VIRUSPERBLOCK 1020 /* Number of virus triangles allocated at once. */
99 
100  // This customized version carves holes defined by the regions.
101  // It'll remove any element which doesn't receive any region.
102  int custom_carveholes(struct mesh *m, struct behavior *b, const int *outside, const int *inside)
103  {
104  struct otri neighbortri;
105  struct otri triangleloop;
106  struct osub subsegloop;
107  triangle **temptri;
108  triangle ptr; /* Temporary variable used by sym(). */
109  double area;
110  double attribute;
111  double triangle_attribute;
112  int sane_mesh;
113 
114  poolinit(& m->viri, sizeof( triangle * ), VIRUSPERBLOCK, VIRUSPERBLOCK, 0);
115 
116  /* Assigns every triangle a regional attribute of -1 */
117  traversalinit(& m->triangles);
118  triangleloop.orient = 0;
119  triangleloop.tri = triangletraverse(m);
120  while ( triangleloop.tri != ( triangle * ) NULL ) {
121  setelemattribute(triangleloop, m->eextras, -1.0);
122  triangleloop.tri = triangletraverse(m);
123  }
124 
125  sane_mesh = 1;
126 
127  /* Loop over all segments */
128  traversalinit(& m->subsegs);
129  subsegloop.ss = subsegtraverse(m);
130  subsegloop.ssorient = 0;
131  while ( subsegloop.ss != ( subseg * ) NULL ) {
132  if ( subsegloop.ss != m->dummysub ) {
133  /* First neighbor. */
134  ssymself(subsegloop);
135  stpivot(subsegloop, neighbortri);
136  if ( neighbortri.tri != m->dummytri ) {
137  area = 0.0;
138  attribute = outside [ mark(subsegloop) - 1 ];
139  triangle_attribute = elemattribute(neighbortri, m->eextras);
140 
141  /* The region no number yet. */
142  if ( triangle_attribute < 0 && attribute >= 0 ) {
143  infect(neighbortri);
144  temptri = ( triangle ** ) poolalloc(& m->viri);
145  * temptri = neighbortri.tri;
146  /* Apply one region's attribute and/or area constraint. */
147  regionplague(m, b, attribute, area);
148  /* The virus pool should be empty now. */
149  } else if ( attribute >= 0 && triangle_attribute >= 0 && attribute != triangle_attribute ) {
150  /* Check for problems. */
151  vertex v1, v2, v3;
152  org(neighbortri, v1);
153  dest(neighbortri, v2);
154  apex(neighbortri, v3);
155  fprintf(stdout, "Error: inconsistent region information (new %d, old %d from %d) (outside) at (%e, %e)\n",
156  ( int ) attribute, ( int ) triangle_attribute, ( int ) mark(subsegloop), ( v1 [ 0 ] + v2 [ 0 ] + v3 [ 0 ] ) / 3, ( v1 [ 1 ] + v2 [ 1 ] + v3 [ 1 ] ) / 3);
157  sane_mesh = 0;
158  }
159  }
160 
161  /* Second neighbor (same procedure). */
162  ssymself(subsegloop);
163  stpivot(subsegloop, neighbortri);
164  if ( neighbortri.tri != m->dummytri ) {
165  area = 0.0;
166  attribute = inside [ mark(subsegloop) - 1 ];
167  triangle_attribute = elemattribute(neighbortri, m->eextras);
168 
169  if ( triangle_attribute < 0 && attribute >= 0 ) {
170  infect(neighbortri);
171  temptri = ( triangle ** ) poolalloc(& m->viri);
172  * temptri = neighbortri.tri;
173  regionplague(m, b, attribute, area);
174  } else if ( attribute >= 0 && triangle_attribute >= 0 && attribute != triangle_attribute ) {
175  vertex v1, v2, v3;
176  org(neighbortri, v1);
177  dest(neighbortri, v2);
178  apex(neighbortri, v3);
179  fprintf(stdout, "Error: inconsistent region information (new %d, old %d from %d) (inside) at (%e, %e)\n",
180  ( int ) attribute, ( int ) triangle_attribute, ( int ) mark(subsegloop), ( v1 [ 0 ] + v2 [ 0 ] + v3 [ 0 ] ) / 3, ( v1 [ 1 ] + v2 [ 1 ] + v3 [ 1 ] ) / 3);
181  sane_mesh = 0;
182  }
183  }
184 
185  subsegloop.ss = subsegtraverse(m);
186  }
187  }
188 
189  /* Remove all triangles with marker 0.0 */
190  traversalinit(& m->triangles);
191  triangleloop.tri = triangletraverse(m);
192  int triangle_number = 0;
193  while ( triangleloop.tri != ( triangle * ) NULL ) {
194  if ( triangleloop.tri != m->dummytri ) {
195  triangle_number++;
196  attribute = elemattribute(triangleloop, m->eextras);
197  if ( attribute == -1.0 ) {
198  fprintf(stderr, "Broken mesh at triangle %d\n", triangle_number);
199  sane_mesh = 0;
200  } else if ( attribute == 0.0 ) {
201  infect(triangleloop);
202  temptri = ( triangle ** ) poolalloc(& m->viri);
203  * temptri = triangleloop.tri;
204  }
205  }
206  triangleloop.tri = triangletraverse(m);
207  }
208  /* Remove the marked elements */
209  plague(m, b);
210 
211  if ( b->regionattrib && !b->refine ) {
212  /* Note the fact that each triangle has an additional attribute. */
213  m->eextras++;
214  }
215 
216  /* Free up memory. */
217  pooldeinit(& m->viri);
218 
219  return sane_mesh;
220  }
221 
222  // Customized triangulation call.
223  void custom_triangulate(char *triswitches, struct triangulateio *in, struct triangulateio *out, struct triangulateio *vorout, const int *outside, const int *inside)
224  {
225  struct mesh m;
226  struct behavior b;
227  //REAL *holearray;
228  //REAL *regionarray;
229  triangleinit(& m);
230  parsecommandline(1, & triswitches, & b);
231  //b.verbose=2;
232  m.steinerleft = b.steiner;
233  transfernodes(& m, & b, in->pointlist, in->pointattributelist,
234  in->pointmarkerlist, in->numberofpoints,
235  in->numberofpointattributes);
236  if ( b.refine ) {
237  m.hullsize = reconstruct(& m, & b, in->trianglelist,
238  in->triangleattributelist, in->trianglearealist,
239  in->numberoftriangles, in->numberofcorners,
240  in->numberoftriangleattributes,
241  in->segmentlist, in->segmentmarkerlist,
242  in->numberofsegments);
243  } else {
244  m.hullsize = delaunay(& m, & b);
245  }
246  m.infvertex1 = ( vertex ) NULL;
247  m.infvertex2 = ( vertex ) NULL;
248  m.infvertex3 = ( vertex ) NULL;
249  if ( b.usesegments ) {
250  m.checksegments = 1;
251  if ( !b.refine ) {
252  formskeleton(& m, & b, in->segmentlist,
253  in->segmentmarkerlist, in->numberofsegments);
254  }
255  }
256  #if 0
257  struct osub subsegloop;
258  traversalinit(& m.subsegs);
259  subsegloop.ss = subsegtraverse(& m);
260  subsegloop.ssorient = 0;
261  while ( subsegloop.ss != ( subseg * ) NULL ) {
262  if ( subsegloop.ss != m.dummysub ) {
263  REAL *p1, *p2;
264  sorg(subsegloop, p1);
265  sdest(subsegloop, p2);
266  printf(" Connected (%f,%f) to (%f,%f)\n", p1 [ 0 ], p1 [ 1 ], p2 [ 0 ], p2 [ 1 ]);
267  subsegloop.ss = subsegtraverse(& m);
268  }
269  }
270  #endif
271  if ( b.poly && ( m.triangles.items > 0 ) ) {
272  //holearray = in->holelist;
273  m.holes = in->numberofholes;
274  //regionarray = in->regionlist;
275  m.regions = in->numberofregions;
276  if ( !b.refine ) {
277  /* Only increase quality if the regions are properly defined. */
278  int sane = custom_carveholes(& m, & b, outside, inside);
279  b.quality *= sane;
280  if ( sane == 0 ) {
281  printf("Probably bad PSLG\n");
282  exit(-1);
283  }
284  }
285  } else {
286  m.holes = 0;
287  m.regions = 0;
288  }
289  if ( b.quality && ( m.triangles.items > 0 ) ) {
290  enforcequality(& m, & b);
291  }
292  m.edges = ( 3l * m.triangles.items + m.hullsize ) / 2l;
293  if ( b.order > 1 ) {
294  highorder(& m, & b);
295  }
296  if ( !b.quiet ) {
297  printf("\n");
298  }
299  if ( b.jettison ) {
300  out->numberofpoints = m.vertices.items - m.undeads;
301  } else {
302  out->numberofpoints = m.vertices.items;
303  }
304  out->numberofpointattributes = m.nextras;
305  out->numberoftriangles = m.triangles.items;
306  out->numberofcorners = ( b.order + 1 ) * ( b.order + 2 ) / 2;
307  out->numberoftriangleattributes = m.eextras;
308  out->numberofedges = m.edges;
309  if ( b.usesegments ) {
310  out->numberofsegments = m.subsegs.items;
311  } else {
312  out->numberofsegments = m.hullsize;
313  }
314  if ( vorout != ( struct triangulateio * ) NULL ) {
315  vorout->numberofpoints = m.triangles.items;
316  vorout->numberofpointattributes = m.nextras;
317  vorout->numberofedges = m.edges;
318  }
319  if ( b.nonodewritten || ( b.noiterationnum && m.readnodefile ) ) {
320  if ( !b.quiet ) {
321  printf("NOT writing vertices.\n");
322  }
323  numbernodes(& m, & b);
324  } else {
325  writenodes(& m, & b, & out->pointlist, & out->pointattributelist,
326  & out->pointmarkerlist);
327  }
328 
329  // Simp. always write the triangles.
330  writeelements(& m, & b, & out->trianglelist, & out->triangleattributelist);
331 
332  if ( b.poly || b.convex ) {
333  writepoly(& m, & b, & out->segmentlist, & out->segmentmarkerlist);
334  out->numberofholes = m.holes;
335  out->numberofregions = m.regions;
336  if ( b.poly ) {
337  out->holelist = in->holelist;
338  out->regionlist = in->regionlist;
339  } else {
340  out->holelist = ( REAL * ) NULL;
341  out->regionlist = ( REAL * ) NULL;
342  }
343  }
344  if ( b.edgesout ) {
345  writeedges(& m, & b, & out->edgelist, & out->edgemarkerlist);
346  }
347  // Simp. no voronoi
348  if ( b.neighbors ) {
349  writeneighbors(& m, & b, & out->neighborlist);
350  }
351  // Simp. No statistics.
352  if ( b.docheck ) {
353  checkmesh(& m, & b);
354  checkdelaunay(& m, & b);
355  }
356  triangledeinit(& m, & b);
357  }
358 }
359 
360 // Convenience function.
361 void clearTriangulateIO(struct triangulateio &t)
362 {
363  t.pointlist = NULL;
364  t.pointattributelist = NULL;
365  t.pointmarkerlist = NULL;
366  t.numberofpoints = 0;
367  t.numberofpointattributes = 0;
368 
369  t.trianglelist = NULL;
370  t.triangleattributelist = NULL;
371  t.trianglearealist = NULL;
372  t.numberoftriangles = 0;
373  t.numberofcorners = 0;
374  t.numberoftriangleattributes = 0;
375 
376  t.segmentlist = NULL;
377  t.segmentmarkerlist = NULL;
378  t.numberofsegments = 0;
379 
380  t.holelist = NULL;
381  t.numberofholes = 0;
382 
383  t.regionlist = NULL;
384  t.numberofregions = 0;
385 }
386 #endif
387 
388 namespace oofem {
390  const IntArray &outside, const IntArray &inside,
391  std :: vector< FloatArray > &nodes, std :: vector< IntArray > &n_markers,
392  std :: vector< IntArray > &triangles, IntArray &t_markers,
393  std :: vector< IntArray > &segments, IntArray &s_markers) const
394 {
395 #ifdef __TRIANGLE_MODULE
396  // 1. Fill the struct for triangle;
397  struct triangulateio mesh;
398  clearTriangulateIO(mesh);
399 
400  // 1.a. Copy over the node data.
401  mesh.numberofpoints = pslg.nx.giveSize();
402  mesh.pointlist = new REAL [ mesh.numberofpoints * 2 ];
403  //mesh.pointmarkerlist = new REAL[mesh.numberofpoints];
404  for ( int i = 0; i < mesh.numberofpoints; ++i ) {
405  mesh.pointlist [ i * 2 ] = pslg.nx(i);
406  mesh.pointlist [ i * 2 + 1 ] = pslg.ny(i);
407  //mesh.pointmarkerlist[i] = pslg.n_marker(i);
408  }
409 
410  // 1.b. Copy over the segment data
411  printf("Copying segment data\n");
412  mesh.numberofsegments = pslg.segment_a.giveSize();
413  mesh.segmentlist = new int [ mesh.numberofsegments * 2 ];
414  for ( int i = 0; i < mesh.numberofsegments; ++i ) {
415  mesh.segmentlist [ i * 2 ] = pslg.segment_a(i);
416  mesh.segmentlist [ i * 2 + 1 ] = pslg.segment_b(i);
417  }
418 
419  if ( pslg.segment_marker.giveSize() > 0 ) {
420  mesh.segmentmarkerlist = new int [ mesh.numberofsegments ];
421  for ( int i = 0; i < mesh.numberofsegments; ++i ) {
422  mesh.segmentmarkerlist [ i ] = pslg.segment_marker(i);
423  }
424  }
425 
426  // 2. Triangulate
427  char options [ 100 ];
428  // Note: Not sure if -A is necessary when using the library interface.
429  sprintf(options, "-p -q %f -a%f %s -A", this->minAngle, this->maxArea, this->quadratic ? "-o2" : "");
430  struct triangulateio output;
431  clearTriangulateIO(output);
432  custom_triangulate( options, & mesh, & output, NULL, outside.givePointer(), inside.givePointer() );
433 
434  // 3. Copy back
435  nodes.resize(output.numberofpoints);
436  //n_markers.resize(output.numberofpoints);
437  for ( int i = 0; i < output.numberofpoints; ++i ) {
438  nodes [ i ].resize(2);
439  nodes [ i ].at(1) = output.pointlist [ i * 2 ];
440  nodes [ i ].at(2) = output.pointlist [ i * 2 + 1 ];
441  //n_markers(i) = output.pointmarkerlist[i]; // Not enough.
442  }
443 
444  triangles.resize(output.numberoftriangles);
445  t_markers.resize(output.numberoftriangles);
446  for ( int i = 0; i < output.numberoftriangles; ++i ) {
447  IntArray &triangle = triangles [ i ];
448  triangle.resize(output.numberofcorners);
449  for ( int j = 0; j < 3; j++ ) { // First three
450  triangle(j) = output.trianglelist [ i * output.numberofcorners + j ];
451  }
452  // Rearrange the strange ordering of the edge nodes.
453  if ( output.numberofcorners == 6 ) {
454  triangle(3) = output.trianglelist [ i * output.numberofcorners + 5 ];
455  triangle(4) = output.trianglelist [ i * output.numberofcorners + 3 ];
456  triangle(5) = output.trianglelist [ i * output.numberofcorners + 4 ];
457  }
458  t_markers.at(i + 1) = ( int ) round(output.triangleattributelist [ i ]);
459  }
460 
461  // A somewhat annoying missing feature of triangle, it won't make the segments quadratic.
462  std :: set< std :: size_t > *node_triangle = NULL;
463  if ( this->quadratic ) {
464  node_triangle = new std :: set< std :: size_t > [ output.numberofpoints ];
465  for ( std :: size_t i = 0; i < triangles.size(); ++i ) {
466  IntArray &triangle = triangles [ i ];
467  for ( int j = 1; j <= 3; ++j ) {
468  node_triangle [ triangle.at(j) - 1 ].insert(i);
469  }
470  }
471  }
472 
473  segments.resize(output.numberofsegments);
474  s_markers.resize(output.numberofsegments);
475  for ( int i = 0; i < output.numberofsegments; ++i ) {
476  IntArray &segment = segments [ i ];
477  segment.resize(this->quadratic ? 3 : 2);
478  segment.at(1) = output.segmentlist [ i * 2 + 0 ];
479  segment.at(2) = output.segmentlist [ i * 2 + 1 ];
480  //segment->at(3) = output.segmentlist[i*3 + 2]; // Quadratic meshes only, not for segments.
481  if ( this->quadratic ) {
482  int a, b, c;
483  std :: set< std :: size_t >tris = node_triangle [ segment.at(1) - 1 ];
484  // Now look up any triangle with the other point included.
485  for ( auto tri: tris ) {
486  IntArray &triangle = triangles [ tri ];
487  if ( ( b = triangle.findFirstIndexOf( segment.at(2) ) ) > 0 ) {
488  a = triangle.findFirstIndexOf( segment.at(1) );
489  if ( a + b == 3 ) {
490  c = 4;
491  } else if ( a + b == 5 ) {
492  c = 5;
493  } else {
494  c = 6;
495  }
496  segment.at(3) = triangle.at(c);
497  break;
498  }
499  }
500  }
501  s_markers.at(i + 1) = output.segmentmarkerlist [ i ];
502  }
503  if ( this->quadratic ) {
504  delete[] node_triangle;
505  }
506 
507  this->fixNodeMarkers(nodes, n_markers, triangles, t_markers, segments, s_markers);
508 
509  // Deleting old memory
510  delete[] mesh.pointlist;
511  //delete[] mesh.pointattributelist;
512  //delete[] mesh.pointmarkerlist;
513  //delete[] mesh.trianglelist;
514  //delete[] mesh.triangleattributelist;
515  //delete[] mesh.trianglearealist;
516  delete[] mesh.segmentlist;
517  delete[] mesh.segmentmarkerlist;
518  //if (mesh.holelist) { delete[] mesh.holelist; }
519  //if (mesh.regionlist) { delete[] mesh.regionlist; }
520 
521  // Holes and regions are referenced in both.
522  free(output.pointlist);
523  //free(output.pointattributelist);
524  free(output.pointmarkerlist);
525  free(output.trianglelist);
526  free(output.triangleattributelist);
527  //free(output.trianglearealist);
528  free(output.segmentlist);
529  free(output.segmentmarkerlist);
530  //free(output.holelist);
531  //free(output.regionlist);
532 #else
533  OOFEM_ERROR("OOFEM is not compiled with support for triangle.");
534 #endif
535 
536  return true;
537 }
538 
539 void TriangleMesherInterface :: fixNodeMarkers(const std :: vector< FloatArray > &nodes, std :: vector< IntArray > &n_markers,
540  const std :: vector< IntArray > &triangles, const IntArray &t_markers,
541  const std :: vector< IntArray > &segments, const IntArray &s_markers)
542 {
543  n_markers.resize( nodes.size() );
544 
545  for ( std :: size_t i = 0; i < segments.size(); ++i ) {
546  for ( int j = 1; j <= segments [ i ].giveSize(); ++j ) {
547  n_markers [ segments [ i ].at(j) - 1 ].insertSortedOnce( s_markers(i) );
548  }
549  }
550  for ( std :: size_t i = 0; i < triangles.size(); ++i ) {
551  for ( int j = 1; j <= triangles [ i ].giveSize(); ++j ) {
552  n_markers [ triangles [ i ].at(j) - 1 ].insertSortedOnce( t_markers(i) );
553  }
554  }
555 }
556 
557 void TriangleMesherInterface :: simplifyPSLG(Triangle_PSLG &coarse, const Triangle_PSLG &pslg, double limit, double minlen)
558 {
559  int segments = pslg.segment_a.giveSize();
560  int nodes = pslg.nx.giveSize();
561 
562  // Calculate the inverted connection node->element
563  std :: set< int > *connectivity = new std :: set< int > [ nodes ];
564  for ( int i = 0; i < segments; i++ ) {
565  connectivity [ pslg.segment_a(i) - 1 ].insert(i + 1);
566  connectivity [ pslg.segment_b(i) - 1 ].insert(i + 1);
567  }
568 
569  // Some conservative error measure
570  IntArray nodeRemoval(nodes);
571  IntArray edgeRemoval(segments);
572  edgeRemoval.zero();
573  IntArray seg_a = pslg.segment_a;
574  IntArray seg_b = pslg.segment_b;
575 
576  nodeRemoval.zero();
577  FloatArray error_x(nodes), error_y(nodes), ab(2), ac(2), bc(2), err_vec(2);
578  error_x.zero();
579  error_y.zero();
580 
581  for ( int j = 0; j < 2; ++j ) {
582  double allowed = j * limit / 2;
583  for ( int i = 1; i <= nodes; i++ ) {
584  std :: set< int > &elems = connectivity [ i - 1 ];
585 #if 1
586  if ( elems.size() < 2 ) {
587  // Lone nodes or edges are removed. It is up for discussion whether or not this is a wanted.
588  nodeRemoval.at(i) = true;
589  if ( elems.size() == 1 ) {
590  int e0 = * elems.begin();
591  edgeRemoval.at(e0) = true;
592  // Find the other node and remove segment e0 from it;
593  int n1 = seg_a.at(e0) == i ? seg_b.at(e0) : seg_a.at(e0);
594  connectivity [ n1 - 1 ].erase(e0);
595  elems.clear();
596  /*printf("Removing edge %d and node = %d\n", e0, n1);
597  * while (connectivity[n1-1].size() < 2) {
598  * nodeRemoval.at(n1) = true;
599  * if (connectivity[n1-1].size() == 1) {
600  * e0 = *connectivity[n1-1].begin();
601  * edgeRemoval.at(e0) = true;
602  * connectivity[n1-1].clear();
603  * n1 = seg_a.at(e0) == i ? seg_b.at(e0) : seg_a.at(e0);
604  * connectivity[n1-1].erase(e0);
605  * } else {
606  * OOFEM_WARNING("HAPPENS ONCE!");
607  * break;
608  * }
609  * }*/
610  }
611  } else
612 #endif
613  if ( elems.size() == 2 ) {
614  int e0 = * elems.begin();
615  int e1 = * ( ++elems.begin() );
616  if ( pslg.segment_marker.at(e0) == pslg.segment_marker.at(e1) ) {
617  double abac, acac;
618  int n0, n1, n2;
619 
620  n0 = seg_a.at(e0) == i ? seg_b.at(e0) : seg_a.at(e0);
621  n1 = i;
622  n2 = seg_a.at(e1) == i ? seg_b.at(e1) : seg_a.at(e1);
623 
624  ab(0) = pslg.nx.at(n1) - pslg.nx.at(n0);
625  ab(1) = pslg.ny.at(n1) - pslg.ny.at(n0);
626  ac(0) = pslg.nx.at(n2) - pslg.nx.at(n0);
627  ac(1) = pslg.ny.at(n2) - pslg.ny.at(n0);
628 
629  abac = ab.dotProduct(ac);
630  acac = ac.computeSquaredNorm();
631 
632  // Find the error (how far the point would be moved to obtain the following line without it).
633  if ( abac <= 0 || acac == 0 ) { // then -ab
634  err_vec(0) = -ab(0);
635  err_vec(1) = -ab(1);
636  } else if ( abac >= acac ) { // then bc
637  err_vec(0) = pslg.nx.at(n2) - pslg.nx.at(n1);
638  err_vec(1) = pslg.ny.at(n2) - pslg.ny.at(n1);
639  } else {
640  err_vec = ac;
641  err_vec.times(abac / acac);
642  err_vec.subtract(ab);
643  }
644 
645  double ev_norm = err_vec.computeNorm();
646  // Max of new or old error;
647  double real_error = 0.0;
648  if ( ev_norm == 0 ) {
649  real_error = 0.0;
650  } else {
651  error_x.at(n1) += err_vec(0);
652  error_y.at(n1) += err_vec(1);
653  real_error = sqrt( error_x.at(n1) * error_x.at(n1) + error_y.at(n1) * error_y.at(n1) );
654  }
655 
656  if ( real_error <= allowed ) {
657  // Mark node for removal, remove second edge, more first edge connection to next node;
658  nodeRemoval.at(i) = true;
659  edgeRemoval.at(e1) = true;
660  connectivity [ n2 - 1 ].erase(e1);
661  connectivity [ n2 - 1 ].insert(e0);
662  if ( seg_a.at(e0) == n1 ) { // Doing the bothersome way to preserve direction of segments.
663  seg_a.at(e0) = n2;
664  } else {
665  seg_b.at(e0) = n2;
666  }
667  elems.clear();
668  // Accumulate the error vector
669  error_x.at(n0) += error_x.at(n1);
670  error_y.at(n0) += error_y.at(n1);
671  error_x.at(n2) += error_x.at(n1);
672  error_y.at(n2) += error_y.at(n1);
673  }
674  }
675  }
676  }
677  }
678 
679  // Deleting the elements which are too short.
680  bool edgeRemoved = true;
681  while ( edgeRemoved ) {
682  edgeRemoved = false;
683  for ( int i = 1; i <= nodes; i++ ) {
684  std :: set< int > &elems = connectivity [ i - 1 ];
685  if ( elems.size() == 2 ) {
686  int e0 = * elems.begin();
687  int e1 = * ( ++elems.begin() );
688  if ( pslg.segment_marker.at(e0) == pslg.segment_marker.at(e1) ) {
689  int n0, n1, n2;
690 
691  n0 = seg_a.at(e0) == i ? seg_b.at(e0) : seg_a.at(e0);
692  n1 = i;
693  n2 = seg_a.at(e1) == i ? seg_b.at(e1) : seg_a.at(e1);
694 
695  ab(0) = pslg.nx.at(n1) - pslg.nx.at(n0);
696  ab(1) = pslg.ny.at(n1) - pslg.ny.at(n0);
697  bc(0) = pslg.nx.at(n2) - pslg.nx.at(n1);
698  bc(1) = pslg.ny.at(n2) - pslg.ny.at(n1);
699 
700  if ( ab.computeSquaredNorm() < minlen * minlen || bc.computeSquaredNorm() < minlen * minlen ) {
701  // Mark node for removal, remove second edge, more first edge connection to next node;
702  nodeRemoval.at(i) = true;
703  edgeRemoval.at(e1) = true;
704  connectivity [ n2 - 1 ].erase(e1);
705  connectivity [ n2 - 1 ].insert(e0);
706  if ( seg_a.at(e0) == n1 ) { // Doing the bothersome way to preserve direction of segments.
707  seg_a.at(e0) = n2;
708  } else {
709  seg_b.at(e0) = n2;
710  }
711  elems.clear();
712  edgeRemoved = true;
713  }
714  }
715  }
716  }
717  }
718 
719 
720  // Remove simple double-connections like o<--->o
721  for ( int i = 1; i <= nodes; i++ ) {
722  std :: set< int > &elems = connectivity [ i - 1 ];
723  if ( elems.size() == 2 ) {
724  int e0 = * elems.begin();
725  int e1 = * ( ++elems.begin() );
726  int n0 = seg_a.at(e0) == i ? seg_b.at(e0) : seg_a.at(e0);
727  int n1 = i;
728  int n2 = seg_a.at(e1) == i ? seg_b.at(e1) : seg_a.at(e1);
729 
730  if ( n0 == n2 ) {
731  // Then we have a simple double connection. We get rid of this node and the two edges.
732  nodeRemoval.at(n1) = true;
733  edgeRemoval.at(e0) = true;
734  edgeRemoval.at(e1) = true;
735  connectivity [ n1 - 1 ].erase(e0);
736  connectivity [ n1 - 1 ].erase(e1);
737  connectivity [ n2 - 1 ].erase(e0);
738  connectivity [ n2 - 1 ].erase(e1);
739  if ( connectivity [ n1 - 1 ].size() == 0 ) {
740  nodeRemoval.at(n2) = true;
741  }
742  }
743  }
744  }
745 
746  delete[] connectivity;
747 
748  // Cleanup
749  int newNodes = 0;
750  IntArray newNumber(nodes);
751  newNumber.zero();
752  coarse.nx.resize(nodes);
753  coarse.ny.resize(nodes);
754  for ( int i = 1; i <= nodes; i++ ) {
755  if ( !nodeRemoval.at(i) ) {
756  newNodes++;
757  coarse.nx.at(newNodes) = pslg.nx.at(i);
758  coarse.ny.at(newNodes) = pslg.ny.at(i);
759  newNumber.at(i) = newNodes;
760  }
761  }
762  coarse.nx.resize(newNodes);
763  coarse.ny.resize(newNodes);
764 
765  int newSegments = 0;
766  coarse.segment_a.resize(segments);
767  coarse.segment_b.resize(segments);
768  coarse.segment_marker.resize(segments);
769 
770  for ( int i = 1; i <= segments; i++ ) {
771  if ( !edgeRemoval.at(i) ) {
772  newSegments++;
773  coarse.segment_a.at(newSegments) = newNumber.at( seg_a.at(i) );
774  coarse.segment_b.at(newSegments) = newNumber.at( seg_b.at(i) );
775  coarse.segment_marker.at(newSegments) = pslg.segment_marker.at(i);
776  }
777  }
778  coarse.segment_a.resize(newSegments);
779  coarse.segment_b.resize(newSegments);
780  coarse.segment_marker.resize(newSegments);
781 }
782 };
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
IntArray segment_b
Second segment connection.
void zero()
Sets all component to zero.
Definition: intarray.C:52
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
static void simplifyPSLG(Triangle_PSLG &coarse, const Triangle_PSLG &pslg, double limit, double minlen=0.0)
Simplifies a PSLG while respecting topology, running in linear time.
bool meshPSLG(const Triangle_PSLG &pslg, const IntArray &outside, const IntArray &inside, std::vector< FloatArray > &nodes, std::vector< IntArray > &n_markers, std::vector< IntArray > &triangles, IntArray &t_markers, std::vector< IntArray > &segments, IntArray &s_markers) const
Constructs a mesh from a PSLG.
const int * givePointer() const
Breaks encapsulation.
Definition: intarray.h:336
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
FloatArray ny
Nodes y coordinates.
static void fixNodeMarkers(const std::vector< FloatArray > &nodes, std::vector< IntArray > &n_markers, const std::vector< IntArray > &triangles, const IntArray &t_markers, const std::vector< IntArray > &segments, const IntArray &s_markers)
Adds all neighboring regions to every node region.
#define OOFEM_ERROR(...)
Definition: error.h:61
Plane straight line graph used as input for meshing with triangle.
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
Class representing vector of real numbers.
Definition: floatarray.h:82
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
int giveSize() const
Definition: intarray.h:203
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
FloatArray nx
Nodes x coordinates.
the oofem namespace is to define a context or scope in which all oofem names are defined.
IntArray segment_a
First segment connection.
int findFirstIndexOf(int value) const
Finds index of first occurrence of given value in array.
Definition: intarray.C:331
IntArray segment_marker
Segment markers.
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:32 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011