1 #include <mystdlib.h>
2
3 #include <meshing.hpp>
4 #include <csg.hpp>
5
6
7 #ifdef SOCKETS
8 #include "../sockets/sockets.hpp"
9 #endif
10
11 #ifndef NOTCL
12 #include <visual.hpp>
13 #endif
14
15
16 #include "nginterface.h"
17 #include "../visualization/soldata.hpp"
18
19
20 #ifdef _MSC_VER
21 // Philippose - 30/01/2009
22 // MSVC Express Edition Support
23 #ifdef MSVC_EXPRESS
24
25 // #include <pthread.h>
26
27 static pthread_t meshingthread;
RunParallel(void * (* fun)(void *),void * in)28 void RunParallel ( void * (*fun)(void *), void * in)
29 {
30 if (netgen::mparam.parthread) // && (ntasks == 1) )
31 {
32 pthread_attr_t attr;
33 pthread_attr_init (&attr);
34 // the following call can be removed if not available:
35 pthread_attr_setstacksize(&attr, 1000000);
36 //pthread_create (&meshingthread, &attr, fun, NULL);
37 pthread_create (&meshingthread, &attr, fun, in);
38 }
39 else
40 fun (in);
41 }
42
43 #else // Using MS VC++ Standard / Enterprise / Professional edition
44
45 // Afx - Threads need different return - value:
46
47 static void* (*sfun)(void *);
fun2(void * val)48 unsigned int fun2 (void * val)
49 {
50 sfun (val);
51 return 0;
52 }
53
RunParallel(void * (* fun)(void *),void * in)54 void RunParallel ( void* (*fun)(void *), void * in)
55 {
56 sfun = fun;
57 if (netgen::mparam.parthread)
58 AfxBeginThread (fun2, in);
59 //AfxBeginThread (fun2, NULL);
60 else
61 fun (in);
62 }
63
64 #endif // #ifdef MSVC_EXPRESS
65
66 #else // For #ifdef _MSC_VER
67
68 // #include <pthread.h>
69
70 static pthread_t meshingthread;
RunParallel(void * (* fun)(void *),void * in)71 void RunParallel ( void * (*fun)(void *), void * in)
72 {
73 bool parthread = netgen::mparam.parthread;
74
75 #ifdef PARALLEL
76 int provided;
77 MPI_Query_thread(&provided);
78 if (provided < 3)
79 if (netgen::ntasks > 1) parthread = false;
80 // cout << "runparallel = " << parthread << endl;
81 #endif
82
83 if (parthread)
84 {
85 pthread_attr_t attr;
86 pthread_attr_init (&attr);
87 // the following call can be removed if not available:
88 pthread_attr_setstacksize(&attr, 1000000);
89 //pthread_create (&meshingthread, &attr, fun, NULL);
90 pthread_create (&meshingthread, &attr, fun, in);
91 }
92 else
93 fun (in);
94 }
95
96 #endif // #ifdef _MSC_VER
97
98
99
100
101
102
103 namespace netgen
104 {
105 #include "writeuser.hpp"
106
107 MeshingParameters mparam;
108
109 // global variable mesh (should not be used in libraries)
110 AutoPtr<Mesh> mesh;
111 NetgenGeometry * ng_geometry = new NetgenGeometry;
112
113 // extern NetgenGeometry * ng_geometry;
114 // extern AutoPtr<Mesh> mesh;
115
116 #ifndef NOTCL
117 extern Tcl_Interp * tcl_interp;
118 #endif
119
120
121 #ifdef OPENGL
122 extern VisualSceneSolution vssolution;
123 #endif
124 extern CSGeometry * ParseCSG (istream & istr);
125
126 #ifdef SOCKETS
127 extern AutoPtr<ClientSocket> clientsocket;
128 //extern Array< AutoPtr < ServerInfo > > servers;
129 extern Array< ServerInfo* > servers;
130 #endif
131
132
133 }
134
135
136 using namespace netgen;
137
138
Ng_LoadGeometry(const char * filename)139 void Ng_LoadGeometry (const char * filename)
140 {
141
142 for (int i = 0; i < geometryregister.Size(); i++)
143 {
144 NetgenGeometry * hgeom = geometryregister[i]->Load (filename);
145 if (hgeom)
146 {
147 delete ng_geometry;
148 ng_geometry = hgeom;
149
150 mesh.Reset();
151 return;
152 }
153 }
154
155 // he: if filename is empty, return
156 // can be used to reset geometry
157 if (strcmp(filename,"")==0)
158 {
159 delete ng_geometry;
160 ng_geometry = new NetgenGeometry();
161 return;
162 }
163
164 cerr << "cannot load geometry '" << filename << "'" << endl;
165 }
166
167
Ng_LoadMeshFromStream(istream & input)168 void Ng_LoadMeshFromStream ( istream & input )
169 {
170 mesh.Reset (new Mesh());
171 mesh -> Load(input);
172
173 for (int i = 0; i < geometryregister.Size(); i++)
174 {
175 NetgenGeometry * hgeom = geometryregister[i]->LoadFromMeshFile (input);
176 if (hgeom)
177 {
178 delete ng_geometry;
179 ng_geometry = hgeom;
180 break;
181 }
182 }
183
184
185 #ifdef LOADOLD
186 if(input.good())
187 {
188 string auxstring;
189 input >> auxstring;
190 if(auxstring == "csgsurfaces")
191 {
192 /*
193 if (geometry)
194 {
195 geometry.Reset (new CSGeometry (""));
196 }
197 if (stlgeometry)
198 {
199 delete stlgeometry;
200 stlgeometry = NULL;
201 }
202 #ifdef OCCGEOMETRY
203 if (occgeometry)
204 {
205 delete occgeometry;
206 occgeometry = NULL;
207 }
208 #endif
209 #ifdef ACIS
210 if (acisgeometry)
211 {
212 delete acisgeometry;
213 acisgeometry = NULL;
214 }
215 #endif
216 geometry2d.Reset (0);
217 */
218 // geometry -> LoadSurfaces(input);
219 CSGeometry * geometry = new CSGeometry ("");
220 geometry -> LoadSurfaces(input);
221
222 delete ng_geometry;
223 ng_geometry = geometry;
224 }
225 }
226 #endif
227 }
228
229
230
231
Ng_LoadMesh(const char * filename)232 void Ng_LoadMesh (const char * filename)
233 {
234 if ( (strlen (filename) > 4) &&
235 strcmp (filename + (strlen (filename)-4), ".vol") != 0 )
236 {
237 mesh.Reset (new Mesh());
238 ReadFile(*mesh,filename);
239
240 //mesh->SetGlobalH (mparam.maxh);
241 //mesh->CalcLocalH();
242 return;
243 }
244
245 ifstream infile(filename);
246 Ng_LoadMeshFromStream(infile);
247 }
248
Ng_LoadMeshFromString(const char * mesh_as_string)249 void Ng_LoadMeshFromString (const char * mesh_as_string)
250 {
251 istringstream instream(mesh_as_string);
252 Ng_LoadMeshFromStream(instream);
253 }
254
255
256
257
Ng_GetDimension()258 int Ng_GetDimension ()
259 {
260 return (mesh) ? mesh->GetDimension() : -1;
261 }
262
Ng_GetNP()263 int Ng_GetNP ()
264 {
265 return (mesh) ? mesh->GetNP() : 0;
266 }
267
Ng_GetNV()268 int Ng_GetNV ()
269 {
270 return (mesh) ? mesh->GetNV() : 0;
271 }
272
Ng_GetNE()273 int Ng_GetNE ()
274 {
275 if(!mesh) return 0;
276 if (mesh->GetDimension() == 3)
277 return mesh->GetNE();
278 else
279 return mesh->GetNSE();
280 }
281
Ng_GetNSE()282 int Ng_GetNSE ()
283 {
284 if(!mesh) return 0;
285 if (mesh->GetDimension() == 3)
286 return mesh->GetNSE();
287 else
288 return mesh->GetNSeg();
289 }
290
Ng_GetPoint(int pi,double * p)291 void Ng_GetPoint (int pi, double * p)
292 {
293 if (pi < 1 || pi > mesh->GetNP())
294 {
295 if (printmessage_importance>0)
296 cout << "Ng_GetPoint: illegal point " << pi << endl;
297 return;
298 }
299
300 const Point3d & hp = mesh->Point (pi);
301 p[0] = hp.X();
302 p[1] = hp.Y();
303 if (mesh->GetDimension() == 3)
304 p[2] = hp.Z();
305 }
306
307
Ng_GetElement(int ei,int * epi,int * np)308 NG_ELEMENT_TYPE Ng_GetElement (int ei, int * epi, int * np)
309 {
310 if (mesh->GetDimension() == 3)
311 {
312 int i;
313 const Element & el = mesh->VolumeElement (ei);
314 for (i = 0; i < el.GetNP(); i++)
315 epi[i] = el.PNum(i+1);
316
317 if (np)
318 *np = el.GetNP();
319
320 if (el.GetType() == PRISM)
321 {
322 // degenerated prism, (should be obsolete)
323 const int map1[] = { 3, 2, 5, 6, 1 };
324 const int map2[] = { 1, 3, 6, 4, 2 };
325 const int map3[] = { 2, 1, 4, 5, 3 };
326
327 const int * map = NULL;
328 int deg1 = 0, deg2 = 0, deg3 = 0;
329 //int deg = 0;
330 if (el.PNum(1) == el.PNum(4)) { map = map1; deg1 = 1; }
331 if (el.PNum(2) == el.PNum(5)) { map = map2; deg2 = 1; }
332 if (el.PNum(3) == el.PNum(6)) { map = map3; deg3 = 1; }
333
334 switch (deg1+deg2+deg3)
335 {
336 {
337 case 1:
338 if (printmessage_importance>0)
339 cout << "degenerated prism found, deg = 1" << endl;
340 for (i = 0; i < 5; i++)
341 epi[i] = el.PNum (map[i]);
342
343 if (np) *np = 5;
344 return NG_PYRAMID;
345 break;
346 }
347 case 2:
348 {
349 if (printmessage_importance>0)
350 cout << "degenerated prism found, deg = 2" << endl;
351 if (!deg1) epi[3] = el.PNum(4);
352 if (!deg2) epi[3] = el.PNum(5);
353 if (!deg3) epi[3] = el.PNum(6);
354
355 if (np) *np = 4;
356 return NG_TET;
357 break;
358 }
359 default:
360 ;
361 }
362
363 }
364
365 return NG_ELEMENT_TYPE (el.GetType());
366 }
367 else
368 {
369 int i;
370 const Element2d & el = mesh->SurfaceElement (ei);
371 for (i = 0; i < el.GetNP(); i++)
372 epi[i] = el.PNum(i+1);
373
374 if (np) *np = el.GetNP();
375 return NG_ELEMENT_TYPE (el.GetType());
376 /*
377 switch (el.GetNP())
378 {
379 case 3: return NG_TRIG;
380 case 4: return NG_QUAD;
381 case 6: return NG_TRIG6;
382 }
383 */
384 }
385
386 // should not occur
387 return NG_TET;
388 }
389
390
Ng_GetElementType(int ei)391 NG_ELEMENT_TYPE Ng_GetElementType (int ei)
392 {
393 if (mesh->GetDimension() == 3)
394 {
395 return NG_ELEMENT_TYPE (mesh->VolumeElement (ei).GetType());
396 }
397 else
398 {
399 const Element2d & el = mesh->SurfaceElement (ei);
400 switch (el.GetNP())
401 {
402 case 3: return NG_TRIG;
403 case 4: return NG_QUAD;
404 case 6: return NG_TRIG6;
405 }
406 }
407
408 // should not occur
409 return NG_TET;
410 }
411
412
413
Ng_GetElementIndex(int ei)414 int Ng_GetElementIndex (int ei)
415 {
416 if (mesh->GetDimension() == 3)
417 return mesh->VolumeElement(ei).GetIndex();
418 else
419 {
420 int ind = mesh->SurfaceElement(ei).GetIndex();
421 ind = mesh->GetFaceDescriptor(ind).BCProperty();
422 return ind;
423 }
424 }
425
Ng_SetElementIndex(const int ei,const int index)426 void Ng_SetElementIndex(const int ei, const int index)
427 {
428 mesh->VolumeElement(ei).SetIndex(index);
429 }
430
Ng_GetElementMaterial(int ei)431 char * Ng_GetElementMaterial (int ei)
432 {
433 static char empty[] = "";
434 if (mesh->GetDimension() == 3)
435 {
436 int ind = mesh->VolumeElement(ei).GetIndex();
437 // cout << "ind = " << ind << endl;
438 const char * mat = mesh->GetMaterial (ind);
439 if (mat)
440 return const_cast<char*> (mat);
441 else
442 return empty;
443 }
444 // add astrid
445 else
446 {
447 int ind = mesh->SurfaceElement(ei).GetIndex();
448 ind = mesh->GetFaceDescriptor(ind).BCProperty();
449 const char * mat = mesh->GetMaterial ( ind );
450 if (mat)
451 return const_cast<char*> (mat);
452 else
453 return empty;
454 }
455 return 0;
456 }
457
Ng_GetDomainMaterial(int dom)458 char * Ng_GetDomainMaterial (int dom)
459 {
460 static char empty[] = "";
461 // astrid
462 if ( 1 ) // mesh->GetDimension() == 3)
463 {
464 const char * mat = mesh->GetMaterial(dom);
465 if (mat)
466 return const_cast<char*> (mat);
467 else
468 return empty;
469 }
470
471 return 0;
472 }
473
Ng_GetUserDataSize(char * id)474 int Ng_GetUserDataSize (char * id)
475 {
476 Array<double> da;
477 mesh->GetUserData (id, da);
478 return da.Size();
479 }
480
Ng_GetUserData(char * id,double * data)481 void Ng_GetUserData (char * id, double * data)
482 {
483 Array<double> da;
484 mesh->GetUserData (id, da);
485 for (int i = 0; i < da.Size(); i++)
486 data[i] = da[i];
487 }
488
489
Ng_GetSurfaceElement(int ei,int * epi,int * np)490 NG_ELEMENT_TYPE Ng_GetSurfaceElement (int ei, int * epi, int * np)
491 {
492 if (mesh->GetDimension() == 3)
493 {
494 const Element2d & el = mesh->SurfaceElement (ei);
495 for (int i = 0; i < el.GetNP(); i++)
496 epi[i] = el[i];
497
498 if (np) *np = el.GetNP();
499
500 return NG_ELEMENT_TYPE (el.GetType());
501 }
502 else
503 {
504 const Segment & seg = mesh->LineSegment (ei);
505
506 if (seg[2] < 0)
507 {
508 epi[0] = seg[0];
509 epi[1] = seg[1];
510
511 if (np) *np = 2;
512 return NG_SEGM;
513 }
514 else
515 {
516 epi[0] = seg[0];
517 epi[1] = seg[1];
518 epi[2] = seg[2];
519
520 if (np) *np = 3;
521 return NG_SEGM3;
522 }
523 }
524
525 return NG_TRIG;
526 }
527
Ng_GetSurfaceElementIndex(int ei)528 int Ng_GetSurfaceElementIndex (int ei)
529 {
530 if (mesh->GetDimension() == 3)
531 return mesh->GetFaceDescriptor(mesh->SurfaceElement(ei).GetIndex()).BCProperty();
532 else
533 return mesh->LineSegment(ei).si;
534 }
535
Ng_GetSurfaceElementSurfaceNumber(int ei)536 int Ng_GetSurfaceElementSurfaceNumber (int ei)
537 {
538 if (mesh->GetDimension() == 3)
539 return mesh->GetFaceDescriptor(mesh->SurfaceElement(ei).GetIndex()).SurfNr();
540 else
541 return mesh->LineSegment(ei).si;
542 }
Ng_GetSurfaceElementFDNumber(int ei)543 int Ng_GetSurfaceElementFDNumber (int ei)
544 {
545 if (mesh->GetDimension() == 3)
546 return mesh->SurfaceElement(ei).GetIndex();
547 else
548 return -1;
549 }
550
551
Ng_GetSurfaceElementBCName(int ei)552 char * Ng_GetSurfaceElementBCName (int ei)
553 {
554 if ( mesh->GetDimension() == 3 )
555 return const_cast<char *>(mesh->GetFaceDescriptor(mesh->SurfaceElement(ei).GetIndex()).GetBCName().c_str());
556 else
557 return const_cast<char *>(mesh->LineSegment(ei).GetBCName().c_str());
558 }
559
560
561 // Inefficient (but maybe safer) version:
562 //void Ng_GetSurfaceElementBCName (int ei, char * name)
563 //{
564 // if ( mesh->GetDimension() == 3 )
565 // strcpy(name,mesh->GetFaceDescriptor(mesh->SurfaceElement(ei).GetIndex()).GetBCName().c_str());
566 // else
567 // strcpy(name,mesh->LineSegment(ei).GetBCName().c_str());
568 //}
569
Ng_GetBCNumBCName(int bcnr)570 char * Ng_GetBCNumBCName (int bcnr)
571 {
572 return const_cast<char *>(mesh->GetBCName(bcnr).c_str());
573 }
574
575
576 // Inefficient (but maybe safer) version:
577 //void Ng_GetBCNumBCName (int bcnr, char * name)
578 //{
579 // strcpy(name,mesh->GetBCName(bcnr).c_str());
580 //}
581
Ng_GetNormalVector(int sei,int locpi,double * nv)582 void Ng_GetNormalVector (int sei, int locpi, double * nv)
583 {
584 nv[0] = 0;
585 nv[1] = 0;
586 nv[2] = 1;
587
588 if (mesh->GetDimension() == 3)
589 {
590 Vec<3> n;
591 Point<3> p;
592 p = mesh->Point (mesh->SurfaceElement(sei).PNum(locpi));
593
594 int surfi = mesh->GetFaceDescriptor(mesh->SurfaceElement(sei).GetIndex()).SurfNr();
595
596 (*testout) << "surfi = " << surfi << endl;
597 #ifdef OCCGEOMETRYxxx
598 OCCGeometry * occgeometry = dynamic_cast<OCCGeometry*> (ng_geometry);
599 if (occgeometry)
600 {
601 PointGeomInfo gi = mesh->SurfaceElement(sei).GeomInfoPi(locpi);
602 occgeometry->GetSurface (surfi).GetNormalVector(p, gi, n);
603 nv[0] = n(0);
604 nv[1] = n(1);
605 nv[2] = n(2);
606 }
607 #endif
608 CSGeometry * geometry = dynamic_cast<CSGeometry*> (ng_geometry);
609 if (geometry)
610 {
611 n = geometry->GetSurface (surfi) -> GetNormalVector(p);
612 nv[0] = n(0);
613 nv[1] = n(1);
614 nv[2] = n(2);
615 }
616 }
617 }
618
619
620
Ng_SetPointSearchStartElement(const int el)621 void Ng_SetPointSearchStartElement(const int el)
622 {
623 mesh->SetPointSearchStartElement(el);
624 }
625
626
Ng_FindElementOfPoint(double * p,double * lami,int build_searchtree,const int * const indices,const int numind)627 int Ng_FindElementOfPoint (double * p, double * lami, int build_searchtree,
628 const int * const indices, const int numind)
629
630 {
631 Array<int> * dummy(NULL);
632 int ind = -1;
633
634 if(indices != NULL)
635 {
636 dummy = new Array<int>(numind);
637 for(int i=0; i<numind; i++) (*dummy)[i] = indices[i];
638 }
639
640 if (mesh->GetDimension() == 3)
641 {
642 Point3d p3d(p[0], p[1], p[2]);
643 ind =
644 mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0);
645 }
646 else
647 {
648 double lam3[3];
649 Point3d p2d(p[0], p[1], 0);
650 ind =
651 mesh->GetElementOfPoint(p2d, lam3, dummy, build_searchtree != 0);
652
653 if (ind > 0)
654 {
655 if(mesh->SurfaceElement(ind).GetType()==QUAD)
656 {
657 lami[0] = lam3[0];
658 lami[1] = lam3[1];
659 }
660 else
661 {
662 lami[0] = 1-lam3[0]-lam3[1];
663 lami[1] = lam3[0];
664 }
665 }
666 }
667
668 delete dummy;
669
670 return ind;
671 }
672
Ng_FindSurfaceElementOfPoint(double * p,double * lami,int build_searchtree,const int * const indices,const int numind)673 int Ng_FindSurfaceElementOfPoint (double * p, double * lami, int build_searchtree,
674 const int * const indices, const int numind)
675
676 {
677 Array<int> * dummy(NULL);
678 int ind = -1;
679
680 if(indices != NULL)
681 {
682 dummy = new Array<int>(numind);
683 for(int i=0; i<numind; i++) (*dummy)[i] = indices[i];
684 }
685
686 if (mesh->GetDimension() == 3)
687 {
688 Point3d p3d(p[0], p[1], p[2]);
689 ind =
690 mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0);
691 }
692 else
693 {
694 //throw NgException("FindSurfaceElementOfPoint for 2D meshes not yet implemented");
695 cerr << "FindSurfaceElementOfPoint for 2D meshes not yet implemented" << endl;
696 }
697
698 delete dummy;
699
700 return ind;
701 }
702
703
Ng_IsElementCurved(int ei)704 int Ng_IsElementCurved (int ei)
705 {
706 if (mesh->GetDimension() == 2)
707 return mesh->GetCurvedElements().IsSurfaceElementCurved (ei-1);
708 else
709 return mesh->GetCurvedElements().IsElementCurved (ei-1);
710 }
711
712
Ng_IsSurfaceElementCurved(int sei)713 int Ng_IsSurfaceElementCurved (int sei)
714 {
715 if (mesh->GetDimension() == 2)
716 return mesh->GetCurvedElements().IsSegmentCurved (sei-1);
717 else
718 return mesh->GetCurvedElements().IsSurfaceElementCurved (sei-1);
719 }
720
721
722
723
Ng_GetElementTransformation(int ei,const double * xi,double * x,double * dxdxi)724 void Ng_GetElementTransformation (int ei, const double * xi,
725 double * x, double * dxdxi)
726 {
727 if (mesh->GetDimension() == 2)
728 {
729 Point<2> xl(xi[0], xi[1]);
730 Point<3> xg;
731 Mat<3,2> dx;
732
733 mesh->GetCurvedElements().CalcSurfaceTransformation (xl, ei-1, xg, dx);
734
735 if (x)
736 {
737 for (int i = 0; i < 2; i++)
738 x[i] = xg(i);
739 }
740
741 if (dxdxi)
742 {
743 for (int i=0; i<2; i++)
744 {
745 dxdxi[2*i] = dx(i,0);
746 dxdxi[2*i+1] = dx(i,1);
747 }
748 }
749 }
750 else
751 {
752 Point<3> xl(xi[0], xi[1], xi[2]);
753 Point<3> xg;
754 Mat<3,3> dx;
755
756 mesh->GetCurvedElements().CalcElementTransformation (xl, ei-1, xg, dx);
757
758 if (x)
759 {
760 for (int i = 0; i < 3; i++)
761 x[i] = xg(i);
762 }
763
764 if (dxdxi)
765 {
766 for (int i=0; i<3; i++)
767 {
768 dxdxi[3*i] = dx(i,0);
769 dxdxi[3*i+1] = dx(i,1);
770 dxdxi[3*i+2] = dx(i,2);
771 }
772 }
773 }
774 }
775
776
777 #ifdef OLD
Ng_GetBufferedElementTransformation(int ei,const double * xi,double * x,double * dxdxi,void * buffer,int buffervalid)778 void Ng_GetBufferedElementTransformation (int ei, const double * xi,
779 double * x, double * dxdxi,
780 void * buffer, int buffervalid)
781 {
782 // buffer = 0;
783 // buffervalid = 0;
784 if (mesh->GetDimension() == 2)
785 {
786 return Ng_GetElementTransformation (ei, xi, x, dxdxi);
787 }
788 else
789 {
790 mesh->GetCurvedElements().CalcElementTransformation (reinterpret_cast<const Point<3> &> (*xi),
791 ei-1,
792 reinterpret_cast<Point<3> &> (*x),
793 reinterpret_cast<Mat<3,3> &> (*dxdxi),
794 buffer, (buffervalid != 0));
795
796 /*
797 Point<3> xl(xi[0], xi[1], xi[2]);
798 Point<3> xg;
799 Mat<3,3> dx;
800 // buffervalid = 0;
801 mesh->GetCurvedElements().CalcElementTransformation (xl, ei-1, xg, dx, buffer, buffervalid);
802
803 // still 1-based arrays
804 if (x)
805 {
806 for (int i = 0; i < 3; i++)
807 x[i] = xg(i);
808 }
809
810 if (dxdxi)
811 {
812 for (int i=0; i<3; i++)
813 {
814 dxdxi[3*i] = dx(i,0);
815 dxdxi[3*i+1] = dx(i,1);
816 dxdxi[3*i+2] = dx(i,2);
817 }
818 }
819 */
820 }
821 }
822 #endif
823
824
825
826
827
828
Ng_GetMultiElementTransformation(int ei,int n,const double * xi,size_t sxi,double * x,size_t sx,double * dxdxi,size_t sdxdxi)829 void Ng_GetMultiElementTransformation (int ei, int n,
830 const double * xi, size_t sxi,
831 double * x, size_t sx,
832 double * dxdxi, size_t sdxdxi)
833 {
834 if (mesh->GetDimension() == 2)
835 mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2> (ei-1, n, xi, sxi, x, sx, dxdxi, sdxdxi);
836 else
837 mesh->GetCurvedElements().CalcMultiPointElementTransformation (ei-1, n, xi, sxi, x, sx, dxdxi, sdxdxi);
838 }
839
840
841
Ng_GetSurfaceElementTransformation(int sei,const double * xi,double * x,double * dxdxi)842 void Ng_GetSurfaceElementTransformation (int sei, const double * xi,
843 double * x, double * dxdxi)
844 {
845 if (mesh->GetDimension() == 2)
846 {
847 Point<3> xg;
848 Vec<3> dx;
849
850 mesh->GetCurvedElements().CalcSegmentTransformation (xi[0], sei-1, xg, dx);
851
852 if (x)
853 for (int i = 0; i < 2; i++)
854 x[i] = xg(i);
855
856 if (dxdxi)
857 for (int i=0; i<2; i++)
858 dxdxi[i] = dx(i);
859
860 }
861 else
862 {
863 Point<2> xl(xi[0], xi[1]);
864 Point<3> xg;
865 Mat<3,2> dx;
866
867 mesh->GetCurvedElements().CalcSurfaceTransformation (xl, sei-1, xg, dx);
868
869 for (int i=0; i<3; i++)
870 {
871 if (x)
872 x[i] = xg(i);
873 if (dxdxi)
874 {
875 dxdxi[2*i] = dx(i,0);
876 dxdxi[2*i+1] = dx(i,1);
877 }
878 }
879 }
880 }
881
882
883
884
885
Ng_GetSegmentIndex(int ei)886 int Ng_GetSegmentIndex (int ei)
887 {
888 const Segment & seg = mesh->LineSegment (ei);
889 return seg.edgenr;
890 }
891
892
Ng_GetSegment(int ei,int * epi,int * np)893 NG_ELEMENT_TYPE Ng_GetSegment (int ei, int * epi, int * np)
894 {
895 const Segment & seg = mesh->LineSegment (ei);
896
897 epi[0] = seg[0];
898 epi[1] = seg[1];
899
900 if (seg[2] < 0)
901 {
902 if (np) *np = 2;
903 return NG_SEGM;
904 }
905 else
906 {
907 epi[2] = seg[2];
908 if (np) *np = 3;
909 return NG_SEGM3;
910 }
911 }
912
913
914
915
916
917
Ng_GetSurfaceElementNeighbouringDomains(const int selnr,int & in,int & out)918 void Ng_GetSurfaceElementNeighbouringDomains(const int selnr, int & in, int & out)
919 {
920 if ( mesh->GetDimension() == 3 )
921 {
922 in = mesh->GetFaceDescriptor(mesh->SurfaceElement(selnr).GetIndex()).DomainIn();
923 out = mesh->GetFaceDescriptor(mesh->SurfaceElement(selnr).GetIndex()).DomainOut();
924 }
925 else
926 {
927 in = mesh -> LineSegment(selnr) . domin;
928 out = mesh -> LineSegment(selnr) . domout;
929 }
930 }
931
932
933 #ifdef PARALLEL
934 // Is Element ei an element of this processor ??
Ng_IsGhostEl(int ei)935 bool Ng_IsGhostEl (int ei)
936 {
937 return false;
938 /*
939 if ( mesh->GetDimension() == 3 )
940 return mesh->VolumeElement(ei).IsGhost();
941 else
942 return false;
943 */
944 }
945
Ng_SetGhostEl(const int ei,const bool aisghost)946 void Ng_SetGhostEl(const int ei, const bool aisghost )
947 {
948 ;
949 /*
950 if ( mesh -> GetDimension () == 3 )
951 mesh -> VolumeElement(ei).SetGhost (aisghost);
952 */
953 }
954
Ng_IsGhostSEl(int ei)955 bool Ng_IsGhostSEl (int ei)
956 {
957 return false;
958 /*
959 if ( mesh -> GetDimension () == 3 )
960 return mesh->SurfaceElement(ei).IsGhost();
961 else
962 return false;
963 */
964 }
965
Ng_SetGhostSEl(const int ei,const bool aisghost)966 void Ng_SetGhostSEl(const int ei, const bool aisghost )
967 {
968 ;
969 /*
970 if ( mesh -> GetDimension () == 3 )
971 mesh -> SurfaceElement(ei).SetGhost (aisghost);
972 */
973 }
974
975
Ng_IsGhostVert(int pnum)976 bool Ng_IsGhostVert ( int pnum )
977 {
978 return false;
979 // return mesh -> Point ( pnum ).IsGhost() ;
980 }
Ng_IsGhostEdge(int ednum)981 bool Ng_IsGhostEdge ( int ednum )
982 {
983 return false;
984 // return mesh -> GetParallelTopology() . IsGhostEdge ( ednum );
985 }
986
Ng_IsGhostFace(int fanum)987 bool Ng_IsGhostFace ( int fanum )
988 {
989 return false;
990 // return mesh -> GetParallelTopology() . IsGhostFace ( fanum );
991 }
992
993 // void Ng_SetGhostVert ( const int pnum, const bool aisghost );
994 // void Ng_SetGhostEdge ( const int ednum, const bool aisghost );
995 // void Ng_SetGhostFace ( const int fanum, const bool aisghost );
996
997
Ng_IsExchangeEl(int elnum)998 bool Ng_IsExchangeEl ( int elnum )
999 { return mesh -> GetParallelTopology() . IsExchangeElement ( elnum ); }
1000
Ng_IsExchangeSEl(int selnum)1001 bool Ng_IsExchangeSEl ( int selnum )
1002 { return mesh -> GetParallelTopology() . IsExchangeSEl ( selnum ); }
1003
Ng_UpdateOverlap()1004 void Ng_UpdateOverlap()
1005 {
1006 ; // mesh->UpdateOverlap();
1007 }
1008
Ng_Overlap()1009 int Ng_Overlap ()
1010 {
1011 return 0;
1012 // return mesh->GetParallelTopology() . Overlap();
1013 }
1014
1015
1016
NgPar_GetLoc2Glob_VolEl(int locnum)1017 int NgPar_GetLoc2Glob_VolEl ( int locnum )
1018 {
1019 return mesh -> GetParallelTopology().GetLoc2Glob_VolEl ( locnum+1) -1;
1020 }
1021
1022 // gibt anzahl an distant pnums zurueck
1023 // * pnums entspricht ARRAY<int[2] >
NgPar_GetDistantNodeNums(int nodetype,int locnum,int * distnums)1024 int NgPar_GetDistantNodeNums ( int nodetype, int locnum, int * distnums )
1025 {
1026 int size;
1027 switch ( nodetype )
1028 {
1029 case 0:
1030 size = mesh->GetParallelTopology().GetDistantPNums( locnum+1, distnums );
1031 break;
1032 case 1:
1033 size = mesh->GetParallelTopology().GetDistantEdgeNums( locnum+1, distnums );
1034 break;
1035 case 2:
1036 size = mesh->GetParallelTopology().GetDistantFaceNums( locnum+1, distnums );
1037 break;
1038 case 3:
1039 size = mesh->GetParallelTopology().GetDistantElNums( locnum+1, distnums );
1040 break;
1041 default:
1042 cerr << "NgPar_GetDistantNodeNums() Unknown nodetype " << nodetype << endl;
1043 size = -1;
1044 }
1045 // 0 - based
1046 for ( int i = 0; i < size; i++ )
1047 distnums[2*i+1]--;
1048
1049 return size;
1050 }
1051
NgPar_GetNDistantNodeNums(int nodetype,int locnum)1052 int NgPar_GetNDistantNodeNums ( int nodetype, int locnum )
1053 {
1054 switch ( nodetype )
1055 {
1056 case 0:
1057 return mesh->GetParallelTopology().GetNDistantPNums( locnum+1 );
1058 case 1:
1059 return mesh->GetParallelTopology().GetNDistantEdgeNums( locnum+1 );
1060 case 2:
1061 return mesh->GetParallelTopology().GetNDistantFaceNums( locnum+1 );
1062 case 3:
1063 return mesh->GetParallelTopology().GetNDistantElNums( locnum+1 );
1064 }
1065 return -1;
1066 }
1067
NgPar_GetDistantPNum(int proc,int locpnum)1068 int NgPar_GetDistantPNum ( int proc, int locpnum )
1069 {
1070 return mesh->GetParallelTopology().GetDistantPNum( proc, locpnum+1) - 1;
1071 }
1072
NgPar_GetDistantEdgeNum(int proc,int locpnum)1073 int NgPar_GetDistantEdgeNum ( int proc, int locpnum )
1074 {
1075 return mesh->GetParallelTopology().GetDistantEdgeNum( proc, locpnum+1) - 1;
1076 }
1077
NgPar_GetDistantFaceNum(int proc,int locpnum)1078 int NgPar_GetDistantFaceNum ( int proc, int locpnum )
1079 {
1080 return mesh->GetParallelTopology().GetDistantFaceNum (proc, locpnum+1 ) - 1;
1081 }
1082
NgPar_GetDistantElNum(int proc,int locelnum)1083 int NgPar_GetDistantElNum ( int proc, int locelnum )
1084 {
1085 return mesh->GetParallelTopology().GetDistantElNum (proc, locelnum+1 ) - 1;
1086 }
1087
NgPar_IsExchangeFace(int fnr)1088 bool NgPar_IsExchangeFace ( int fnr )
1089 {
1090 return (mesh->GetParallelTopology().GetNDistantFaceNums( fnr+1 ) > 0);
1091 // return mesh->GetParallelTopology().IsExchangeFace ( fnr+1 );
1092 }
1093
NgPar_IsExchangeVert(int vnum)1094 bool NgPar_IsExchangeVert ( int vnum )
1095 {
1096 return (mesh->GetParallelTopology().GetNDistantPNums( vnum+1 ) > 0);
1097 // return mesh->GetParallelTopology().IsExchangeVert ( vnum+1 );
1098 }
1099
NgPar_IsExchangeEdge(int ednum)1100 bool NgPar_IsExchangeEdge ( int ednum )
1101 {
1102 return (mesh->GetParallelTopology().GetNDistantEdgeNums( ednum+1 ) > 0);
1103 // return mesh->GetParallelTopology().IsExchangeEdge ( ednum+1 );
1104 }
1105
NgPar_IsExchangeElement(int elnum)1106 bool NgPar_IsExchangeElement ( int elnum )
1107 {
1108 return (mesh->GetParallelTopology().GetNDistantElNums( elnum+1 ) > 0);
1109 // return mesh->GetParallelTopology().IsExchangeElement ( elnum+1 );
1110 }
1111
1112
NgPar_PrintParallelMeshTopology()1113 void NgPar_PrintParallelMeshTopology ()
1114 {
1115 mesh -> GetParallelTopology().Print ();
1116 }
1117
1118
1119 #endif
1120
Ng_SetRefinementFlag(int ei,int flag)1121 void Ng_SetRefinementFlag (int ei, int flag)
1122 {
1123 if (mesh->GetDimension() == 3)
1124 {
1125 mesh->VolumeElement(ei).SetRefinementFlag (flag != 0);
1126 mesh->VolumeElement(ei).SetStrongRefinementFlag (flag >= 10);
1127 }
1128 else
1129 {
1130 mesh->SurfaceElement(ei).SetRefinementFlag (flag != 0);
1131 mesh->SurfaceElement(ei).SetStrongRefinementFlag (flag >= 10);
1132 }
1133 }
1134
Ng_SetSurfaceRefinementFlag(int ei,int flag)1135 void Ng_SetSurfaceRefinementFlag (int ei, int flag)
1136 {
1137 if (mesh->GetDimension() == 3)
1138 {
1139 mesh->SurfaceElement(ei).SetRefinementFlag (flag != 0);
1140 mesh->SurfaceElement(ei).SetStrongRefinementFlag (flag >= 10);
1141 }
1142 }
1143
1144
Ng_Refine(NG_REFINEMENT_TYPE reftype)1145 void Ng_Refine (NG_REFINEMENT_TYPE reftype)
1146 {
1147 NgLock meshlock (mesh->MajorMutex(), 1);
1148
1149 BisectionOptions biopt;
1150 biopt.usemarkedelements = 1;
1151 biopt.refine_p = 0;
1152 biopt.refine_hp = 0;
1153 if (reftype == NG_REFINE_P)
1154 biopt.refine_p = 1;
1155 if (reftype == NG_REFINE_HP)
1156 biopt.refine_hp = 1;
1157
1158 const Refinement & ref = ng_geometry->GetRefinement();
1159
1160 // Refinement * ref;
1161 MeshOptimize2d * opt = NULL;
1162
1163 /*
1164 if (geometry2d)
1165 ref = new Refinement2d(*geometry2d);
1166 else if (stlgeometry)
1167 ref = new RefinementSTLGeometry(*stlgeometry);
1168 #ifdef OCCGEOMETRY
1169 else if (occgeometry)
1170 ref = new OCCRefinementSurfaces (*occgeometry);
1171 #endif
1172 #ifdef ACIS
1173 else if (acisgeometry)
1174 {
1175 ref = new ACISRefinementSurfaces (*acisgeometry);
1176 opt = new ACISMeshOptimize2dSurfaces(*acisgeometry);
1177 ref->Set2dOptimizer(opt);
1178 }
1179 #endif
1180 else if (geometry && mesh->GetDimension() == 3)
1181 {
1182 ref = new RefinementSurfaces(*geometry);
1183 opt = new MeshOptimize2dSurfaces(*geometry);
1184 ref->Set2dOptimizer(opt);
1185 }
1186 else
1187 {
1188 ref = new Refinement();
1189 }
1190 */
1191
1192 ref.Bisect (*mesh, biopt);
1193
1194 mesh -> UpdateTopology();
1195 mesh -> GetCurvedElements().SetIsHighOrder (false);
1196
1197 // mesh -> GetCurvedElements().BuildCurvedElements (ref, mparam.elementorder);
1198 // delete ref;
1199 delete opt;
1200 }
1201
Ng_SecondOrder()1202 void Ng_SecondOrder ()
1203 {
1204 const_cast<Refinement&> (ng_geometry->GetRefinement()).MakeSecondOrder(*mesh);
1205 /*
1206 if (stlgeometry)
1207 {
1208 RefinementSTLGeometry ref (*stlgeometry);
1209 ref.MakeSecondOrder (*mesh);
1210 }
1211
1212 else if (geometry2d)
1213 {
1214 Refinement2d ref (*geometry2d);
1215 ref.MakeSecondOrder (*mesh);
1216 }
1217
1218 else if (geometry && mesh->GetDimension() == 3)
1219
1220 {
1221 RefinementSurfaces ref (*geometry);
1222 ref.MakeSecondOrder (*mesh);
1223 }
1224 else
1225 {
1226 if (printmessage_importance>0)
1227 cout << "no geom" << endl;
1228 Refinement ref;
1229 ref.MakeSecondOrder (*mesh);
1230 }
1231 */
1232 mesh -> UpdateTopology();
1233 }
1234
1235 /*
1236 void Ng_HPRefinement (int levels)
1237 {
1238 Refinement * ref;
1239
1240 if (stlgeometry)
1241 ref = new RefinementSTLGeometry (*stlgeometry);
1242 else if (geometry2d)
1243 ref = new Refinement2d (*geometry2d);
1244 else
1245 ref = new RefinementSurfaces (*geometry);
1246
1247
1248 HPRefinement (*mesh, ref, levels);
1249 }
1250
1251 void Ng_HPRefinement (int levels, double parameter)
1252 {
1253 Refinement * ref;
1254
1255 if (stlgeometry)
1256 ref = new RefinementSTLGeometry (*stlgeometry);
1257 else if (geometry2d)
1258 ref = new Refinement2d (*geometry2d);
1259 else
1260 ref = new RefinementSurfaces (*geometry);
1261
1262
1263 HPRefinement (*mesh, ref, levels, parameter);
1264 }
1265 */
1266
Ng_HPRefinement(int levels,double parameter,bool setorders,bool ref_level)1267 void Ng_HPRefinement (int levels, double parameter, bool setorders,
1268 bool ref_level)
1269 {
1270 NgLock meshlock (mesh->MajorMutex(), true);
1271 Refinement & ref = const_cast<Refinement&> (ng_geometry -> GetRefinement());
1272 HPRefinement (*mesh, &ref, levels);
1273 /*
1274 Refinement * ref;
1275
1276 if (stlgeometry)
1277 ref = new RefinementSTLGeometry (*stlgeometry);
1278 else if (geometry2d)
1279 ref = new Refinement2d (*geometry2d);
1280 else
1281 ref = new RefinementSurfaces (*geometry);
1282
1283 HPRefinement (*mesh, ref, levels, parameter, setorders, ref_level);
1284 */
1285 }
1286
1287
Ng_HighOrder(int order,bool rational)1288 void Ng_HighOrder (int order, bool rational)
1289 {
1290 NgLock meshlock (mesh->MajorMutex(), true);
1291 /*
1292 Refinement * ref;
1293
1294 if (stlgeometry)
1295 ref = new RefinementSTLGeometry (*stlgeometry);
1296 #ifdef OCCGEOMETRY
1297 else if (occgeometry)
1298 ref = new OCCRefinementSurfaces (*occgeometry);
1299 #endif
1300 #ifdef ACIS
1301 else if (acisgeometry)
1302 {
1303 ref = new ACISRefinementSurfaces (*acisgeometry);
1304 }
1305 #endif
1306 else if (geometry2d)
1307 ref = new Refinement2d (*geometry2d);
1308 else
1309 {
1310 ref = new RefinementSurfaces (*geometry);
1311 }
1312 */
1313 // cout << "parameter 1: " << argv[1] << " (conversion to int = " << atoi(argv[1]) << ")" << endl;
1314
1315
1316 mesh -> GetCurvedElements().BuildCurvedElements (&const_cast<Refinement&> (ng_geometry -> GetRefinement()),
1317 order, rational);
1318 mesh -> SetNextMajorTimeStamp();
1319
1320 /*
1321 if(mesh)
1322 mesh -> GetCurvedElements().BuildCurvedElements (ref, order, rational);
1323 */
1324
1325 // delete ref;
1326 }
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
Ng_ME_GetNVertices(NG_ELEMENT_TYPE et)1339 int Ng_ME_GetNVertices (NG_ELEMENT_TYPE et)
1340 {
1341 switch (et)
1342 {
1343 case NG_SEGM:
1344 case NG_SEGM3:
1345 return 2;
1346
1347 case NG_TRIG:
1348 case NG_TRIG6:
1349 return 3;
1350
1351 case NG_QUAD:
1352 return 4;
1353
1354 case NG_TET:
1355 case NG_TET10:
1356 return 4;
1357
1358 case NG_PYRAMID:
1359 return 5;
1360
1361 case NG_PRISM:
1362 case NG_PRISM12:
1363 return 6;
1364
1365 case NG_HEX:
1366 return 8;
1367
1368 default:
1369 cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
1370 }
1371 return 0;
1372 }
1373
Ng_ME_GetNEdges(NG_ELEMENT_TYPE et)1374 int Ng_ME_GetNEdges (NG_ELEMENT_TYPE et)
1375 {
1376 switch (et)
1377 {
1378 case NG_SEGM:
1379 case NG_SEGM3:
1380 return 1;
1381
1382 case NG_TRIG:
1383 case NG_TRIG6:
1384 return 3;
1385
1386 case NG_QUAD:
1387 return 4;
1388
1389 case NG_TET:
1390 case NG_TET10:
1391 return 6;
1392
1393 case NG_PYRAMID:
1394 return 8;
1395
1396 case NG_PRISM:
1397 case NG_PRISM12:
1398 return 9;
1399
1400 case NG_HEX:
1401 return 12;
1402
1403 default:
1404 cerr << "Ng_ME_GetNEdges, illegal element type " << et << endl;
1405 }
1406 return 0;
1407 }
1408
1409
Ng_ME_GetNFaces(NG_ELEMENT_TYPE et)1410 int Ng_ME_GetNFaces (NG_ELEMENT_TYPE et)
1411 {
1412 switch (et)
1413 {
1414 case NG_SEGM:
1415 case NG_SEGM3:
1416 return 0;
1417
1418 case NG_TRIG:
1419 case NG_TRIG6:
1420 return 1;
1421
1422 case NG_QUAD:
1423 case NG_QUAD6:
1424 return 1;
1425
1426 case NG_TET:
1427 case NG_TET10:
1428 return 4;
1429
1430 case NG_PYRAMID:
1431 return 5;
1432
1433 case NG_PRISM:
1434 case NG_PRISM12:
1435 return 5;
1436
1437 case NG_HEX:
1438 return 6;
1439
1440 default:
1441 cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl;
1442 }
1443 return 0;
1444 }
1445
1446
Ng_ME_GetVertices(NG_ELEMENT_TYPE et)1447 const NG_POINT * Ng_ME_GetVertices (NG_ELEMENT_TYPE et)
1448 {
1449 static double segm_points [][3] =
1450 { { 1, 0, 0 },
1451 { 0, 0, 0 } };
1452
1453 static double trig_points [][3] =
1454 { { 1, 0, 0 },
1455 { 0, 1, 0 },
1456 { 0, 0, 0 } };
1457
1458 static double quad_points [][3] =
1459 { { 0, 0, 0 },
1460 { 1, 0, 0 },
1461 { 1, 1, 0 },
1462 { 0, 1, 0 } };
1463
1464 static double tet_points [][3] =
1465 { { 1, 0, 0 },
1466 { 0, 1, 0 },
1467 { 0, 0, 1 },
1468 { 0, 0, 0 } };
1469
1470 static double pyramid_points [][3] =
1471 {
1472 { 0, 0, 0 },
1473 { 1, 0, 0 },
1474 { 1, 1, 0 },
1475 { 0, 1, 0 },
1476 { 0, 0, 1-1e-7 },
1477 };
1478
1479 static double prism_points[][3] =
1480 {
1481 { 1, 0, 0 },
1482 { 0, 1, 0 },
1483 { 0, 0, 0 },
1484 { 1, 0, 1 },
1485 { 0, 1, 1 },
1486 { 0, 0, 1 }
1487 };
1488
1489 switch (et)
1490 {
1491 case NG_SEGM:
1492 case NG_SEGM3:
1493 return segm_points;
1494
1495 case NG_TRIG:
1496 case NG_TRIG6:
1497 return trig_points;
1498
1499 case NG_QUAD:
1500 case NG_QUAD6:
1501 return quad_points;
1502
1503 case NG_TET:
1504 case NG_TET10:
1505 return tet_points;
1506
1507 case NG_PYRAMID:
1508 return pyramid_points;
1509
1510 case NG_PRISM:
1511 case NG_PRISM12:
1512 return prism_points;
1513
1514 case NG_HEX:
1515 default:
1516 cerr << "Ng_ME_GetVertices, illegal element type " << et << endl;
1517 }
1518 return 0;
1519 }
1520
1521
1522
Ng_ME_GetEdges(NG_ELEMENT_TYPE et)1523 const NG_EDGE * Ng_ME_GetEdges (NG_ELEMENT_TYPE et)
1524 {
1525 static int segm_edges[1][2] =
1526 { { 1, 2 }};
1527
1528 static int trig_edges[3][2] =
1529 { { 3, 1 },
1530 { 3, 2 },
1531 { 1, 2 }};
1532
1533 static int quad_edges[4][2] =
1534 { { 1, 2 },
1535 { 4, 3 },
1536 { 1, 4 },
1537 { 2, 3 }};
1538
1539
1540 static int tet_edges[6][2] =
1541 { { 4, 1 },
1542 { 4, 2 },
1543 { 4, 3 },
1544 { 1, 2 },
1545 { 1, 3 },
1546 { 2, 3 }};
1547
1548 static int prism_edges[9][2] =
1549 { { 3, 1 },
1550 { 1, 2 },
1551 { 3, 2 },
1552 { 6, 4 },
1553 { 4, 5 },
1554 { 6, 5 },
1555 { 3, 6 },
1556 { 1, 4 },
1557 { 2, 5 }};
1558
1559 static int pyramid_edges[8][2] =
1560 { { 1, 2 },
1561 { 2, 3 },
1562 { 1, 4 },
1563 { 4, 3 },
1564 { 1, 5 },
1565 { 2, 5 },
1566 { 3, 5 },
1567 { 4, 5 }};
1568
1569
1570
1571 switch (et)
1572 {
1573 case NG_SEGM:
1574 case NG_SEGM3:
1575 return segm_edges;
1576
1577 case NG_TRIG:
1578 case NG_TRIG6:
1579 return trig_edges;
1580
1581 case NG_QUAD:
1582 case NG_QUAD6:
1583 return quad_edges;
1584
1585 case NG_TET:
1586 case NG_TET10:
1587 return tet_edges;
1588
1589 case NG_PYRAMID:
1590 return pyramid_edges;
1591
1592 case NG_PRISM:
1593 case NG_PRISM12:
1594 return prism_edges;
1595
1596 case NG_HEX:
1597 default:
1598 cerr << "Ng_ME_GetEdges, illegal element type " << et << endl;
1599 }
1600 return 0;
1601 }
1602
1603
Ng_ME_GetFaces(NG_ELEMENT_TYPE et)1604 const NG_FACE * Ng_ME_GetFaces (NG_ELEMENT_TYPE et)
1605 {
1606 static int tet_faces[4][4] =
1607 { { 4, 2, 3, 0 },
1608 { 4, 1, 3, 0 },
1609 { 4, 1, 2, 0 },
1610 { 1, 2, 3, 0 } };
1611
1612 static int prism_faces[5][4] =
1613 {
1614 { 1, 2, 3, 0 },
1615 { 4, 5, 6, 0 },
1616 { 3, 1, 4, 6 },
1617 { 1, 2, 5, 4 },
1618 { 2, 3, 6, 5 }
1619 };
1620
1621 static int pyramid_faces[5][4] =
1622 {
1623 { 1, 2, 5, 0 },
1624 { 2, 3, 5, 0 },
1625 { 3, 4, 5, 0 },
1626 { 4, 1, 5, 0 },
1627 { 1, 2, 3, 4 }
1628 };
1629
1630 static int trig_faces[1][4] =
1631 {
1632 { 1, 2, 3, 0 },
1633 };
1634
1635 switch (et)
1636 {
1637 case NG_TET:
1638 case NG_TET10:
1639 return tet_faces;
1640
1641 case NG_PRISM:
1642 case NG_PRISM12:
1643 return prism_faces;
1644
1645 case NG_PYRAMID:
1646 return pyramid_faces;
1647
1648
1649 case NG_SEGM:
1650 case NG_SEGM3:
1651
1652 case NG_TRIG:
1653 case NG_TRIG6:
1654 return trig_faces;
1655 case NG_QUAD:
1656
1657
1658 case NG_HEX:
1659
1660 default:
1661 cerr << "Ng_ME_GetFaces, illegal element type " << et << endl;
1662 }
1663 return 0;
1664 }
1665
1666
1667
Ng_UpdateTopology()1668 void Ng_UpdateTopology()
1669 {
1670 if (mesh)
1671 mesh -> UpdateTopology();
1672 }
1673
Ng_SelectMesh(Ng_Mesh newmesh)1674 Ng_Mesh Ng_SelectMesh (Ng_Mesh newmesh)
1675 {
1676 Mesh * hmesh = mesh.Ptr();
1677 mesh.Ptr() = (Mesh*)newmesh;
1678 return hmesh;
1679 }
1680
1681
1682
Ng_GetNEdges()1683 int Ng_GetNEdges()
1684 {
1685 return mesh->GetTopology().GetNEdges();
1686 }
Ng_GetNFaces()1687 int Ng_GetNFaces()
1688 {
1689 return mesh->GetTopology().GetNFaces();
1690 }
1691
1692
1693
Ng_GetElement_Edges(int elnr,int * edges,int * orient)1694 int Ng_GetElement_Edges (int elnr, int * edges, int * orient)
1695 {
1696 const MeshTopology & topology = mesh->GetTopology();
1697 if (mesh->GetDimension() == 3)
1698 return topology.GetElementEdges (elnr, edges, orient);
1699 else
1700 return topology.GetSurfaceElementEdges (elnr, edges, orient);
1701 }
1702
Ng_GetElement_Faces(int elnr,int * faces,int * orient)1703 int Ng_GetElement_Faces (int elnr, int * faces, int * orient)
1704 {
1705 const MeshTopology & topology = mesh->GetTopology();
1706 if (mesh->GetDimension() == 3)
1707 return topology.GetElementFaces (elnr, faces, orient);
1708 else
1709 {
1710 faces[0] = elnr;
1711 if (orient) orient[0] = 0;
1712 return 1;
1713 }
1714 }
1715
Ng_GetSurfaceElement_Edges(int elnr,int * edges,int * orient)1716 int Ng_GetSurfaceElement_Edges (int elnr, int * edges, int * orient)
1717 {
1718 const MeshTopology & topology = mesh->GetTopology();
1719 if (mesh->GetDimension() == 3)
1720 return topology.GetSurfaceElementEdges (elnr, edges, orient);
1721 else
1722 {
1723 if (orient)
1724 topology.GetSegmentEdge(elnr, edges[0], orient[0]);
1725 else
1726 edges[0] = topology.GetSegmentEdge(elnr);
1727 }
1728 return 1;
1729 /*
1730 int i, ned;
1731 const MeshTopology & topology = mesh->GetTopology();
1732 Array<int> ia;
1733 topology.GetSurfaceElementEdges (elnr, ia);
1734 ned = ia.Size();
1735 for (i = 1; i <= ned; i++)
1736 edges[i-1] = ia.Get(i);
1737
1738 if (orient)
1739 {
1740 topology.GetSurfaceElementEdgeOrientations (elnr, ia);
1741 for (i = 1; i <= ned; i++)
1742 orient[i-1] = ia.Get(i);
1743 }
1744 return ned;
1745 */
1746 }
1747
Ng_GetSurfaceElement_Face(int selnr,int * orient)1748 int Ng_GetSurfaceElement_Face (int selnr, int * orient)
1749 {
1750 if (mesh->GetDimension() == 3)
1751 {
1752 const MeshTopology & topology = mesh->GetTopology();
1753 if (orient)
1754 *orient = topology.GetSurfaceElementFaceOrientation (selnr);
1755 return topology.GetSurfaceElementFace (selnr);
1756 }
1757 return -1;
1758 }
1759
Ng_GetFace_Vertices(int fnr,int * vert)1760 int Ng_GetFace_Vertices (int fnr, int * vert)
1761 {
1762 const MeshTopology & topology = mesh->GetTopology();
1763 ArrayMem<int,4> ia;
1764 topology.GetFaceVertices (fnr, ia);
1765 for (int i = 0; i < ia.Size(); i++)
1766 vert[i] = ia[i];
1767 // cout << "face verts = " << ia << endl;
1768 return ia.Size();
1769 }
1770
1771
Ng_GetFace_Edges(int fnr,int * edge)1772 int Ng_GetFace_Edges (int fnr, int * edge)
1773 {
1774 const MeshTopology & topology = mesh->GetTopology();
1775 ArrayMem<int,4> ia;
1776 topology.GetFaceEdges (fnr, ia);
1777 for (int i = 0; i < ia.Size(); i++)
1778 edge[i] = ia[i];
1779 return ia.Size();
1780 }
1781
Ng_GetEdge_Vertices(int ednr,int * vert)1782 void Ng_GetEdge_Vertices (int ednr, int * vert)
1783 {
1784 const MeshTopology & topology = mesh->GetTopology();
1785 topology.GetEdgeVertices (ednr, vert[0], vert[1]);
1786 }
1787
1788
Ng_GetNVertexElements(int vnr)1789 int Ng_GetNVertexElements (int vnr)
1790 {
1791 if (mesh->GetDimension() == 3)
1792 return mesh->GetTopology().GetVertexElements(vnr).Size();
1793 else
1794 return mesh->GetTopology().GetVertexSurfaceElements(vnr).Size();
1795 }
1796
Ng_GetVertexElements(int vnr,int * els)1797 void Ng_GetVertexElements (int vnr, int * els)
1798 {
1799 FlatArray<int> ia(0,0);
1800 if (mesh->GetDimension() == 3)
1801 ia = mesh->GetTopology().GetVertexElements(vnr);
1802 else
1803 ia = mesh->GetTopology().GetVertexSurfaceElements(vnr);
1804 for (int i = 0; i < ia.Size(); i++)
1805 els[i] = ia[i];
1806 }
1807
1808
Ng_GetElementOrder(int enr)1809 int Ng_GetElementOrder (int enr)
1810 {
1811 if (mesh->GetDimension() == 3)
1812 return mesh->VolumeElement(enr).GetOrder();
1813 else
1814 return mesh->SurfaceElement(enr).GetOrder();
1815 }
1816
Ng_GetElementOrders(int enr,int * ox,int * oy,int * oz)1817 void Ng_GetElementOrders (int enr, int * ox, int * oy, int * oz)
1818 {
1819 if (mesh->GetDimension() == 3)
1820 mesh->VolumeElement(enr).GetOrder(*ox, *oy, *oz);
1821 else
1822 mesh->SurfaceElement(enr).GetOrder(*ox, *oy, *oz);
1823 }
1824
Ng_SetElementOrder(int enr,int order)1825 void Ng_SetElementOrder (int enr, int order)
1826 {
1827 if (mesh->GetDimension() == 3)
1828 return mesh->VolumeElement(enr).SetOrder(order);
1829 else
1830 return mesh->SurfaceElement(enr).SetOrder(order);
1831 }
1832
Ng_SetElementOrders(int enr,int ox,int oy,int oz)1833 void Ng_SetElementOrders (int enr, int ox, int oy, int oz)
1834 {
1835 if (mesh->GetDimension() == 3)
1836 mesh->VolumeElement(enr).SetOrder(ox, oy, oz);
1837 else
1838 mesh->SurfaceElement(enr).SetOrder(ox, oy);
1839 }
1840
1841
Ng_GetSurfaceElementOrder(int enr)1842 int Ng_GetSurfaceElementOrder (int enr)
1843 {
1844 return mesh->SurfaceElement(enr).GetOrder();
1845 }
1846
1847 //HERBERT: falsche Anzahl von Argumenten
1848 //void Ng_GetSurfaceElementOrders (int enr, int * ox, int * oy, int * oz)
Ng_GetSurfaceElementOrders(int enr,int * ox,int * oy)1849 void Ng_GetSurfaceElementOrders (int enr, int * ox, int * oy)
1850 {
1851 int d;
1852 mesh->SurfaceElement(enr).GetOrder(*ox, *oy, d);
1853 }
1854
Ng_SetSurfaceElementOrder(int enr,int order)1855 void Ng_SetSurfaceElementOrder (int enr, int order)
1856 {
1857 return mesh->SurfaceElement(enr).SetOrder(order);
1858 }
1859
Ng_SetSurfaceElementOrders(int enr,int ox,int oy)1860 void Ng_SetSurfaceElementOrders (int enr, int ox, int oy)
1861 {
1862 mesh->SurfaceElement(enr).SetOrder(ox, oy);
1863 }
1864
1865
Ng_GetNLevels()1866 int Ng_GetNLevels ()
1867 {
1868 return (mesh) ? mesh->mglevels : 0;
1869 }
1870
1871
Ng_GetParentNodes(int ni,int * parents)1872 void Ng_GetParentNodes (int ni, int * parents)
1873 {
1874 if (ni <= mesh->mlbetweennodes.Size())
1875 {
1876 parents[0] = mesh->mlbetweennodes.Get(ni).I1();
1877 parents[1] = mesh->mlbetweennodes.Get(ni).I2();
1878 }
1879 else
1880 parents[0] = parents[1] = 0;
1881 }
1882
1883
Ng_GetParentElement(int ei)1884 int Ng_GetParentElement (int ei)
1885 {
1886 if (mesh->GetDimension() == 3)
1887 {
1888 if (ei <= mesh->mlparentelement.Size())
1889 return mesh->mlparentelement.Get(ei);
1890 }
1891 else
1892 {
1893 if (ei <= mesh->mlparentsurfaceelement.Size())
1894 return mesh->mlparentsurfaceelement.Get(ei);
1895 }
1896 return 0;
1897 }
1898
1899
Ng_GetParentSElement(int ei)1900 int Ng_GetParentSElement (int ei)
1901 {
1902 if (mesh->GetDimension() == 3)
1903 {
1904 if (ei <= mesh->mlparentsurfaceelement.Size())
1905 return mesh->mlparentsurfaceelement.Get(ei);
1906 }
1907 else
1908 {
1909 return 0;
1910 }
1911 return 0;
1912 }
1913
1914
1915
1916
1917
Ng_GetClusterRepVertex(int pi)1918 int Ng_GetClusterRepVertex (int pi)
1919 {
1920 return mesh->GetClusters().GetVertexRepresentant(pi);
1921 }
1922
Ng_GetClusterRepEdge(int pi)1923 int Ng_GetClusterRepEdge (int pi)
1924 {
1925 return mesh->GetClusters().GetEdgeRepresentant(pi);
1926 }
1927
Ng_GetClusterRepFace(int pi)1928 int Ng_GetClusterRepFace (int pi)
1929 {
1930 return mesh->GetClusters().GetFaceRepresentant(pi);
1931 }
1932
Ng_GetClusterRepElement(int pi)1933 int Ng_GetClusterRepElement (int pi)
1934 {
1935 return mesh->GetClusters().GetElementRepresentant(pi);
1936 }
1937
1938
1939
1940
1941
Ng_GetNPeriodicVertices(int idnr)1942 int Ng_GetNPeriodicVertices (int idnr)
1943 {
1944 Array<INDEX_2> apairs;
1945 mesh->GetIdentifications().GetPairs (idnr, apairs);
1946 return apairs.Size();
1947 }
1948
1949
1950 // pairs should be an integer array of 2*npairs
Ng_GetPeriodicVertices(int idnr,int * pairs)1951 void Ng_GetPeriodicVertices (int idnr, int * pairs)
1952 {
1953 Array<INDEX_2> apairs;
1954 mesh->GetIdentifications().GetPairs (idnr, apairs);
1955 for (int i = 0; i < apairs.Size(); i++)
1956 {
1957 pairs[2*i] = apairs[i].I1();
1958 pairs[2*i+1] = apairs[i].I2();
1959 }
1960
1961 }
1962
1963
1964
Ng_GetNPeriodicEdges(int idnr)1965 int Ng_GetNPeriodicEdges (int idnr)
1966 {
1967 Array<INDEX,PointIndex::BASE> map;
1968 //const MeshTopology & top = mesh->GetTopology();
1969 int nse = mesh->GetNSeg();
1970
1971 int cnt = 0;
1972 // for (int id = 1; id <= mesh->GetIdentifications().GetMaxNr(); id++)
1973 {
1974 mesh->GetIdentifications().GetMap(idnr, map);
1975 //(*testout) << "ident-map " << id << ":" << endl << map << endl;
1976
1977 for (SegmentIndex si = 0; si < nse; si++)
1978 {
1979 PointIndex other1 = map[(*mesh)[si][0]];
1980 PointIndex other2 = map[(*mesh)[si][1]];
1981 // (*testout) << "seg = " << (*mesh)[si] << "; other = "
1982 // << other1 << "-" << other2 << endl;
1983 if (other1 && other2 && mesh->IsSegment (other1, other2))
1984 {
1985 cnt++;
1986 }
1987 }
1988 }
1989 return cnt;
1990 }
1991
Ng_GetPeriodicEdges(int idnr,int * pairs)1992 void Ng_GetPeriodicEdges (int idnr, int * pairs)
1993 {
1994 Array<INDEX,PointIndex::BASE> map;
1995 const MeshTopology & top = mesh->GetTopology();
1996 int nse = mesh->GetNSeg();
1997
1998 int cnt = 0;
1999 // for (int id = 1; id <= mesh->GetIdentifications().GetMaxNr(); id++)
2000 {
2001 mesh->GetIdentifications().GetMap(idnr, map);
2002
2003 //(*testout) << "map = " << map << endl;
2004
2005 for (SegmentIndex si = 0; si < nse; si++)
2006 {
2007 PointIndex other1 = map[(*mesh)[si][0]];
2008 PointIndex other2 = map[(*mesh)[si][1]];
2009 if (other1 && other2 && mesh->IsSegment (other1, other2))
2010 {
2011 SegmentIndex otherseg = mesh->SegmentNr (other1, other2);
2012 pairs[cnt++] = top.GetSegmentEdge (si+1);
2013 pairs[cnt++] = top.GetSegmentEdge (otherseg+1);
2014 }
2015 }
2016 }
2017 }
2018
2019
2020
Ng_PushStatus(const char * str)2021 void Ng_PushStatus (const char * str)
2022 {
2023 PushStatus (MyStr (str));
2024 }
2025
Ng_PopStatus()2026 void Ng_PopStatus ()
2027 {
2028 PopStatus ();
2029 }
2030
Ng_SetThreadPercentage(double percent)2031 void Ng_SetThreadPercentage (double percent)
2032 {
2033 SetThreadPercent (percent);
2034 }
2035
Ng_GetStatus(char ** str,double & percent)2036 void Ng_GetStatus (char ** str, double & percent)
2037 {
2038 MyStr s;
2039 GetStatus(s,percent);
2040 *str = new char[s.Length()+1];
2041 strcpy(*str,s.c_str());
2042 }
2043
2044
Ng_SetTerminate(void)2045 void Ng_SetTerminate(void)
2046 {
2047 multithread.terminate = 1;
2048 }
Ng_UnSetTerminate(void)2049 void Ng_UnSetTerminate(void)
2050 {
2051 multithread.terminate = 0;
2052 }
2053
Ng_ShouldTerminate(void)2054 int Ng_ShouldTerminate(void)
2055 {
2056 return multithread.terminate;
2057 }
2058
Ng_SetRunning(int flag)2059 void Ng_SetRunning(int flag)
2060 {
2061 multithread.running = flag;
2062 }
Ng_IsRunning()2063 int Ng_IsRunning()
2064 {
2065 return multithread.running;
2066 }
2067
2068 ///// Added by Roman Stainko ....
Ng_GetVertex_Elements(int vnr,int * elems)2069 int Ng_GetVertex_Elements( int vnr, int* elems )
2070 {
2071 const MeshTopology& topology = mesh->GetTopology();
2072 ArrayMem<int,4> indexArray;
2073 topology.GetVertexElements( vnr, indexArray );
2074
2075 for( int i=0; i<indexArray.Size(); i++ )
2076 elems[i] = indexArray[i];
2077
2078 return indexArray.Size();
2079 }
2080
2081 ///// Added by Roman Stainko ....
Ng_GetVertex_SurfaceElements(int vnr,int * elems)2082 int Ng_GetVertex_SurfaceElements( int vnr, int* elems )
2083 {
2084 const MeshTopology& topology = mesh->GetTopology();
2085 ArrayMem<int,4> indexArray;
2086 topology.GetVertexSurfaceElements( vnr, indexArray );
2087
2088 for( int i=0; i<indexArray.Size(); i++ )
2089 elems[i] = indexArray[i];
2090
2091 return indexArray.Size();
2092 }
2093
2094 ///// Added by Roman Stainko ....
Ng_GetVertex_NElements(int vnr)2095 int Ng_GetVertex_NElements( int vnr )
2096 {
2097 const MeshTopology& topology = mesh->GetTopology();
2098 ArrayMem<int,4> indexArray;
2099 topology.GetVertexElements( vnr, indexArray );
2100
2101 return indexArray.Size();
2102 }
2103
2104 ///// Added by Roman Stainko ....
Ng_GetVertex_NSurfaceElements(int vnr)2105 int Ng_GetVertex_NSurfaceElements( int vnr )
2106 {
2107 const MeshTopology& topology = mesh->GetTopology();
2108 ArrayMem<int,4> indexArray;
2109 topology.GetVertexSurfaceElements( vnr, indexArray );
2110
2111 return indexArray.Size();
2112 }
2113
2114
2115
2116 #ifdef SOCKETS
Ng_SocketClientOpen(const int port,const char * host)2117 int Ng_SocketClientOpen( const int port, const char * host )
2118 {
2119 try
2120 {
2121 if(host)
2122 clientsocket.Reset(new ClientSocket(port,host));
2123 else
2124 clientsocket.Reset(new ClientSocket(port));
2125 }
2126 catch( SocketException e)
2127 {
2128 cerr << e.Description() << endl;
2129 return 0;
2130 }
2131 return 1;
2132 }
2133
Ng_SocketClientWrite(const char * write,char ** reply)2134 void Ng_SocketClientWrite( const char * write, char** reply)
2135 {
2136 string output = write;
2137 (*clientsocket) << output;
2138 string sreply;
2139 (*clientsocket) >> sreply;
2140 *reply = new char[sreply.size()+1];
2141 strcpy(*reply,sreply.c_str());
2142 }
2143
2144
Ng_SocketClientClose(void)2145 void Ng_SocketClientClose ( void )
2146 {
2147 clientsocket.Reset(NULL);
2148 }
2149
2150
Ng_SocketClientGetServerHost(const int number,char ** host)2151 void Ng_SocketClientGetServerHost ( const int number, char ** host )
2152 {
2153 *host = new char[servers[number]->host.size()+1];
2154 strcpy(*host,servers[number]->host.c_str());
2155 }
2156
Ng_SocketClientGetServerPort(const int number,int * port)2157 void Ng_SocketClientGetServerPort ( const int number, int * port )
2158 {
2159 *port = servers[number]->port;
2160 }
2161
Ng_SocketClientGetServerClientID(const int number,int * id)2162 void Ng_SocketClientGetServerClientID ( const int number, int * id )
2163 {
2164 *id = servers[number]->clientid;
2165 }
2166
2167 #endif // SOCKETS
2168
2169
2170
2171
2172 #ifdef PARALLEL
Ng_SetElementPartition(const int elnr,const int part)2173 void Ng_SetElementPartition ( const int elnr, const int part )
2174 {
2175 mesh->VolumeElement(elnr+1).SetPartition(part);
2176
2177 }
Ng_GetElementPartition(const int elnr)2178 int Ng_GetElementPartition ( const int elnr )
2179 {
2180 return mesh->VolumeElement(elnr+1).GetPartition();
2181 }
2182 #endif
2183
2184
Ng_InitPointCurve(double red,double green,double blue)2185 void Ng_InitPointCurve(double red, double green, double blue)
2186 {
2187 mesh->InitPointCurve(red, green, blue);
2188 }
2189
Ng_AddPointCurvePoint(const double * point)2190 void Ng_AddPointCurvePoint(const double * point)
2191 {
2192 Point3d pt;
2193 pt.X() = point[0];
2194 pt.Y() = point[1];
2195 pt.Z() = point[2];
2196 mesh->AddPointCurvePoint(pt);
2197 }
2198
2199
Ng_SaveMesh(const char * meshfile)2200 void Ng_SaveMesh ( const char * meshfile )
2201 {
2202 mesh -> Save(string(meshfile));
2203 }
2204
2205
Ng_Bisect_WithInfo(const char * refinementfile,double ** qualityloss,int * qualityloss_size)2206 int Ng_Bisect_WithInfo ( const char * refinementfile, double ** qualityloss, int * qualityloss_size )
2207 {
2208 BisectionOptions biopt;
2209 biopt.outfilename = NULL; // "ngfepp.vol";
2210 biopt.femcode = "fepp";
2211 biopt.refinementfilename = refinementfile;
2212
2213 Refinement * ref = const_cast<Refinement*> (&ng_geometry -> GetRefinement());
2214 MeshOptimize2d * opt = NULL;
2215 /*
2216 if (stlgeometry)
2217 ref = new RefinementSTLGeometry(*stlgeometry);
2218 #ifdef OCCGEOMETRY
2219 else if (occgeometry)
2220 ref = new OCCRefinementSurfaces (*occgeometry);
2221 #endif
2222 #ifdef ACIS
2223 else if (acisgeometry)
2224 {
2225 ref = new ACISRefinementSurfaces(*acisgeometry);
2226 opt = new ACISMeshOptimize2dSurfaces(*acisgeometry);
2227 ref->Set2dOptimizer(opt);
2228 }
2229 #endif
2230 else
2231 {
2232 ref = new RefinementSurfaces(*geometry);
2233 opt = new MeshOptimize2dSurfaces(*geometry);
2234 ref->Set2dOptimizer(opt);
2235 }
2236 */
2237 #ifdef ACIS
2238 if (acisgeometry)
2239 {
2240 // ref = new ACISRefinementSurfaces(*acisgeometry);
2241 opt = new ACISMeshOptimize2dSurfaces(*acisgeometry);
2242 ref->Set2dOptimizer(opt);
2243 }
2244 else
2245 #endif
2246 {
2247 // ref = new RefinementSurfaces(*geometry);
2248 CSGeometry * geometry = dynamic_cast<CSGeometry*> (ng_geometry);
2249 if (geometry)
2250 {
2251 opt = new MeshOptimize2dSurfaces(*geometry);
2252 ref->Set2dOptimizer(opt);
2253 }
2254 }
2255
2256 if(!mesh->LocalHFunctionGenerated())
2257 mesh->CalcLocalH(mparam.grading);
2258
2259 mesh->LocalHFunction().SetGrading (mparam.grading);
2260
2261 Array<double> * qualityloss_arr = NULL;
2262 if(qualityloss != NULL)
2263 qualityloss_arr = new Array<double>;
2264
2265 ref -> Bisect (*mesh, biopt, qualityloss_arr);
2266
2267 int retval = 0;
2268
2269 if(qualityloss != NULL)
2270 {
2271 *qualityloss = new double[qualityloss_arr->Size()+1];
2272
2273 for(int i = 0; i<qualityloss_arr->Size(); i++)
2274 (*qualityloss)[i+1] = (*qualityloss_arr)[i];
2275
2276 retval = qualityloss_arr->Size();
2277
2278 delete qualityloss_arr;
2279 }
2280
2281 mesh -> UpdateTopology();
2282 mesh -> GetCurvedElements().BuildCurvedElements (ref, mparam.elementorder);
2283
2284 multithread.running = 0;
2285 delete ref;
2286 delete opt;
2287
2288 return retval;
2289 }
2290
Ng_Bisect(const char * refinementfile)2291 void Ng_Bisect ( const char * refinementfile )
2292 {
2293 Ng_Bisect_WithInfo( refinementfile, NULL, NULL );
2294 }
2295
2296
2297
2298
2299
2300 /*
2301 number of nodes of type nt
2302 nt = 0 is Vertex
2303 nt = 1 is Edge
2304 nt = 2 is Face
2305 nt = 3 is Cell
2306 */
Ng_GetNNodes(int nt)2307 int Ng_GetNNodes (int nt)
2308 {
2309 switch (nt)
2310 {
2311 case 0: return mesh -> GetNV();
2312 case 1: return mesh->GetTopology().GetNEdges();
2313 case 2: return mesh->GetTopology().GetNFaces();
2314 case 3: return mesh -> GetNE();
2315 }
2316 return -1;
2317 }
2318
2319
Ng_GetClosureNodes(int nt,int nodenr,int nodeset,int * nodes)2320 int Ng_GetClosureNodes (int nt, int nodenr, int nodeset, int * nodes)
2321 {
2322 switch (nt)
2323 {
2324 case 3: // The closure of a cell
2325 {
2326 int cnt = 0;
2327 if (nodeset & 1) // Vertices
2328 {
2329 const Element & el = (*mesh)[ElementIndex(nodenr)];
2330 for (int i = 0; i < el.GetNP(); i++)
2331 {
2332 nodes[cnt++] = 0;
2333 nodes[cnt++] = el[i] - PointIndex::BASE;
2334 }
2335 }
2336
2337 if (nodeset & 2) // Edges
2338 {
2339 int edges[12];
2340 int ned;
2341 ned = mesh->GetTopology().GetElementEdges (nodenr+1, edges, 0);
2342 for (int i = 0; i < ned; i++)
2343 {
2344 nodes[cnt++] = 1;
2345 nodes[cnt++] = edges[i]-1;
2346 }
2347 }
2348
2349 if (nodeset & 4) // Faces
2350 {
2351 int faces[12];
2352 int nfa;
2353 nfa = mesh->GetTopology().GetElementFaces (nodenr+1, faces, 0);
2354 for (int i = 0; i < nfa; i++)
2355 {
2356 nodes[cnt++] = 2;
2357 nodes[cnt++] = faces[i]-1;
2358 }
2359 }
2360
2361 if (nodeset & 8) // Cell
2362 {
2363 nodes[cnt++] = 3;
2364 nodes[cnt++] = nodenr;
2365 }
2366
2367 return cnt/2;
2368 }
2369 default:
2370 {
2371 cerr << "GetClosureNodes not implemented for Nodetype " << nt << endl;
2372 }
2373 }
2374 return 0;
2375 }
2376
2377
2378
Ng_GetNElements(int dim)2379 int Ng_GetNElements (int dim)
2380 {
2381 switch (dim)
2382 {
2383 case 0: return mesh -> GetNV();
2384 case 1: return mesh -> GetNSeg();
2385 case 2: return mesh -> GetNSE();
2386 case 3: return mesh -> GetNE();
2387 }
2388 return -1;
2389 }
2390
2391
2392
2393 /*
2394 closure nodes of element
2395 nodeset is bit-coded, bit 0 includes Vertices, bit 1 edges, etc
2396 E.g., nodeset = 6 includes edge and face nodes
2397 nodes is pair of integers (nodetype, nodenr)
2398 return value is number of nodes
2399 */
Ng_GetElementClosureNodes(int dim,int elementnr,int nodeset,int * nodes)2400 int Ng_GetElementClosureNodes (int dim, int elementnr, int nodeset, int * nodes)
2401 {
2402 switch (dim)
2403 {
2404 case 3: // The closure of a volume element = CELL
2405 {
2406 return Ng_GetClosureNodes (3, elementnr, nodeset, nodes);
2407 }
2408 case 2:
2409 {
2410 int cnt = 0;
2411 if (nodeset & 1) // Vertices
2412 {
2413 const Element2d & el = (*mesh)[SurfaceElementIndex(elementnr)];
2414 for (int i = 0; i < el.GetNP(); i++)
2415 {
2416 nodes[cnt++] = 0;
2417 nodes[cnt++] = el[i] - PointIndex::BASE;
2418 }
2419 }
2420
2421 if (nodeset & 2) // Edges
2422 {
2423 int edges[12];
2424 int ned;
2425 ned = mesh->GetTopology().GetSurfaceElementEdges (elementnr+1, edges, 0);
2426 for (int i = 0; i < ned; i++)
2427 {
2428 nodes[cnt++] = 1;
2429 nodes[cnt++] = edges[i]-1;
2430 }
2431 }
2432
2433 if (nodeset & 4) // Faces
2434 {
2435 int face = mesh->GetTopology().GetSurfaceElementFace (elementnr+1);
2436 nodes[cnt++] = 2;
2437 nodes[cnt++] = face-1;
2438 }
2439
2440 return cnt/2;
2441 }
2442 default:
2443 {
2444 cerr << "GetClosureNodes not implemented for Element of dimension " << dim << endl;
2445 }
2446 }
2447 return 0;
2448 }
2449