"Fossies" - the Fresh Open Source Software Archive

Member "gmsh-4.3.0-source/contrib/Netgen/libsrc/meshing/adfront3.cpp" (23 Feb 2019, 18223 Bytes) of package /linux/privat/gmsh-4.3.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "adfront3.cpp" see the Fossies "Dox" file reference documentation.

    1 #include <mystdlib.h>
    2 #include "meshing.hpp"
    3 
    4 
    5 /* ********************** FrontPoint ********************** */
    6 
    7 namespace netgen
    8 {
    9 
   10 FrontPoint3 :: FrontPoint3 () 
   11 { 
   12   globalindex = -1;
   13   nfacetopoint = 0; 
   14   frontnr = 1000; 
   15   cluster = 0;
   16 }
   17 
   18 
   19 FrontPoint3 :: FrontPoint3 (const Point<3> & ap, PointIndex agi)
   20 { 
   21   p = ap; 
   22   globalindex = agi;
   23   nfacetopoint = 0; 
   24   frontnr = 1000; 
   25   cluster = 0;
   26 }
   27 
   28 
   29 
   30 /* ********************** FrontFace ********************** */
   31 
   32 FrontFace :: FrontFace () 
   33 { 
   34   qualclass = 1; 
   35   oldfront = 0; 
   36   hashvalue = 0;
   37   cluster = 0;
   38 }
   39 
   40 FrontFace :: FrontFace (const MiniElement2d & af)
   41 { 
   42   f = af; 
   43   oldfront = 0; 
   44   qualclass = 1; 
   45   hashvalue = 0;
   46 }
   47 
   48 void FrontFace :: Invalidate ()
   49 { 
   50   f.Delete(); 
   51   oldfront = 0; 
   52   qualclass = 1000; 
   53 }
   54 
   55 
   56 
   57 
   58 /* ********************** AddFront ********************** */
   59  
   60 
   61 AdFront3 :: AdFront3 ()
   62 {
   63   nff = 0;
   64   nff4 = 0;
   65   vol = 0;
   66 
   67   hashon = 1;
   68   hashcreated = 0;
   69   if (hashon) 
   70     hashtable.Init(&points, &faces);
   71 
   72   facetree = NULL;
   73   connectedpairs = NULL;
   74 
   75   rebuildcounter = -1;
   76   lasti = 0;
   77   minval = -1;
   78 }
   79 
   80 
   81 AdFront3 :: ~AdFront3 ()
   82 {
   83   delete facetree;
   84   delete connectedpairs;
   85 }
   86 
   87 void AdFront3 :: GetPoints (Array<Point<3> > & apoints) const
   88 {
   89   for (PointIndex pi = PointIndex::BASE; 
   90        pi < points.Size()+PointIndex::BASE; pi++)
   91     
   92     apoints.Append (points[pi].P());
   93 }
   94 
   95 
   96 PointIndex AdFront3 :: AddPoint (const Point<3> & p, PointIndex globind)
   97 {
   98   if (delpointl.Size())
   99     {
  100       PointIndex pi = delpointl.Last();
  101       delpointl.DeleteLast ();
  102       
  103       points[pi] = FrontPoint3 (p, globind);
  104       return pi;
  105     }
  106   else
  107     {
  108       points.Append (FrontPoint3 (p, globind));
  109       return points.Size()-1+PointIndex::BASE;
  110     }
  111 }
  112 
  113 
  114 INDEX AdFront3 :: AddFace (const MiniElement2d & aface)
  115 {
  116   int i, minfn;
  117 
  118   nff++;
  119 
  120   for (i = 0; i < aface.GetNP(); i++)
  121     points[aface[i]].AddFace();
  122 
  123   const Point3d & p1 = points[aface[0]].P();
  124   const Point3d & p2 = points[aface[1]].P();
  125   const Point3d & p3 = points[aface[2]].P();
  126 
  127   vol += 1.0/6.0 * (p1.X() + p2.X() + p3.X()) *
  128     ( (p2.Y() - p1.Y()) * (p3.Z() - p1.Z()) -
  129       (p2.Z() - p1.Z()) * (p3.Y() - p1.Y()) );
  130 
  131   if (aface.GetNP() == 4)
  132     {
  133       nff4++;
  134       const Point3d & p4 = points[aface[3]].P();      
  135       vol += 1.0/6.0 * (p1.X() + p3.X() + p4.X()) *
  136     ( (p3.Y() - p1.Y()) * (p4.Z() - p1.Z()) -
  137       (p3.Z() - p1.Z()) * (p4.Y() - p1.Y()) );
  138     }
  139 
  140 
  141   minfn = 1000;
  142   for (i = 0; i < aface.GetNP(); i++)
  143     {
  144       int fpn = points[aface[i]].FrontNr();
  145       if (i == 0 || fpn < minfn)
  146     minfn = fpn;
  147     }
  148 
  149 
  150   int cluster = 0;
  151   for (i = 1; i <= aface.GetNP(); i++)
  152     {
  153       if (points[aface.PNum(i)].cluster)
  154     cluster = points[aface.PNum(i)].cluster;
  155     }
  156   for (i = 1; i <= aface.GetNP(); i++)
  157     points[aface.PNum(i)].cluster = cluster;
  158 
  159 
  160   for (i = 1; i <= aface.GetNP(); i++)
  161     points[aface.PNum(i)].DecFrontNr (minfn+1);
  162 
  163   int nfn = faces.Append(FrontFace (aface));
  164   faces.Elem(nfn).cluster = cluster;
  165 
  166   if (hashon && hashcreated) 
  167     hashtable.AddElem(aface, nfn);
  168 
  169   return nfn;
  170 }
  171 
  172 
  173 
  174 void AdFront3 :: DeleteFace (INDEX fi)
  175 {
  176   nff--;
  177 
  178   for (int i = 1; i <= faces.Get(fi).Face().GetNP(); i++)
  179     {
  180       PointIndex pi = faces.Get(fi).Face().PNum(i);
  181       points[pi].RemoveFace();
  182       if (!points[pi].Valid())
  183     delpointl.Append (pi);
  184     }
  185 
  186   const MiniElement2d & face = faces.Get(fi).Face();
  187   const Point3d & p1 = points[face.PNum(1)].P();
  188   const Point3d & p2 = points[face.PNum(2)].P();
  189   const Point3d & p3 = points[face.PNum(3)].P();
  190 
  191   vol -= 1.0/6.0 * (p1.X() + p2.X() + p3.X()) *
  192     ( (p2.Y() - p1.Y()) * (p3.Z() - p1.Z()) -
  193       (p2.Z() - p1.Z()) * (p3.Y() - p1.Y()) );
  194 
  195   if (face.GetNP() == 4)
  196     {
  197       const Point3d & p4 = points[face.PNum(4)].P();      
  198       vol -= 1.0/6.0 * (p1.X() + p3.X() + p4.X()) *
  199     ( (p3.Y() - p1.Y()) * (p4.Z() - p1.Z()) -
  200       (p3.Z() - p1.Z()) * (p4.Y() - p1.Y()) );
  201 
  202       nff4--;
  203     }
  204 
  205   faces.Elem(fi).Invalidate();
  206 }
  207 
  208 
  209 INDEX AdFront3 :: AddConnectedPair (const INDEX_2 & apair)
  210 {
  211   if (!connectedpairs)
  212     connectedpairs = new TABLE<int, PointIndex::BASE> (GetNP());
  213 
  214   connectedpairs->Add (apair.I1(), apair.I2());
  215   connectedpairs->Add (apair.I2(), apair.I1());
  216 
  217   return 0;
  218 }
  219 
  220 
  221 void AdFront3 :: CreateTrees ()
  222 {
  223   int i, j;
  224   PointIndex pi;
  225   Point3d pmin, pmax;
  226 
  227   for (pi = PointIndex::BASE; 
  228        pi < GetNP()+PointIndex::BASE; pi++)
  229     {
  230       const Point<3> & p = GetPoint(pi);
  231       if (pi == PointIndex::BASE)
  232     {
  233       pmin = p;
  234       pmax = p;
  235     }
  236       else
  237     {
  238       pmin.SetToMin (p);
  239       pmax.SetToMax (p);
  240     }
  241     }
  242 
  243   pmax = pmax + 0.5 * (pmax - pmin);
  244   pmin = pmin + 0.5 * (pmin - pmax);
  245 
  246   delete facetree;
  247   facetree = new Box3dTree (pmin, pmax);
  248   
  249   for (i = 1; i <= GetNF(); i++)
  250     {
  251       const MiniElement2d & el = GetFace(i);
  252       pmin = GetPoint (el[0]);
  253       pmax = pmin;
  254       for (j = 1; j < 3; j++)
  255     {
  256       const Point<3> & p = GetPoint (el[j]);
  257       pmin.SetToMin (p);
  258       pmax.SetToMax (p);
  259     }
  260       pmax = pmax + 0.01 * (pmax - pmin);
  261       pmin = pmin + 0.01 * (pmin - pmax);
  262       //      (*testout) << "insert " << i << ": " << pmin << " - " << pmax << "\n";
  263       facetree -> Insert (pmin, pmax, i);
  264     }
  265 }
  266 
  267 
  268 void AdFront3 :: GetIntersectingFaces (const Point<3> & pmin, const Point<3> & pmax, 
  269                        Array<int> & ifaces) const
  270 {
  271   facetree -> GetIntersecting (pmin, pmax, ifaces);
  272 }
  273 
  274 void AdFront3 :: GetFaceBoundingBox (int i, Box3d & box) const
  275 {
  276   const FrontFace & face = faces.Get(i);
  277   box.SetPoint (points[face.f[0]].p);
  278   box.AddPoint (points[face.f[1]].p);
  279   box.AddPoint (points[face.f[2]].p);
  280 }
  281 
  282 void AdFront3 :: RebuildInternalTables ()
  283 {
  284   static int timer_a = NgProfiler::CreateTimer ("Adfront3::RebuildInternal A");
  285   static int timer_b = NgProfiler::CreateTimer ("Adfront3::RebuildInternal B");
  286   static int timer_c = NgProfiler::CreateTimer ("Adfront3::RebuildInternal C");
  287   static int timer_d = NgProfiler::CreateTimer ("Adfront3::RebuildInternal D");
  288 
  289 
  290   NgProfiler::StartTimer (timer_a);   
  291   int hi = 0;
  292   for (int i = 1; i <= faces.Size(); i++)
  293     if (faces.Get(i).Valid())
  294       {
  295     hi++;
  296     if (hi < i)
  297       faces.Elem(hi) = faces.Get(i);
  298       }
  299   
  300   faces.SetSize (nff);
  301 
  302   int np = points.Size();
  303 
  304   for (int i = PointIndex::BASE; 
  305        i < np+PointIndex::BASE; i++)
  306     points[i].cluster = i;
  307   
  308   NgProfiler::StopTimer (timer_a);    
  309   NgProfiler::StartTimer (timer_b);   
  310 
  311   int change;
  312   do
  313     {
  314       change = 0;
  315       for (int i = 1; i <= faces.Size(); i++)
  316     {
  317       const MiniElement2d & el = faces.Get(i).Face();
  318 
  319       int mini = points[el.PNum(1)].cluster;
  320       int maxi = mini;
  321       
  322       for (int j = 2; j <= 3; j++)
  323         {
  324           int ci = points[el.PNum(j)].cluster;
  325           if (ci < mini) mini = ci;
  326           if (ci > maxi) maxi = ci;
  327         }
  328 
  329       if (mini < maxi)
  330         {
  331           change = 1;
  332           for (int j = 1; j <= 3; j++)
  333         points[el.PNum(j)].cluster = mini;
  334         }
  335     }
  336     }
  337   while (change);
  338 
  339 
  340   NgProfiler::StopTimer (timer_b);    
  341   NgProfiler::StartTimer (timer_c);   
  342 
  343 
  344 
  345 
  346   BitArrayChar<PointIndex::BASE> usecl(np);
  347   usecl.Clear();
  348   for (int i = 1; i <= faces.Size(); i++)
  349     {
  350       usecl.Set (points[faces.Get(i).Face().PNum(1)].cluster);
  351       faces.Elem(i).cluster =
  352     points[faces.Get(i).Face().PNum(1)].cluster;
  353     }
  354   int cntcl = 0;
  355   for (int i = PointIndex::BASE; 
  356        i < np+PointIndex::BASE; i++)
  357     if (usecl.Test(i))
  358       cntcl++;
  359 
  360   Array<double, PointIndex::BASE> clvol (np);
  361   clvol = 0.0;
  362 
  363   for (int i = 1; i <= faces.Size(); i++)
  364     {
  365       const MiniElement2d & face = faces.Get(i).Face();
  366 
  367       const Point3d p1 = points[face.PNum(1)].P();      
  368       const Point3d p2 = points[face.PNum(2)].P();      
  369       const Point3d p3 = points[face.PNum(3)].P();      
  370       
  371       double vi = 1.0/6.0 * (p1.X() + p2.X() + p3.X()) *
  372     ( (p2.Y() - p1.Y()) * (p3.Z() - p1.Z()) -
  373       (p2.Z() - p1.Z()) * (p3.Y() - p1.Y()) );
  374       
  375       if (face.GetNP() == 4)
  376     {
  377       const Point3d p4 = points[face.PNum(4)].P();      
  378       vi += 1.0/6.0 * (p1.X() + p3.X() + p4.X()) *
  379         ( (p3.Y() - p1.Y()) * (p4.Z() - p1.Z()) -
  380           (p3.Z() - p1.Z()) * (p4.Y() - p1.Y()) );
  381     }
  382      
  383       clvol[faces.Get(i).cluster] += vi;
  384     }
  385 
  386   NgProfiler::StopTimer (timer_c);    
  387   NgProfiler::StartTimer (timer_d);   
  388 
  389 
  390 
  391   int negvol = 0;
  392   for (int i = PointIndex::BASE; 
  393        i < clvol.Size()+PointIndex::BASE; i++)
  394     {
  395       if (clvol[i] < 0)
  396     negvol = 1;
  397     }
  398 
  399   if (negvol)
  400     {
  401       for (int i = 1; i <= faces.Size(); i++)
  402     faces.Elem(i).cluster = 1;
  403       for (int i = PointIndex::BASE; 
  404        i < points.Size()+PointIndex::BASE; i++)
  405     points[i].cluster = 1;
  406     }
  407 
  408   if (hashon) 
  409     hashtable.Create();
  410 
  411   NgProfiler::StopTimer (timer_d);    
  412 }
  413 
  414 
  415 
  416 int AdFront3 :: SelectBaseElement ()
  417 {
  418   int i, hi, fstind;
  419 
  420   /*
  421   static int minval = -1;
  422   static int lasti = 0;
  423   static int counter = 0;
  424   */
  425   if (rebuildcounter <= 0)
  426     {
  427       RebuildInternalTables();
  428       rebuildcounter = nff / 10 + 1;
  429       
  430       lasti = 0;
  431     }
  432   rebuildcounter--;
  433 
  434   /*
  435   if (faces.Size() > 2 * nff)
  436     {
  437       // compress facelist
  438 
  439       RebuildInternalTables ();
  440       lasti = 0;
  441     }
  442     */
  443   
  444   fstind = 0;
  445 
  446   for (i = lasti+1; i <= faces.Size() && !fstind; i++)
  447     if (faces.Elem(i).Valid())
  448       {
  449     hi = faces.Get(i).QualClass() +
  450       points[faces.Get(i).Face().PNum(1)].FrontNr() +
  451       points[faces.Get(i).Face().PNum(2)].FrontNr() +
  452       points[faces.Get(i).Face().PNum(3)].FrontNr();
  453     
  454     if (hi <= minval)
  455       {
  456         minval = hi;
  457         fstind = i;
  458         lasti = fstind;
  459       }
  460       }
  461   
  462   if (!fstind)
  463     {
  464       minval = INT_MAX;
  465       for (i = 1; i <= faces.Size(); i++)
  466     if (faces.Elem(i).Valid())
  467       {
  468         hi = faces.Get(i).QualClass() +
  469           points[faces.Get(i).Face().PNum(1)].FrontNr() +
  470           points[faces.Get(i).Face().PNum(2)].FrontNr() +
  471           points[faces.Get(i).Face().PNum(3)].FrontNr();
  472         
  473         if (hi <= minval)
  474           {
  475         minval = hi;
  476         fstind = i;
  477         lasti = 0;
  478           }
  479       }
  480     }
  481 
  482 
  483   return fstind;
  484 }
  485 
  486 
  487 
  488 int AdFront3 :: GetLocals (int fstind,
  489                Array<Point3d > & locpoints,
  490                Array<MiniElement2d> & locfaces,   // local index
  491                Array<PointIndex> & pindex,
  492                Array<INDEX> & findex,
  493                INDEX_2_HASHTABLE<int> & getconnectedpairs,
  494                float xh,
  495                float relh,
  496                INDEX& facesplit)
  497 {
  498   static int timer = NgProfiler::CreateTimer ("AdFront3::GetLocals");
  499   NgProfiler::RegionTimer reg (timer);
  500 
  501 
  502   if (hashon && faces.Size() < 500) { hashon=0; }
  503   if (hashon && !hashcreated) 
  504     {
  505       hashtable.Create(); 
  506       hashcreated=1;
  507     }
  508 
  509   INDEX i, j;
  510   PointIndex pstind;
  511   INDEX pi;
  512   Point3d midp, p0;
  513 
  514   //  static Array<int, PointIndex::BASE> invpindex;
  515   
  516   Array<MiniElement2d> locfaces2;           //all local faces in radius xh
  517   Array<int> locfaces3;           // all faces in outer radius relh
  518   Array<INDEX> findex2;
  519 
  520   locfaces2.SetSize(0);
  521   locfaces3.SetSize(0);
  522   findex2.SetSize(0);
  523 
  524   int cluster = faces.Get(fstind).cluster;
  525 
  526   pstind = faces.Get(fstind).Face().PNum(1);
  527   p0 = points[pstind].P();
  528   
  529   locfaces2.Append(faces.Get(fstind).Face());
  530   findex2.Append(fstind);
  531 
  532 
  533   Box3d b1 (p0 - Vec3d(xh, xh, xh), p0 + Vec3d (xh, xh, xh));
  534 
  535   if (hashon)
  536     {
  537       hashtable.GetLocals(locfaces2, findex2, fstind, p0, xh);
  538     }
  539   else
  540     {
  541       for (i = 1; i <= faces.Size(); i++)
  542     {
  543       const MiniElement2d & face = faces.Get(i).Face();
  544       if (faces.Get(i).cluster == cluster && faces.Get(i).Valid() && i != fstind)
  545         {
  546           Box3d b2;
  547           b2.SetPoint (points[face[0]].P());
  548           b2.AddPoint (points[face[1]].P());
  549           b2.AddPoint (points[face[2]].P());
  550 
  551           if (b1.Intersect (b2))
  552         {
  553           locfaces2.Append(faces.Get(i).Face());
  554           findex2.Append(i);
  555         }
  556         }
  557     }
  558     }
  559 
  560   //local faces for inner radius:
  561   for (i = 1; i <= locfaces2.Size(); i++)
  562     {
  563       const MiniElement2d & face = locfaces2.Get(i);
  564       const Point3d & p1 = points[face[0]].P();
  565       const Point3d & p2 = points[face[1]].P();
  566       const Point3d & p3 = points[face[2]].P();
  567 
  568       midp = Center (p1, p2, p3);
  569 
  570       if (Dist2 (midp, p0) <= relh * relh || i == 1)
  571     {
  572           locfaces.Append(locfaces2.Get(i));
  573       findex.Append(findex2.Get(i));
  574     }
  575       else
  576     locfaces3.Append (i);
  577     }
  578   
  579   facesplit=locfaces.Size();
  580   
  581   
  582   //local faces for outer radius:
  583   for (i = 1; i <= locfaces3.Size(); i++)
  584     {
  585       locfaces.Append (locfaces2.Get(locfaces3.Get(i)));
  586       findex.Append (findex2.Get(locfaces3.Get(i)));
  587     }
  588 
  589 
  590   invpindex.SetSize (points.Size());
  591   for (i = 1; i <= locfaces.Size(); i++)
  592     for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
  593       {
  594     pi = locfaces.Get(i).PNum(j);
  595     invpindex[pi] = -1;
  596       }
  597 
  598   for (i = 1; i <= locfaces.Size(); i++)
  599     {
  600       for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
  601     {
  602       pi = locfaces.Get(i).PNum(j);
  603       if (invpindex[pi] == -1)
  604         {
  605           pindex.Append (pi);
  606           invpindex[pi] = pindex.Size();  // -1+PointIndex::BASE;
  607           locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P());
  608         }
  609       else
  610         locfaces.Elem(i).PNum(j) = invpindex[pi];
  611 
  612     }
  613     }
  614 
  615 
  616 
  617   if (connectedpairs)
  618     {
  619       for (i = 1; i <= locpoints.Size(); i++)
  620     {
  621       int pind = pindex.Get(i);
  622       if (pind >= 1 && pind <= connectedpairs->Size ())
  623         {
  624           for (j = 1; j <= connectedpairs->EntrySize(pind); j++)
  625         {
  626           int oi = connectedpairs->Get(pind, j);
  627           int other = invpindex.Get(oi);
  628           if (other >= 1 && other <= pindex.Size() && 
  629               pindex.Get(other) == oi)
  630             {
  631               // INDEX_2 coned(i, other);
  632               // coned.Sort();
  633               // (*testout) << "connected: " << locpoints.Get(i) << "-" << locpoints.Get(other) << endl;
  634               getconnectedpairs.Set (INDEX_2::Sort (i, other), 1);
  635             }
  636         }
  637         }
  638     }
  639     }
  640   
  641 
  642   /*
  643     // add isolated points
  644   for (i = 1; i <= points.Size(); i++)
  645     if (points.Elem(i).Valid() && Dist (points.Elem(i).P(), p0) <= xh)
  646       {
  647     if (!invpindex.Get(i))
  648       {
  649         locpoints.Append (points.Get(i).P());
  650         pindex.Append (i);
  651         invpindex.Elem(i) = pindex.Size();
  652       }
  653       }
  654       */
  655   return faces.Get(fstind).QualClass();
  656 }
  657 
  658 
  659 // returns all points connected with fi
  660 void AdFront3 :: GetGroup (int fi,
  661                Array<MeshPoint> & grouppoints,
  662                Array<MiniElement2d> & groupelements,
  663                Array<PointIndex> & pindex,
  664                Array<INDEX> & findex) 
  665 {
  666   // static Array<char> pingroup;
  667   int i, j, changed;
  668 
  669   pingroup.SetSize(points.Size());
  670 
  671   pingroup = 0;
  672   for (j = 1; j <= 3; j++)
  673     pingroup.Elem (faces.Get(fi).Face().PNum(j)) = 1;
  674 
  675   do
  676     {
  677       changed = 0;
  678 
  679       for (i = 1; i <= faces.Size(); i++)
  680     if (faces.Get(i).Valid())
  681       {
  682         const MiniElement2d & face = faces.Get(i).Face();
  683 
  684         int fused = 0;
  685         for (j = 1; j <= 3; j++)
  686           if (pingroup.Elem(face.PNum(j))) 
  687         fused++;
  688             
  689         if (fused >= 2)
  690           for (j = 1; j <= 3; j++)
  691         if (!pingroup.Elem(face.PNum(j)))
  692           {
  693             pingroup.Elem(face.PNum(j)) = 1;
  694             changed = 1;
  695           }
  696       }
  697 
  698     }
  699   while (changed);
  700 
  701 
  702   //  static Array<int> invpindex;
  703   invpindex.SetSize (points.Size());
  704   
  705 
  706   for (i = 1; i <= points.Size(); i++)
  707     if (points.Get(i).Valid())
  708       {
  709     grouppoints.Append (points.Get(i).P());
  710     pindex.Append (i);
  711     invpindex.Elem(i) = pindex.Size();
  712       }
  713 
  714   for (i = 1; i <= faces.Size(); i++)
  715     if (faces.Get(i).Valid())
  716       {
  717     int fused = 0;
  718     for (j = 1; j <= 3; j++)
  719       if (pingroup.Get(faces.Get(i).Face().PNum(j)))
  720         fused++;
  721 
  722     if (fused >= 2)
  723       {
  724         groupelements.Append (faces.Get(i).Face());
  725         findex.Append (i);
  726       }
  727       }
  728       
  729   for (i = 1; i <= groupelements.Size(); i++)
  730     for (j = 1; j <= 3; j++)
  731       {
  732     groupelements.Elem(i).PNum(j) =
  733       invpindex.Get(groupelements.Elem(i).PNum(j));
  734       }
  735 
  736   /*
  737   for (i = 1; i <= groupelements.Size(); i++)
  738     for (j = 1; j <= 3; j++)
  739       for (k = 1; k <= grouppoints.Size(); k++)
  740         if (pindex.Get(k) == groupelements.Get(i).PNum(j))
  741           {
  742         groupelements.Elem(i).PNum(j) = k;
  743         break;
  744           }
  745   */          
  746 }
  747 
  748 
  749 void AdFront3 :: SetStartFront (int /* baseelnp */)
  750 {
  751   INDEX i;
  752   int j;
  753 
  754   for (i = 1; i <= faces.Size(); i++)
  755     if (faces.Get(i).Valid())
  756       {
  757     const MiniElement2d & face = faces.Get(i).Face();
  758     for (j = 1; j <= 3; j++)
  759       points[face.PNum(j)].DecFrontNr(0);
  760       }
  761 
  762   /*
  763   if (baseelnp)
  764     {
  765       for (i = 1; i <= faces.Size(); i++)
  766     if (faces.Get(i).Valid() && faces.Get(i).Face().GetNP() != baseelnp)
  767       faces.Elem(i).qualclass = 1000;
  768     }
  769     */
  770 }
  771 
  772 
  773 bool AdFront3 :: Inside (const Point<3> & p) const
  774 {
  775   int cnt;
  776   Vec3d n, v1, v2;
  777   DenseMatrix a(3), ainv(3);
  778   Vector b(3), u(3);
  779 
  780   // random numbers:
  781   n.X() = 0.123871;
  782   n.Y() = 0.15432;
  783   n.Z() = -0.43989;
  784 
  785   cnt = 0;
  786   for (int i = 1; i <= faces.Size(); i++)
  787     if (faces.Get(i).Valid())
  788       {
  789     const Point<3> & p1 = points[faces.Get(i).Face().PNum(1)].P();
  790     const Point<3> & p2 = points[faces.Get(i).Face().PNum(2)].P();
  791     const Point<3> & p3 = points[faces.Get(i).Face().PNum(3)].P();
  792 
  793     v1 = p2 - p1;
  794     v2 = p3 - p1;
  795 
  796     a(0, 0) = v1.X();
  797     a(1, 0) = v1.Y();
  798     a(2, 0) = v1.Z();
  799     a(0, 1) = v2.X();
  800     a(1, 1) = v2.Y();
  801     a(2, 1) = v2.Z();
  802     a(0, 2) = -n.X();
  803     a(1, 2) = -n.Y();
  804     a(2, 2) = -n.Z();
  805 
  806     b(0) = p(0) - p1(0);
  807     b(1) = p(1) - p1(1);
  808     b(2) = p(2) - p1(2);
  809 
  810     CalcInverse (a, ainv);
  811     ainv.Mult (b, u);
  812 
  813     if (u(0) >= 0 && u(1) >= 0 && u(0)+u(1) <= 1 &&
  814         u(2) > 0)
  815       {
  816         cnt++;
  817       }
  818       }
  819 
  820   return ((cnt % 2) != 0);
  821 }
  822 
  823 
  824 
  825 
  826 
  827 int AdFront3 :: SameSide (const Point<3> & lp1, const Point<3> & lp2,
  828               const Array<int> * testfaces) const
  829 {
  830   const Point<3> *line[2];
  831   line[0] = &lp1;
  832   line[1] = &lp2;
  833 
  834 
  835   Point3d pmin(lp1);
  836   Point3d pmax(lp1);
  837   pmin.SetToMin (lp2);
  838   pmax.SetToMax (lp2);
  839   
  840   ArrayMem<int, 100> aprif;
  841   aprif.SetSize(0);
  842   
  843   if (!testfaces)
  844     facetree->GetIntersecting (pmin, pmax, aprif);
  845   else
  846     for (int i = 1; i <= testfaces->Size(); i++)
  847       aprif.Append (testfaces->Get(i));
  848 
  849   int cnt = 0;
  850   for (int ii = 1; ii <= aprif.Size(); ii++)
  851     {
  852       int i = aprif.Get(ii);
  853       
  854       if (faces.Get(i).Valid())
  855     {
  856       const Point<3> *tri[3];
  857       tri[0] = &points[faces.Get(i).Face().PNum(1)].P();
  858       tri[1] = &points[faces.Get(i).Face().PNum(2)].P();
  859       tri[2] = &points[faces.Get(i).Face().PNum(3)].P();
  860           
  861       if (IntersectTriangleLine (&tri[0], &line[0]))
  862         cnt++;
  863     }
  864     }
  865 
  866   return ((cnt+1) % 2);
  867 }
  868 }