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

double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,double h1,SVector3 direction){
  double x,y,z;
  double average;
  double h2;
  double coeffA,coeffB;

  x = point.x() + h1*direction.x();
  y = point.y() + h1*direction.y();
  z = point.z() + h1*direction.z();
  
  if(inside_domain(octree,x,y,z)){
    h2 = get_size(x,y,z);
  }
  else h2 = h1;

  coeffA = 1.0;
  coeffB = 0.16;
  
  if(h2>h1){
    average = coeffA*h1 + (1.0-coeffA)*h2;
  }
  else{
    average = coeffB*h1 + (1.0-coeffB)*h2;
  }
	
  return average;
}

static void _relocateVertexOfPyramid(MVertex *ver,
                                     const std::vector<MElement *> &lt,
                                     double relax)
{
  if(ver->onWhat()->dim() != 3) return;
  double x = 0.0, y = 0.0, z = 0.0;
  int N = 0;
  MElement *pyramid = NULL;

  for(std::size_t i = 0; i < lt.size(); i++) {
    double XCG = 0.0, YCG = 0.0, ZCG = 0.0;
    if(lt[i]->getNumVertices() == 5)
      pyramid = lt[i];
    else {
      for(std::size_t j = 0; j < lt[i]->getNumVertices(); j++) {
        XCG += lt[i]->getVertex(j)->x();
        YCG += lt[i]->getVertex(j)->y();
        ZCG += lt[i]->getVertex(j)->z();
      }
      x += XCG;
      y += YCG;
      z += ZCG;
      N += lt[i]->getNumVertices();
    }
  }
  x /= N;
  y /= N;
  z /= N;

  if(pyramid) {
    MFace q = pyramid->getFace(4);
    double A = q.approximateArea();
    SVector3 n = q.normal();
    n.normalize();
    SPoint3 c = q.barycenter();
    SVector3 d(x - c.x(), y - c.y(), z - c.z());
    if(dot(n, d) < 0) n = n * (-1.0);
    double H = .5 * sqrt(fabs(A));
    double XOPT = c.x() + relax * H * n.x();
    double YOPT = c.y() + relax * H * n.y();
    double ZOPT = c.z() + relax * H * n.z();
    double FULL_MOVE_OBJ =
      objective_function(1.0, ver, XOPT, YOPT, ZOPT, lt, true);
    // printf("relax %g obj %g\n",relax,FULL_MOVE_OBJ);
    if(FULL_MOVE_OBJ > 0.1) {
      ver->x() = XOPT;
      ver->y() = YOPT;
      ver->z() = ZOPT;
      return;
    }
  }
}

void Centerline::buildKdTree()
{
  FILE * f = Fopen("myPOINTS.pos","w");
  fprintf(f, "View \"\"{\n");

  int nbPL = 3;  //10 points per line
  //int nbNodes  = (lines.size()+1) + (nbPL*lines.size());
  int nbNodes  = (colorp.size()) + (nbPL*lines.size());

  ANNpointArray nodes = annAllocPts(nbNodes, 3);
  int ind = 0;
  std::map<MVertex*, int>::iterator itp = colorp.begin();
  while (itp != colorp.end()){
     MVertex *v = itp->first;
     nodes[ind][0] = v->x();
     nodes[ind][1] = v->y();
     nodes[ind][2] = v->z();
     itp++; ind++;
  }
  for(unsigned int k = 0; k < lines.size(); ++k){
   MVertex *v0 = lines[k]->getVertex(0);
   MVertex *v1 = lines[k]->getVertex(1);
   SVector3 P0(v0->x(),v0->y(), v0->z());
   SVector3 P1(v1->x(),v1->y(), v1->z());
   for (int j= 1; j < nbPL+1; j++){
     double inc = (double)j/(double)(nbPL+1);
     SVector3 Pj = P0+inc*(P1-P0);
     nodes[ind][0] = Pj.x();
     nodes[ind][1] = Pj.y();
     nodes[ind][2] = Pj.z();
     ind++;
   }
 }

 kdtree = new ANNkd_tree(nodes, nbNodes, 3);

 for(int i = 0; i < nbNodes; ++i){
   fprintf(f, "SP(%g,%g,%g){%g};\n",
	   nodes[i][0], nodes[i][1],nodes[i][2],1.0);
 }
 fprintf(f,"};\n");
 fclose(f);
}

void highOrderTools::computeMetricInfo(GFace *gf,
                                       MElement *e,
                                       fullMatrix<double> &J,
                                       fullMatrix<double> &JT,
                                       fullVector<double> &D)
{
  int nbNodes = e->getNumVertices();
  //  printf("ELEMENT --\n");
  for (int j = 0; j < nbNodes; j++){
    SPoint2 param;
    reparamMeshVertexOnFace(e->getVertex(j), gf, param);
    //    printf("%g %g vs %g %g %g\n",param.x(),param.y(),
    //	   e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());

    Pair<SVector3,SVector3> der = gf->firstDer(param);

    int XJ = j;
    int YJ = j + nbNodes;
    int ZJ = j + 2 * nbNodes;
    int UJ = j;
    int VJ = j + nbNodes;
    J(XJ,UJ) = der.first().x();
    J(YJ,UJ) = der.first().y();
    J(ZJ,UJ) = der.first().z();
    J(XJ,VJ) = der.second().x();
    J(YJ,VJ) = der.second().y();
    J(ZJ,VJ) = der.second().z();

    JT(UJ,XJ) = der.first().x();
    JT(UJ,YJ) = der.first().y();
    JT(UJ,ZJ) = der.first().z();
    JT(VJ,XJ) = der.second().x();
    JT(VJ,YJ) = der.second().y();
    JT(VJ,ZJ) = der.second().z();

    SVector3 ss = getSSL(e->getVertex(j));
    GPoint gp = gf->point(param);
    D(XJ) = (gp.x() - ss.x());
    D(YJ) = (gp.y() - ss.y());
    D(ZJ) = (gp.z() - ss.z());
  }
}

void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int order) const
{
  if(!_orig)
    return;

  if (_orig->getDim()==getDim())
    return _orig->getGradShapeFunctions(u, v, w, s, order);

  int nsf = _orig->getNumShapeFunctions();
  double gradsuvw[1256][3];
  _orig->getGradShapeFunctions(u, v, w, gradsuvw, order);

  double jac[3][3];
  double invjac[3][3];
  _orig->getJacobian(u, v, w, jac);
  inv3x3(jac, invjac);
  MEdge edge = getBaseElement()->getEdge(0);
  SVector3 tang = edge.tangent();

  double gradxyz[3];
  double projgradxyz[3];
  for (int i=0; i<nsf; ++i)
  {
    // (i) get the cartesian coordinates of the gradient
    gradxyz[0] = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
    gradxyz[1] = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
    gradxyz[2] = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];

    // (ii) projection of the gradient on edges in the cartesian space
    SVector3 grad(&gradxyz[0]);
    double prodscal = dot(tang,grad);
    projgradxyz[0] = prodscal * tang.x();
    projgradxyz[1] = prodscal * tang.y();
    projgradxyz[2] = prodscal * tang.z();

    // (iii) get the parametric coordinates of the projection in the parametric space of the parent element
    s[i][0] = jac[0][0] * projgradxyz[0] + jac[0][1] * projgradxyz[1] + jac[0][2] * projgradxyz[2];
    s[i][1] = jac[1][0] * projgradxyz[0] + jac[1][1] * projgradxyz[1] + jac[1][2] * projgradxyz[2];
    s[i][2] = jac[2][0] * projgradxyz[0] + jac[2][1] * projgradxyz[1] + jac[2][2] * projgradxyz[2];
  }
}

int main(int argc,char *argv[])
{
  if(argc < 6){
    printf("Usage: %s file lx ly lz rmax [levels=1] [refcs=1]\n", argv[0]);
    printf("where\n");
    printf("  'file' contains a CAD model\n");
    printf("  'lx', 'ly' and 'lz' are the sizes of the elements along the"
           " x-, y- and z-axis at the coarsest level\n");
    printf("  'rmax' is the radius of the largest sphere that can be inscribed"
           " in the structure\n");
    printf("  'levels' sets the number of levels in the grid\n");
    printf("  'refcs' selects if curved surfaces should be refined\n");
    return -1;
  }

  GmshInitialize();
  GmshSetOption("General", "Terminal", 1.);
  GmshMergeFile(argv[1]);
  double lx = atof(argv[2]), ly = atof(argv[3]), lz = atof(argv[4]);
  double rmax = atof(argv[5]);
  int levels = (argc > 6) ? atof(argv[6]) : 1;
  int refineCurvedSurfaces = (argc > 7) ? atof(argv[7]) : 1;

  // minimum distance between points in the cloud at the coarsest
  // level
  double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));

  // radius of the "tube" created around parts to refine at the
  // coarsest level
  double rtube = std::max(lx, std::max(ly, lz)) * 2.;

  GModel *gm = GModel::current();

  std::vector<SPoint3> points;
  Msg::Info("Filling coarse point cloud on surfaces");
  for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
    (*fit)->fillPointCloud(sampling, &points);
  Msg::Info("  %d points in the surface cloud", (int)points.size());

  std::vector<SPoint3> refinePoints;
  if(levels > 1){
    double s = sampling / pow(2., levels - 1);
    Msg::Info("Filling refined point cloud on curves and curved surfaces");
    for (GModel::eiter eit = gm->firstEdge(); eit != gm->lastEdge(); eit++)
      fillPointCloud(*eit, s, refinePoints);

    // FIXME: refine this by computing e.g. "mean" curvature
    if(refineCurvedSurfaces){
      for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
        if((*fit)->geomType() != GEntity::Plane)
          (*fit)->fillPointCloud(2 * s, &refinePoints);
    }
    Msg::Info("  %d points in the refined cloud", (int)refinePoints.size());
  }

  SBoundingBox3d bb;
  for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
  for(unsigned int i = 0; i < refinePoints.size(); i++) bb += refinePoints[i];
  bb.scale(1.21, 1.21, 1.21);
  SVector3 range = bb.max() - bb.min();
  int NX = range.x() / lx;
  int NY = range.y() / ly;
  int NZ = range.z() / lz;
  if(NX < 2) NX = 2;
  if(NY < 2) NY = 2;
  if(NZ < 2) NZ = 2;

  Msg::Info("  bounding box min: %g %g %g -- max: %g %g %g",
            bb.min().x(), bb.min().y(), bb.min().z(),
            bb.max().x(), bb.max().y(), bb.max().z());
  Msg::Info("  Nx=%d Ny=%d Nz=%d", NX, NY, NZ);

  cartesianBox<double> box(bb.min().x(), bb.min().y(), bb.min().z(),
                           SVector3(range.x(), 0, 0),
                           SVector3(0, range.y(), 0),
                           SVector3(0, 0, range.z()),
                           NX, NY, NZ, levels);

  Msg::Info("Inserting active cells in the cartesian grid");
  Msg::Info("  level %d", box.getLevel());
  for (unsigned int i = 0; i < points.size(); i++)
    insertActiveCells(points[i].x(), points[i].y(), points[i].z(), rmax, box);

  cartesianBox<double> *parent = &box, *child;
  while((child = parent->getChildBox())){
    Msg::Info("  level %d", child->getLevel());
    for(unsigned int i = 0; i < refinePoints.size(); i++)
      insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
                        rtube / pow(2., (levels - child->getLevel())), *child);
    parent = child;
  }

  // remove child cells that do not entirely fill parent cell or for
  // which there is no parent neighbor; then remove parent cells that
  // have children
  Msg::Info("Removing cells to match X-FEM mesh topology constraints");
  removeBadChildCells(&box);
  removeParentCellsWithChildren(&box);

  // we generate duplicate nodes at this point so we can easily access
//.........这里部分代码省略.........

bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad)
{
  double a = NORM.x();
  double b = NORM.y();
  double c = NORM.z();
  double d = -a * PT.x() - b * PT.y() - c * PT.z();

  int maxStep = 20;
  const double EPS = 0.007;

  //std::set<MEdge,Less_Edge> allEdges;
  std::vector<MEdge> allEdges;
  for(unsigned int i = 0; i < triangles.size(); i++){
    for ( unsigned int j= 0; j <  3; j++){
      allEdges.push_back(triangles[i]->getEdge(j));
      //allEdges.insert(triangles[i]->getEdge(j));
    }
  }
  std::unique(allEdges.begin(), allEdges.end());

  bool closedCut = false;
  int step = 0;
  while (!closedCut && step < maxStep){
    double rad = 1.1*maxRad+0.05*step*maxRad;
    std::map<MEdge,MVertex*,Less_Edge> cutEdges;
    std::vector<MVertex*> cutVertices;
    std::vector<MTriangle*> newTris;
    std::set<MEdge,Less_Edge> newCut;
    cutEdges.clear();
    cutVertices.clear();
    newTris.clear();
    newCut.clear();
    // for (std::set<MEdge,Less_Edge>::iterator it = allEdges.begin();
    // 	 it != allEdges.end() ; ++it){
    // MEdge me = *it;
    for (unsigned int j = 0; j < allEdges.size(); j++){
      MEdge me = allEdges[j];
      SVector3 P1(me.getVertex(0)->x(),me.getVertex(0)->y(), me.getVertex(0)->z());
      SVector3 P2(me.getVertex(1)->x(),me.getVertex(1)->y(), me.getVertex(1)->z());
      double V1 = a * P1.x() + b * P1.y() + c * P1.z() + d;
      double V2 = a * P2.x() + b * P2.y() + c * P2.z() + d;
      bool inters = (V1*V2<=0.0) ? true: false;
      bool inDisk = ((norm(P1-PT) < rad ) || (norm(P2-PT) < rad)) ? true : false;
      double rdist = -V1/(V2-V1);
      if (inters && rdist > EPS && rdist < 1.-EPS){
	SVector3 PZ = P1+rdist*(P2-P1);
	if (inDisk){
          MVertex *newv = new MVertex (PZ.x(), PZ.y(), PZ.z());
          cutEdges.insert(std::make_pair(me,newv));
        }
      }
      else if (inters && rdist <= EPS && inDisk )
	cutVertices.push_back(me.getVertex(0));
      else if (inters && rdist >= 1.-EPS && inDisk)
	cutVertices.push_back(me.getVertex(1));
    }
    for(unsigned int i = 0; i < triangles.size(); i++){
      cutTriangle(triangles[i], cutEdges,cutVertices, newTris, newCut);
    }
    if (isClosed(newCut)) {
      triangles.clear();
      triangles = newTris;
      theCut.insert(newCut.begin(),newCut.end());
      break;
    }
    else {
      step++;
      //if (step == maxStep) {printf("no closed cut %d \n", (int)newCut.size()); };
      // // triangles = newTris;
      // // theCut.insert(newCut.begin(),newCut.end());
      // char name[256];
      // sprintf(name, "myCUT-%d.pos", step);
      // FILE * f2 = Fopen(name,"w");
      // fprintf(f2, "View \"\"{\n");
      // std::set<MEdge,Less_Edge>::iterator itp =  newCut.begin();
      // while (itp != newCut.end()){
      // 	MEdge l = *itp;
      // 	fprintf(f2, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
      // 		  l.getVertex(0)->x(), l.getVertex(0)->y(), l.getVertex(0)->z(),
      // 		  l.getVertex(1)->x(), l.getVertex(1)->y(), l.getVertex(1)->z(),
      // 		  1.0,1.0);
      // 	  itp++;
      // 	}
      // 	fprintf(f2,"};\n");
      // 	fclose(f2);
    }
  }


  if (step < maxStep){
    //printf("cutByDisk OK step =%d  \n", step);
    return true;
  }
  else {
    //printf("cutByDisk not succeeded \n");
    return false;
  }

}

void highOrderTools::applySmoothingTo(std::vector<MElement*> &all, GFace *gf)
{
#ifdef HAVE_TAUCS
  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
  linearSystemPETSc<double> *lsys = new  linearSystemPETSc<double>;
#endif
  // compute the straight sided positions of high order nodes that are
  // on the edges of the face in the UV plane
  dofManager<double> myAssembler(lsys);
  elasticityTerm El(0, 1.0, CTX::instance()->mesh.hoPoissonRatio, _tag);
  std::vector<MElement*> layer, v;
  double minD;
  getDistordedElements(all, CTX::instance()->mesh.hoThresholdMin, v, minD);
  int numBad = v.size();
  const int nbLayers = CTX::instance()->mesh.hoNLayers;
  for (int i = 0; i < nbLayers; i++){
    addOneLayer(all, v, layer);
    v.insert(v.end(), layer.begin(), layer.end());
  }

  if (!v.size()) return;
  Msg::Info("Smoothing high order mesh : model face %d (%d elements considered in "
            "the elastic analogy, worst mapping %12.5E, %3d bad elements)", gf->tag(),
            v.size(),minD,numBad);

  addOneLayer(all, v, layer);
  std::set<MVertex*>::iterator it;
  std::set<MVertex*> verticesToMove;

  // on the last layer, fix displacement to 0
  for (unsigned int i = 0; i < layer.size(); i++){
    for (int j = 0; j < layer[i]->getNumVertices(); j++){
      MVertex *vert = layer[i]->getVertex(j);
      myAssembler.fixVertex(vert, 0, _tag, 0);
      myAssembler.fixVertex(vert, 1, _tag, 0);
    }
  }

  // fix all vertices that cannot move
  for (unsigned int i = 0; i < v.size(); i++){
    moveToStraightSidedLocation(v[i]);
    for (int j = 0; j < v[i]->getNumVertices(); j++){
      MVertex *vert = v[i]->getVertex(j);
      if (vert->onWhat()->dim() < 2){
	double du = 0, dv = 0;
        myAssembler.fixVertex(vert, 0, _tag, du);
        myAssembler.fixVertex(vert, 1, _tag, dv);
      }
    }
  }

  // number the other DOFs
  for (unsigned int i = 0; i < v.size(); i++){
    for (int j = 0; j < v[i]->getNumVertices(); j++){
      MVertex *vert = v[i]->getVertex(j);
      myAssembler.numberVertex(vert, 0, _tag);
      myAssembler.numberVertex(vert, 1, _tag);
      verticesToMove.insert(vert);
    }
  }

  double dx0 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
  double dx = dx0;
  Msg::Debug(" dx0 = %12.5E", dx0);
  int iter = 0;
  while(0){
    double dx2 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
    Msg::Debug(" dx2  = %12.5E", dx2);
    if (fabs(dx2 - dx) < 1.e-4 * dx0)break;
    if (iter++ > 2)break;
    dx = dx2;
  }

  for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
    SPoint2 param;
    if ((*it)->onWhat()->dim() == 2){
      reparamMeshVertexOnFace(*it, gf, param);
      GPoint gp = gf->point(param);
      (*it)->x() = gp.x();
      (*it)->y() = gp.y();
      (*it)->z() = gp.z();
      _targetLocation[*it] = SVector3(gp.x(), gp.y(), gp.z());
    }
    else{
      SVector3 p =  getTL(*it);
      (*it)->x() = p.x();
      (*it)->y() = p.y();
      (*it)->z() = p.z();
    }
  }
  delete lsys;
}