← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to contrib/DiscreteIntegration/Integration3D.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-27 17:36:40 UTC
  • mfrom: (1.2.9 upstream) (8.1.2 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090927173640-oxyhzt0eadjfrlwz
[Christophe Prud'homme]
* New upstream release
  + solver code refactoring
  + better IDE integration.

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/DiscreteIntegration/Integration3D.cpp
498
498
 
499
499
// DI_Point methods ----------------------------------------------------------------------------------------------
500
500
DI_Point::DI_Point (double x, double y, double z, gLevelset &ls) : x_(x), y_(y), z_(z) {
501
 
  Ls.push_back(adjustLs(ls(x, y, z)));
 
501
  addLs(ls(x, y, z));
502
502
}
503
503
DI_Point & DI_Point::operator= (const DI_Point &rhs) {
504
504
  if(this != &rhs){
506
506
  }
507
507
  return *this;
508
508
}
 
509
void DI_Point::addLs (const double ls) {
 
510
  Ls.push_back(adjustLs(ls));
 
511
}
509
512
void DI_Point::addLs (const DI_Element *e) {
510
 
  Ls.push_back(e->evalLs(x_, y_, z_));
 
513
  addLs(e->evalLs(x_, y_, z_));
511
514
}
512
515
void DI_Point::chooseLs (const gLevelset *Lsi) {
513
516
  if(Ls.size() < 2) return;
528
531
}
529
532
void DI_Point::computeLs (const gLevelset &ls) {
530
533
  Ls.clear();
531
 
  double l = ls(x_, y_, z_);
532
 
  Ls.push_back(adjustLs(l));
 
534
  addLs(ls(x_, y_, z_));
533
535
}
534
536
bool DI_Point::equal(const DI_Point &p) const {
535
537
  return (fabs(x() - p.x()) < EQUALITY_TOL && fabs(y() - p.y()) < EQUALITY_TOL && fabs(z() - p.z()) < EQUALITY_TOL);
660
662
  for(int i = 0; i < nbMid(); i++)
661
663
    mid_[i]->addLs(ls[nbVert()+i]);
662
664
}
 
665
void DI_Element::addLs (int primTag, std::map<int, double> nodeLs[8]) {
 
666
  for(int i = 0; i < nbVert(); i++)
 
667
    pts_[i]->addLs(nodeLs[i][primTag]);
 
668
}
663
669
void DI_Element::addLs (const DI_Element *e) {
664
670
  if(e->sizeLs() < 1) return;
665
671
  for(int i = 0; i < nbVert(); i++)
673
679
  }
674
680
  for(int j = 0; j < nbVert(); ++j) {
675
681
    double ls = Ls(e->x(j), e->y(j), e->z(j));
676
 
    pts_[j]->addLs(adjustLs(ls));
 
682
    pts_[j]->addLs(ls);
677
683
  }
678
684
  for(int j = 0; j < nbMid(); ++j) {
679
685
    std::vector<int> s(nbVert()); int n;
682
688
    for(int k = 0; k < n; k++){
683
689
      xc += e->x(s[k]); yc += e->y(s[k]); zc += e->z(s[k]); }
684
690
    double ls = Ls(xc/n, yc/n, zc/n);
685
 
    mid_[j]->addLs(adjustLs(ls));
 
691
    mid_[j]->addLs(ls);
686
692
  }
687
693
}
688
694
void DI_Element::chooseLs (const gLevelset *Lsi) {
1485
1491
// ----------------------------------------------------------------------------
1486
1492
 
1487
1493
// Split a reference line cut by a level set into sublines
1488
 
void DI_Line::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
1494
void DI_Line::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
1489
1495
                          std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cuttingPoints) const
1490
1496
{
1491
1497
  if (!isCrossed(pt(0), pt(1))) {
1514
1520
  subTriangles.push_back(DI_Triangle(pt(2), p02, p12)); // 2-02-12
1515
1521
}
1516
1522
// Split a reference triangle cut by a level set into subtriangles
1517
 
void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
1523
void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
1518
1524
                             std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
1519
1525
                             std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cuttingPoints) const
1520
1526
{
1628
1634
}
1629
1635
 
1630
1636
// Split a reference tetrahedron cut by a level set (defined in a hex) into sub tetrahedra and create triangles on the surface of the level set
1631
 
void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
1637
void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
1632
1638
                       std::vector<DI_Tetra> &subTetras, std::vector<DI_Triangle> &surfTriangles,
1633
1639
                       std::vector<DI_CuttingPoint> &cuttingPoints, std::vector<DI_QualError> &QError) const
1634
1640
{
2043
2049
// -----------------------------------------------------------------------------
2044
2050
 
2045
2051
// cut a line into sublines along the levelset curve
2046
 
void DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
 
2052
bool DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
2047
2053
                   std::vector<DI_CuttingPoint> &cp, const int polynomialOrder,
2048
 
                   std::vector<DI_Line> &lines, int recurLevel) const
 
2054
                   std::vector<DI_Line> &lines, int recurLevel, std::map<int, double> nodeLs[2]) const
2049
2055
{
2050
2056
  //printf(" "); print();
2051
2057
  std::vector<const gLevelset *> RPN, RPNi;
2057
2063
  std::vector<DI_CuttingPoint> ll_cp;
2058
2064
 
2059
2065
  RecurElement re(&ll);
2060
 
  bool signChange = re.cut(recurLevel, this, Ls);
2061
 
  re.pushSubLines(ll_subLines);
 
2066
  bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
 
2067
  pushSubElements(&re, ll_subLines);
2062
2068
 
2063
2069
  if(signChange){
2064
2070
    for(int l = 0; l < (int)RPN.size(); l++) {
2065
 
      const gLevelset *Lsi = RPN[l]; //printf("LS(0,0)=%g LS(1,1)=%g\n",(*Lsi)(0,0,0),(*Lsi)(1,1,0));
 
2071
      const gLevelset *Lsi = RPN[l];
2066
2072
      RPNi.push_back(Lsi);
2067
2073
      if(Lsi->isPrimitive()) {
2068
 
        ll.addLs(this, *Lsi);
 
2074
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) ll.addLs(Lsi->getTag(), nodeLs);
 
2075
        else ll.addLs(this, *Lsi);
2069
2076
        int nbSubLn = ll_subLines.size();
2070
2077
        int nbCP = ll_cp.size();
2071
2078
        for(int i = 0; i < nbSubLn; i++) ll_subLines[i].addLs(&ll);
2090
2097
    for(int l = 0; l < (int)RPN.size(); l++) {
2091
2098
      const gLevelset *Lsi = RPN[l];
2092
2099
      RPNi.push_back(Lsi);
2093
 
      if(Lsi->isPrimitive())
2094
 
        ll.addLs(this, *Lsi);
 
2100
      if(Lsi->isPrimitive()) {
 
2101
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) ll.addLs(Lsi->getTag(), nodeLs);
 
2102
        else ll.addLs(this, *Lsi);
 
2103
      }
2095
2104
    }
2096
2105
  }
2097
2106
 
2110
2119
      if(cp[i].equal(ll_cp[p])) {isIn = true; break;}
2111
2120
    if(!isIn) cp.push_back(ll_cp[p]);
2112
2121
  }
 
2122
  return signChange;
2113
2123
}
2114
2124
 
2115
2125
// cut a line into sublines along one levelset primitive
2116
 
void DI_Line::cut(const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
2126
void DI_Line::cut(const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
2117
2127
           std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cp) const
2118
2128
{
2119
2129
  // check if the line is cut by the level set
2133
2143
}
2134
2144
 
2135
2145
// cut a triangle into subtriangles along the levelset curve
2136
 
void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
 
2146
bool DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
2137
2147
                       std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
2138
2148
                       const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
2139
2149
                       std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
2140
 
                       std::vector<DI_Line> &surfLines, int recurLevel) const
 
2150
                       std::vector<DI_Line> &surfLines, int recurLevel, std::map<int, double> nodeLs[3]) const
2141
2151
{
2142
2152
  //printf(" ");print();
2143
2153
  std::vector<const gLevelset *> RPN, RPNi;
2151
2161
  std::vector<DI_CuttingPoint> tt_cp;
2152
2162
 
2153
2163
  RecurElement re(&tt);
2154
 
  bool signChange = re.cut(recurLevel, this, Ls);
2155
 
  re.pushSubTriangles(tt_subTriangles); //for(int i=0;i<(int)tt_subTriangles.size(); i++) tt_subTriangles[i].print();
 
2164
  bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
 
2165
  pushSubElements(&re, tt_subTriangles);
2156
2166
 
2157
2167
  if(signChange){
2158
2168
    for(int l = 0; l < (int)RPN.size(); l++) {
2159
2169
      const gLevelset *Lsi = RPN[l]; //printf("LS(0,0)=%g LS(1,1)=%g\n",(*Lsi)(0,0,0),(*Lsi)(1,1,0));
2160
2170
      RPNi.push_back(Lsi);
2161
2171
      if(Lsi->isPrimitive()) {
2162
 
        tt.addLs(this, *Lsi);
 
2172
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
 
2173
        else tt.addLs(this, *Lsi);
2163
2174
        int nbSubQ = tt_subQuads.size(), nbSubQ1 = nbSubQ;
2164
2175
        int nbSubTr = tt_subTriangles.size(), nbSubTr1 = nbSubTr;
2165
2176
        int nbSurfLn = tt_surfLines.size(), nbSurfLn1 = nbSurfLn;
2236
2247
    for(int l = 0; l < (int)RPN.size(); l++) {
2237
2248
      const gLevelset *Lsi = RPN[l];
2238
2249
      RPNi.push_back(Lsi);
2239
 
      if(Lsi->isPrimitive())
2240
 
        tt.addLs(this, *Lsi);
 
2250
      if(Lsi->isPrimitive()) {
 
2251
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
 
2252
        else tt.addLs(this, *Lsi);
 
2253
      }
2241
2254
    }
2242
2255
  }
2243
2256
 
2258
2271
      DI_Point vP[3];
2259
2272
      for(int i = 0; i < 3; i++) {
2260
2273
        double ls = Ls(x(i), y(i), z(i));
2261
 
        vP[i] = DI_Point(tt.x(i), tt.y(i), tt.z(i), adjustLs(ls));
 
2274
        vP[i] = DI_Point(tt.x(i), tt.y(i), tt.z(i), ls);
2262
2275
      }
2263
2276
      DI_Triangle trt(vP[0], vP[1], vP[2]); // reference triangle
2264
2277
      tt_subTriangles.push_back(trt);
2265
2278
    }
2266
2279
  }*/
2267
2280
 
2268
 
  //printf("tt = "); tt.printls();
2269
 
 
2270
2281
  for(int q = 0; q < (int)tt_subQuads.size(); q++) {
2271
2282
    tt_subQuads[q].computeLsTagDom(&tt, RPN);
2272
2283
    DI_Quad tt_subQ = tt_subQuads[q];
2282
2293
    subTriangles.push_back(tt_subTr);
2283
2294
  }
2284
2295
  for(int l = 0; l < (int)tt_surfLines.size(); l++) {
2285
 
    tt_surfLines[l].computeLsTagBound(tt_subQuads, tt_subTriangles);  //tt_surfLines[l].printls();
 
2296
    tt_surfLines[l].computeLsTagBound(tt_subQuads, tt_subTriangles);
2286
2297
    if(tt_surfLines[l].lsTag() == -1) continue;
2287
2298
    DI_Line tt_surfLn = tt_surfLines[l];
2288
 
    mappingEl(&tt_surfLn);                      //tt_surfLn.printls();
 
2299
    mappingEl(&tt_surfLn);
2289
2300
    tt_surfLn.integrationPoints(polynomialOrderL, &tt_surfLines[l], &tt, RPN, ipS);
2290
2301
    surfLines.push_back(tt_surfLn);
2291
2302
  }
2297
2308
      if(cp[i].equal(tt_cp[p])) {isIn = true; break;}
2298
2309
    if(!isIn) cp.push_back(tt_cp[p]);
2299
2310
  }
 
2311
  return signChange;
2300
2312
}
2301
2313
 
2302
2314
// cut a triangle into subtriangles along one levelset primitive
2303
 
void DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
2315
void DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
2304
2316
                    std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
2305
2317
                    std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cp) const
2306
2318
{
2325
2337
}
2326
2338
 
2327
2339
// cut a quadrangle into subtriangles along the levelset curve
2328
 
void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
 
2340
bool DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
2329
2341
                   std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
2330
2342
                   const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
2331
2343
                   std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
2332
 
                   std::vector<DI_Line> &surfLines, int recurLevel) const
 
2344
                   std::vector<DI_Line> &surfLines, int recurLevel, std::map<int, double> nodeLs[4]) const
2333
2345
{
2334
2346
  printf(" "); printls();
2335
2347
  std::vector<const gLevelset *> RPN, RPNi;
2343
2355
  std::vector<DI_CuttingPoint> qq_cp;
2344
2356
 
2345
2357
  RecurElement re(&qq);
2346
 
  bool signChange = re.cut(recurLevel, this, Ls, 1.e-5);
2347
 
  re.pushSubQuads(qq_subQuads);
2348
 
 
2349
 
  qq.addLs(this, Ls); printf("qq:"); qq.printls();
2350
 
  for(int i = 0; i < (int)qq_subQuads.size(); i++) {
2351
 
    qq_subQuads[i].addLs(&qq);
2352
 
    printf("qi:"); qq_subQuads[i].printls();
2353
 
    qq_subQuads[i].clearLs();
2354
 
  }
2355
 
  qq.clearLs();
 
2358
  bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
 
2359
  pushSubElements(&re, qq_subQuads);
2356
2360
 
2357
2361
  if(signChange) {
2358
2362
    for(int l = 0; l < (int)RPN.size(); l++) {
2359
2363
      const gLevelset *Lsi = RPN[l];
2360
2364
      RPNi.push_back(Lsi);
2361
2365
      if(Lsi->isPrimitive()) {
2362
 
        qq.addLs(this, *Lsi);
 
2366
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) qq.addLs(Lsi->getTag(), nodeLs);
 
2367
        else qq.addLs(this, *Lsi);
2363
2368
        int nbSubQ = qq_subQuads.size(), nbSubQ1 = nbSubQ;
2364
2369
        int nbSubTr = qq_subTriangles.size(), nbSubTr1 = nbSubTr;
2365
2370
        int nbSurfLn = qq_surfLines.size(), nbSurfLn1 = nbSurfLn;
2440
2445
    for(int l = 0; l < (int)RPN.size(); l++) {
2441
2446
      const gLevelset *Lsi = RPN[l];
2442
2447
      RPNi.push_back(Lsi);
2443
 
      if(Lsi->isPrimitive())
2444
 
        qq.addLs(this, *Lsi);
 
2448
      if(Lsi->isPrimitive()) {
 
2449
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) qq.addLs(Lsi->getTag(), nodeLs);
 
2450
        else qq.addLs(this, *Lsi);
 
2451
      }
2445
2452
    }
2446
2453
  }
2447
2454
 
2459
2466
      DI_Point vP[4];
2460
2467
      for(int i = 0; i < 4; i++) {
2461
2468
        double ls = Ls(x(i), y(i), z(i));
2462
 
        vP[i] = DI_Point(qq.x(i), qq.y(i), qq.z(i), adjustLs(ls));
 
2469
        vP[i] = DI_Point(qq.x(i), qq.y(i), qq.z(i), ls);
2463
2470
      }
2464
2471
      DI_Quad qt(vP[0], vP[1], vP[2], vP[3]); // reference quad
2465
2472
      qq_subQuads.push_back(qt);
2496
2503
      if(cp[i].equal(qq_cp[p])) {isIn = true; break;}
2497
2504
    if(!isIn) cp.push_back(qq_cp[p]);
2498
2505
  }
 
2506
  return signChange;
2499
2507
}
2500
2508
 
2501
2509
// cut a quadrangle into subtriangles along the levelset curve
2502
 
void DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
2510
void DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
2503
2511
                std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
2504
2512
                std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cp) const
2505
2513
{
2535
2543
}
2536
2544
 
2537
2545
// cut a tetrahedron into subtetrahedra along the levelset surface
2538
 
void DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
 
2546
bool DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
2539
2547
                    std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
2540
2548
                    const int polynomialOrderT, const int polynomialOrderQ, const int polynomialOrderTr,
2541
2549
                    std::vector<DI_Tetra> &subTetras, std::vector<DI_Quad> &surfQuads,
2542
 
                    std::vector<DI_Triangle> &surfTriangles, int recurLevel) const
 
2550
                    std::vector<DI_Triangle> &surfTriangles, int recurLevel, std::map<int, double> nodeLs[4]) const
2543
2551
{
2544
2552
  //printf(" ");print();
2545
2553
  std::vector<const gLevelset *> RPN, RPNi;
2556
2564
  std::vector<DI_QualError> QError;
2557
2565
 
2558
2566
  RecurElement re(&tt);
2559
 
  bool signChange = re.cut(recurLevel, this, Ls);
2560
 
  re.pushSubTetras(tt_subTetras);
 
2567
  bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
 
2568
  pushSubElements(&re, tt_subTetras);
2561
2569
 
2562
2570
  if(signChange) {
2563
2571
    for(int l = 0; l < (int)RPN.size(); l++) {
2564
2572
      const gLevelset *Lsi = RPN[l];
2565
2573
      RPNi.push_back(Lsi);
2566
2574
      if(Lsi->isPrimitive()) {
2567
 
        tt.addLs(this, *Lsi);
 
2575
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
 
2576
        else tt.addLs(this, *Lsi);
2568
2577
        int nbSubT = tt_subTetras.size(), nbSubT1 = nbSubT;
2569
2578
        int nbSurfQ = tt_surfQuads.size();
2570
2579
        int nbSurfTr = tt_surfTriangles.size(), nbSurfTr1 = nbSurfTr;
2624
2633
    for(int l = 0; l < (int)RPN.size(); l++) {
2625
2634
      const gLevelset *Lsi = RPN[l];
2626
2635
      RPNi.push_back(Lsi);
2627
 
      if(Lsi->isPrimitive())
2628
 
        tt.addLs(this, *Lsi);
 
2636
      if(Lsi->isPrimitive()) {
 
2637
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
 
2638
        else tt.addLs(this, *Lsi);
 
2639
      }
2629
2640
    }
2630
2641
  }
2631
2642
 
2644
2655
      DI_Point vP[4];
2645
2656
      for(int i = 0; i < 4; i++) {
2646
2657
        double ls = Ls(x(i), y(i), z(i));
2647
 
        vP[i] = DI_Point(tt.x(i), tt.y(i), tt.z(i), adjustLs(ls));
 
2658
        vP[i] = DI_Point(tt.x(i), tt.y(i), tt.z(i), ls);
2648
2659
      }
2649
2660
      DI_Tetra tet(vP[0], vP[1], vP[2], vP[3]); // reference tetra
2650
2661
      tt_subTetras.push_back(tet);
2685
2696
      if(cp[i].equal(tt_cp[p])) {isIn = true; break;}
2686
2697
    if(!isIn) cp.push_back(tt_cp[p]);
2687
2698
  }
 
2699
  return signChange;
2688
2700
}
2689
2701
 
2690
2702
// cut a tetrahedron into subtetrahedra along the levelset surface
2691
 
void DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
2703
void DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
2692
2704
                 std::vector<DI_Tetra> &subTetras, std::vector<DI_Quad> &surfQuads,
2693
2705
                 std::vector<DI_Triangle> &surfTriangles, std::vector<DI_CuttingPoint> &cp,
2694
2706
                 std::vector<DI_QualError> &QError) const
2713
2725
}
2714
2726
 
2715
2727
// cut a hex into subtetras along the levelset surface
2716
 
void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
 
2728
bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
2717
2729
                   std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
2718
2730
                   const int polynomialOrderH, const int polynomialOrderT,
2719
2731
                   const int polynomialOrderQ, const int polynomialOrderTr,
2720
2732
                   std::vector<DI_Hexa> &subHexas, std::vector<DI_Tetra> &subTetras,
2721
2733
                   std::vector<DI_Quad> &surfQuads, std::vector<DI_Triangle> &surfTriangles,
2722
 
                   std::vector<DI_Line> &frontLines, int recurLevel) const
 
2734
                   std::vector<DI_Line> &frontLines, int recurLevel, std::map<int, double> nodeLs[8]) const
2723
2735
{
2724
2736
  printf(" "); print();
2725
2737
 
2736
2748
  std::vector<DI_QualError> QError;
2737
2749
 
2738
2750
  RecurElement re(&hh);
2739
 
  bool signChange = re.cut(recurLevel, this, Ls);
2740
 
  re.pushSubHexas(hh_subHexas);
 
2751
  bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
 
2752
  pushSubElements(&re, hh_subHexas);
2741
2753
 
2742
2754
  if(signChange){
2743
2755
    for(int l = 0; l < (int)RPN.size(); l++) {
2744
2756
      const gLevelset *Lsi = RPN[l];
2745
2757
      RPNi.push_back(Lsi);
2746
2758
      if(Lsi->isPrimitive()) {
2747
 
        hh.addLs(this, *Lsi);
 
2759
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) hh.addLs(Lsi->getTag(), nodeLs);
 
2760
        else hh.addLs(this, *Lsi);
2748
2761
        int nbSubH = hh_subHexas.size();
2749
2762
        int nbSubT = hh_subTetras.size(), nbSubT1 = nbSubT;
2750
2763
        int nbSurfQ = hh_surfQuads.size(), nbSurfQ1 = nbSurfQ;
2858
2871
    for(int l = 0; l < (int)RPN.size(); l++) {
2859
2872
      const gLevelset *Lsi = RPN[l];
2860
2873
      RPNi.push_back(Lsi);
2861
 
      if(Lsi->isPrimitive())
2862
 
        hh.addLs(this, *Lsi);
 
2874
      if(Lsi->isPrimitive()) {
 
2875
        if(nodeLs && nodeLs[0].count(Lsi->getTag())) hh.addLs(Lsi->getTag(), nodeLs);
 
2876
        else hh.addLs(this, *Lsi);
 
2877
      }
2863
2878
    }
2864
2879
  }
2865
2880
 
2881
2896
      DI_Point vP[8];
2882
2897
      for(int i = 0; i < 8; i++) {
2883
2898
        double ls = Ls(x(i), y(i), z(i));
2884
 
        vP[i] = DI_Point(hh.x(i), hh.y(i), hh.z(i), adjustLs(ls));
 
2899
        vP[i] = DI_Point(hh.x(i), hh.y(i), hh.z(i), ls);
2885
2900
      }
2886
2901
      DI_Hexa hex(vP[0], vP[1], vP[2], vP[3], vP[4], vP[5], vP[6], vP[7]); // reference hexa
2887
2902
      hh_subHexas.push_back(hex);
2937
2952
      if(cp[i].equal(hh_cp[p])) {isIn = true; break;}
2938
2953
    if(!isIn) cp.push_back(hh_cp[p]);
2939
2954
  }
 
2955
  return signChange;
2940
2956
}
2941
2957
 
2942
2958
// cut a hex into subtetras along the levelset surface
2943
 
void DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
 
2959
void DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
2944
2960
                   std::vector<DI_Hexa> &Hexas, std::vector<DI_Tetra> &subTetras,
2945
2961
                   std::vector<DI_Quad> &surfQuads, std::vector<DI_Triangle> &surfTriangles,
2946
2962
                   std::vector<DI_CuttingPoint> &cp, std::vector<DI_QualError> &QError) const