netgen/libsrc/meshing/improve2.hpp

#ifndef FILE_IMPROVE2
#define FILE_IMPROVE2

template<typename TINDEX>
void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges )
{
  static Timer tbuild_edges("Build edges"); RegionTimer reg(tbuild_edges);

  static constexpr int tetedges[6][2] =
    { { 0, 1 }, { 0, 2 }, { 0, 3 },
        { 1, 2 }, { 1, 3 }, { 2, 3 } };

  int ntasks = 2*ngcore::TaskManager::GetMaxThreads();
  Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);

  ParallelFor(IntRange(ntasks), [&] (int ti)
    {
      auto myrange = mesh.Points().Range().Split(ti, ntasks);
      ArrayMem<std::tuple<PointIndex,PointIndex>, 100> local_edges;
      for (auto pi : myrange)
      {
        local_edges.SetSize(0);

        for(auto ei : elementsonnode[pi])
        {
            const auto & elem = mesh[ei];
            if (elem.IsDeleted()) continue;

            for (int j = 0; j < 6; j++)
            {
                PointIndex pi0 = elem[tetedges[j][0]];
                PointIndex pi1 = elem[tetedges[j][1]];
                if (pi1 < pi0) Swap(pi0, pi1);
                if(pi0==pi)
                    local_edges.Append(std::make_tuple(pi0, pi1));
            }
        }
        QuickSort(local_edges);

        auto edge_prev = std::make_tuple<PointIndex, PointIndex>(-1,-1);

        for(auto edge : local_edges)
            if(edge != edge_prev)
            {
                task_edges[ti].Append(edge);
                edge_prev = edge;
            }
      }
    }, ntasks);

  int num_edges = 0;
  for (auto & edg : task_edges)
      num_edges += edg.Size();
  edges.SetAllocSize(num_edges);
  for (auto & edg : task_edges)
      edges.Append(edg);
}


class Neighbour
{
  int nr[3];
  int orient[3];

public:
  Neighbour () { ; }

  void SetNr (int side, int anr) { nr[side] = anr; }
  int GetNr (int side) { return nr[side]; }

  void SetOrientation (int side, int aorient) { orient[side] = aorient; }
  int GetOrientation (int side) { return orient[side]; }
};

///
class MeshOptimize2d
{
  int faceindex;
  int improveedges;
  double metricweight;
  int writestatus;

public:
  ///
  MeshOptimize2d ();
  virtual ~MeshOptimize2d() { ; }
  ///
  void ImproveMesh (Mesh & mesh2d, const MeshingParameters & mp);
  void ImproveMeshJacobian (Mesh & mesh2d, const MeshingParameters & mp);
  void ImproveVolumeMesh (Mesh & mesh);
  void ProjectBoundaryPoints(NgArray<int> & surfaceindex, 
			     const NgArray<Point<3>* > & from, NgArray<Point<3>* > & dest);

  bool EdgeSwapping (Mesh & mesh, const int usemetric, Array<Neighbour> &neighbors, Array<bool> &swapped,
    const SurfaceElementIndex t1, const int edge, const int t, Array<int,PointIndex> &pdef, const bool check_only=false);
  void EdgeSwapping (Mesh & mesh, int usemetric);
  void CombineImprove (Mesh & mesh);
  void SplitImprove (Mesh & mesh);

  void GenericImprove (Mesh & mesh);


  void SetFaceIndex (int fi) { faceindex = fi; }
  void SetImproveEdges (int ie) { improveedges = ie; }
  void SetMetricWeight (double mw) { metricweight = mw; }
  void SetWriteStatus (int ws) { writestatus = ws; }



  ///
  virtual void ProjectPoint (INDEX /* surfind */, Point<3> & /* p */) const { };

  /// project point, use gi as initial value, and compute new gi
  virtual int ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const 
  { ProjectPoint (surfind, p); return CalcPointGeomInfo (surfind, gi, p); }

  ///
  virtual void ProjectPoint2 (INDEX /* surfind */, INDEX /* surfind2 */, Point<3> & /* p */) const { };

  /// liefert zu einem 3d-Punkt die geominfo (Dreieck) und liefert 1, wenn erfolgreich, 
  /// 0, wenn nicht (Punkt ausserhalb von chart)
  virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point<3> & /*p3*/) const
    { gi.trignum = 1; return 1;};

  virtual int CalcPointGeomInfo(int /* surfind */, PointGeomInfo& gi, const Point<3> & p3) const
    { return CalcPointGeomInfo (gi, p3); }

  ///
  virtual void GetNormalVector(INDEX surfind, const Point<3>  & p, PointGeomInfo & gi, Vec<3> & n) const;
  virtual void GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const;

  void CheckMeshApproximation (Mesh & mesh);


  ///
  friend class Opti2SurfaceMinFunction;
  ///
  friend class Opti2EdgeMinFunction;
  ///
  friend double Opti2FunctionValueGrad (const Vector & x, Vector & grad);
  ///
  friend double Opti2EdgeFunctionValueGrad (const Vector & x, Vector & grad);



};


extern void CalcTriangleBadness (double x2, double x3, double y3, 
				 double metricweight,
				 double h, double & badness, 
				 double & g1x, double & g1y);




extern double CalcTriangleBadness (const Point<3> & p1, 
				   const Point<3> & p2, 
				   const Point<3> & p3,
				   double metricweight,
				   double h);

extern double CalcTriangleBadness (const Point<3> & p1, 
				   const Point<3> & p2, 
				   const Point<3> & p3,
				   const Vec<3> & n,
				   double metricweight,
				   double h);

#endif