meshGFaceBDS.cpp

// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.

#include <stdlib.h>
#include "GmshMessage.h"
#include "robustPredicates.h"
#include "meshGFace.h"
#include "meshGFaceOptimize.h"
#include "BackgroundMesh.h"
#include "GVertex.h"
#include "GEdge.h"
#include "GFace.h"
#include "discreteFace.h"
#include "MVertex.h"
#include "MElement.h"
#include "Context.h"
#include "GPoint.h"
#include "GModel.h"
#include "Numeric.h"
#include "BDS.h"
#include "qualityMeasures.h"
#include "Field.h"
#include "OS.h"

static void getDegeneratedVertices(
  BDS_Mesh &m, std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap,
  std::set<MVertex *, MVertexLessThanNum> &degenerated,
  std::vector<BDS_Edge *> &degenerated_edges)
{
  degenerated.clear();

  std::vector<BDS_Edge *>::iterator it = m.edges.begin();

  while(it != m.edges.end()) {
    BDS_Edge *e = *it;
    if(!e->deleted && e->numfaces() == 1) {
      std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp1 =
        recoverMap->find(e->p1);
      std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp2 =
        recoverMap->find(e->p2);
      if(itp1 != recoverMap->end() && itp2 != recoverMap->end() &&
         itp1->second == itp2->second) {
        degenerated.insert(itp1->second);
        //	itp1->first->_degeneratedTo = itp2->first;
        //	itp2->first->_degeneratedTo = itp1->first;
        degenerated_edges.push_back(e);
      }
    }
    ++it;
  }
}

double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
{
  const double dx = p1->X - p2->X;
  const double dy = p1->Y - p2->Y;
  const double dz = p1->Z - p2->Z;
  return std::sqrt(dx * dx + dy * dy + dz * dz);
}

inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
  GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u), 0.5 * (p1->v + p2->v)));

  if(!GP.succeeded()) return computeEdgeLinearLength(p1, p2);

  const double dx1 = p1->X - GP.x();
  const double dy1 = p1->Y - GP.y();