1 /**************************************************************************/
2 /* File: nglib.cc */
3 /* Author: Joachim Schoeberl */
4 /* Date: 7. May. 2000 */
5 /**************************************************************************/
6
7 /*
8
9 Interface to the netgen meshing kernel
10
11 */
12
13
14 #include <mystdlib.h>
15 #include <myadt.hpp>
16
17 #include <linalg.hpp>
18 #include <csg.hpp>
19 #include <stlgeom.hpp>
20 #include <geometry2d.hpp>
21 #include <meshing.hpp>
22
23
24
25 // #include <FlexLexer.h>
26
27 namespace netgen {
28 extern void MeshFromSpline2D (SplineGeometry2d & geometry,
29 Mesh *& mesh,
30 MeshingParameters & mp);
31 }
32
33
34
35
36
37
38
39 namespace nglib {
40 #include "nglib.h"
41 }
42
43 using namespace netgen;
44
45 // constants and types:
46
47 namespace nglib
48 {
49 // initialize, deconstruct Netgen library:
Ng_Init()50 void Ng_Init ()
51 {
52 mycout = &cout;
53 myerr = &cerr;
54 testout = new ofstream ("test.out");
55 }
56
Ng_Exit()57 void Ng_Exit ()
58 {
59 ;
60 }
61
62
Ng_NewMesh()63 Ng_Mesh * Ng_NewMesh ()
64 {
65 Mesh * mesh = new Mesh;
66 mesh->AddFaceDescriptor (FaceDescriptor (1, 1, 0, 1));
67 return (Ng_Mesh*)(void*)mesh;
68 }
69
Ng_DeleteMesh(Ng_Mesh * mesh)70 void Ng_DeleteMesh (Ng_Mesh * mesh)
71 {
72 delete (Mesh*)mesh;
73 }
74
75 // return bc property for surface element i
76 int
EG_GetSurfaceElementBCProperty(Ng_Mesh * mesh,int i)77 EG_GetSurfaceElementBCProperty(Ng_Mesh * mesh, int i)
78 {
79 int j = ((Mesh*)mesh)->SurfaceElement(i).GetIndex();
80 int k = ((Mesh*)mesh)->GetFaceDescriptor(j).BCProperty();
81 return k;
82 }
83
84 // return surface and volume element in pi
85 Ng_Surface_Element_Type
EG_GetSurfaceElement(Ng_Mesh * mesh,int num,int * pi,int * ptrignum)86 EG_GetSurfaceElement (Ng_Mesh * mesh, int num, int * pi, int * ptrignum)
87 {
88 const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
89 for (int i = 1; i <= el.GetNP(); i++) {
90 pi[i-1] = el.PNum(i);
91 ptrignum[i-1] = el.GeomInfoPi(i).trignum;
92 }
93 Ng_Surface_Element_Type et;
94 switch (el.GetNP())
95 {
96 case 3: et = NG_TRIG; break;
97 case 4: et = NG_QUAD; break;
98 case 6: et = NG_TRIG6; break;
99 }
100 return et;
101 }
102
103 // return segment bcnum
104 void
EG_GetSegmentBCProperty(Ng_Mesh * mesh,Ng_Geometry_2D * geom,int num,int * bcnum)105 EG_GetSegmentBCProperty (Ng_Mesh *mesh, Ng_Geometry_2D * geom, int num, int * bcnum)
106 {
107 const Segment & seg = ((Mesh*)mesh)->LineSegment(num);
108
109 int edgenum = seg.edgenr;
110
111 SplineGeometry2d *geom2d = (SplineGeometry2d*)geom;
112
113 SplineSegment &spline = geom2d->GetSpline(num);
114
115 if(bcnum)
116 *bcnum = spline.bc;
117 }
118
119 // feeds points, surface elements and volume elements to the mesh
Ng_AddPoint(Ng_Mesh * mesh,double * x)120 void Ng_AddPoint (Ng_Mesh * mesh, double * x)
121 {
122 Mesh * m = (Mesh*)mesh;
123 m->AddPoint (Point3d (x[0], x[1], x[2]));
124 }
125
Ng_AddSurfaceElement(Ng_Mesh * mesh,Ng_Surface_Element_Type et,int * pi)126 void Ng_AddSurfaceElement (Ng_Mesh * mesh, Ng_Surface_Element_Type et,
127 int * pi)
128 {
129 Mesh * m = (Mesh*)mesh;
130 Element2d el (3);
131 el.SetIndex (1);
132 el.PNum(1) = pi[0];
133 el.PNum(2) = pi[1];
134 el.PNum(3) = pi[2];
135 m->AddSurfaceElement (el);
136 }
137
Ng_AddVolumeElement(Ng_Mesh * mesh,Ng_Volume_Element_Type et,int * pi)138 void Ng_AddVolumeElement (Ng_Mesh * mesh, Ng_Volume_Element_Type et,
139 int * pi)
140 {
141 Mesh * m = (Mesh*)mesh;
142 Element el (4);
143 el.SetIndex (1);
144 el.PNum(1) = pi[0];
145 el.PNum(2) = pi[1];
146 el.PNum(3) = pi[2];
147 el.PNum(4) = pi[3];
148 m->AddVolumeElement (el);
149 }
150
151 // ask for number of points, surface and volume elements
Ng_GetNP(Ng_Mesh * mesh)152 int Ng_GetNP (Ng_Mesh * mesh)
153 {
154 return ((Mesh*)mesh) -> GetNP();
155 }
156
Ng_GetNSE(Ng_Mesh * mesh)157 int Ng_GetNSE (Ng_Mesh * mesh)
158 {
159 return ((Mesh*)mesh) -> GetNSE();
160 }
161
Ng_GetNE(Ng_Mesh * mesh)162 int Ng_GetNE (Ng_Mesh * mesh)
163 {
164 return ((Mesh*)mesh) -> GetNE();
165 }
166
167
168 // return point coordinates
Ng_GetPoint(Ng_Mesh * mesh,int num,double * x)169 void Ng_GetPoint (Ng_Mesh * mesh, int num, double * x)
170 {
171 const Point3d & p = ((Mesh*)mesh)->Point(num);
172 x[0] = p.X();
173 x[1] = p.Y();
174 x[2] = p.Z();
175 }
176
177 // return surface and volume element in pi
178 Ng_Surface_Element_Type
Ng_GetSurfaceElement(Ng_Mesh * mesh,int num,int * pi)179 Ng_GetSurfaceElement (Ng_Mesh * mesh, int num, int * pi)
180 {
181 const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
182 for (int i = 1; i <= el.GetNP(); i++)
183 pi[i-1] = el.PNum(i);
184 Ng_Surface_Element_Type et;
185 switch (el.GetNP())
186 {
187 case 3: et = NG_TRIG; break;
188 case 4: et = NG_QUAD; break;
189 case 6: et = NG_TRIG6; break;
190 }
191 return et;
192 }
193
194 Ng_Volume_Element_Type
Ng_GetVolumeElement(Ng_Mesh * mesh,int num,int * pi)195 Ng_GetVolumeElement (Ng_Mesh * mesh, int num, int * pi)
196 {
197 const Element & el = ((Mesh*)mesh)->VolumeElement(num);
198 for (int i = 1; i <= el.GetNP(); i++)
199 pi[i-1] = el.PNum(i);
200 Ng_Volume_Element_Type et;
201 switch (el.GetNP())
202 {
203 case 4: et = NG_TET; break;
204 case 5: et = NG_PYRAMID; break;
205 case 6: et = NG_PRISM; break;
206 case 10: et = NG_TET10; break;
207 }
208 return et;
209 }
210
Ng_RestrictMeshSizeGlobal(Ng_Mesh * mesh,double h)211 void Ng_RestrictMeshSizeGlobal (Ng_Mesh * mesh, double h)
212 {
213 ((Mesh*)mesh) -> SetGlobalH (h);
214 }
215
Ng_RestrictMeshSizePoint(Ng_Mesh * mesh,double * p,double h)216 void Ng_RestrictMeshSizePoint (Ng_Mesh * mesh, double * p, double h)
217 {
218 ((Mesh*)mesh) -> RestrictLocalH (Point3d (p[0], p[1], p[2]), h);
219 }
220
Ng_RestrictMeshSizeBox(Ng_Mesh * mesh,double * pmin,double * pmax,double h)221 void Ng_RestrictMeshSizeBox (Ng_Mesh * mesh, double * pmin, double * pmax, double h)
222 {
223 for (double x = pmin[0]; x < pmax[0]; x += h)
224 for (double y = pmin[1]; y < pmax[1]; y += h)
225 for (double z = pmin[2]; z < pmax[2]; z += h)
226 ((Mesh*)mesh) -> RestrictLocalH (Point3d (x, y, z), h);
227 }
228
229
230 // generates volume mesh from surface mesh
Ng_GenerateVolumeMesh(Ng_Mesh * mesh,Ng_Meshing_Parameters * mp)231 Ng_Result Ng_GenerateVolumeMesh (Ng_Mesh * mesh, Ng_Meshing_Parameters * mp)
232 {
233 Mesh * m = (Mesh*)mesh;
234
235
236 MeshingParameters mparam;
237 mparam.maxh = mp->maxh;
238 mparam.meshsizefilename = mp->meshsize_filename;
239
240 m->CalcLocalH();
241
242 MeshVolume (mparam, *m);
243 RemoveIllegalElements (*m);
244 OptimizeVolume (mparam, *m);
245
246 return NG_OK;
247 }
248
249
250
251 // 2D Meshing Functions:
252
253
254
Ng_AddPoint_2D(Ng_Mesh * mesh,double * x)255 void Ng_AddPoint_2D (Ng_Mesh * mesh, double * x)
256 {
257 Mesh * m = (Mesh*)mesh;
258
259 m->AddPoint (Point3d (x[0], x[1], 0));
260 }
261
Ng_AddBoundarySeg_2D(Ng_Mesh * mesh,int pi1,int pi2)262 void Ng_AddBoundarySeg_2D (Ng_Mesh * mesh, int pi1, int pi2)
263 {
264 Mesh * m = (Mesh*)mesh;
265
266 Segment seg;
267 seg.p1 = pi1;
268 seg.p2 = pi2;
269 m->AddSegment (seg);
270 }
271
272
Ng_GetNP_2D(Ng_Mesh * mesh)273 int Ng_GetNP_2D (Ng_Mesh * mesh)
274 {
275 Mesh * m = (Mesh*)mesh;
276 return m->GetNP();
277 }
278
Ng_GetNE_2D(Ng_Mesh * mesh)279 int Ng_GetNE_2D (Ng_Mesh * mesh)
280 {
281 Mesh * m = (Mesh*)mesh;
282 return m->GetNSE();
283 }
284
Ng_GetNSeg_2D(Ng_Mesh * mesh)285 int Ng_GetNSeg_2D (Ng_Mesh * mesh)
286 {
287 Mesh * m = (Mesh*)mesh;
288 return m->GetNSeg();
289 }
290
Ng_GetPoint_2D(Ng_Mesh * mesh,int num,double * x)291 void Ng_GetPoint_2D (Ng_Mesh * mesh, int num, double * x)
292 {
293 Mesh * m = (Mesh*)mesh;
294
295 Point<3> & p = m->Point(num);
296 x[0] = p(0);
297 x[1] = p(1);
298 }
299
Ng_GetElement_2D(Ng_Mesh * mesh,int num,int * pi,int * matnum)300 void Ng_GetElement_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum)
301 {
302 const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num);
303 for (int i = 1; i <= 3; i++)
304 pi[i-1] = el.PNum(i);
305 if (matnum)
306 *matnum = el.GetIndex();
307 }
308
309
Ng_GetSegment_2D(Ng_Mesh * mesh,int num,int * pi,int * matnum)310 void Ng_GetSegment_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum)
311 {
312 const Segment & seg = ((Mesh*)mesh)->LineSegment(num);
313 pi[0] = seg.p1;
314 pi[1] = seg.p2;
315
316 if (matnum)
317 *matnum = seg.edgenr;
318 }
319
320
321
322
Ng_LoadGeometry_2D(const char * filename)323 Ng_Geometry_2D * Ng_LoadGeometry_2D (const char * filename)
324 {
325 SplineGeometry2d * geom = new SplineGeometry2d();
326 geom -> Load (filename);
327 return (Ng_Geometry_2D *)geom;
328 }
329
Ng_GenerateMesh_2D(Ng_Geometry_2D * geom,Ng_Mesh ** mesh,Ng_Meshing_Parameters * mp)330 Ng_Result Ng_GenerateMesh_2D (Ng_Geometry_2D * geom,
331 Ng_Mesh ** mesh,
332 Ng_Meshing_Parameters * mp)
333 {
334 // use global variable mparam
335 // MeshingParameters mparam;
336 mparam.maxh = mp->maxh;
337 mparam.meshsizefilename = mp->meshsize_filename;
338 mparam.quad = mp->quad_dominated;
339
340 Mesh * m;
341 MeshFromSpline2D (*(SplineGeometry2d*)geom, m, mparam);
342
343 cout << m->GetNSE() << " elements, " << m->GetNP() << " points" << endl;
344
345 *mesh = (Ng_Mesh*)m;
346 return NG_OK;
347 }
348
Ng_HP_Refinement(Ng_Geometry_2D * geom,Ng_Mesh * mesh,int levels)349 void Ng_HP_Refinement (Ng_Geometry_2D * geom,
350 Ng_Mesh * mesh,
351 int levels)
352 {
353 Refinement2d ref(*(SplineGeometry2d*)geom);
354 HPRefinement (*(Mesh*)mesh, &ref, levels);
355 }
356
Ng_HP_Refinement(Ng_Geometry_2D * geom,Ng_Mesh * mesh,int levels,double parameter)357 void Ng_HP_Refinement (Ng_Geometry_2D * geom,
358 Ng_Mesh * mesh,
359 int levels, double parameter)
360 {
361 Refinement2d ref(*(SplineGeometry2d*)geom);
362 HPRefinement (*(Mesh*)mesh, &ref, levels, parameter);
363 }
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379 ARRAY<STLReadTriangle> readtrias; //only before initstlgeometry
380 ARRAY<Point<3> > readedges; //only before init stlgeometry
381
Ng_SaveMesh(Ng_Mesh * mesh,const char * filename)382 void Ng_SaveMesh(Ng_Mesh * mesh, const char* filename)
383 {
384 ((Mesh*)mesh)->Save(filename);
385 }
386
Ng_LoadMesh(const char * filename)387 Ng_Mesh * Ng_LoadMesh(const char* filename)
388 {
389 Mesh * mesh = new Mesh;
390 mesh->Load(filename);
391 return ( (Ng_Mesh*)mesh );
392 }
393
394 // loads geometry from STL file
Ng_STL_LoadGeometry(const char * filename,int binary)395 Ng_STL_Geometry * Ng_STL_LoadGeometry (const char * filename, int binary)
396 {
397 int i;
398 STLGeometry geom;
399 STLGeometry* geo;
400 ifstream ist(filename);
401
402 if (binary)
403 {
404 geo = geom.LoadBinary(ist);
405 }
406 else
407 {
408 geo = geom.Load(ist);
409 }
410
411 readtrias.SetSize(0);
412 readedges.SetSize(0);
413
414 Point3d p;
415 Vec3d normal;
416 double p1[3];
417 double p2[3];
418 double p3[3];
419 double n[3];
420
421 Ng_STL_Geometry * geo2 = Ng_STL_NewGeometry();
422
423 for (i = 1; i <= geo->GetNT(); i++)
424 {
425 const STLTriangle& t = geo->GetTriangle(i);
426 p = geo->GetPoint(t.PNum(1));
427 p1[0] = p.X(); p1[1] = p.Y(); p1[2] = p.Z();
428 p = geo->GetPoint(t.PNum(2));
429 p2[0] = p.X(); p2[1] = p.Y(); p2[2] = p.Z();
430 p = geo->GetPoint(t.PNum(3));
431 p3[0] = p.X(); p3[1] = p.Y(); p3[2] = p.Z();
432 normal = t.Normal();
433 n[0] = normal.X(); n[1] = normal.Y(); n[2] = normal.Z();
434
435 Ng_STL_AddTriangle(geo2, p1, p2, p3, n);
436 }
437
438 return geo2;
439 }
440
441 // generate new STL Geometry
Ng_STL_NewGeometry()442 Ng_STL_Geometry * Ng_STL_NewGeometry ()
443 {
444 return (Ng_STL_Geometry*)(void*)new STLGeometry;
445 }
446
447 // after adding triangles (and edges) initialize
Ng_STL_InitSTLGeometry(Ng_STL_Geometry * geom)448 Ng_Result Ng_STL_InitSTLGeometry (Ng_STL_Geometry * geom)
449 {
450 STLGeometry* geo = (STLGeometry*)geom;
451 geo->InitSTLGeometry(readtrias);
452 readtrias.SetSize(0);
453
454 if (readedges.Size() != 0)
455 {
456 int i;
457 /*
458 for (i = 1; i <= readedges.Size(); i+=2)
459 {
460 cout << "e(" << readedges.Get(i) << "," << readedges.Get(i+1) << ")" << endl;
461 }
462 */
463 geo->AddEdges(readedges);
464 }
465
466 if (geo->GetStatus() == STLTopology::STL_GOOD || geo->GetStatus() == STLTopology::STL_WARNING) return NG_OK;
467 return NG_SURFACE_INPUT_ERROR;
468 }
469
470 // automatically generates edges:
Ng_STL_MakeEdges(Ng_STL_Geometry * geom,Ng_Mesh * mesh,Ng_Meshing_Parameters * mp)471 Ng_Result Ng_STL_MakeEdges (Ng_STL_Geometry * geom,
472 Ng_Mesh* mesh,
473 Ng_Meshing_Parameters * mp)
474 {
475 STLGeometry* stlgeometry = (STLGeometry*)geom;
476 Mesh* me = (Mesh*)mesh;
477
478 MeshingParameters mparam;
479
480 mparam.maxh = mp->maxh;
481 mparam.meshsizefilename = mp->meshsize_filename;
482
483 me -> SetGlobalH (mparam.maxh);
484 me -> SetLocalH (stlgeometry->GetBoundingBox().PMin() - Vec3d(10, 10, 10),
485 stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10),
486 0.3);
487
488 me -> LoadLocalMeshSize (mp->meshsize_filename);
489 /*
490 if (mp->meshsize_filename)
491 {
492 ifstream infile (mp->meshsize_filename);
493 if (!infile.good()) return NG_FILE_NOT_FOUND;
494 me -> LoadLocalMeshSize (infile);
495 }
496 */
497
498 STLMeshing (*stlgeometry, *me);
499
500 stlgeometry->edgesfound = 1;
501 stlgeometry->surfacemeshed = 0;
502 stlgeometry->surfaceoptimized = 0;
503 stlgeometry->volumemeshed = 0;
504
505 return NG_OK;
506 }
507
508
509 // generates mesh, empty mesh be already created.
Ng_STL_GenerateSurfaceMesh(Ng_STL_Geometry * geom,Ng_Mesh * mesh,Ng_Meshing_Parameters * mp)510 Ng_Result Ng_STL_GenerateSurfaceMesh (Ng_STL_Geometry * geom,
511 Ng_Mesh* mesh,
512 Ng_Meshing_Parameters * mp)
513 {
514 STLGeometry* stlgeometry = (STLGeometry*)geom;
515 Mesh* me = (Mesh*)mesh;
516
517 MeshingParameters mparam;
518
519 mparam.maxh = mp->maxh;
520 mparam.meshsizefilename = mp->meshsize_filename;
521
522 /*
523 me -> SetGlobalH (mparam.maxh);
524 me -> SetLocalH (stlgeometry->GetBoundingBox().PMin() - Vec3d(10, 10, 10),
525 stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10),
526 0.3);
527 */
528 /*
529 STLMeshing (*stlgeometry, *me);
530
531 stlgeometry->edgesfound = 1;
532 stlgeometry->surfacemeshed = 0;
533 stlgeometry->surfaceoptimized = 0;
534 stlgeometry->volumemeshed = 0;
535 */
536 int retval = STLSurfaceMeshing (*stlgeometry, *me);
537 if (retval == MESHING3_OK)
538 {
539 (*mycout) << "Success !!!!" << endl;
540 stlgeometry->surfacemeshed = 1;
541 stlgeometry->surfaceoptimized = 0;
542 stlgeometry->volumemeshed = 0;
543 }
544 else if (retval == MESHING3_OUTERSTEPSEXCEEDED)
545 {
546 (*mycout) << "ERROR: Give up because of too many trials. Meshing aborted!" << endl;
547 }
548 else if (retval == MESHING3_TERMINATE)
549 {
550 (*mycout) << "Meshing Stopped!" << endl;
551 }
552 else
553 {
554 (*mycout) << "ERROR: Surface meshing not successful. Meshing aborted!" << endl;
555 }
556
557
558 STLSurfaceOptimization (*stlgeometry, *me, mparam);
559
560 return NG_OK;
561 }
562
563
564 // fills STL Geometry
565 // positive orientation
566 // normal vector may be null-pointer
Ng_STL_AddTriangle(Ng_STL_Geometry * geom,double * p1,double * p2,double * p3,double * nv)567 void Ng_STL_AddTriangle (Ng_STL_Geometry * geom,
568 double * p1, double * p2, double * p3, double * nv)
569 {
570 Point<3> apts[3];
571 apts[0] = Point<3>(p1[0],p1[1],p1[2]);
572 apts[1] = Point<3>(p2[0],p2[1],p2[2]);
573 apts[2] = Point<3>(p3[0],p3[1],p3[2]);
574
575 Vec<3> n;
576 if (!nv)
577 n = Cross (apts[0]-apts[1], apts[0]-apts[2]);
578 else
579 n = Vec<3>(nv[0],nv[1],nv[2]);
580
581 readtrias.Append(STLReadTriangle(apts,n));
582 }
583
584 // add (optional) edges:
Ng_STL_AddEdge(Ng_STL_Geometry * geom,double * p1,double * p2)585 void Ng_STL_AddEdge (Ng_STL_Geometry * geom,
586 double * p1, double * p2)
587 {
588 readedges.Append(Point3d(p1[0],p1[1],p1[2]));
589 readedges.Append(Point3d(p2[0],p2[1],p2[2]));
590 }
591
592
593
Ng_Meshing_Parameters()594 Ng_Meshing_Parameters :: Ng_Meshing_Parameters()
595 {
596 maxh = 1000;
597 fineness = 0.5;
598 secondorder = 0;
599 meshsize_filename = 0;
600 quad_dominated = 0;
601 }
602
603
604 }
605
606
607 // compatibility functions:
608
609 namespace netgen
610 {
611
612 char geomfilename[255];
613
MyError(const char * ch)614 void MyError (const char * ch)
615 {
616 cerr << ch;
617 }
618
619 //Destination for messages, errors, ...
Ng_PrintDest(const char * s)620 void Ng_PrintDest(const char * s)
621 {
622 (*mycout) << s << flush;
623 }
624
GetTime()625 double GetTime ()
626 {
627 return 0;
628 }
629
ResetTime()630 void ResetTime ()
631 {
632 ;
633 }
634
MyBeep(int i)635 void MyBeep (int i)
636 {
637 ;
638 }
639
640 }
641