C++ SVector3::norm方法代码示例

bool GEdge::XYZToU(const double X, const double Y, const double Z,
                   double &u, const double relax) const
{
  const int MaxIter = 25;
  const int NumInitGuess = 11;

  double err;//, err2;
  int iter;

  Range<double> uu = parBounds(0);
  double uMin = uu.low();
  double uMax = uu.high();

  SVector3 Q(X, Y, Z), P;

  double init[NumInitGuess];

  for (int i = 0; i < NumInitGuess; i++)
    init[i] = uMin + (uMax - uMin) / (NumInitGuess - 1) * i;

  for(int i = 0; i < NumInitGuess; i++){
    u = init[i];
    double uNew = u;
    //err2 = 1.0;
    iter = 1;

    SVector3 dPQ = P - Q;
    err = dPQ.norm();

    if (err < 1.e-8 * CTX::instance()->lc) return true;

    while(iter++ < MaxIter && err > 1e-8 * CTX::instance()->lc) {
      SVector3 der = firstDer(u);
      uNew = u - relax * dot(dPQ,der) / dot(der,der);
      uNew = std::min(uMax,std::max(uMin,uNew));
      P = position(uNew);

      dPQ = P - Q;
      err = dPQ.norm();
      //err2 = fabs(uNew - u);
      u = uNew;
    }

    if (err < 1e-8 * CTX::instance()->lc) return true;
  }

  if(relax > 1.e-2) {
    //    Msg::Info("point %g %g %g on edge %d : Relaxation factor = %g",
    //              Q.x(), Q.y(), Q.z(), 0.75 * relax);
    return XYZToU(Q.x(), Q.y(), Q.z(), u, 0.75 * relax);
  }

  //  Msg::Error("Could not converge reparametrisation of point (%e,%e,%e) on edge %d",
  //             Q.x(), Q.y(), Q.z(), tag());
  return false;
}

double distMaxStraight(MElement *el)
{
  const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
  const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
  int nV = lagrange->points.size1();
  int nV1 = lagrange1->points.size1();
  int dim = lagrange1->dimension;
  SPoint3 sxyz[256];
  for (int i = 0; i < nV1; ++i) {
    sxyz[i] = el->getVertex(i)->point();
  }
  for (int i = nV1; i < nV; ++i) {
    double f[256];
    lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1),
                 dim < 3 ? 0 : lagrange->points(i, 2), f);
    for (int j = 0; j < nV1; ++j)
      sxyz[i] += sxyz[j] * f[j];
  }

  double maxdx = 0.0;
  for (int iV = nV1; iV < nV; iV++) {
    SVector3 d = el->getVertex(iV)->point()-sxyz[iV];
    double dx = d.norm();
    if (dx > maxdx) maxdx = dx;
  }

  return maxdx;
}

double angle3Vertices(const MVertex *p1, const MVertex *p2, const MVertex *p3)
{
  SVector3 a(p1->x() - p2->x(), p1->y() - p2->y(), p1->z() - p2->z());
  SVector3 b(p3->x() - p2->x(), p3->y() - p2->y(), p3->z() - p2->z());
  SVector3 c = crossprod(a, b);
  double sinA = c.norm();
  double cosA = dot(a, b);
  return atan2 (sinA, cosA);
}

double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1,
                           double x2, double x3, double tau)
{
  // Create a new possible center in the area between x2 and x3, closer to x2
  double x4 = x2 + GOLDEN2 * (x3 - x2);

  // Evaluate termination criterion
  if (fabs(x3 - x1) < tau * (fabs(x2) + fabs(x4)))
    return (x3 + x1) / 2;

  const SVector3 dp4 = q - ge->position(x4);
  const SVector3 dp2 = q - ge->position(x2);

  const double d4 = dp4.norm();
  const double d2 = dp2.norm();

  if (d4 < d2)
    return goldenSectionSearch(ge, q, x2, x4, x3, tau);
  else
    return goldenSectionSearch(ge,q , x4, x2, x1, tau);
}

GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
{
  // printf("looking for closest point in curve %d to point %g %g\n",tag(),q.x(),q.y());

  const int nbSamples = 100;

  Range<double> interval = parBounds(0);

  double tMin = std::min(interval.high(), interval.low());
  double tMax = std::max(interval.high(), interval.low());

  double DMIN = 1.e22;
  double topt = tMin;
  const double DT = (tMax - tMin) / (nbSamples - 1.) ;
  for (int i = 0; i < nbSamples; i++){
    t = tMin + i * DT;
    const SVector3 dp = q - position(t);
    const double D = dp.norm();
    if (D < DMIN) {
      topt = t;
      DMIN = D;
    }
  }

  // printf("parameter %g as an initial guess (dist = %g)\n",topt,DMIN);

  if (topt == tMin)
    t = goldenSectionSearch (this, q, topt, topt + DT/2, topt + DT,  1.e-9);
  else if (topt == tMax)
    t = goldenSectionSearch (this, q, topt - DT, topt - DT/2 , topt, 1.e-9);
  else
    t = goldenSectionSearch (this, q, topt - DT, topt, topt + DT, 1.e-9);

  const SVector3 dp = q - position(t);
  // const double D = dp.norm();
  // printf("after golden section parameter %g  (dist = %g)\n",t,D);

  return point(t);
}

void meshGEdge::operator() (GEdge *ge)
{
#if defined(HAVE_ANN)
  FieldManager *fields = ge->model()->getFields();
  BoundaryLayerField *blf = 0;
  Field *bl_field = fields->get(fields->getBoundaryLayerField());
  blf = dynamic_cast<BoundaryLayerField*> (bl_field);
#else
  bool blf = false;
#endif

  ge->model()->setCurrentMeshEntity(ge);

  if(ge->geomType() == GEntity::DiscreteCurve) return;
  if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
  if(ge->meshAttributes.method == MESH_NONE) return;
  if(CTX::instance()->mesh.meshOnlyVisible && !ge->getVisibility()) return;

  // look if we are doing the STL triangulation
  std::vector<MVertex*> &mesh_vertices = ge->mesh_vertices ;
  std::vector<MLine*> &lines = ge->lines ;

  deMeshGEdge dem;
  dem(ge);

  if(MeshExtrudedCurve(ge)) return;

  if (ge->meshMaster() != ge){
    GEdge *gef = dynamic_cast<GEdge*> (ge->meshMaster());
    if (gef->meshStatistics.status == GEdge::PENDING) return;
    Msg::Info("Meshing curve %d (%s) as a copy of %d", ge->tag(),
              ge->getTypeString().c_str(), ge->meshMaster()->tag());
    copyMesh(gef, ge, ge->masterOrientation);
    ge->meshStatistics.status = GEdge::DONE;
    return;
  }

  Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());

  // compute bounds
  Range<double> bounds = ge->parBounds(0);
  double t_begin = bounds.low();
  double t_end = bounds.high();

  // first compute the length of the curve by integrating one
  double length;
  std::vector<IntPoint> Points;
  if(ge->geomType() == GEntity::Line &&
      ge->getBeginVertex() == ge->getEndVertex() &&
      //do not consider closed lines as degenerated
      (ge->position(0.5) - ge->getBeginVertex()->xyz()).norm() < CTX::instance()->geom.tolerance)
    length = 0.; // special case t avoid infinite loop in integration
  else
    length = Integration(ge, t_begin, t_end, F_One, Points, 1.e-8 * CTX::instance()->lc);
  ge->setLength(length);
  Points.clear();

  if(length < CTX::instance()->mesh.toleranceEdgeLength){
    ge->setTooSmall(true);
  }

  // Integrate detJ/lc du
  double a;
  int N;
  if(length == 0. && CTX::instance()->mesh.toleranceEdgeLength == 0.){
    Msg::Warning("Curve %d has a zero length", ge->tag());
    a = 0.;
    N = 1;
  }
  else if(ge->degenerate(0)){
    a = 0.;
    N = 1;
  }
  else if(ge->meshAttributes.method == MESH_TRANSFINITE){
    a = Integration(ge, t_begin, t_end, F_Transfinite, Points,
                    CTX::instance()->mesh.lcIntegrationPrecision);
    N = ge->meshAttributes.nbPointsTransfinite;
    if(CTX::instance()->mesh.flexibleTransfinite && CTX::instance()->mesh.lcFactor)
      N /= CTX::instance()->mesh.lcFactor;
  }
  else{
    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf){
      a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
                      CTX::instance()->mesh.lcIntegrationPrecision);
    }
    else{
       a = Integration(ge, t_begin, t_end, F_Lc, Points,
                      CTX::instance()->mesh.lcIntegrationPrecision);
    }

    // we should maybe provide an option to disable the smoothing
    for (unsigned int i = 0; i < Points.size(); i++){
      IntPoint &pt = Points[i];
      SVector3 der = ge->firstDer(pt.t);
      pt.xp = der.norm();
    }
    a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
    N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.99));
  }

//.........这里部分代码省略.........