BackgroundMesh.cpp

// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.

#include "GmshMessage.h"
#include "BackgroundMesh.h"
#include "Numeric.h"
#include "Context.h"
#include "GVertex.h"
#include "GEdge.h"
#include "GEdgeCompound.h"
#include "GFace.h"
#include "GFaceCompound.h"
#include "GModel.h"
#include "Field.h"
#include "MElement.h"
#include "MElementOctree.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MVertex.h"

#if defined(HAVE_SOLVER)
#include "dofManager.h"
#include "laplaceTerm.h"
#include "linearSystemGMM.h"
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
#include "linearSystemPETSc.h"
#endif

// computes the characteristic length of the mesh at a vertex in order
// to have the geometry captured with accuracy. A parameter called
// CTX::instance()->mesh.minCircPoints tells the minimum number of points per
// radius of curvature

#if defined(HAVE_ANN)
static int _NBANN = 2;
#endif

SMetric3 buildMetricTangentToCurve(SVector3 &t, double l_t, double l_n)
{
  if (l_t == 0.0) return SMetric3(1.e-22);
  SVector3 a;
  if (fabs(t(0)) <= fabs(t(1)) && fabs(t(0)) <= fabs(t(2))){
    a = SVector3(1,0,0);
  }
  else if (fabs(t(1)) <= fabs(t(0)) && fabs(t(1)) <= fabs(t(2))){
    a = SVector3(0,1,0);
  }
  else{
    a = SVector3(0,0,1);
  }
  SVector3 b = crossprod (t,a);
  SVector3 c = crossprod (b,t);
  b.normalize();
  c.normalize();
  t.normalize();
  SMetric3 Metric (1./(l_t*l_t),1./(l_n*l_n),1./(l_n*l_n),t,b,c);
  //  printf("bmttc %g %g %g %g %g\n",l_t,l_n,Metric(0,0),Metric(0,1),Metric(1,1));
  return Metric;
}

SMetric3 buildMetricTangentToSurface(SVector3 &t1, SVector3 &t2,
                                     double l_t1, double l_t2, double l_n)
{
  t1.normalize();
  t2.normalize();
  SVector3 n = crossprod (t1,t2);