← Back to branch summary

~ubuntu-branches/ubuntu/dapper/gmsh/dapper

« back to all changes in this revision

Viewing changes to Netgen/libsrc/interface/nginterface.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2005-09-23 23:24:14 UTC
  • mfrom: (2.1.1 sarge)
  • Revision ID: james.westby@ubuntu.com-20050923232414-iao03ey38rd8pg0m
Tags: 1.60.1-1ubuntu2
gcc4 build failure - meshtype.hpp. (Closes BTS #324310)

Show diffs side-by-side

added added

removed removed

Lines of Context:
Netgen/libsrc/interface/nginterface.cpp
6
6
#include <geometry2d.hpp>
7
7
#include <stlgeom.hpp>
8
8
 
 
9
#ifdef OCCGEOMETRY
 
10
#include <occgeom.hpp>
 
11
#endif
 
12
 
 
13
 
9
14
#include <visual.hpp>
10
15
 
11
16
#include "nginterface.h"
12
17
// #include <FlexLexer.h>
13
18
 
14
19
 
 
20
// #include <mystdlib.h>
15
21
 
16
22
 
17
23
namespace netgen
20
26
  extern VisualSceneMesh vsmesh;
21
27
  extern Tcl_Interp * tcl_interp;
22
28
 
23
 
  extern SplineGeometry2d * geometry2d;
24
 
  extern CSGeometry * parsegeom;
25
 
  extern CSGeometry * geometry;
 
29
  extern AutoPtr<SplineGeometry2d> geometry2d;
 
30
  extern AutoPtr<CSGeometry> geometry;
26
31
  extern STLGeometry * stlgeometry;
 
32
#ifdef OCCGEOMETRY
 
33
  extern OCCGeometry * occgeometry;
 
34
#endif
27
35
 
 
36
#ifdef OPENGL
28
37
  extern VisualSceneSolution vssolution;
 
38
#endif
29
39
  extern CSGeometry * ParseCSG (istream & istr);
30
40
}
31
41
 
47
57
{
48
58
  ifstream infile (filename);
49
59
 
50
 
  if (geometry) 
51
 
    delete geometry;
52
 
 
53
 
  if (strcmp (&filename[strlen(filename)-3], "geo") == 0)
 
60
  geometry.Reset ();
 
61
  geometry2d.Reset ();
 
62
 
 
63
#ifdef OCCGEOMETRY
 
64
  delete occgeometry;
 
65
  occgeometry = 0;
 
66
#endif
 
67
 
 
68
  if ((strcmp (&filename[strlen(filename)-3], "geo") == 0) ||
 
69
      (strcmp (&filename[strlen(filename)-3], "GEO") == 0) ||
 
70
      (strcmp (&filename[strlen(filename)-3], "Geo") == 0))
54
71
    {
55
 
      /*
56
 
      geometry = new CSGeometry(filename);
57
 
      
58
 
      parsegeom = geometry;
59
 
      
60
 
      lexer = new yyFlexLexer (&infile);
61
 
      
62
 
      extern int yyparse ();
63
 
      yyparse ();
64
 
      delete lexer; 
65
 
      */
66
 
 
67
 
      geometry = netgen::ParseCSG (infile);
 
72
      geometry.Reset (netgen::ParseCSG (infile));
68
73
 
69
74
      if (!geometry)
70
 
        throw NgException ("input file not found");
 
75
        {
 
76
          geometry.Reset (new CSGeometry ());
 
77
          throw NgException ("input file not found");
 
78
        }
71
79
 
72
80
      geometry -> FindIdenticSurfaces(1e-6);
 
81
 
 
82
      double detail = atof (Tcl_GetVar (tcl_interp, "geooptions.detail", 0));
 
83
      double facets = atof (Tcl_GetVar (tcl_interp, "geooptions.facets", 0));
73
84
      Box<3> box (geometry->BoundingBox());
74
 
 
75
 
      geometry->CalcTriangleApproximation (box, 0.01, 10);
76
 
    }
 
85
      
 
86
      if (atoi (Tcl_GetVar (tcl_interp, "geooptions.drawcsg", 0)))
 
87
        geometry->CalcTriangleApproximation(box, detail, facets);
 
88
 
 
89
      //      geometry->CalcTriangleApproximation (box, 0.01, 10);
 
90
    }
 
91
 
 
92
  else if (strcmp (&filename[strlen(filename)-4], "in2d") == 0)
 
93
    {
 
94
      geometry2d.Reset (new SplineGeometry2d());
 
95
      geometry2d -> Load (filename);
 
96
    }
 
97
 
 
98
  else if ((strcmp (&filename[strlen(filename)-3], "stl") == 0) ||
 
99
           (strcmp (&filename[strlen(filename)-3], "STL") == 0) ||
 
100
           (strcmp (&filename[strlen(filename)-3], "Stl") == 0))
 
101
    {
 
102
      ifstream infile(filename);
 
103
      stlgeometry = STLGeometry :: Load (infile);
 
104
      stlgeometry->edgesfound = 0;
 
105
      Mesh meshdummy;
 
106
      stlgeometry->Clear();
 
107
      stlgeometry->BuildEdges();
 
108
      stlgeometry->MakeAtlas(meshdummy);
 
109
      stlgeometry->CalcFaceNums();
 
110
      stlgeometry->AddFaceEdges();
 
111
      stlgeometry->LinkEdges();
 
112
    }
 
113
 
 
114
#ifdef OCCGEOMETRY
 
115
  else if ((strcmp (&filename[strlen(filename)-4], "iges") == 0) ||
 
116
           (strcmp (&filename[strlen(filename)-3], "igs") == 0) ||
 
117
           (strcmp (&filename[strlen(filename)-3], "IGS") == 0) ||
 
118
           (strcmp (&filename[strlen(filename)-4], "IGES") == 0))
 
119
    {
 
120
      PrintMessage (1, "Load IGES geometry file ", filename);
 
121
      occgeometry = LoadOCC_IGES (filename);
 
122
    }
 
123
  else if ((strcmp (&filename[strlen(filename)-4], "step") == 0) ||
 
124
           (strcmp (&filename[strlen(filename)-3], "stp") == 0) ||
 
125
           (strcmp (&filename[strlen(filename)-3], "STP") == 0) ||
 
126
           (strcmp (&filename[strlen(filename)-4], "STEP") == 0))
 
127
    {
 
128
      PrintMessage (1, "Load STEP geometry file ", filename);
 
129
      occgeometry = LoadOCC_STEP (filename);
 
130
    }
 
131
#endif
77
132
  else
78
133
    {
79
 
      geometry = new CSGeometry("");
80
 
 
81
 
      if (strcmp (&filename[strlen(filename)-4], "in2d") == 0)
82
 
        {
83
 
 
84
 
          if (geometry2d)
85
 
            delete geometry2d;
86
 
          geometry2d = new SplineGeometry2d();
87
 
          geometry2d -> Load (filename);
88
 
        }
89
 
 
90
 
      else
91
 
        {
92
 
          cerr << "Unknown geometry extension!!" << endl;
93
 
        }
 
134
      cerr << "Unknown geometry extension!!" << endl;
94
135
    }
95
136
}                          
96
137
 
293
334
    {
294
335
      const Segment & seg = mesh->LineSegment (ei);
295
336
 
296
 
      if (!seg.pmid)
 
337
      if (seg.pmid < 0)
297
338
        {
298
339
          epi[0] = seg.p1;
299
340
          epi[1] = seg.p2;
329
370
  nv[0] = 0; 
330
371
  nv[1] = 0;
331
372
  nv[2] = 1;
 
373
 
 
374
  (*testout) << "Ng_GetNormalVector (sei = " << sei << ", locpi = " << locpi << ")" << endl;
332
375
  
333
376
  if (mesh->GetDimension() == 3)
334
377
    {
337
380
      p = mesh->Point (mesh->SurfaceElement(sei).PNum(locpi));
338
381
 
339
382
      int surfi = mesh->GetFaceDescriptor(mesh->SurfaceElement(sei).GetIndex()).SurfNr();
 
383
 
 
384
      (*testout) << "surfi = " << surfi << endl;
 
385
#ifdef OCCGEOMETRY
 
386
      if (occgeometry)
 
387
        {
 
388
          PointGeomInfo gi = mesh->SurfaceElement(sei).GeomInfoPi(locpi);
 
389
          occgeometry->GetSurface (surfi).GetNormalVector(p, gi, n);
 
390
          nv[0] = n(0);
 
391
          nv[1] = n(1);
 
392
          nv[2] = n(2);
 
393
        }
 
394
      else
 
395
#endif
340
396
      if (geometry)
341
397
        {
342
 
          geometry->GetSurface (surfi) -> GetNormalVector(p, n);
 
398
          (*testout) << "geometry defined" << endl;
 
399
          n = geometry->GetSurface (surfi) -> GetNormalVector(p);
 
400
          (*testout) << "aus is" << endl;
343
401
          nv[0] = n(0);
344
402
          nv[1] = n(1);
345
403
          nv[2] = n(2);
348
406
}
349
407
 
350
408
 
351
 
int Ng_FindElementOfPoint (double * p, double * lami)
 
409
int Ng_FindElementOfPoint (double * p, double * lami, int build_searchtree, int index)
352
410
{
353
411
  if (mesh->GetDimension() == 3)
354
412
    {
355
413
      Point3d p3d(p[0], p[1], p[2]);
356
414
      int ind = 
357
 
        mesh->GetElementOfPoint(p3d, lami);
 
415
        mesh->GetElementOfPoint(p3d, lami, build_searchtree != 0, index);
358
416
      return ind;
359
417
    }
360
418
  else
362
420
      double lam3[3];
363
421
      Point3d p2d(p[0], p[1], 0);
364
422
      int ind = 
365
 
        mesh->GetElementOfPoint(p2d, lam3);
 
423
        mesh->GetElementOfPoint(p2d, lam3, build_searchtree != 0, index);
366
424
      lami[0] = lam3[0];
367
425
      lami[1] = lam3[1];
368
426
      return ind;
382
440
 
383
441
      mesh->GetCurvedElements().CalcSurfaceTransformation (xl, ei-1, xg, dx);
384
442
 
385
 
      // still 1-based arrays
386
443
      if (x)
387
444
        {
388
445
          for (int i = 0; i < 2; i++)
471
528
 
472
529
 
473
530
 
 
531
void Ng_GetSurfaceElementNeighbouringDomains(const int selnr, int & in, int & out)
 
532
{
 
533
  in = mesh->GetFaceDescriptor((*mesh)[static_cast<SurfaceElementIndex>(selnr)].GetIndex()).DomainIn();
 
534
  out = mesh->GetFaceDescriptor((*mesh)[static_cast<SurfaceElementIndex>(selnr)].GetIndex()).DomainOut();
 
535
}
 
536
 
474
537
 
475
538
void Ng_SetRefinementFlag (int ei, int flag)
476
539
{
491
554
{
492
555
  BisectionOptions biopt;
493
556
  biopt.usemarkedelements = 1;
 
557
  biopt.refine_p = 0;
494
558
  biopt.refine_hp = 0;
495
 
  if (reftype == NG_REFINE_P || reftype == NG_REFINE_HP)
 
559
  if (reftype == NG_REFINE_P)
 
560
    biopt.refine_p = 1;
 
561
  if (reftype == NG_REFINE_HP)
496
562
    biopt.refine_hp = 1;
 
563
  Refinement * ref;
497
564
 
498
 
  if (geometry && mesh->GetDimension() == 3)
499
 
    {
500
 
      RefinementSurfaces ref (*geometry);
501
 
      ref.Bisect (*mesh, biopt);
502
 
    }
 
565
  if (geometry2d)
 
566
    ref = new Refinement2d(*geometry2d);
503
567
  else if (stlgeometry)
504
 
    {
505
 
      RefinementSTLGeometry ref (*stlgeometry);
506
 
      ref.Bisect (*mesh, biopt);
507
 
    }
508
 
  else if (geometry2d)
509
 
    {
510
 
      Refinement2d ref (*geometry2d);
511
 
      ref.Bisect (*mesh, biopt);
512
 
    }
 
568
    ref = new RefinementSTLGeometry(*stlgeometry);
 
569
#ifdef OCCGEOMETRY
 
570
  else if (occgeometry)
 
571
    ref = new OCCRefinementSurfaces (*occgeometry);
 
572
#endif
 
573
  else if (geometry && mesh->GetDimension() == 3)
 
574
    ref = new RefinementSurfaces(*geometry);
513
575
  else
514
576
    {
515
 
      cout << "No geometry available" << endl;
516
 
      Refinement ref;
517
 
      ref.Bisect (*mesh, biopt);
 
577
      ref = new Refinement();
518
578
    }
519
579
 
 
580
  ref -> Bisect (*mesh, biopt);
 
581
 
520
582
  mesh -> UpdateTopology();
 
583
  // mesh -> GetCurvedElements().BuildCurvedElements (ref, mparam.elementorder);
 
584
  delete ref;
521
585
}
522
586
 
523
587
void Ng_SecondOrder ()
552
616
 
553
617
void Ng_HPRefinement (int levels)
554
618
{
555
 
  HPRefinement (*mesh, levels);
 
619
  Refinement * ref;
 
620
 
 
621
  if (stlgeometry)
 
622
    ref = new RefinementSTLGeometry (*stlgeometry);
 
623
  else if (geometry2d)
 
624
    ref = new Refinement2d (*geometry2d);
 
625
  else
 
626
    ref = new RefinementSurfaces (*geometry);
 
627
 
 
628
 
 
629
  HPRefinement (*mesh, ref, levels);
556
630
}
557
631
 
558
632
 
562
636
 
563
637
  if (stlgeometry)
564
638
    ref = new RefinementSTLGeometry (*stlgeometry);
 
639
#ifdef OCCGEOMETRY
 
640
  else if (occgeometry)
 
641
    ref = new OCCRefinementSurfaces (*occgeometry);
 
642
#endif
565
643
  else if (geometry2d)
566
644
    ref = new Refinement2d (*geometry2d);
567
645
  else
1019
1097
}
1020
1098
 
1021
1099
 
 
1100
int Ng_GetNVertexElements (int vnr)
 
1101
{
 
1102
  if (mesh->GetDimension() == 3)
 
1103
    return mesh->GetTopology().GetVertexElements(vnr).Size();
 
1104
  else
 
1105
    return mesh->GetTopology().GetVertexSurfaceElements(vnr).Size();
 
1106
}
 
1107
 
 
1108
void Ng_GetVertexElements (int vnr, int * els)
 
1109
{
 
1110
  FlatArray<int> ia(0,0);
 
1111
  if (mesh->GetDimension() == 3)
 
1112
    ia = mesh->GetTopology().GetVertexElements(vnr);
 
1113
  else
 
1114
    ia = mesh->GetTopology().GetVertexSurfaceElements(vnr);
 
1115
  for (int i = 0; i < ia.Size(); i++)
 
1116
    els[i] = ia[i];
 
1117
}
 
1118
 
 
1119
 
1022
1120
int Ng_GetElementOrder (int enr)
1023
1121
{
1024
1122
  if (mesh->GetDimension() == 3)
1122
1220
 
1123
1221
void Ng_SetSolutionData (Ng_SolutionData * soldata)
1124
1222
{
 
1223
#ifdef OPENGL
1125
1224
  //   vssolution.ClearSolutionData ();
1126
1225
  VisualSceneSolution::SolData * vss = new VisualSceneSolution::SolData;
1127
1226
 
1140
1239
  vss->soltype = VisualSceneSolution::SolType (soldata->soltype);
1141
1240
  vss->solclass = soldata->solclass;
1142
1241
  vssolution.AddSolutionData (vss);
1143
 
}
 
1242
#endif
 
1243
}
 
1244
 
 
1245
void Ng_ClearSolutionData ()
 
1246
{
 
1247
  vssolution.ClearSolutionData();
 
1248
}
 
1249
 
1144
1250
 
1145
1251
 
1146
1252
void Ng_Redraw ()
1147
1253
{
 
1254
#ifdef OPENGL
1148
1255
  vssolution.UpdateSolutionTimeStamp();
1149
1256
  Render();
 
1257
#endif
1150
1258
}
1151
1259
 
1152
1260
 
1153
1261
void Ng_SetVisualizationParameter (const char * name, const char * value)
1154
1262
{
 
1263
#ifdef OPENGL
1155
1264
  char buf[100];
1156
1265
  sprintf (buf, "visoptions.%s", name);
1157
1266
  cout << "name = " << name << ", value = " << value << endl;
1158
1267
  cout << "set tcl-variable " << buf << " to " << value << endl;
1159
1268
  Tcl_SetVar (tcl_interp, buf, const_cast<char*> (value), 0);
1160
1269
  Tcl_Eval (tcl_interp, "Ng_Vis_Set parameters;");
 
1270
#endif
1161
1271
}
1162
1272
 
1163
1273
 
1198
1308
  for (i = 1; i <= ne; i++)
1199
1309
    {
1200
1310
      int j;
1201
 
      Element2d tri(3);
 
1311
      Element2d tri(TRIG);
1202
1312
      tri.SetIndex(1); //faceind
1203
1313
      
1204
1314
      for (j = 1; j <= 3; j++)
1240
1350
      pairs[2*i+1] = apairs[i].I2();
1241
1351
    }
1242
1352
      
1243
 
        /*<<<<<< nginterface.cpp
1244
 
      infile >> np;
1245
 
      for (i = 1; i <= np; i++)
1246
 
        {
1247
 
          Point3d p;
1248
 
          infile >> p.X() >> p.Y() >> p.Z();
1249
 
          if (firsttime)
1250
 
            mesh->AddPoint (p);
1251
 
          else
1252
 
            mesh->Point(i)=p;
1253
 
        }
1254
 
 
1255
 
      //firsttime = 0;
1256
 
      Ng_Redraw();
1257
 
      //}
1258
 
}
1259
 
 
1260
 
                
1261
 
int Ng_GetNPeriodicVertices ()
1262
 
{
1263
 
  ARRAY<INDEX_2> apairs;
1264
 
  mesh->GetIdentifications().GetPairs (0, apairs);
1265
 
  return apairs.Size();
1266
 
}
1267
 
 
1268
 
 
1269
 
// pairs should be an integer array of 2*npairs
1270
 
void Ng_GetPeriodicVertices (int * pairs)
1271
 
{
1272
 
  ARRAY<INDEX_2> apairs;
1273
 
  mesh->GetIdentifications().GetPairs (0, apairs);
1274
 
  for (int i = 0; i < apairs.Size(); i++)
1275
 
    {
1276
 
      pairs[2*i] = apairs[i].I1();
1277
 
      pairs[2*i+1] = apairs[i].I2();
1278
 
    }
1279
 
=======
1280
 
      if (firsttime)
1281
 
        mesh->AddSurfaceElement (tri);
1282
 
    }
1283
 
  
1284
 
  infile >> np;
1285
 
  for (i = 1; i <= np; i++)
1286
 
    {
1287
 
      Point3d p;
1288
 
      infile >> p.X() >> p.Y() >> p.Z();
1289
 
      if (firsttime)
1290
 
        mesh->AddPoint (p);
1291
 
      else
1292
 
        mesh->Point(i)=p;
1293
 
    }
1294
 
  
1295
 
  //firsttime = 0;
1296
 
  Ng_Redraw();
1297
 
>>>>>>> 1.9 
1298
 
*/
1299
 
}
 
1353
}
 
1354
 
 
1355
 
 
1356
 
 
1357
int Ng_GetNPeriodicEdges ()
 
1358
{
 
1359
  ARRAY<INDEX,PointIndex::BASE> map;
 
1360
  const MeshTopology & top = mesh->GetTopology();
 
1361
  int nse = mesh->GetNSeg();
 
1362
 
 
1363
  int cnt = 0;
 
1364
  //  for (int id = 1; id <= mesh->GetIdentifications().GetMaxNr(); id++)
 
1365
    {
 
1366
      mesh->GetIdentifications().GetMap(0, map);
 
1367
      //(*testout) << "ident-map " << id << ":" << endl << map << endl;
 
1368
 
 
1369
      for (SegmentIndex si = 0; si < nse; si++)
 
1370
        {
 
1371
          PointIndex other1 = map[(*mesh)[si].p1];
 
1372
          PointIndex other2 = map[(*mesh)[si].p2];
 
1373
          //  (*testout) << "seg = " << (*mesh)[si] << "; other = " 
 
1374
          //     << other1 << "-" << other2 << endl;
 
1375
          if (other1 && other2 && mesh->IsSegment (other1, other2))
 
1376
            {
 
1377
              cnt++;
 
1378
            }
 
1379
        }
 
1380
    }
 
1381
  return cnt;
 
1382
}
 
1383
 
 
1384
void Ng_GetPeriodicEdges (int * pairs)
 
1385
{
 
1386
  ARRAY<INDEX,PointIndex::BASE> map;
 
1387
  const MeshTopology & top = mesh->GetTopology();
 
1388
  int nse = mesh->GetNSeg();
 
1389
 
 
1390
  int cnt = 0;
 
1391
  //  for (int id = 1; id <= mesh->GetIdentifications().GetMaxNr(); id++)
 
1392
    {
 
1393
      mesh->GetIdentifications().GetMap(0, map);
 
1394
      
 
1395
      //(*testout) << "map = " << map << endl;
 
1396
 
 
1397
      for (SegmentIndex si = 0; si < nse; si++)
 
1398
        {
 
1399
          PointIndex other1 = map[(*mesh)[si].p1];
 
1400
          PointIndex other2 = map[(*mesh)[si].p2];
 
1401
          if (other1 && other2 && mesh->IsSegment (other1, other2))
 
1402
            {
 
1403
              SegmentIndex otherseg = mesh->SegmentNr (other1, other2);
 
1404
              pairs[cnt++] = top.GetSegmentEdge (si+1);
 
1405
              pairs[cnt++] = top.GetSegmentEdge (otherseg+1);
 
1406
            }
 
1407
        }
 
1408
    }
 
1409
}
 
1410
 
 
1411
 
 
1412
 
 
1413
void Ng_PushStatus (const char * str)
 
1414
{
 
1415
  PushStatus (MyStr (str));
 
1416
}
 
1417
 
 
1418
void Ng_PopStatus ()
 
1419
{
 
1420
  PopStatus ();
 
1421
}
 
1422
 
 
1423
void Ng_SetThreadPercentage (double percent)
 
1424
{
 
1425
  SetThreadPercent (percent);
 
1426
}
 
1427
 
 
1428
 
 
1429
///// Added by Roman Stainko ....
 
1430
int Ng_GetVertex_Elements( int vnr, int* elems )
 
1431
{
 
1432
  const MeshTopology& topology = mesh->GetTopology();
 
1433
  ArrayMem<int,4> indexArray;
 
1434
  topology.GetVertexElements( vnr, indexArray );
 
1435
  
 
1436
  for( int i=0; i<indexArray.Size(); i++ )
 
1437
    elems[i] = indexArray[i];
 
1438
  
 
1439
  return indexArray.Size();
 
1440
}
 
1441
 
 
1442
///// Added by Roman Stainko ....
 
1443
int Ng_GetVertex_SurfaceElements( int vnr, int* elems )
 
1444
{
 
1445
  const MeshTopology& topology = mesh->GetTopology();
 
1446
  ArrayMem<int,4> indexArray;
 
1447
  topology.GetVertexSurfaceElements( vnr, indexArray );
 
1448
  
 
1449
  for( int i=0; i<indexArray.Size(); i++ )
 
1450
    elems[i] = indexArray[i];
 
1451
  
 
1452
  return indexArray.Size();
 
1453
}
 
1454
 
 
1455
///// Added by Roman Stainko ....
 
1456
int Ng_GetVertex_NElements( int vnr )
 
1457
{
 
1458
  const MeshTopology& topology = mesh->GetTopology();
 
1459
  ArrayMem<int,4> indexArray;
 
1460
  topology.GetVertexElements( vnr, indexArray );
 
1461
  
 
1462
  return indexArray.Size();
 
1463
}
 
1464
 
 
1465
///// Added by Roman Stainko ....
 
1466
int Ng_GetVertex_NSurfaceElements( int vnr )
 
1467
{
 
1468
  const MeshTopology& topology = mesh->GetTopology();
 
1469
  ArrayMem<int,4> indexArray;
 
1470
  topology.GetVertexSurfaceElements( vnr, indexArray );
 
1471
 
 
1472
  return indexArray.Size();
 
1473
}
 
1474
 
1300
1475
 
1301
1476