1 /*****************************************************************************
2 * *
3 * Elmer, A Finite Element Software for Multiphysical Problems *
4 * *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland *
6 * *
7 * This program is free software; you can redistribute it and/or *
8 * modify it under the terms of the GNU General Public License *
9 * as published by the Free Software Foundation; either version 2 *
10 * of the License, or (at your option) any later version. *
11 * *
12 * This program is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with this program (in file fem/GPL-2); if not, write to the *
19 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, *
20 * Boston, MA 02110-1301, USA. *
21 * *
22 *****************************************************************************/
23
24 /*****************************************************************************
25 * *
26 * ElmerGUI meshutils *
27 * *
28 *****************************************************************************
29 * *
30 * Authors: Mikko Lyly, Juha Ruokolainen and Peter R�back *
31 * Email: Juha.Ruokolainen@csc.fi *
32 * Web: http://www.csc.fi/elmer *
33 * Address: CSC - IT Center for Science Ltd. *
34 * Keilaranta 14 *
35 * 02101 Espoo, Finland *
36 * *
37 * Original Date: 15 Mar 2008 *
38 * *
39 *****************************************************************************/
40
41 #include <QtCore>
42 #include <iostream>
43 #include <stdlib.h>
44 #include <cmath>
45 #include "meshutils.h"
46 using namespace std;
47
48 template <typename T>
AllocateDynamicArray(int nRows,int nCols)49 T **AllocateDynamicArray(int nRows, int nCols) {
50 T **dynamicArray;
51 dynamicArray = new T*[nRows];
52 for( int i = 0 ; i < nRows ; i++ )
53 dynamicArray[i] = new T [nCols];
54 return dynamicArray;
55 }
56
57 template <typename T>
FreeDynamicArray(T ** dArray)58 void FreeDynamicArray(T** dArray) {
59 delete [] *dArray;
60 delete [] dArray;
61 }
62
63 template <typename T>
AllocateDynamicVector(int nRows)64 T *AllocateDynamicVector(int nRows) {
65 T *dynamicArray;
66 dynamicArray = new T[nRows];
67 return dynamicArray;
68 }
69
70 template <typename T>
FreeDynamicVector(T * dArray)71 void FreeDynamicVector(T* dArray) {
72 delete [] dArray;
73 }
74
Meshutils()75 Meshutils::Meshutils()
76 {
77 }
78
79
~Meshutils()80 Meshutils::~Meshutils()
81 {
82 }
83
84
85 // Clear mesh...
86 //-----------------------------------------------------------------------------
clearMesh(mesh_t * mesh)87 void Meshutils::clearMesh(mesh_t *mesh)
88 {
89 if(mesh == NULL)
90 return;
91
92 mesh->clear();
93 delete mesh;
94
95 cout << "Old mesh cleared" << endl;
96 cout.flush();
97 }
98
99
100
101 // Find surface elements for 3D elements when they are not provided
102 //----------------------------------------------------------------------------
findSurfaceElements(mesh_t * mesh)103 void Meshutils::findSurfaceElements(mesh_t *mesh)
104 {
105 #define RESETENTRY1 \
106 h->node[0] = UNKNOWN; \
107 h->node[1] = UNKNOWN; \
108 h->element[0] = UNKNOWN; \
109 h->element[1] = UNKNOWN; \
110 h->index[0] = UNKNOWN; \
111 h->index[1] = UNKNOWN; \
112 h->next = NULL;
113
114 class hashEntry {
115 public:
116 int node[2];
117 int element[2];
118 int index[2];
119 int face;
120 hashEntry *next;
121 };
122
123
124 if(mesh->getElements() == 0) return;
125
126 int keys = mesh->getNodes();
127 hashEntry *hash = new hashEntry[keys];
128
129 bool found;
130 hashEntry *h;
131
132 for(int i=0; i<keys; i++) {
133 h = &hash[i];
134 RESETENTRY1;
135 }
136
137 // TODO: only 1st and 2nd order elements
138
139 static int familyfaces[9] = {0, 0, 0, 0, 0, 4, 5, 5, 6};
140
141 static int faceedges8[] = {4, 4, 4, 4, 4, 4};
142 static int facemap8[][8] = {{0,1,2,3,8,9,10,11}, {4,5,6,7,16,17,18,19}, {0,1,5,4,8,13,16,12}, {1,2,6,5,9,14,17,13}, {2,3,7,6,10,15,18,14}, {3,0,4,7,11,12,19,15}};
143
144 static int faceedges7[] = {3, 3, 4, 4, 4};
145 static int facemap7[][8] = {{0,1,2,6,7,8}, {3,4,5,12,13,14}, {0,1,4,3,6,10,12,9}, {1,2,5,4,7,11,13,10}, {2,0,3,5,8,9,14,11}};
146
147 static int faceedges6[] = {4, 3, 3, 3, 3};
148 static int facemap6[][8] = {{0,1,2,3,5,6,7,8}, {0,1,4,5,10,9}, {1,2,4,6,11,10}, {2,3,4,7,12,11}, {3,0,4,8,9,12}};
149
150 static int faceedges5[4] = {3, 3, 3, 3};
151 static int facemap5[][6] = {{0,1,2,4,5,6}, {0,1,3,4,8,7}, {1,2,3,5,9,8}, {2,0,3,6,7,9}};
152
153
154 for(int i=0; i < mesh->getElements(); i++) {
155 element_t *e = mesh->getElement(i);
156
157 int facenodes = 0;
158 int *facemap = NULL;
159 int n[4];
160
161 int family = e->getCode() / 100;
162 int faces = familyfaces[family];
163
164 for(int f=0; f<faces; f++) {
165 if(family == 5) {
166 facenodes = faceedges5[f];
167 facemap = &facemap5[f][0];
168 }
169 else if(family == 6) {
170 facenodes = faceedges6[f];
171 facemap = &facemap6[f][0];
172 }
173 else if(family == 7) {
174 facenodes = faceedges7[f];
175 facemap = &facemap7[f][0];
176 }
177 else if(family == 8) {
178 facenodes = faceedges8[f];
179 facemap = &facemap8[f][0];
180 }
181
182 for(int j=0; j < facenodes; j++)
183 n[j] = e->getNodeIndex(facemap[j]);
184
185 // Order indexes in an increasing order
186 for(int k=facenodes-1;k>0;k--) {
187 for(int j=0;j<k;j++) {
188 if(n[j] > n[j+1]) {
189 int tmp = n[j+1];
190 n[j+1] = n[j];
191 n[j] = tmp;
192 }
193 }
194 }
195
196 // three nodes define a face uniquely also for rectangles
197 h = &hash[n[0]];
198 found = false;
199 while(h->next) {
200 if((h->node[0] == n[1]) && (h->node[1] == n[2])) {
201 found = true;
202 break;
203 }
204 h = h->next;
205 }
206
207 if(!found) {
208 h->node[0] = n[1];
209 h->node[1] = n[2];
210 h->element[0] = i;
211 h->index[0] = e->getIndex();
212 h->face = f;
213
214 h->next = new hashEntry;
215 h = h->next;
216 RESETENTRY1;
217 } else {
218 h->index[1] = e->getIndex();
219 h->element[1] = i;
220 }
221 }
222 }
223
224 // count faces that have different materials at either sides:
225 int allsurfaces = 0;
226 int surfaces = 0;
227 int maxindex1 = 0;
228 int maxindex2 = 0;
229 for(int i=0; i<keys; i++) {
230 h = &hash[i];
231 while(h->next){
232 if(h->element[0] > UNKNOWN) allsurfaces++;
233 if(h->index[0] != h->index[1]) surfaces++;
234 if(h->index[0] > maxindex1) maxindex1 = h->index[0];
235 if(h->index[1] > maxindex2) maxindex2 = h->index[1];
236 h = h->next;
237 }
238 }
239 cout << "Found " << surfaces << " interface faces of " << allsurfaces << endl;
240
241
242 // Create a index table such that all combinations of materials
243 // get different BC index
244 //int indextable[maxindex1+1][maxindex2+1];
245 int **indextable = AllocateDynamicArray<int>(maxindex1+1, maxindex2+1);
246 int index1,index2;
247 for(int i=0;i<=maxindex1;i++)
248 for(int j=0;j<=maxindex2;j++)
249 indextable[i][j] = 0;
250
251 for(int i=0; i<keys; i++) {
252 h = &hash[i];
253 while(h->next){
254 if(h->index[0] != h->index[1]) {
255 index1 = h->index[0];
256 index2 = h->index[1];
257 if(index2 == -1) index2 = 0;
258 indextable[index1][index2] = 1;
259 }
260 h = h->next;
261 }
262 }
263 index1=0;
264 for(int i=0;i<=maxindex1;i++)
265 for(int j=0;j<=maxindex2;j++)
266 if(indextable[i][j]) indextable[i][j] = ++index1;
267
268 cout << "Boundaries were numbered up to index " << index1 << endl;
269
270 // Finally set the surfaces:
271 mesh->setSurfaces(surfaces);
272 mesh->newSurfaceArray(mesh->getSurfaces());
273
274 surfaces = 0;
275 for(int i=0; i<keys; i++) {
276 h = &hash[i];
277
278 while(h->next) {
279 if(h->index[0] != h->index[1]) {
280 surface_t *s = mesh->getSurface(surfaces);
281
282 s->setElements(1);
283 s->newElementIndexes(2);
284 s->setElementIndex(0, h->element[0]);
285 s->setElementIndex(1, h->element[1]);
286 if(s->getElementIndex(1) >= 0) s->setElements(2); // ?????
287
288 element_t *e = mesh->getElement(h->element[0]);
289 int code = e->getCode();
290 int family = code / 100;
291 int f = h->face;
292
293 int faceedges = 0;
294 int *facemap = NULL;
295 int degree = 1;
296
297 if(family == 5) {
298 faceedges = faceedges5[f];
299 if( code == 510) degree = 2;
300 facemap = &facemap5[f][0];
301 }
302 else if(family == 6) {
303 faceedges = faceedges6[f];
304 if( code == 613) degree = 2;
305 facemap = &facemap6[f][0];
306 }
307 else if(family == 7) {
308 faceedges = faceedges7[f];
309 if( code == 715) degree = 2;
310 facemap = &facemap7[f][0];
311 }
312 else if(family == 8) {
313 faceedges = faceedges8[f];
314 if( code == 820 ) degree = 2;
315 facemap = &facemap8[f][0];
316 }
317
318 int facenodes = degree * faceedges;
319
320 s->setNodes(facenodes);
321 s->setCode(100 * faceedges + facenodes);
322 s->newNodeIndexes(s->getNodes());
323 for(int j=0; j < s->getNodes(); j++)
324 s->setNodeIndex(j, e->getNodeIndex(facemap[j]));
325
326 index1 = h->index[0];
327 index2 = h->index[1];
328 if(index2 < 0) index2 = 0;
329
330 s->setIndex(indextable[index1][index2]);
331
332 s->setEdges(s->getNodes());
333 s->newEdgeIndexes(s->getEdges());
334 for(int j=0; j < s->getEdges(); j++)
335 s->setEdgeIndex(j, UNKNOWN);
336
337 s->setNature(PDE_BOUNDARY);
338 surfaces++;
339 }
340 h = h->next;
341 }
342 }
343
344 delete [] hash;
345 }
346
347
348
349 // Find parent elements for existing surfaces...
350 //----------------------------------------------------------------------------
findEdgeElementParents(mesh_t * mesh)351 void Meshutils::findEdgeElementParents(mesh_t *mesh)
352 {
353 #define RESETENTRY0 \
354 h->node[0] = UNKNOWN; \
355 h->node[1] = UNKNOWN; \
356 h->element[0] = UNKNOWN; \
357 h->element[1] = UNKNOWN; \
358 h->next = NULL;
359
360 class hashEntry {
361 public:
362 int node[2];
363 int element[2];
364 hashEntry *next;
365 };
366
367 int keys = mesh->getNodes();
368 hashEntry *hash = new hashEntry[keys];
369
370 bool found;
371 hashEntry *h;
372
373 for(int i=0; i<keys; i++) {
374 h = &hash[i];
375 RESETENTRY0;
376 }
377
378 // TODO: only tetrahedron at the moment
379
380 static int edgemap[][2] = {{0,1}, {1,2}, {2,0}};
381
382 for(int i = 0; i < mesh->getSurfaces(); i++) {
383 surface_t *s = mesh->getSurface(i);
384
385 for(int f = 0; f < 3; f++) {
386 int n0 = s->getNodeIndex(edgemap[f][0]);
387 int n1 = s->getNodeIndex(edgemap[f][1]);
388
389 if(n1 < n0) {
390 int tmp = n1;
391 n1 = n0;
392 n0 = tmp;
393 }
394
395 h = &hash[n0];
396 found = false;
397 while(h->next) {
398 if(h->node[0] == n1) {
399 found = true;
400 break;
401 }
402 h = h->next;
403 }
404
405 if(!found) {
406 h->node[0] = n1;
407 h->element[0] = i;
408 h->next = new hashEntry;
409 h = h->next;
410 RESETENTRY0;
411 } else {
412 h->element[1] = i;
413 }
414 }
415 }
416
417 // count faces:
418 int edges = 0;
419 for(int i = 0; i < keys; i++) {
420 h = &hash[i];
421 while((h = h->next) != NULL)
422 edges++;
423 }
424
425 cout << "Found total of " << edges << " edges" << endl;
426
427 // Finally find parents:
428 for(int i = 0; i < mesh->getEdges(); i++) {
429 edge_t *e = mesh->getEdge(i);
430
431 int n0 = e->getNodeIndex(0);
432 int n1 = e->getNodeIndex(1);
433
434 if(n1 < n0) {
435 int tmp = n1;
436 n1 = n0;
437 n0 = tmp;
438 }
439
440 h = &hash[n0];
441 while(h->next) {
442 if(h->node[0] == n1) {
443
444 // should we deallocate s->element if it exists?
445 e->setSurfaces(2);
446 e->newSurfaceIndexes(2);
447
448 e->setSurfaceIndex(0, h->element[0]);
449 e->setSurfaceIndex(1, h->element[1]);
450 }
451 h = h->next;
452 }
453 }
454
455 delete [] hash;
456 }
457
458
459 // Find parent elements for existing surfaces...
460 //----------------------------------------------------------------------------
findSurfaceElementParents(mesh_t * mesh)461 void Meshutils::findSurfaceElementParents(mesh_t *mesh)
462 {
463 #define RESETENTRY0 \
464 h->node[0] = UNKNOWN; \
465 h->node[1] = UNKNOWN; \
466 h->element[0] = UNKNOWN; \
467 h->element[1] = UNKNOWN; \
468 h->next = NULL;
469
470 class hashEntry {
471 public:
472 int node[2];
473 int element[2];
474 hashEntry *next;
475 };
476
477 int keys = mesh->getNodes();
478 hashEntry *hash = new hashEntry[keys];
479
480 bool found;
481 hashEntry *h;
482
483 for(int i=0; i<keys; i++) {
484 h = &hash[i];
485 RESETENTRY0;
486 }
487
488 // TODO: only tetrahedron at the moment
489
490 static int facemap[][3] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,2,3}};
491
492 for(int i=0; i < mesh->getElements(); i++) {
493 element_t *e = mesh->getElement(i);
494
495 for(int f=0; f<4; f++) {
496 int n0 = e->getNodeIndex(facemap[f][0]);
497 int n1 = e->getNodeIndex(facemap[f][1]);
498 int n2 = e->getNodeIndex(facemap[f][2]);
499
500 if(n2 < n1) {
501 int tmp = n2;
502 n2 = n1;
503 n1 = tmp;
504 }
505
506 if(n2 < n0) {
507 int tmp = n2;
508 n2 = n0;
509 n0 = tmp;
510 }
511
512 if(n1 < n0) {
513 int tmp = n1;
514 n1 = n0;
515 n0 = tmp;
516 }
517
518 h = &hash[n0];
519 found = false;
520 while(h->next) {
521 if((h->node[0] == n1) && (h->node[1] == n2)) {
522 found = true;
523 break;
524 }
525 h = h->next;
526 }
527
528 if(!found) {
529 h->node[0] = n1;
530 h->node[1] = n2;
531 h->element[0] = i;
532 h->next = new hashEntry;
533 h = h->next;
534 RESETENTRY0;
535 } else {
536 h->element[1] = i;
537 }
538 }
539 }
540
541 // count faces:
542 int faces = 0;
543 for(int i=0; i<keys; i++) {
544 h = &hash[i];
545 while((h = h->next) != NULL)
546 faces++;
547 }
548
549 cout << "Found total of " << faces << " faces" << endl;
550
551 // Finally find parents:
552 for(int i=0; i < mesh->getSurfaces(); i++) {
553 surface_t *s = mesh->getSurface(i);
554
555 int n0 = s->getNodeIndex(0);
556 int n1 = s->getNodeIndex(1);
557 int n2 = s->getNodeIndex(2);
558
559 if(n2 < n1) {
560 int tmp = n2;
561 n2 = n1;
562 n1 = tmp;
563 }
564
565 if(n2 < n0) {
566 int tmp = n2;
567 n2 = n0;
568 n0 = tmp;
569 }
570
571 if(n1 < n0) {
572 int tmp = n1;
573 n1 = n0;
574 n0 = tmp;
575 }
576
577 h = &hash[n0];
578 while(h->next) {
579 if((h->node[0] == n1) && (h->node[1] == n2)) {
580
581 // should we deallocate s->element if it exists?
582 s->setElements(2);
583 s->newElementIndexes(2);
584
585 s->setElementIndex(0, h->element[0]);
586 s->setElementIndex(1, h->element[1]);
587 }
588 h = h->next;
589 }
590 }
591
592 delete [] hash;
593 }
594
595
596
597 // Find points for edge elements...
598 //-----------------------------------------------------------------------------
findEdgeElementPoints(mesh_t * mesh)599 void Meshutils::findEdgeElementPoints(mesh_t *mesh)
600 {
601 class hashEntry {
602 public:
603 int edges;
604 int *edge;
605 };
606
607 int keys = mesh->getNodes();
608
609 hashEntry *hash = new hashEntry[keys];
610
611 for(int i = 0; i < keys; i++) {
612 hashEntry *h = &hash[i];
613 h->edges = 0;
614 h->edge = NULL;
615 }
616
617 for(int i = 0; i < mesh->getEdges(); i++) {
618 edge_t *e = mesh->getEdge(i);
619
620 if(e->getNature() == PDE_BOUNDARY) {
621 for(int k = 0; k < 2; k++) {
622 int n = e->getNodeIndex(k);
623
624 hashEntry *h = &hash[n];
625
626 bool found = false;
627 for(int j = 0; j < h->edges; j++) {
628 if(h->edge[j] == i) {
629 found = true;
630 break;
631 }
632 }
633
634 if(!found) {
635 int *tmp = new int[h->edges+1];
636 for(int j = 0; j < h->edges; j++)
637 tmp[j] = h->edge[j];
638 tmp[h->edges] = i;
639 delete [] h->edge;
640
641 h->edges++;
642 h->edge = tmp;
643 }
644 }
645 }
646 }
647
648 // count points:
649 int count = 0;
650 for(int i = 0; i < keys; i++) {
651 hashEntry *h = &hash[i];
652 if(h->edges > 0)
653 count++;
654 }
655
656 cout << "Found " << count << " points on boundary edges" << endl;
657 cout.flush();
658
659 // delete old points, if any:
660 if(mesh->getPoints() > 0) {
661 cout << "Deleting old points and creating new" << endl;
662 cout.flush();
663 for(int i = 0; i < mesh->getPoints(); i++) {
664 mesh->getPoint(i)->deleteNodeIndexes();
665 mesh->getPoint(i)->deleteEdgeIndexes();
666 }
667 mesh->deletePointArray();
668 }
669
670 mesh->setPoints(count);
671 mesh->newPointArray(mesh->getPoints());
672
673 count = 0;
674 for(int i = 0; i < keys; i++) {
675 hashEntry *h = &hash[i];
676
677 if(h->edges > 0) {
678 point_t *p = mesh->getPoint(count++);
679 p->setNodes(1);
680 p->newNodeIndexes(1);
681 p->setNodeIndex(0, i);
682 p->setEdges(h->edges);
683 p->newEdgeIndexes(p->getEdges());
684 for(int j = 0; j < p->getEdges(); j++) {
685 p->setEdgeIndex(j, h->edge[j]);
686 }
687 p->setSharp(false);
688 }
689 }
690
691 // delete temp stuff
692 for(int i = 0; i < keys; i++) {
693 hashEntry *h = &hash[i];
694 if(h->edges > 0)
695 delete [] h->edge;
696 }
697
698 delete [] hash;
699
700 // Inverse map
701 cout << "Constructing inverse map from edges to points" << endl;
702 cout.flush();
703
704 for(int i=0; i < mesh->getPoints(); i++) {
705 point_t *p = mesh->getPoint(i);
706
707 for(int j=0; j < p->getEdges(); j++) {
708 int k = p->getEdgeIndex(j);
709 if ( k<0 ) continue;
710
711 edge_t *e = mesh->getEdge(k);
712
713 // allocate space for two points, if not yet done:
714 if(e->getPoints() < 2) {
715 e->setPoints(2);
716 e->newPointIndexes(2);
717 e->setPointIndex(0, -1);
718 e->setPointIndex(1, -1);
719 }
720
721 for(int r=0; r < e->getPoints(); r++) {
722 if(e->getPointIndex(r) < 0) {
723 e->setPointIndex(r, i);
724 break;
725 }
726 }
727 }
728 }
729 }
730
731
732
733 // Find edges for surface elements...
734 //-----------------------------------------------------------------------------
findSurfaceElementEdges(mesh_t * mesh)735 void Meshutils::findSurfaceElementEdges(mesh_t *mesh)
736 {
737 #define RESETENTRY \
738 h->surfaces = 0; \
739 h->surface = NULL; \
740 h->parentindex[0] = UNKNOWN; \
741 h->parentindex[1] = UNKNOWN; \
742 h->next = NULL;
743
744 int keys = mesh->getNodes();
745
746 class hashEntry {
747 public:
748 int nodes;
749 int *node;
750 int nature;
751 int index;
752 int surfaces;
753 int *surface;
754 int parentindex[2];
755 hashEntry *next;
756 };
757
758 hashEntry *hash = new hashEntry[keys];
759
760 bool found;
761 hashEntry *h;
762
763 for(int i=0; i<keys; i++) {
764 h = &hash[i];
765 RESETENTRY;
766 }
767
768 bool createindexes = ((!mesh->getEdges()) && (!mesh->getElements()));
769
770 // ????? should we test for mesh->edges ?????
771 // if ( mesh->edge && (mesh->getEdges() > 0) ) {
772 if ( mesh->getEdges() > 0 ) {
773 // add existing edges first:
774 for( int i=0; i<mesh->getEdges(); i++ )
775 {
776 edge_t *edge = mesh->getEdge(i);
777 int n0 = edge->getNodeIndex(0);
778 int n1 = edge->getNodeIndex(1);
779
780 int m = (n0<n1) ? n0 : n1;
781 int n = (n0<n1) ? n1 : n0;
782
783 h = &hash[m];
784 found = false;
785 while(h->next) {
786 if(h->node[0] == n) {
787 found = true;
788 break;
789 }
790 h = h->next;
791 }
792
793 if(!found) {
794 h->nodes = edge->getNodes() - 1;
795 h->node = new int[h->nodes];
796 h->node[0] = n;
797 for( int j=1; j<h->nodes; j++ )
798 {
799 h->node[j] = edge->getNodeIndex(j+1);
800 }
801 h->surfaces = edge->getSurfaces();
802 h->surface = new int[edge->getSurfaces()];
803 for( int j=0; j < edge->getSurfaces(); j++ ) {
804 h->surface[j] = edge->getSurfaceIndex(j);
805 }
806 h->index = edge->getIndex();
807 h->nature = edge->getNature();
808 h->next = new hashEntry;
809 h = h->next;
810 RESETENTRY;
811 }
812 }
813
814 mesh->setEdges(0);
815 mesh->deleteEdgeArray(); // delete [] mesh->edge;
816 }
817
818
819 static int triedgemap[][4] = { {0,1,3,6}, {1,2,4,7}, {2,0,5,8} };
820 static int quadedgemap[][4] = {{0,1,4,8}, {1,2,5,9}, {2,3,6,10}, {3,0,7,11}};
821
822
823 for(int i=0; i < mesh->getSurfaces(); i++) {
824 surface_t *s = mesh->getSurface(i);
825
826 // loop over edges
827 for(int e=0; e < s->getEdges(); e++) {
828 int n0, n1,nrest[2] = { -1,-1 };
829 if((int)(s->getCode() / 100) == 3) {
830 n0 = s->getNodeIndex(triedgemap[e][0]);
831 n1 = s->getNodeIndex(triedgemap[e][1]);
832 if ( s->getCode() >= 306) nrest[0] = s->getNodeIndex(triedgemap[e][2]);
833 if ( s->getCode() >= 309) nrest[1] = s->getNodeIndex(triedgemap[e][3]);
834 } else if((int)(s->getCode() / 100) == 4) {
835 n0 = s->getNodeIndex(quadedgemap[e][0]);
836 n1 = s->getNodeIndex(quadedgemap[e][1]);
837 // corrected 4.11.2008:
838 if ( s->getCode() >= 408) nrest[0] = s->getNodeIndex(quadedgemap[e][2]);
839 if ( s->getCode() >= 412) nrest[1] = s->getNodeIndex(quadedgemap[e][3]);
840 } else {
841 cout << "findBoundaryElementEdges: error: unknown element code" << endl;
842 exit(0);
843 }
844
845 int m = (n0<n1) ? n0 : n1;
846 int n = (n0<n1) ? n1 : n0;
847
848 h = &hash[m];
849 found = false;
850 while(h->next) {
851 if(h->node[0] == n) {
852 found = true;
853 break;
854 }
855 h = h->next;
856 }
857
858 if(!found) {
859 h->nodes = 1;
860 h->nodes += nrest[0] >=0 ? 1:0;
861 h->nodes += nrest[1] >=0 ? 1:0;
862 h->node = new int[h->nodes];
863 h->node[0] = n;
864 for( int j=1; j<h->nodes; j++ )
865 h->node[j] = nrest[j-1];
866
867 h->surfaces = 1;
868 h->surface = new int[1];
869 h->surface[0] = i;
870 h->parentindex[0] = s->getIndex();
871 h->index = UNKNOWN;
872 h->nature = PDE_UNKNOWN;
873 h->next = new hashEntry;
874 h = h->next;
875 RESETENTRY;
876 } else {
877 int *tmp = new int[h->surfaces];
878
879 found = false;
880 for(int j=0; j<h->surfaces; j++)
881 {
882 tmp[j] = h->surface[j];
883 if ( tmp[j] == i ) found = true;
884 }
885 if ( found ) {
886 delete [] tmp;
887 } else {
888 delete [] h->surface;
889 h->surface = new int[h->surfaces+1];
890 for(int j=0; j<h->surfaces; j++)
891 h->surface[j] = tmp[j];
892 h->surface[h->surfaces++] = i;
893 if( s->getIndex() < h->parentindex[0]) {
894 h->parentindex[1] = h->parentindex[0];
895 h->parentindex[0] = s->getIndex();
896 }
897 else {
898 h->parentindex[1] = s->getIndex();
899 }
900 delete [] tmp;
901 }
902 }
903 }
904 }
905
906 // count edges:
907 int edges = 0;
908 for(int i=0; i<keys; i++) {
909 h = &hash[i];
910 while((h = h->next) != NULL)
911 edges++;
912 }
913
914
915
916 cout << "Found " << edges << " edges on boundary" << endl;
917
918 mesh->setEdges(edges);
919 mesh->newEdgeArray(edges); // edge = new edge_t[edges];
920
921 if(createindexes) {
922 cout << "Creating edge indexes " << endl;
923
924 int maxindex1 = UNKNOWN;
925 int maxindex2 = UNKNOWN;
926
927 for(int i=0; i<keys; i++) {
928 h = &hash[i];
929 while(h->next){
930 if(h->parentindex[0] > maxindex1) maxindex1 = h->parentindex[0];
931 if(h->parentindex[1] > maxindex2) maxindex2 = h->parentindex[1];
932 h = h->next;
933 }
934 }
935
936 // Create a index table such that all combinations of materials
937 // get different BC index
938 //int indextable[maxindex1+2][maxindex2+2];
939 int **indextable = AllocateDynamicArray<int>(maxindex1+2, maxindex2+2);
940 int index1,index2;
941 for(int i=-1;i<=maxindex1;i++)
942 for(int j=-1;j<=maxindex2;j++)
943 indextable[i+1][j+1] = 0;
944
945 int edgebcs=0;
946 for(int i=0; i<keys; i++) {
947 h = &hash[i];
948 while(h->next){
949 if(h->parentindex[0] != h->parentindex[1]) {
950 index1 = h->parentindex[0];
951 index2 = h->parentindex[1];
952 indextable[index1+1][index2+1] = 1;
953 edgebcs += 1;
954 }
955 h = h->next;
956 }
957 }
958
959 index1=0;
960 for(int i=-1;i<=maxindex1;i++)
961 for(int j=-1;j<=maxindex2;j++)
962 if(indextable[i+1][j+1]) indextable[i+1][j+1] = ++index1;
963
964 cout << edgebcs << " boundary edges were numbered up to index " << index1 << endl;
965
966 for(int i=0; i<keys; i++) {
967 h = &hash[i];
968 while(h->next){
969 if(h->parentindex[0] != h->parentindex[1]) {
970 index1 = h->parentindex[0];
971 index2 = h->parentindex[1];
972 h->index = indextable[index1+1][index2+1];
973 h->nature = PDE_BOUNDARY;
974 }
975 h = h->next;
976 }
977 }
978 }
979
980
981 // Create edges:
982 edges = 0;
983 for(int i=0; i<keys; i++) {
984 h = &hash[i];
985 while(h->next) {
986 edge_t *e = mesh->getEdge(edges++);
987
988 e->setNature(h->nature);
989 e->setNodes(h->nodes + 1);
990 e->newNodeIndexes(e->getNodes());
991 e->setNodeIndex(0, i); // ?????
992 for( int j=1; j < e->getNodes(); j++ )
993 e->setNodeIndex(j, h->node[j-1]);
994
995 e->setCode(200 + e->getNodes());
996
997 e->setSurfaces(h->surfaces);
998 e->newSurfaceIndexes(max(e->getSurfaces(), 2));
999 e->setSurfaceIndex(0, -1);
1000 e->setSurfaceIndex(1, -1);
1001
1002 for(int j=0; j < e->getSurfaces(); j++)
1003 e->setSurfaceIndex(j, h->surface[j]);
1004
1005 e->setSharp(false);
1006
1007 e->setIndex(h->index);
1008 e->setPoints(0);
1009 h = h->next;
1010 }
1011 }
1012
1013 delete [] hash;
1014
1015 // Inverse map
1016 for(int i=0; i < mesh->getEdges(); i++) {
1017 edge_t *e = mesh->getEdge(i);
1018
1019 for(int j=0; j < e->getSurfaces(); j++) {
1020 int k = e->getSurfaceIndex(j);
1021 if ( k< 0 ) continue;
1022
1023 surface_t *s = mesh->getSurface(k);
1024
1025 for(int r=0; r < s->getEdges(); r++) {
1026 if(s->getEdgeIndex(r) < 0) {
1027 s->setEdgeIndex(r, i);
1028 break;
1029 }
1030 }
1031 }
1032 }
1033
1034 #if 0
1035 cout << "*********************" << endl;
1036 for(int i=0; i<mesh->getEdges(); i++)
1037 cout << "Edge " << i << " nodes " << mesh->edge[i].node[0] << " "<< mesh->edge[i].node[0] << endl;
1038
1039 for(int i=0; i < mesh->getSurfaces(); i++)
1040 cout << "Surface " << i << " nodes "
1041 << mesh->surface[i].node[0] << " "
1042 << mesh->surface[i].node[1] << " "
1043 << mesh->surface[i].node[2] << " "
1044 << " Edges "
1045 << mesh->surface[i].edge[0] << " "
1046 << mesh->surface[i].edge[1] << " "
1047 << mesh->surface[i].edge[2] << " "
1048 << " Parents "
1049 << mesh->surface[i].element[0] << " "
1050 << mesh->surface[i].element[1] << " "
1051 << endl;
1052
1053 cout.flush();
1054 #endif
1055
1056 }
1057
1058
1059
1060
1061 // Find sharp points for edge elements...
1062 //-----------------------------------------------------------------------------
findSharpPoints(mesh_t * mesh,double limit)1063 void Meshutils::findSharpPoints(mesh_t *mesh, double limit)
1064 {
1065 qreal t0[3], t1[3];
1066
1067 cout << "Limit: " << limit << " degrees" << endl;
1068 cout.flush();
1069
1070 double angle;
1071 int count = 0;
1072 point_t *point = NULL;
1073 edge_t *edge = NULL;
1074 Helpers *helpers = new Helpers;
1075
1076 for(int i=0; i < mesh->getPoints(); i++) {
1077 point = mesh->getPoint(i);
1078
1079 if(point->getEdges() == 2) {
1080 int n = point->getNodeIndex(0);
1081
1082 int e0 = point->getEdgeIndex(0);
1083 int e1 = point->getEdgeIndex(1);
1084
1085 edge = mesh->getEdge(e0);
1086 int n0 = edge->getNodeIndex(0);
1087 if(edge->getNodeIndex(1) != n)
1088 n0 = edge->getNodeIndex(1);
1089
1090 edge = mesh->getEdge(e1);
1091 int n1 = edge->getNodeIndex(0);
1092 if(edge->getNodeIndex(1) != n)
1093 n1 = edge->getNodeIndex(1);
1094
1095 // unit tangent from node to node0
1096 t0[0] = mesh->getNode(n0)->getX(0) - mesh->getNode(n)->getX(0);
1097 t0[1] = mesh->getNode(n0)->getX(1) - mesh->getNode(n)->getX(1);
1098 t0[2] = mesh->getNode(n0)->getX(2) - mesh->getNode(n)->getX(2);
1099
1100 // unit tangent from node to node1
1101 t1[0] = mesh->getNode(n1)->getX(0) - mesh->getNode(n)->getX(0);
1102 t1[1] = mesh->getNode(n1)->getX(1) - mesh->getNode(n)->getX(1);
1103 t1[2] = mesh->getNode(n1)->getX(2) - mesh->getNode(n)->getX(2);
1104
1105 helpers->normalize(t0);
1106 helpers->normalize(t1);
1107
1108 double cosofangle = t0[0]*t1[0] + t0[1]*t1[1] + t0[2]*t1[2];
1109 angle = acos(cosofangle) / MYPI * 180.0;
1110 } else {
1111 angle = 0.0;
1112 }
1113
1114 point->setSharp(false);
1115 if(sqrt(angle*angle) < (180.0-limit) ) {
1116 point->setSharp(true);
1117 count++;
1118 }
1119 }
1120
1121 cout << "Found " << count << " sharp points" << endl;
1122 delete helpers;
1123 }
1124
1125
1126
1127 // Find sharp edges for surface elements...
1128 //-----------------------------------------------------------------------------
findSharpEdges(mesh_t * mesh,double limit)1129 void Meshutils::findSharpEdges(mesh_t *mesh, double limit)
1130 {
1131
1132 cout << "Limit: " << limit << " degrees" << endl;
1133 cout.flush();
1134
1135 double angle;
1136 int count = 0;
1137
1138 for(int i=0; i<mesh->getEdges(); i++) {
1139 edge_t *edge = mesh->getEdge(i);
1140
1141 if(edge->getSurfaces() == 2) {
1142 int s0 = edge->getSurfaceIndex(0);
1143 int s1 = edge->getSurfaceIndex(1);
1144 double *n0 = mesh->getSurface(s0)->getNormalVec();
1145 double *n1 = mesh->getSurface(s1)->getNormalVec();
1146 double cosofangle = n0[0]*n1[0] + n0[1]*n1[1] + n0[2]*n1[2];
1147 cosofangle = std::abs(cosofangle);
1148 angle = acos(cosofangle) / MYPI * 180.0;
1149 } else {
1150 angle = 180.0;
1151 }
1152
1153 edge->setSharp(false);
1154 if(sqrt(angle*angle) > limit) {
1155 edge->setSharp(true);
1156 count++;
1157 }
1158 }
1159
1160 cout << "Found " << count << " sharp edges" << endl;
1161 }
1162
1163
1164
1165 // Divide edge by sharp points...
1166 //-----------------------------------------------------------------------------
divideEdgeBySharpPoints(mesh_t * mesh)1167 int Meshutils::divideEdgeBySharpPoints(mesh_t *mesh)
1168 {
1169 class Bc {
1170 public:
1171 void propagateIndex(mesh_t* mesh, int index, int i) {
1172 edge_t *edge = mesh->getEdge(i);
1173
1174 // index is ok
1175 if(!edge->getSelected() || (edge->getIndex() != UNKNOWN)
1176 || (edge->getNature() != PDE_BOUNDARY)) return;
1177
1178 // set index
1179 edge->setIndex(index);
1180
1181 // propagate index
1182 for(int j=0; j < edge->getPoints(); j++) {
1183 int k = edge->getPointIndex(j);
1184 point_t *point = mesh->getPoint(k);
1185
1186 // skip sharp points
1187 if(!point->isSharp()) {
1188 for(int m = 0; m < point->getEdges(); m++) {
1189 int n = point->getEdgeIndex(m);
1190 propagateIndex(mesh, index, n);
1191 }
1192 }
1193
1194 }
1195 }
1196 };
1197
1198
1199 // reset bc-indices on edges:
1200 int index = 0;
1201 int count = 0;
1202
1203 for(int i=0; i < mesh->getEdges(); i++)
1204 {
1205 edge_t *edge = mesh->getEdge(i);
1206 if((edge->getNature() == PDE_BOUNDARY) && !edge->getSelected())
1207 index = max(index, edge->getIndex());
1208 }
1209 index++;
1210
1211 // int edge_index[index];
1212 int *edge_index = new int[index];
1213
1214 for( int i=0; i<index; i++ )
1215 edge_index[i] = UNKNOWN;
1216
1217 index = 0;
1218 for(int i=0; i < mesh->getEdges(); i++)
1219 {
1220 edge_t *edge = mesh->getEdge(i);
1221 if (edge->getNature() == PDE_BOUNDARY)
1222 {
1223 if ( edge->getSelected() ) {
1224 count++;
1225 edge->setIndex(UNKNOWN);
1226 } else {
1227 if ( edge_index[edge->getIndex()] == UNKNOWN )
1228 edge_index[edge->getIndex()] = ++index;
1229 edge->setIndex(edge_index[edge->getIndex()]);
1230 }
1231 }
1232 }
1233
1234 if ( count==0 ) {
1235 cout << "No boundary edges to divde." << endl;
1236 return 0;
1237 }
1238
1239 Bc *bc = new Bc;
1240
1241 // recursively determine boundary parts:
1242 for(int i=0; i < mesh->getEdges(); i++)
1243 {
1244 edge_t *edge = mesh->getEdge(i);
1245 if(edge->getSelected() && (edge->getIndex() == UNKNOWN) && (edge->getNature() == PDE_BOUNDARY))
1246 bc->propagateIndex(mesh, ++index, i);
1247 }
1248 index++;
1249
1250 // Create a hopefully mesh independent indexing of groupings to enable
1251 // reapplying merge/division operations after remeshing. The indices
1252 // are given based on group bounding box corner distances from a given
1253 // point. Fails if two groups have same bbox, which should not happen
1254 // often though (?)
1255
1256 //double xmin[index], ymin[index], zmin[index];
1257 //double xmax[index], ymax[index], zmax[index];
1258 //double xc = 0, yc = 0, zc = 0, dist[index];
1259 //int cc[index], order[index], sorder[index];
1260 //double g_xmin,g_xmax,g_ymin,g_ymax,g_zmin,g_zmax;
1261
1262 double *xmin = new double[index];
1263 double *ymin = new double[index];
1264 double *zmin = new double[index];
1265 double *xmax = new double[index];
1266 double *ymax = new double[index];
1267 double *zmax = new double[index];
1268 double xc = 0, yc = 0, zc = 0;
1269 double *dist = new double[index];
1270 int *cc = new int[index];
1271 int *order = new int[index];
1272 int *sorder = new int[index];
1273 double g_xmin,g_xmax,g_ymin,g_ymax,g_zmin,g_zmax;
1274
1275 for( int i=0; i<index; i++ )
1276 {
1277 cc[i] = 0;
1278 order[i] = i;
1279 xmin[i] = ymin[i] = zmin[i] = 1e20;
1280 xmax[i] = ymax[i] = zmax[i] = -1e20;
1281 }
1282
1283 for( int i=0; i<mesh->getEdges(); i++ )
1284 {
1285 edge_t *edge = mesh->getEdge(i);
1286 if (edge->getNature() == PDE_BOUNDARY) {
1287 int k = edge->getIndex();
1288 for( int j=0; j < edge->getNodes(); j++ ) {
1289 int n = edge->getNodeIndex(j);
1290 cc[k]++;
1291 xmin[k] = min(xmin[k], mesh->getNode(n)->getX(0));
1292 ymin[k] = min(ymin[k], mesh->getNode(n)->getX(1));
1293 zmin[k] = min(zmin[k], mesh->getNode(n)->getX(2));
1294
1295 xmax[k] = max(xmax[k], mesh->getNode(n)->getX(0));
1296 ymax[k] = max(ymax[k], mesh->getNode(n)->getX(1));
1297 zmax[k] = max(zmax[k], mesh->getNode(n)->getX(2));
1298 }
1299 }
1300 }
1301
1302 g_xmin = g_ymin = g_zmin = 1e20;
1303 g_xmax = g_ymax = g_zmax = -1e20;
1304 for( int i=0; i<index; i++)
1305 {
1306 g_xmin = min(xmin[i],g_xmin);
1307 g_ymin = min(ymin[i],g_ymin);
1308 g_zmin = min(zmin[i],g_zmin);
1309
1310 g_xmax = max(xmax[i],g_xmax);
1311 g_ymax = max(ymax[i],g_ymax);
1312 g_zmax = max(zmax[i],g_zmax);
1313 }
1314
1315 double g_scale = max(max(g_xmax-g_xmin,g_ymax-g_ymin),g_zmax-g_zmin);
1316 double g_xp = g_xmax + 32.1345 * g_scale;
1317 double g_yp = g_ymin - 5.3*MYPI * g_scale;
1318 double g_zp = g_zmax + 8.1234 * g_scale;
1319
1320 for( int i=0; i<index; i++ )
1321 {
1322 dist[i] = 0;
1323 if ( cc[i]>0 ) {
1324 for( int j=0; j<8; j++ ) {
1325 switch(j) {
1326 case 0: xc=xmin[i]; yc=ymin[i]; zc=zmin[i]; break;
1327 case 1: xc=xmax[i]; break;
1328 case 2: yc=xmax[i]; break;
1329 case 3: xc=xmin[i]; break;
1330 case 4: zc=zmax[i]; break;
1331 case 5: yc=ymin[i]; break;
1332 case 6: xc=xmax[i]; break;
1333 case 7: yc=ymax[i]; break;
1334 }
1335 dist[i] += (xc-g_xp)*(xc-g_xp);
1336 dist[i] += (yc-g_yp)*(yc-g_yp);
1337 dist[i] += (zc-g_zp)*(zc-g_zp);
1338 }
1339 }
1340 }
1341
1342 sort_index( index, dist, order );
1343 for( int i=0; i<index; i++ )
1344 sorder[order[i]] = i;
1345
1346 for( int i=0; i<mesh->getEdges(); i++ )
1347 if ( mesh->getEdge(i)->getNature() == PDE_BOUNDARY )
1348 mesh->getEdge(i)->setIndex(sorder[mesh->getEdge(i)->getIndex()]);
1349
1350 --index;
1351 cout << "Edge divided into " << index << " parts" << endl;
1352
1353 delete bc;
1354
1355 return index;
1356 }
1357
1358
sort_index(int n,double * a,int * b)1359 void Meshutils::sort_index(int n, double *a, int *b)
1360 {
1361 QList<QPair<double, int> > list;
1362
1363 for(int i = 0; i < n; i++)
1364 list << qMakePair(a[i], b[i]);
1365
1366 qSort(list);
1367
1368 for(int i = 0; i < n; i++) {
1369 a[i] = list[i].first;
1370 b[i] = list[i].second;
1371 }
1372 }
1373
1374
1375 // Divide surface by sharp edges...
1376 //-----------------------------------------------------------------------------
divideSurfaceBySharpEdges(mesh_t * mesh)1377 int Meshutils::divideSurfaceBySharpEdges(mesh_t *mesh)
1378 {
1379 class Bc {
1380 public:
1381 void propagateIndex(mesh_t* mesh, int index, int i) {
1382 surface_t *surf = mesh->getSurface(i);
1383
1384 // index is ok
1385 if(!surf->getSelected() || (surf->getIndex() != UNKNOWN)
1386 || (surf->getNature() != PDE_BOUNDARY) ) return;
1387
1388 // set index
1389 surf->setIndex(index);
1390
1391 // propagate index
1392 for(int j=0; j < surf->getEdges(); j++) {
1393 int k = surf->getEdgeIndex(j);
1394 edge_t *edge = mesh->getEdge(k);
1395
1396 // skip sharp edges
1397 if(!edge->isSharp()) {
1398 for(int m=0; m < edge->getSurfaces(); m++) {
1399 int n = edge->getSurfaceIndex(m);
1400 propagateIndex(mesh, index, n);
1401 }
1402 }
1403 }
1404 }
1405 };
1406
1407 // reset bc-indices:
1408 int count = 0;
1409 int index = 0;
1410
1411 for( int i=0; i<mesh->getSurfaces(); i++ )
1412 {
1413 surface_t *surf = mesh->getSurface(i);
1414 if( (surf->getNature() == PDE_BOUNDARY) && !surf->getSelected() )
1415 index = max(index, surf->getIndex());
1416 }
1417 index++;
1418
1419 // int surf_index[index];
1420 int *surf_index = new int[index];
1421
1422 for( int i=0; i<index; i++ )
1423 surf_index[i] = UNKNOWN;
1424
1425 index = 0;
1426 for(int i=0; i < mesh->getSurfaces(); i++)
1427 {
1428 surface_t *surf = mesh->getSurface(i);
1429 if (surf->getNature() == PDE_BOUNDARY) {
1430 if ( surf->getSelected() ) {
1431 count++;
1432 surf->setIndex(UNKNOWN);
1433 } else {
1434 if ( surf_index[surf->getIndex()] == UNKNOWN )
1435 surf_index[surf->getIndex()] = ++index;
1436 surf->setIndex(surf_index[surf->getIndex()]);
1437 }
1438 }
1439 }
1440
1441 if ( count==0 ) {
1442 cout << "No boundary surfaces to divde." << endl;
1443 return 0;
1444 }
1445
1446
1447 // recursively determine boundary parts:
1448 Bc *bc = new Bc;
1449
1450 for(int i=0; i < mesh->getSurfaces(); i++)
1451 {
1452 surface_t *surf = mesh->getSurface(i);
1453 if(surf->getSelected() && (surf->getIndex() == UNKNOWN) && (surf->getNature() == PDE_BOUNDARY))
1454 bc->propagateIndex(mesh, ++index, i);
1455 }
1456 index++;
1457
1458 // Create a hopefully mesh independent indexing of groupings to enable
1459 // reapplying merge/division operations after remeshing. The indices
1460 // are given based on group bounding box corner distances from a given
1461 // point. Fails if two groups have same bbox, which should not happen
1462 // often though (?)
1463
1464 //double xmin[index], ymin[index], zmin[index];
1465 //double xmax[index], ymax[index], zmax[index];
1466 //double xc = 0, yc = 0, zc = 0, dist[index];
1467 //int cc[index], order[index], sorder[index];
1468 //double g_xmin,g_xmax,g_ymin,g_ymax,g_zmin,g_zmax;
1469
1470 double *xmin = new double[index];
1471 double *ymin = new double[index];
1472 double *zmin = new double[index];
1473 double *xmax = new double[index];
1474 double *ymax = new double[index];
1475 double *zmax = new double[index];
1476 double xc = 0, yc = 0, zc = 0;
1477 double *dist = new double[index];
1478 int *cc = new int[index];
1479 int *order = new int[index];
1480 int *sorder = new int[index];
1481 double g_xmin,g_xmax,g_ymin,g_ymax,g_zmin,g_zmax;
1482
1483 for( int i=0; i<index; i++ )
1484 {
1485 cc[i] = 0;
1486 order[i] = i;
1487 xmin[i] = ymin[i] = zmin[i] = 1e20;
1488 xmax[i] = ymax[i] = zmax[i] = -1e20;
1489 }
1490
1491 for( int i=0; i < mesh->getSurfaces(); i++ )
1492 {
1493 surface_t *surf = mesh->getSurface(i);
1494 if ( mesh->getSurface(i)->getNature() == PDE_BOUNDARY ) {
1495 int k = surf->getIndex();
1496 for( int j=0; j < surf->getNodes(); j++ ) {
1497 int n = surf->getNodeIndex(j);
1498 cc[k]++;
1499 xmin[k] = min( xmin[k], mesh->getNode(n)->getX(0) );
1500 ymin[k] = min( ymin[k], mesh->getNode(n)->getX(1) );
1501 zmin[k] = min( zmin[k], mesh->getNode(n)->getX(2) );
1502
1503 xmax[k] = max( xmax[k], mesh->getNode(n)->getX(0) );
1504 ymax[k] = max( ymax[k], mesh->getNode(n)->getX(1) );
1505 zmax[k] = max( zmax[k], mesh->getNode(n)->getX(2) );
1506 }
1507 }
1508 }
1509
1510 g_xmin = g_ymin = g_zmin = 1e20;
1511 g_xmax = g_ymax = g_zmax = -1e20;
1512 for( int i=0; i<index; i++)
1513 {
1514 g_xmin = min(xmin[i],g_xmin);
1515 g_ymin = min(ymin[i],g_ymin);
1516 g_zmin = min(zmin[i],g_zmin);
1517
1518 g_xmax = max(xmax[i],g_xmax);
1519 g_ymax = max(ymax[i],g_ymax);
1520 g_zmax = max(zmax[i],g_zmax);
1521 }
1522
1523 double g_scale = max(max(g_xmax-g_xmin,g_ymax-g_ymin),g_zmax-g_zmin);
1524 double g_xp = g_xmax + 32.1345 * g_scale;
1525 double g_yp = g_ymin - 5.3*MYPI * g_scale;
1526 double g_zp = g_zmax + 8.1234 * g_scale;
1527
1528 for( int i=0; i<index; i++ )
1529 {
1530 dist[i] = 0;
1531 if ( cc[i]>0 ) {
1532 for( int j=0; j<8; j++ ) {
1533 switch(j) {
1534 case 0: xc=xmin[i]; yc=ymin[i]; zc=zmin[i]; break;
1535 case 1: xc=xmax[i]; break;
1536 case 2: yc=xmax[i]; break;
1537 case 3: xc=xmin[i]; break;
1538 case 4: zc=zmax[i]; break;
1539 case 5: yc=ymin[i]; break;
1540 case 6: xc=xmax[i]; break;
1541 case 7: yc=ymax[i]; break;
1542 }
1543 dist[i] += (xc-g_xp)*(xc-g_xp);
1544 dist[i] += (yc-g_yp)*(yc-g_yp);
1545 dist[i] += (zc-g_zp)*(zc-g_zp);
1546 }
1547 }
1548 }
1549
1550 sort_index( index, dist, order );
1551 for( int i=0; i<index; i++ )
1552 sorder[order[i]] = i;
1553
1554 for( int i=0; i < mesh->getSurfaces(); i++ )
1555 if ( mesh->getSurface(i)->getNature() == PDE_BOUNDARY )
1556 mesh->getSurface(i)->setIndex(sorder[mesh->getSurface(i)->getIndex()]);
1557 --index;
1558
1559 cout << "Surface divided into " << index << " parts" << endl;
1560
1561 delete bc;
1562
1563 return index;
1564 }
1565
1566
1567
1568 // Find surface element normals...
1569 //-----------------------------------------------------------------------------
findSurfaceElementNormals(mesh_t * mesh)1570 void Meshutils::findSurfaceElementNormals(mesh_t *mesh)
1571 {
1572 static qreal a[3], b[3], c[3];
1573 double center_surface[3], center_element[3], center_difference[3];
1574 Helpers *helpers = new Helpers;
1575 int u, v, w, e0, e1, i0, i1, bigger;
1576
1577 for(int i=0; i < mesh->getSurfaces(); i++) {
1578 surface_t *surface = mesh->getSurface(i);
1579
1580 u = surface->getNodeIndex(0);
1581 v = surface->getNodeIndex(1);
1582
1583 if((int)(surface->getCode() / 100) == 3) {
1584 w = surface->getNodeIndex(2);
1585 } else if((int)(surface->getCode() / 100) == 4) {
1586 w = surface->getNodeIndex(3);
1587 } else {
1588 cout << "findBoundaryElementNormals: error: unknown code" << endl;
1589 cout.flush();
1590 exit(0);
1591 }
1592
1593 // Calculate normal (modulo sign):
1594 a[0] = mesh->getNode(v)->getX(0) - mesh->getNode(u)->getX(0);
1595 a[1] = mesh->getNode(v)->getX(1) - mesh->getNode(u)->getX(1);
1596 a[2] = mesh->getNode(v)->getX(2) - mesh->getNode(u)->getX(2);
1597
1598 b[0] = mesh->getNode(w)->getX(0) - mesh->getNode(u)->getX(0);
1599 b[1] = mesh->getNode(w)->getX(1) - mesh->getNode(u)->getX(1);
1600 b[2] = mesh->getNode(w)->getX(2) - mesh->getNode(u)->getX(2);
1601
1602 helpers->crossProduct(a,b,c);
1603 helpers->normalize(c);
1604
1605 surface->setNormal(0, -c[0]);
1606 surface->setNormal(1, -c[1]);
1607 surface->setNormal(2, -c[2]);
1608
1609 // Determine sign:
1610 //----------------
1611
1612 // a) which parent element has bigger index?
1613
1614 e0 = surface->getElementIndex(0);
1615 e1 = surface->getElementIndex(1);
1616
1617 if( (e0<0) && (e1<0) ) {
1618 // both parents unknown
1619 bigger = -1;
1620 } else if(e1<0) {
1621 // e0 known, e1 unknown
1622 bigger = e0;
1623 } else if (e0<0) {
1624 // e0 unknown, e1 known
1625 bigger = e1;
1626 } else {
1627 // both parents known
1628 bigger = e0;
1629 i0 = mesh->getElement(e0)->getIndex();
1630 i1 = mesh->getElement(e1)->getIndex();
1631 if(i1 > i0)
1632 bigger = e1;
1633 }
1634
1635 // b) normal should point to the parent with smaller index:
1636
1637 if(bigger > -1) {
1638
1639 // Compute center point of the surface element:
1640 center_surface[0] = 0.0;
1641 center_surface[1] = 0.0;
1642 center_surface[2] = 0.0;
1643
1644 for(int i=0; i < surface->getNodes(); i++) {
1645 int j = surface->getNodeIndex(i);
1646 node_t *n = mesh->getNode(j);
1647 center_surface[0] += n->getX(0);
1648 center_surface[1] += n->getX(1);
1649 center_surface[2] += n->getX(2);
1650 }
1651
1652 center_surface[0] /= (double)(surface->getNodes());
1653 center_surface[1] /= (double)(surface->getNodes());
1654 center_surface[2] /= (double)(surface->getNodes());
1655
1656 element_t *e = mesh->getElement(bigger);
1657
1658 // compute center point of the parent element:
1659 center_element[0] = 0.0;
1660 center_element[1] = 0.0;
1661 center_element[2] = 0.0;
1662
1663 for(int i=0; i < e->getNodes(); i++) {
1664 int j = e->getNodeIndex(i);
1665 node_t *n = mesh->getNode(j);
1666 center_element[0] += n->getX(0);
1667 center_element[1] += n->getX(1);
1668 center_element[2] += n->getX(2);
1669 }
1670
1671 center_element[0] /= (double)(e->getNodes());
1672 center_element[1] /= (double)(e->getNodes());
1673 center_element[2] /= (double)(e->getNodes());
1674
1675 // difference of the centers:
1676 center_difference[0] = center_element[0] - center_surface[0];
1677 center_difference[1] = center_element[1] - center_surface[1];
1678 center_difference[2] = center_element[2] - center_surface[2];
1679
1680 // dot product must be negative
1681 double dp = center_difference[0]*c[0]
1682 + center_difference[1]*c[1]
1683 + center_difference[2]*c[2];
1684
1685 if(dp > 0.0) {
1686 surface->setNormal(0, -surface->getNormal(0));
1687 surface->setNormal(1, -surface->getNormal(1));
1688 surface->setNormal(2, -surface->getNormal(2));
1689
1690 // change orientation of the surface element:
1691 if(surface->getCode() == 303) {
1692 int tmp = surface->getNodeIndex(1);
1693 surface->setNodeIndex(1, surface->getNodeIndex(2));
1694 surface->setNodeIndex(2, tmp);
1695
1696 } else if(surface->getCode() == 404) {
1697 int tmp = surface->getNodeIndex(1);
1698 surface->setNodeIndex(1, surface->getNodeIndex(3));
1699 surface->setNodeIndex(3, tmp);
1700
1701 } else if(surface->getCode() == 306) {
1702 int tmp = surface->getNodeIndex(1);
1703 surface->setNodeIndex(1, surface->getNodeIndex(2));
1704 surface->setNodeIndex(2, tmp);
1705 tmp = surface->getNodeIndex(3);
1706 surface->setNodeIndex(3, surface->getNodeIndex(5));
1707 surface->setNodeIndex(5, tmp);
1708
1709 } else if(surface->getCode() == 408) {
1710 int tmp = surface->getNodeIndex(1);
1711 surface->setNodeIndex(1, surface->getNodeIndex(3));
1712 surface->setNodeIndex(3, tmp);
1713 tmp = surface->getNodeIndex(4);
1714 surface->setNodeIndex(4, surface->getNodeIndex(7));
1715 surface->setNodeIndex(7, tmp);
1716 tmp = surface->getNodeIndex(5);
1717 surface->setNodeIndex(5, surface->getNodeIndex(6));
1718 surface->setNodeIndex(6, tmp);
1719
1720 } else {
1721 cout << "findSurfaceElementNormals: error: unable to change element orientation" << endl;
1722 cout.flush();
1723 exit(0);
1724 }
1725 }
1726 }
1727 }
1728
1729 for( int i=0; i < mesh->getSurfaces(); i++ )
1730 {
1731 surface_t *surface = mesh->getSurface(i);
1732 int n = surface->getCode() / 100;
1733 for(int j=0; j<n; j++ )
1734 {
1735 surface->setVertexNormalVec(j, surface->getNormalVec());
1736 }
1737 }
1738
1739
1740 // List surfaces connected to nodes
1741 //class n_s_t {
1742 //public:
1743 //int index;
1744 //n_s_t *next;
1745 //} n_s[mesh->getNodes()];
1746
1747 class n_s_t {
1748 public:
1749 int index;
1750 n_s_t *next;
1751 };
1752
1753 n_s_t *n_s = new n_s_t[mesh->getNodes()];
1754
1755 for( int i=0; i<mesh->getNodes(); i++ )
1756 {
1757 n_s[i].index = -1;
1758 n_s[i].next = NULL;
1759 }
1760
1761 for( int i=0; i < mesh->getSurfaces(); i++ )
1762 {
1763 surface_t *surface = mesh->getSurface(i);
1764 int n = surface->getCode() / 100;
1765
1766 for( int j=0; j<n; j++ )
1767 {
1768 n_s_t *p = &n_s[surface->getNodeIndex(j)];
1769 if ( p->index >= 0 ) {
1770 n_s_t *q = new n_s_t;
1771 q->next = p->next;
1772 p->next = q;
1773 q->index = i;
1774 } else p->index = i;
1775 }
1776 }
1777
1778 // average normals over surfaces connected to vertices if
1779 // normals within the limit_angle:
1780 double limit_angle = cos(50.*3.14159/180.);
1781
1782 for( int i=0; i < mesh->getSurfaces(); i++ )
1783 {
1784 surface_t *surf1 = mesh->getSurface(i);
1785 int n = surf1->getCode() / 100;
1786
1787 for( int j=0; j<n; j++ )
1788 {
1789 n_s_t *p = &n_s[surf1->getNodeIndex(j)];
1790 for( ; p && p->index>=0; p=p->next )
1791 {
1792 if ( p->index == i ) continue;
1793
1794 surface_t *surf2 = mesh->getSurface(p->index);
1795 double s = 0.;
1796
1797 s += surf1->getNormal(0) * surf2->getNormal(0);
1798 s += surf1->getNormal(1) * surf2->getNormal(1);
1799 s += surf1->getNormal(2) * surf2->getNormal(2);
1800 if ( fabs(s) > limit_angle )
1801 {
1802 if ( s > 0 ) {
1803 surf1->addVertexNormalVec(j, surf2->getNormalVec());
1804 } else {
1805 surf1->subVertexNormalVec(j, surf2->getNormalVec());
1806 }
1807 }
1808 }
1809 }
1810 }
1811
1812 // delete lists:
1813 for( int i=0; i<mesh->getNodes(); i++ )
1814 {
1815 n_s_t *p=&n_s[i], *q;
1816 p = p->next;
1817 while( p )
1818 {
1819 q = p->next;
1820 delete p;
1821 p = q;
1822 }
1823 }
1824
1825 // And finally normalize:
1826 for( int i=0; i < mesh->getSurfaces(); i++ )
1827 {
1828 surface_t *surface = mesh->getSurface(i);
1829 int n = surface->getCode() / 100;
1830
1831 for(int j=0; j<n; j++ )
1832 {
1833 double* v = surface->getVertexNormalVec(j);
1834 double s = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
1835 if ( s != 0 ) {
1836 s = sqrt(s);
1837 v[0] /= s;
1838 v[1] /= s;
1839 v[2] /= s;
1840 }
1841 }
1842 }
1843
1844
1845 delete helpers;
1846 }
1847
1848
1849
1850 // Increase elementorder from 1 to 2
1851 //----------------------------------------------------------------------------
increaseElementOrder(mesh_t * mesh)1852 void Meshutils::increaseElementOrder(mesh_t *mesh)
1853 {
1854 class hashEntry {
1855 public:
1856 int node1;
1857 int noedge;
1858 hashEntry *next;
1859 };
1860
1861 if((mesh->getElements() == 0) && (mesh->getSurfaces() == 0)) return;
1862
1863 int keys = mesh->getNodes();
1864 hashEntry *hash = new hashEntry[keys];
1865 hashEntry *h;
1866
1867 for(int i=0; i<keys; i++) {
1868 h = &hash[i];
1869 h->node1 = UNKNOWN;
1870 h->noedge = UNKNOWN;
1871 h->next = NULL;
1872 }
1873
1874 static int familylinnodes[9] = {0, 1, 2, 3, 4, 4, 5, 6, 8};
1875 static int familyedges[9] = {0, 0, 1, 3, 4, 6, 8, 9, 12};
1876
1877 static int edgemap8[][2] = {{0,1}, {1,2}, {2,3}, {3,0}, {0,4}, {1,5}, {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {7,4}};
1878 static int edgemap7[][2] = {{0,1}, {1,2}, {2,0}, {0,3}, {1,4}, {2,5}, {3,4}, {4,5}, {5,6}};
1879 static int edgemap6[][2] = {{0,1}, {1,2}, {2,3}, {3,0}, {0,4}, {1,4}, {2,4}, {3,4}};
1880 static int edgemap5[][2] = {{0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3}};
1881 static int edgemap4[][2] = {{0,1}, {1,2}, {2,3}, {3,0}};
1882 static int edgemap3[][2] = {{0,1}, {1,2}, {2,0}};
1883 static int edgemap2[][2] = {{0,1}};
1884
1885
1886 // First go through elements and find all new edges introduced by the 2nd order
1887 int noedges = 0;
1888 for(int i=0; i < mesh->getElements() + mesh->getSurfaces(); i++) {
1889
1890 element_t *e;
1891 if(i < mesh->getElements())
1892 e = mesh->getElement(i);
1893 else
1894 e = mesh->getSurface(i - mesh->getElements());
1895
1896 int *edgemap = NULL;
1897 int family = e->getCode() / 100;
1898 int edges = familyedges[family];
1899
1900 // Skip undtermined and nonlinear element
1901 if((e->getNodes() < 2) || (e->getNodes() > familylinnodes[family])) continue;
1902
1903 for(int f=0; f<edges; f++) {
1904 if(family == 3)
1905 edgemap = &edgemap3[f][0];
1906 else if(family == 4)
1907 edgemap = &edgemap4[f][0];
1908 if(family == 5)
1909 edgemap = &edgemap5[f][0];
1910 else if(family == 6)
1911 edgemap = &edgemap6[f][0];
1912 else if(family == 7)
1913 edgemap = &edgemap7[f][0];
1914 else if(family == 8)
1915 edgemap = &edgemap8[f][0];
1916
1917 int n0 = e->getNodeIndex(edgemap[0]);
1918 int n1 = e->getNodeIndex(edgemap[1]);
1919
1920 if(n0 < 0 || n1 < 0) continue;
1921
1922 // Order indexes in an increasing order
1923 if(n0 > n1) {
1924 int tmp = n0;
1925 n0 = n1;
1926 n1 = tmp;
1927 }
1928
1929 h = &hash[n0];
1930 bool found = false;
1931 while(h->next) {
1932 if(h->node1 == n1) {
1933 found = true;
1934 break;
1935 }
1936 h = h->next;
1937 }
1938
1939 if(!found) {
1940 h->node1 = n1;
1941 h->noedge = noedges;
1942
1943 h->next = new hashEntry;
1944 h = h->next;
1945 h->node1 = UNKNOWN;
1946 h->noedge = UNKNOWN;
1947 h->next = NULL;
1948 noedges++;
1949 }
1950 }
1951 }
1952
1953 cout << "Found " << noedges << " edges in the mesh " << endl;
1954
1955 if(noedges == 0) {
1956 return;
1957 delete [] hash;
1958 }
1959
1960 // Then redifine the mesh using the additional nodes
1961 int quadnodes = mesh->getNodes() + noedges;
1962 node_t *quadnode = new node_t[quadnodes];
1963
1964 // Copy linear nodes
1965 for(int i = 0; i < mesh->getNodes(); i++) {
1966 quadnode[i].setX(0, mesh->getNode(i)->getX(0));
1967 quadnode[i].setX(1, mesh->getNode(i)->getX(1));
1968 quadnode[i].setX(2, mesh->getNode(i)->getX(2));
1969 quadnode[i].setIndex(mesh->getNode(i)->getIndex());
1970 }
1971
1972 // Quadratic nodes are taken to be means of the linear nodes
1973 for(int i=0; i<keys; i++) {
1974 int j = 0;
1975 h = &hash[i];
1976 while(h->next){
1977 if(h->noedge > UNKNOWN) {
1978 j++;
1979 int n0 = i;
1980 int n1 = h->node1;
1981 int j = mesh->getNodes() + h->noedge;
1982 quadnode[j].setX(0, 0.5*(quadnode[n0].getX(0) + quadnode[n1].getX(0)));
1983 quadnode[j].setX(1, 0.5*(quadnode[n0].getX(1) + quadnode[n1].getX(1)));
1984 quadnode[j].setX(2, 0.5*(quadnode[n0].getX(2) + quadnode[n1].getX(2)));
1985 quadnode[j].setIndex(UNKNOWN);
1986 }
1987 h = h->next;
1988 }
1989 }
1990 mesh->deleteNodeArray();
1991 mesh->setNodeArray(quadnode);
1992
1993 cout << "Set " << quadnodes << " additional nodes" << endl;
1994
1995
1996 int tmpnode[27];
1997
1998 for(int i=0; i < mesh->getElements() + mesh->getSurfaces() + mesh->getEdges(); i++) {
1999
2000 element_t *e;
2001 if(i < mesh->getElements())
2002 e = mesh->getElement(i);
2003 else if (i < mesh->getElements() + mesh->getSurfaces())
2004 e = mesh->getSurface(i - mesh->getElements());
2005 else
2006 e = mesh->getEdge(i - mesh->getElements() - mesh->getSurfaces());
2007
2008 int *edgemap = NULL;
2009 int family = e->getCode() / 100;
2010 int edges = familyedges[family];
2011
2012 /* Not a linear element */
2013 if((e->getNodes() < 2) || (e->getNodes() > familylinnodes[family])) continue;
2014
2015 int linnodes = e->getNodes();
2016 for(int j=0;j<linnodes;j++)
2017 tmpnode[j] = e->getNodeIndex(j);
2018
2019 e->deleteNodeIndexes();
2020 e->setCode(e->getCode() + edges);
2021 e->setNodes(e->getNodes() + edges);
2022 e->newNodeIndexes(e->getNodes());
2023 for(int j=0;j<linnodes;j++)
2024 e->setNodeIndex(j, tmpnode[j]);
2025
2026
2027 for(int f=0; f<edges; f++) {
2028
2029 if(family == 2)
2030 edgemap = &edgemap2[f][0];
2031 if(family == 3)
2032 edgemap = &edgemap3[f][0];
2033 else if(family == 4)
2034 edgemap = &edgemap4[f][0];
2035 if(family == 5)
2036 edgemap = &edgemap5[f][0];
2037 else if(family == 6)
2038 edgemap = &edgemap6[f][0];
2039 else if(family == 7)
2040 edgemap = &edgemap7[f][0];
2041 else if(family == 8)
2042 edgemap = &edgemap8[f][0];
2043
2044 int n0 = e->getNodeIndex(edgemap[0]);
2045 int n1 = e->getNodeIndex(edgemap[1]);
2046 if((n0 < 0) || (n1 < 0)) continue;
2047
2048 // Order indexes in an increasing order
2049 if(n0 > n1) {
2050 int tmp = n0;
2051 n0 = n1;
2052 n1 = tmp;
2053 }
2054
2055 h = &hash[n0];
2056 bool found = false;
2057 while(h->next) {
2058 if(h->node1 == n1) {
2059 found = true;
2060 e->setNodeIndex(linnodes+f, h->noedge + mesh->getNodes());
2061 break;
2062 }
2063 h = h->next;
2064 }
2065 if(!found) {
2066 cout << "All edges should be found!? " << endl;
2067 }
2068 }
2069 }
2070 mesh->setNodes(quadnodes);
2071
2072 delete [] hash;
2073 }
2074
2075
2076
2077 // Decrease elementorder from >1 to 1
2078 //----------------------------------------------------------------------------
decreaseElementOrder(mesh_t * mesh)2079 void Meshutils::decreaseElementOrder(mesh_t *mesh)
2080 {
2081 if((mesh->getElements() == 0) && (mesh->getSurfaces() == 0)) return;
2082
2083 static int familylinnodes[9] = {0, 1, 2, 3, 4, 4, 5, 6, 8};
2084
2085 int *activenodes = new int[mesh->getNodes()];
2086 for(int i=0; i < mesh->getNodes(); i++)
2087 activenodes[i] = -1;
2088
2089 // int noedges = 0;
2090 for(int i=0; i < mesh->getElements() + mesh->getSurfaces(); i++) {
2091
2092 element_t *e;
2093 if(i < mesh->getElements())
2094 e = mesh->getElement(i);
2095 else
2096 e = mesh->getSurface(i - mesh->getElements());
2097
2098 int family = e->getCode() / 100;
2099 int linnodes = familylinnodes[family];
2100
2101 for(int j=0;j<linnodes;j++)
2102 activenodes[e->getNodeIndex(j)] = 0;
2103 }
2104
2105 int linnodes = 0;
2106 for(int i=0; i < mesh->getNodes(); i++)
2107 if(activenodes[i] > -1) activenodes[i] = linnodes++;
2108
2109 cout << "Setting " << linnodes << " linear nodes of the total " << mesh->getNodes() << " nodes" << endl;
2110 if(linnodes == mesh->getNodes()) return;
2111
2112
2113 // Copy linear nodes
2114 node_t *linnode = new node_t[linnodes];
2115 for(int i=0; i < mesh->getNodes(); i++) {
2116 int j = activenodes[i];
2117 if(j > -1) {
2118 linnode[j].setX(0, mesh->getNode(i)->getX(0));
2119 linnode[j].setX(1, mesh->getNode(i)->getX(1));
2120 linnode[j].setX(2, mesh->getNode(i)->getX(2));
2121 }
2122 }
2123 mesh->deleteNodeArray();
2124 mesh->setNodeArray(linnode);
2125 mesh->setNodes(linnodes);
2126
2127 int tmpnode[8];
2128
2129 for(int i=0; i < mesh->getElements() + mesh->getSurfaces() + mesh->getEdges(); i++) {
2130
2131 element_t *e;
2132 if(i < mesh->getElements())
2133 e = mesh->getElement(i);
2134 else if (i < mesh->getElements() + mesh->getSurfaces())
2135 e = mesh->getSurface(i - mesh->getElements());
2136 else
2137 e = mesh->getEdge(i - mesh->getElements() - mesh->getSurfaces());
2138
2139 int family = e->getCode() / 100;
2140 int linnodes = familylinnodes[family];
2141
2142 for(int j=0;j<linnodes;j++)
2143 tmpnode[j] = e->getNodeIndex(j);
2144
2145 e->deleteNodeIndexes();
2146 e->setCode(100*family + linnodes);
2147 e->setNodes(linnodes);
2148 e->newNodeIndexes(e->getNodes());
2149
2150 for(int j=0;j<linnodes;j++)
2151 e->setNodeIndex(j, activenodes[tmpnode[j]]);
2152 }
2153 }
2154
2155 // Clean up hanging sharp edges (for visualization)...
2156 //----------------------------------------------------------------------------
cleanHangingSharpEdges(mesh_t * mesh)2157 int Meshutils::cleanHangingSharpEdges(mesh_t *mesh)
2158 {
2159 int count_total = 0;
2160 int count_removed = 0;
2161
2162 if(mesh->getEdges() == 0)
2163 return 0;
2164
2165 for(int i = 0; i < mesh->getEdges(); i++) {
2166 edge_t *edge = mesh->getEdge(i);
2167 if(edge->isSharp()) {
2168 count_total++;
2169 if(edge->getSurfaces() == 2) {
2170 // has two parents
2171 int i0 = edge->getSurfaceIndex(0);
2172 int i1 = edge->getSurfaceIndex(1);
2173 surface_t *s0 = NULL;
2174 surface_t *s1 = NULL;
2175 if(i0 >= 0)
2176 s0 = mesh->getSurface(i0);
2177 if(i1 >= 0)
2178 s1 = mesh->getSurface(i1);
2179 if((s0 != NULL) && (s1 != NULL)) {
2180 int index0 = s0->getIndex();
2181 int index1 = s1->getIndex();
2182 if(index0 == index1) {
2183 // remove
2184 count_removed++;
2185 edge->setSharp(false);
2186 }
2187 }
2188 }
2189 }
2190 }
2191
2192 return count_removed;
2193 }
2194