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 mesh_t (Elmer mesh structure) *
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 #include <iostream>
41 #include <fstream>
42 #include "meshtype.h"
43 using namespace std;
44
45 // node_t
46 //-----------------------------------------------------------------------------
node_t()47 node_t::node_t()
48 {
49 }
50
~node_t()51 node_t::~node_t()
52 {
53 }
54
setXvec(double * y)55 void node_t::setXvec(double* y)
56 {
57 this->x[0] = y[0];
58 this->x[1] = y[1];
59 this->x[2] = y[2];
60 }
61
getXvec()62 double* node_t::getXvec()
63 {
64 return &this->x[0];
65 }
66
setX(int n,double y)67 void node_t::setX(int n, double y)
68 {
69 this->x[n] = y;
70 }
71
getX(int n) const72 double node_t::getX(int n) const
73 {
74 return this->x[n];
75 }
76
setIndex(int n)77 void node_t::setIndex(int n)
78 {
79 this->index = n;
80
81 }
82
getIndex() const83 int node_t::getIndex() const
84 {
85 return this->index;
86 }
87
88 // element_t
89 //-----------------------------------------------------------------------------
element_t()90 element_t::element_t()
91 {
92 }
93
~element_t()94 element_t::~element_t()
95 {
96 }
97
setNature(int n)98 void element_t::setNature(int n)
99 {
100 this->nature = n;
101 }
102
getNature() const103 int element_t::getNature() const
104 {
105 return this->nature;
106 }
107
setCode(int n)108 void element_t::setCode(int n)
109 {
110 this->code = n;
111 }
112
getCode() const113 int element_t::getCode() const
114 {
115 return this->code;
116 }
117
setNodes(int n)118 void element_t::setNodes(int n)
119 {
120 this->nodes = n;
121 }
122
getNodes() const123 int element_t::getNodes() const
124 {
125 return this->nodes;
126 }
127
setIndex(int n)128 void element_t::setIndex(int n)
129 {
130 this->index = n;
131 }
132
getIndex() const133 int element_t::getIndex() const
134 {
135 return this->index;
136 }
137
setSelected(int n)138 void element_t::setSelected(int n)
139 {
140 this->selected = n;
141 }
142
getSelected() const143 int element_t::getSelected() const
144 {
145 return this->selected;
146 }
147
getNodeIndex(int n) const148 int element_t::getNodeIndex(int n) const
149 {
150 return this->node[n];
151 }
152
setNodeIndex(int m,int n)153 void element_t::setNodeIndex(int m, int n)
154 {
155 this->node[m] = n;
156 }
157
newNodeIndexes(int n)158 void element_t::newNodeIndexes(int n)
159 {
160 this->node = new int[n];
161 }
162
deleteNodeIndexes()163 void element_t::deleteNodeIndexes()
164 {
165 delete [] this->node;
166 }
167
getNodeIndexes() const168 int* element_t::getNodeIndexes() const
169 {
170 return this->node;
171 }
172
173 // point_t
174 //-----------------------------------------------------------------------------
point_t()175 point_t::point_t()
176 {
177 }
178
~point_t()179 point_t::~point_t()
180 {
181 }
182
setSharp(bool b)183 void point_t::setSharp(bool b)
184 {
185 this->sharp_point = b;
186 }
187
isSharp() const188 bool point_t::isSharp() const
189 {
190 return this->sharp_point;
191 }
192
setEdges(int n)193 void point_t::setEdges(int n)
194 {
195 this->edges = n;
196 }
197
getEdges() const198 int point_t::getEdges() const
199 {
200 return this->edges;
201 }
202
getEdgeIndex(int n) const203 int point_t::getEdgeIndex(int n) const
204 {
205 return this->edge[n];
206 }
207
setEdgeIndex(int m,int n)208 void point_t::setEdgeIndex(int m, int n)
209 {
210 this->edge[m] = n;
211 }
212
newEdgeIndexes(int n)213 void point_t::newEdgeIndexes(int n)
214 {
215 this->edge = new int[n];
216 }
217
deleteEdgeIndexes()218 void point_t::deleteEdgeIndexes()
219 {
220 delete [] this->edge;
221 }
222
223 // edge_t
224 //-----------------------------------------------------------------------------
edge_t()225 edge_t::edge_t()
226 {
227 }
228
~edge_t()229 edge_t::~edge_t()
230 {
231 }
232
setSharp(bool b)233 void edge_t::setSharp(bool b)
234 {
235 this->sharp_edge = b;
236 }
237
isSharp() const238 bool edge_t::isSharp() const
239 {
240 return this->sharp_edge;
241 }
242
setPoints(int n)243 void edge_t::setPoints(int n)
244 {
245 this->points = n;
246 }
247
getPoints() const248 int edge_t::getPoints() const
249 {
250 return this->points;
251 }
252
setPointIndex(int m,int n)253 void edge_t::setPointIndex(int m, int n)
254 {
255 this->point[m] = n;
256 }
257
getPointIndex(int n) const258 int edge_t::getPointIndex(int n) const
259 {
260 return this->point[n];
261 }
262
newPointIndexes(int n)263 void edge_t::newPointIndexes(int n)
264 {
265 this->point = new int[n];
266 }
267
deletePointIndexes()268 void edge_t::deletePointIndexes()
269 {
270 delete [] this->point;
271 }
272
setSurfaces(int n)273 void edge_t::setSurfaces(int n)
274 {
275 this->surfaces = n;
276 }
277
getSurfaces() const278 int edge_t::getSurfaces() const
279 {
280 return this->surfaces;
281 }
282
setSurfaceIndex(int m,int n)283 void edge_t::setSurfaceIndex(int m, int n)
284 {
285 this->surface[m] = n;
286 }
287
getSurfaceIndex(int n) const288 int edge_t::getSurfaceIndex(int n) const
289 {
290 return this->surface[n];
291 }
292
newSurfaceIndexes(int n)293 void edge_t::newSurfaceIndexes(int n)
294 {
295 this->surface = new int[n];
296 }
297
deleteSurfaceIndexes()298 void edge_t::deleteSurfaceIndexes()
299 {
300 delete [] this->surface;
301 }
302
303 // surface_t
304 //-----------------------------------------------------------------------------
surface_t()305 surface_t::surface_t()
306 {
307 }
308
~surface_t()309 surface_t::~surface_t()
310 {
311 }
312
setEdges(int n)313 void surface_t::setEdges(int n)
314 {
315 this->edges = n;
316 }
317
getEdges() const318 int surface_t::getEdges() const
319 {
320 return this->edges;
321 }
322
setEdgeIndex(int m,int n)323 void surface_t::setEdgeIndex(int m, int n)
324 {
325 this->edge[m] = n;
326 }
327
getEdgeIndex(int n) const328 int surface_t::getEdgeIndex(int n) const
329 {
330 return this->edge[n];
331 }
332
newEdgeIndexes(int n)333 void surface_t::newEdgeIndexes(int n)
334 {
335 this->edge = new int[n];
336 }
337
deleteEdgeIndexes()338 void surface_t::deleteEdgeIndexes()
339 {
340 delete [] this->edge;
341 }
342
setElements(int n)343 void surface_t::setElements(int n)
344 {
345 this->elements = n;
346 }
347
getElements() const348 int surface_t::getElements() const
349 {
350 return this->elements;
351 }
352
setElementIndex(int m,int n)353 void surface_t::setElementIndex(int m, int n)
354 {
355 this->element[m] = n;
356 }
357
getElementIndex(int n) const358 int surface_t::getElementIndex(int n) const
359 {
360 return this->element[n];
361 }
362
newElementIndexes(int n)363 void surface_t::newElementIndexes(int n)
364 {
365 this->element = new int[n];
366 }
367
deleteElementIndexes()368 void surface_t::deleteElementIndexes()
369 {
370 delete [] this->element;
371 }
372
getNormalVec()373 double* surface_t::getNormalVec()
374 {
375 return &this->normal[0];
376 }
377
setNormalVec(double * d)378 void surface_t::setNormalVec(double* d)
379 {
380 this->normal[0] = d[0];
381 this->normal[1] = d[1];
382 this->normal[2] = d[2];
383 }
384
getNormal(int n) const385 double surface_t::getNormal(int n) const
386 {
387 return this->normal[n];
388 }
389
setNormal(int n,double d)390 void surface_t::setNormal(int n, double d)
391 {
392 this->normal[n] = d;
393 }
394
setVertexNormalVec(int n,double * d)395 void surface_t::setVertexNormalVec(int n, double* d)
396 {
397 this->vertex_normals[n][0] = d[0];
398 this->vertex_normals[n][1] = d[1];
399 this->vertex_normals[n][2] = d[2];
400 }
401
addVertexNormalVec(int n,double * d)402 void surface_t::addVertexNormalVec(int n, double* d)
403 {
404 this->vertex_normals[n][0] += d[0];
405 this->vertex_normals[n][1] += d[1];
406 this->vertex_normals[n][2] += d[2];
407 }
408
subVertexNormalVec(int n,double * d)409 void surface_t::subVertexNormalVec(int n, double* d)
410 {
411 this->vertex_normals[n][0] -= d[0];
412 this->vertex_normals[n][1] -= d[1];
413 this->vertex_normals[n][2] -= d[2];
414 }
415
getVertexNormalVec(int n)416 double* surface_t::getVertexNormalVec(int n)
417 {
418 return &this->vertex_normals[n][0];
419 }
420
421 // mesh_t
422 //-----------------------------------------------------------------------------
mesh_t()423 mesh_t::mesh_t()
424 {
425 this->setDefaults();
426 }
427
~mesh_t()428 mesh_t::~mesh_t()
429 {
430 }
431
isUndefined() const432 bool mesh_t::isUndefined() const
433 {
434 if((cdim < 0) || (dim < 0) || (nodes < 1))
435 return true;
436
437 return false;
438 }
439
clear()440 void mesh_t::clear()
441 {
442 delete [] element;
443 delete [] surface;
444 delete [] edge;
445 delete [] point;
446 delete [] node;
447
448 setDefaults();
449 }
450
setDefaults()451 void mesh_t::setDefaults()
452 {
453 cdim = -1;
454 dim = -1;
455 nodes = 0;
456 node = 0;
457 points = 0;
458 point = 0;
459 edges = 0;
460 edge = 0;
461 surfaces = 0;
462 surface = 0;
463 elements = 0;
464 element = 0;
465 }
466
467 // Load Elmer mesh files and populate mesh structures
468 //---------------------------------------------------------------------------
load(char * dirName)469 bool mesh_t::load(char* dirName)
470 {
471 char fileName[1024];
472 ifstream mesh_header;
473 ifstream mesh_nodes;
474 ifstream mesh_elements;
475 ifstream mesh_boundary;
476
477 // Header:
478 //--------
479 sprintf(fileName, "%s/mesh.header", dirName);
480 mesh_header.open(fileName);
481
482 if(!mesh_header.is_open()) {
483 cout << "Mesh: load: unable to open " << fileName << endl;
484 return false;
485 }
486
487 int nodes, elements, boundaryelements, types, type, ntype;
488
489 mesh_header >> nodes >> elements >> boundaryelements;
490 mesh_header >> types;
491
492 int elements_zero_d = 0;
493 int elements_one_d = 0;
494 int elements_two_d = 0;
495 int elements_three_d = 0;
496
497 for(int i = 0; i < types; i++) {
498 mesh_header >> type >> ntype;
499
500 switch(type/100) {
501 case 1:
502 elements_zero_d += ntype;
503 break;
504 case 2:
505 elements_one_d += ntype;
506 break;
507 case 3:
508 case 4:
509 elements_two_d += ntype;
510 break;
511 case 5:
512 case 6:
513 case 7:
514 case 8:
515 elements_three_d += ntype;
516 break;
517 default:
518 cout << "Unknown element family (possibly not implamented)" << endl;
519 cout.flush();
520 return false;
521 // exit(0);
522 }
523 }
524
525 cout << "Summary:" << endl;
526 cout << "Nodes: " << nodes << endl;
527 cout << "Point elements: " << elements_zero_d << endl;
528 cout << "Edge elements: " << elements_one_d << endl;
529 cout << "Surface elements: " << elements_two_d << endl;
530 cout << "Volume elements: " << elements_three_d << endl;
531 cout.flush();
532
533 // Set mesh dimension:
534 this->dim = 0;
535
536 if(elements_one_d > 0)
537 this->dim = 1;
538
539 if(elements_two_d > 0)
540 this->dim = 2;
541
542 if(elements_three_d > 0)
543 this->dim = 3;
544
545 this->nodes = nodes;
546 node = new node_t[nodes];
547
548 this->points = elements_zero_d;
549 point = new point_t[this->points];
550
551 this->edges = elements_one_d;
552 edge = new edge_t[this->edges];
553
554 this->surfaces = elements_two_d;
555 surface = new surface_t[this->surfaces];
556
557 this->elements = elements_three_d;
558 element = new element_t[this->elements];
559
560 mesh_header.close();
561
562 // Nodes:
563 //-------
564 sprintf(fileName, "%s/mesh.nodes", dirName);
565 mesh_nodes.open(fileName);
566
567 if(!mesh_nodes.is_open()) {
568 cout << "Mesh: load: unable to open " << fileName << endl;
569 return false;
570 }
571
572 int number, index;
573 double x, y, z;
574
575 for(int i = 0; i < nodes; i++) {
576 node_t *node = &this->node[i];
577 mesh_nodes >> number >> index >> x >> y >> z;
578 node->setX(0, x);
579 node->setX(1, y);
580 node->setX(2, z);
581 node->setIndex(index);
582 }
583
584 mesh_nodes.close();
585
586 // Elements:
587 //----------
588 sprintf(fileName, "%s/mesh.elements", dirName);
589 mesh_elements.open(fileName);
590
591 if(!mesh_elements.is_open()) {
592 cout << "Mesh: load: unable to open " << fileName << endl;
593 return false;
594 }
595
596 int current_point = 0;
597 int current_edge = 0;
598 int current_surface = 0;
599 int current_element = 0;
600
601 point_t *point = NULL;
602 edge_t *edge = NULL;
603 surface_t *surface = NULL;
604 element_t *element = NULL;
605
606 for(int i = 0; i < elements; i++) {
607 mesh_elements >> number >> index >> type;
608
609 switch(type/100) {
610 case 1:
611 point = &this->point[current_point++];
612 point->setNature(PDE_BULK);
613 point->setIndex(index);
614 point->setCode(type);
615 point->setNodes(point->getCode() % 100);
616 point->newNodeIndexes(point->getNodes());
617 for(int j = 0; j < point->getNodes(); j++) {
618 int k;
619 mesh_elements >> k;
620 point->setNodeIndex(j, k-1);
621 }
622 point->setEdges(2);
623 point->newEdgeIndexes(point->getEdges());
624 point->setEdgeIndex(0, -1);
625 point->setEdgeIndex(1, -1);
626 break;
627
628 case 2:
629 edge = &this->edge[current_edge++];
630 edge->setNature(PDE_BULK);
631 edge->setIndex(index);
632 edge->setCode(type);
633 edge->setNodes(edge->getCode() % 100);
634 edge->newNodeIndexes(edge->getNodes());
635 for(int j = 0; j < edge->getNodes(); j++) {
636 int k;
637 mesh_elements >> k;
638 edge->setNodeIndex(j, k-1);
639 }
640 edge->setSurfaces(0);
641 edge->newSurfaceIndexes(edge->getSurfaces());
642 edge->setSurfaceIndex(0, -1);
643 edge->setSurfaceIndex(1, -1);
644
645 break;
646
647 case 3:
648 case 4:
649 surface = &this->surface[current_surface++];
650 surface->setNature(PDE_BULK);
651 surface->setIndex(index);
652 surface->setCode(type);
653 surface->setNodes(surface->getCode() % 100);
654 surface->newNodeIndexes(surface->getNodes());
655 for(int j = 0; j < surface->getNodes(); j++) {
656 int k;
657 mesh_elements >> k;
658 surface->setNodeIndex(j, k-1);
659 }
660 surface->setEdges((int)(surface->getCode() / 100));
661 surface->newEdgeIndexes(surface->getEdges());
662 for(int j = 0; j < surface->getEdges(); j++)
663 surface->setEdgeIndex(j, -1);
664 surface->setElements(2);
665 surface->newElementIndexes(surface->getElements());
666 surface->setElementIndex(0, -1);
667 surface->setElementIndex(1, -1);
668
669 break;
670
671 case 5:
672 case 6:
673 case 7:
674 case 8:
675 element = &this->element[current_element++];
676 element->setNature(PDE_BULK);
677 element->setIndex(index);
678 element->setCode(type);
679 element->setNodes(element->getCode() % 100);
680 element->newNodeIndexes(element->getNodes());
681 for(int j = 0; j < element->getNodes(); j++) {
682 int k;
683 mesh_elements >> k;
684 element->setNodeIndex(j, k-1);
685 }
686 break;
687
688 default:
689 cout << "Unknown element type (possibly not implemented" << endl;
690 cout.flush();
691 return false;
692 // exit(0);
693 }
694 }
695
696 mesh_elements.close();
697
698 // Boundary elements:
699 //-------------------
700 sprintf(fileName, "%s/mesh.boundary", dirName);
701 mesh_boundary.open(fileName);
702
703 if(!mesh_boundary.is_open()) {
704 cout << "Mesh: load: unable to open " << fileName << endl;
705 return false;
706 }
707
708 int parent0, parent1;
709
710 for(int i = 0; i < boundaryelements; i++) {
711 mesh_boundary >> number >> index >> parent0 >> parent1 >> type;
712
713 switch(type/100) {
714 case 1:
715 point = &this->point[current_point++];
716 point->setNature(PDE_BOUNDARY);
717 point->setIndex(index);
718 point->setEdges(2);
719 point->newEdgeIndexes(point->getEdges());
720 point->setEdgeIndex(0, parent0 - 1);
721 point->setEdgeIndex(1, parent0 - 1);
722 point->setCode(type);
723 point->setNodes(point->getCode() % 100);
724 point->newNodeIndexes(point->getNodes());
725 for(int j = 0; j < point->getNodes(); j++) {
726 int k;
727 mesh_boundary >> k;
728 point->setNodeIndex(j, k-1);
729 }
730 break;
731
732 case 2:
733 edge = &this->edge[current_edge++];
734 edge->setNature(PDE_BOUNDARY);
735 edge->setIndex(index);
736 edge->setSurfaces(2);
737 edge->newSurfaceIndexes(edge->getSurfaces());
738 edge->setSurfaceIndex(0, parent0 - 1);
739 edge->setSurfaceIndex(1, parent1 - 1);
740 edge->setCode(type);
741 edge->setNodes(edge->getCode() % 100);
742 edge->newNodeIndexes(edge->getNodes());
743 for(int j = 0; j < edge->getNodes(); j++) {
744 int k;
745 mesh_boundary >> k;
746 edge->setNodeIndex(j, k-1);
747 }
748
749 break;
750
751 case 3:
752 case 4:
753 surface = &this->surface[current_surface++];
754 surface->setNature(PDE_BOUNDARY);
755 surface->setIndex(index);
756 surface->setElements(2);
757 surface->newElementIndexes(surface->getElements());
758 surface->setElementIndex(0, parent0 - 1);
759 surface->setElementIndex(1, parent1 - 1);
760 surface->setCode(type);
761 surface->setNodes(surface->getCode() % 100);
762 surface->newNodeIndexes(surface->getNodes());
763 for(int j = 0; j < surface->getNodes(); j++) {
764 int k;
765 mesh_boundary >> k;
766 surface->setNodeIndex(j, k-1);
767 }
768 surface->setEdges((int)(surface->getCode() / 100));
769 surface->newEdgeIndexes(surface->getEdges());
770 for(int j = 0; j < surface->getEdges(); j++)
771 surface->setEdgeIndex(j, -1);
772
773 break;
774
775 case 5:
776 case 6:
777 case 7:
778 case 8:
779 // these can't be boundary elements
780 break;
781
782 default:
783 break;
784 }
785 }
786
787 mesh_boundary.close();
788
789 this->boundingBox();
790
791 return true;
792 }
793
794 // Save Elmer mesh files and populate mesh structures
795 //---------------------------------------------------------------------------
save(char * dirName)796 bool mesh_t::save(char *dirName)
797 {
798 char fileName[1024];
799 ofstream mesh_header;
800 ofstream mesh_nodes;
801 ofstream mesh_elements;
802 ofstream mesh_boundary;
803
804 // Elmer's elements codes are smaller than 1000
805 int maxcode = 1000;
806 int *bulk_by_type = new int[maxcode];
807 int *boundary_by_type = new int[maxcode];
808
809 for(int i = 0; i < maxcode; i++) {
810 bulk_by_type[i] = 0;
811 boundary_by_type[i] = 0;
812 }
813
814 for(int i = 0; i < elements; i++) {
815 element_t *e = &element[i];
816
817 if(e->getNature() == PDE_BULK)
818 bulk_by_type[e->getCode()]++;
819
820 if(e->getNature() == PDE_BOUNDARY)
821 boundary_by_type[e->getCode()]++;
822 }
823
824 for(int i = 0; i < surfaces; i++) {
825 surface_t *s = &surface[i];
826
827 if(s->getNature() == PDE_BULK)
828 bulk_by_type[s->getCode()]++;
829
830 if(s->getNature() == PDE_BOUNDARY)
831 boundary_by_type[s->getCode()]++;
832 }
833
834 for(int i = 0; i < edges; i++) {
835 edge_t *e = &edge[i];
836
837 if(e->getNature() == PDE_BULK)
838 bulk_by_type[e->getCode()]++;
839
840 if(e->getNature() == PDE_BOUNDARY)
841 boundary_by_type[e->getCode()]++;
842 }
843
844 for(int i = 0; i < points; i++) {
845 point_t *p = &point[i];
846
847 if(p->getNature() == PDE_BULK)
848 bulk_by_type[p->getCode()]++;
849
850 if(p->getNature() == PDE_BOUNDARY)
851 boundary_by_type[p->getCode()]++;
852 }
853
854 int bulk_elements = 0;
855 int boundary_elements = 0;
856 int element_types = 0;
857
858 for(int i = 0; i < maxcode; i++) {
859 bulk_elements += bulk_by_type[i];
860 boundary_elements += boundary_by_type[i];
861
862 if((bulk_by_type[i] > 0) || (boundary_by_type[i] > 0))
863 element_types++;
864 }
865
866 // Header:
867 //---------
868 sprintf(fileName, "%s/mesh.header", dirName);
869
870 mesh_header.open(fileName);
871
872 if(!mesh_header.is_open()) {
873 cout << "Unable to open " << fileName << endl;
874 return false;
875 }
876
877 cout << "Saving " << nodes << " nodes" << endl;
878 cout << "Saving " << bulk_elements << " elements" << endl;
879 cout << "Saving " << boundary_elements << " boundary elements" << endl;
880 cout.flush();
881
882 mesh_header << nodes << " ";
883 mesh_header << bulk_elements << " ";
884 mesh_header << boundary_elements << endl;
885 mesh_header << element_types << endl;
886
887 for(int i = 0; i < maxcode; i++) {
888 int j = bulk_by_type[i] + boundary_by_type[i];
889 if(j > 0)
890 mesh_header << i << " " << j << endl;
891 }
892
893 mesh_header.close();
894
895 // Nodes:
896 //--------
897 sprintf(fileName, "%s/mesh.nodes", dirName);
898
899 mesh_nodes.open(fileName);
900
901 if(!mesh_nodes.is_open()) {
902 cout << "Unable to open " << fileName << endl;
903 return false;
904 }
905
906 for(int i = 0; i < this->nodes; i++) {
907 node_t *node = &this->node[i];
908
909 int ind = node->getIndex();
910
911 mesh_nodes << i+1 << " " << ind << " ";
912 mesh_nodes << node->getX(0) << " ";
913 mesh_nodes << node->getX(1) << " ";
914 mesh_nodes << node->getX(2) << endl;
915 }
916
917 mesh_nodes.close();
918
919
920 // Elements:
921 //----------
922 sprintf(fileName, "%s/mesh.elements", dirName);
923
924 mesh_elements.open(fileName);
925
926 if(!mesh_elements.is_open()) {
927 cout << "Unable to open " << fileName << endl;
928 return false;
929 }
930
931 int current = 0;
932
933 for(int i = 0; i < this->elements; i++) {
934 element_t *e = &this->element[i];
935
936 int ind = e->getIndex();
937
938 if(ind < 1)
939 ind = 1;
940
941 if(e->getNature() == PDE_BULK) {
942 mesh_elements << ++current << " ";
943 mesh_elements << ind << " ";
944 mesh_elements << e->getCode() << " ";
945
946 for(int j = 0; j < e->getNodes(); j++)
947 mesh_elements << e->getNodeIndex(j) + 1 << " ";
948
949 mesh_elements << endl;
950 }
951 }
952
953 for(int i = 0; i < this->surfaces; i++) {
954 surface_t *s = &this->surface[i];
955
956 int ind = s->getIndex();
957
958 if(ind < 1)
959 ind = 1;
960
961 if(s->getNature() == PDE_BULK) {
962 mesh_elements << ++current << " ";
963 mesh_elements << ind << " ";
964 mesh_elements << s->getCode() << " ";
965
966 for(int j = 0; j < s->getNodes(); j++)
967 mesh_elements << s->getNodeIndex(j) + 1 << " ";
968
969 mesh_elements << endl;
970 }
971 }
972
973 for(int i = 0; i < this->edges; i++) {
974 edge_t *e = &this->edge[i];
975
976 int ind = e->getIndex();
977
978 if(ind < 1)
979 ind = 1;
980
981 if(e->getNature() == PDE_BULK) {
982 mesh_elements << ++current << " ";
983 mesh_elements << ind << " ";
984 mesh_elements << e->getCode() << " ";
985
986 for(int j = 0; j < e->getNodes(); j++)
987 mesh_elements << e->getNodeIndex(j) + 1 << " ";
988
989 mesh_elements << endl;
990 }
991 }
992
993 for(int i = 0; i < this->points; i++) {
994 point_t *p = &this->point[i];
995
996 int ind = p->getIndex();
997
998 if(ind < 1)
999 ind = 1;
1000
1001 if(p->getNature() == PDE_BULK) {
1002 mesh_elements << ++current << " ";
1003 mesh_elements << ind << " ";
1004 mesh_elements << p->getCode() << " ";
1005
1006 for(int j = 0; j < p->getNodes(); j++)
1007 mesh_elements << p->getNodeIndex(j) + 1 << " ";
1008
1009 mesh_elements << endl;
1010 }
1011 }
1012
1013 mesh_elements.close();
1014
1015 // Boundary elements:
1016 //-------------------
1017 sprintf(fileName, "%s/mesh.boundary", dirName);
1018
1019 mesh_boundary.open(fileName);
1020
1021 if(!mesh_boundary.is_open()) {
1022 cout << "Unable to open " << fileName << endl;
1023 return false;
1024 }
1025
1026 current = 0;
1027
1028 for(int i = 0; i < this->surfaces; i++) {
1029 surface_t *s = &this->surface[i];
1030
1031 if(s->getNature() == PDE_BULK)
1032 continue;
1033
1034 int e0 = s->getElementIndex(0) + 1;
1035 int e1 = s->getElementIndex(1) + 1;
1036
1037 if(e0 < 0)
1038 e0 = 0;
1039
1040 if(e1 < 0)
1041 e1 = 0;
1042
1043 int ind = s->getIndex();
1044
1045 if(ind < 1)
1046 ind = 1;
1047
1048 if(s->getNature() == PDE_BOUNDARY) {
1049 mesh_boundary << ++current << " ";
1050 mesh_boundary << ind << " ";
1051 mesh_boundary << e0 << " " << e1 << " ";
1052 mesh_boundary << s->getCode() << " ";
1053
1054 for(int j = 0; j < s->getNodes(); j++)
1055 mesh_boundary << s->getNodeIndex(j) + 1 << " ";
1056
1057 mesh_boundary << endl;
1058 }
1059 }
1060
1061
1062 for(int i = 0; i < this->edges; i++) {
1063 edge_t *e = &this->edge[i];
1064
1065 if(e->getNature() == PDE_BULK)
1066 continue;
1067
1068 int s0 = e->getSurfaceIndex(0) + 1;
1069 int s1 = e->getSurfaceIndex(1) + 1;
1070
1071 if(s0 < 0)
1072 s0 = 0;
1073
1074 if(s1 < 0)
1075 s1 = 0;
1076
1077 int ind = e->getIndex();
1078
1079 if(ind < 1)
1080 ind = 1;
1081
1082 if(e->getNature() == PDE_BOUNDARY) {
1083 mesh_boundary << ++current << " ";
1084 mesh_boundary << ind << " ";
1085 mesh_boundary << s0 << " " << s1 << " ";
1086 mesh_boundary << e->getCode() << " ";
1087
1088 for(int j = 0; j < e->getNodes(); j++)
1089 mesh_boundary << e->getNodeIndex(j) + 1 << " ";
1090
1091 mesh_boundary << endl;
1092 }
1093 }
1094
1095 for(int i = 0; i < this->points; i++) {
1096 point_t *p = &this->point[i];
1097
1098 int e0 = p->getEdgeIndex(0) + 1;
1099 int e1 = p->getEdgeIndex(1) + 1;
1100
1101 if(e0 < 0)
1102 e0 = 0;
1103
1104 if(e1 < 0)
1105 e1 = 0;
1106
1107 int ind = p->getIndex();
1108
1109 if(ind < 1)
1110 ind = 1;
1111
1112 if(p->getNature() == PDE_BOUNDARY) {
1113 mesh_boundary << ++current << " ";
1114 mesh_boundary << ind << " ";
1115 mesh_boundary << e0 << " " << e1 << " ";
1116 mesh_boundary << p->getCode() << " ";
1117
1118 for(int j = 0; j < p->getNodes(); j++)
1119 mesh_boundary << p->getNodeIndex(j) + 1 << " ";
1120
1121 mesh_boundary << endl;
1122 }
1123 }
1124
1125 mesh_boundary.close();
1126
1127 delete [] bulk_by_type;
1128 delete [] boundary_by_type;
1129
1130 return true;
1131 }
1132
1133 // Bounding box...
1134 //-----------------------------------------------------------------------------
boundingBox()1135 double* mesh_t::boundingBox()
1136 {
1137 double *result = new double[10];
1138
1139 double xmin = +9e9;
1140 double xmax = -9e9;
1141
1142 double ymin = +9e9;
1143 double ymax = -9e9;
1144
1145 double zmin = +9e9;
1146 double zmax = -9e9;
1147
1148 for(int i=0; i < this->nodes; i++) {
1149 node_t *node = &this->node[i];
1150
1151 if(node->getX(0) > xmax)
1152 xmax = node->getX(0);
1153
1154 if(node->getX(0) < xmin)
1155 xmin = node->getX(0);
1156
1157 if(node->getX(1) > ymax)
1158 ymax = node->getX(1);
1159
1160 if(node->getX(1) < ymin)
1161 ymin = node->getX(1);
1162
1163 if(node->getX(2) > zmax)
1164 zmax = node->getX(2);
1165
1166 if(node->getX(2) < zmin)
1167 zmin = node->getX(2);
1168 }
1169
1170 double xmid = (xmin + xmax)/2.0;
1171 double ymid = (ymin + ymax)/2.0;
1172 double zmid = (zmin + zmax)/2.0;
1173
1174 double xlen = (xmax - xmin)/2.0;
1175 double ylen = (ymax - ymin)/2.0;
1176 double zlen = (zmax - zmin)/2.0;
1177
1178 double s = xlen;
1179
1180 if(ylen > s)
1181 s = ylen;
1182
1183 if(zlen > s)
1184 s = zlen;
1185
1186 s *= 1.1;
1187
1188 bool sx = xmin==xmax;
1189 bool sy = ymin==ymax;
1190 bool sz = zmin==zmax;
1191
1192 this->cdim = 3;
1193
1194 if((sz && sy) || (sz && sx) || (sx && sy))
1195 this->cdim = 1;
1196
1197 else if(sz || sy || sx)
1198 this->cdim = 2;
1199
1200 result[0] = xmin;
1201 result[1] = xmax;
1202 result[2] = ymin;
1203 result[3] = ymax;
1204 result[4] = zmin;
1205 result[5] = zmax;
1206
1207 result[6] = xmid;
1208 result[7] = ymid;
1209 result[8] = zmid;
1210
1211 result[9] = s;
1212
1213 return result;
1214 }
1215
setCdim(int n)1216 void mesh_t::setCdim(int n)
1217 {
1218 this->cdim = n;
1219 }
1220
getCdim() const1221 int mesh_t::getCdim() const
1222 {
1223 return this->cdim;
1224 }
1225
setDim(int n)1226 void mesh_t::setDim(int n)
1227 {
1228 this->dim = n;
1229 }
1230
getDim() const1231 int mesh_t::getDim() const
1232 {
1233 return this->dim;
1234 }
1235
setNodes(int n)1236 void mesh_t::setNodes(int n)
1237 {
1238 this->nodes = n;
1239 }
1240
getNodes() const1241 int mesh_t::getNodes() const
1242 {
1243 return this->nodes;
1244 }
1245
setPoints(int n)1246 void mesh_t::setPoints(int n)
1247 {
1248 this->points = n;
1249 }
1250
getPoints() const1251 int mesh_t::getPoints() const
1252 {
1253 return this->points;
1254 }
1255
setEdges(int n)1256 void mesh_t::setEdges(int n)
1257 {
1258 this->edges = n;
1259 }
1260
getEdges() const1261 int mesh_t::getEdges() const
1262 {
1263 return this->edges;
1264 }
1265
setSurfaces(int n)1266 void mesh_t::setSurfaces(int n)
1267 {
1268 this->surfaces = n;
1269 }
1270
getSurfaces() const1271 int mesh_t::getSurfaces() const
1272 {
1273 return this->surfaces;
1274 }
1275
setElements(int n)1276 void mesh_t::setElements(int n)
1277 {
1278 this->elements = n;
1279 }
1280
getElements() const1281 int mesh_t::getElements() const
1282 {
1283 return this->elements;
1284 }
1285
getNode(int n)1286 node_t* mesh_t::getNode(int n)
1287 {
1288 return &this->node[n];
1289 }
1290
setNodeArray(node_t * n)1291 void mesh_t::setNodeArray(node_t* n)
1292 {
1293 this->node = n;
1294 }
1295
newNodeArray(int n)1296 void mesh_t::newNodeArray(int n)
1297 {
1298 this->node = new node_t[n];
1299 }
1300
deleteNodeArray()1301 void mesh_t::deleteNodeArray()
1302 {
1303 delete [] this->node;
1304 }
1305
getPoint(int n)1306 point_t* mesh_t::getPoint(int n)
1307 {
1308 return &this->point[n];
1309 }
1310
setPointArray(point_t * p)1311 void mesh_t::setPointArray(point_t* p)
1312 {
1313 this->point = p;
1314 }
1315
newPointArray(int n)1316 void mesh_t::newPointArray(int n)
1317 {
1318 this->point = new point_t[n];
1319 }
1320
deletePointArray()1321 void mesh_t::deletePointArray()
1322 {
1323 delete [] this->point;
1324 }
1325
getEdge(int n)1326 edge_t* mesh_t::getEdge(int n)
1327 {
1328 return &this->edge[n];
1329 }
1330
setEdgeArray(edge_t * e)1331 void mesh_t::setEdgeArray(edge_t* e)
1332 {
1333 this->edge = e;
1334 }
1335
newEdgeArray(int n)1336 void mesh_t::newEdgeArray(int n)
1337 {
1338 this->edge = new edge_t[n];
1339 }
1340
deleteEdgeArray()1341 void mesh_t::deleteEdgeArray()
1342 {
1343 delete [] this->edge;
1344 }
1345
getSurface(int n)1346 surface_t* mesh_t::getSurface(int n)
1347 {
1348 return &this->surface[n];
1349 }
1350
setSurfaceArray(surface_t * s)1351 void mesh_t::setSurfaceArray(surface_t* s)
1352 {
1353 this->surface = s;
1354 }
1355
newSurfaceArray(int n)1356 void mesh_t::newSurfaceArray(int n)
1357 {
1358 this->surface = new surface_t[n];
1359 }
1360
deleteSurfaceArray()1361 void mesh_t::deleteSurfaceArray()
1362 {
1363 delete [] this->surface;
1364 }
1365
getElement(int n)1366 element_t* mesh_t::getElement(int n)
1367 {
1368 return &this->element[n];
1369 }
1370
setElementArray(element_t * e)1371 void mesh_t::setElementArray(element_t* e)
1372 {
1373 this->element = e;
1374 }
1375
newElementArray(int n)1376 void mesh_t::newElementArray(int n)
1377 {
1378 this->element = new element_t[n];
1379 }
1380
deleteElementArray()1381 void mesh_t::deleteElementArray()
1382 {
1383 delete [] this->element;
1384 }
1385