C++ GEdge::point方法代码示例

bool OptHOM::addBndObjGrad(double factor, double &Obj, alglib::real_1d_array &gradObj)
{
  // set the mesh to its present position
  std::vector<SPoint3> xyz,uvw;
  mesh.getGEntityPositions(xyz,uvw);
  mesh.updateGEntityPositions();

  //could be better (e.g. store the model in the Mesh:: datastrucure)

  GModel *gm = GModel::current();

  // for all model edges, compute the error between the geometry and the mesh

  maxDistCAD = 0.0;
  double distCAD = 0.0;

  for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
    // do not do straight lines
    if ((*it)->geomType() == GEntity::Line)continue;
    // look at all mesh lines

    std::vector<bool> doWeCompute((*it)->lines.size());
    for (unsigned int i=0;i<(*it)->lines.size(); i++){
      doWeCompute[i] = false;
      for (unsigned int j=0;j<(*it)->lines[i]->getNumVertices(); j++){
        int index = mesh.getFreeVertexStartIndex((*it)->lines[i]->getVertex(j));
        if (index >=0){
          doWeCompute[i] = true;
          continue;
        }
      }
    }

    std::vector<double> dist((*it)->lines.size());
    for (unsigned int i=0;i<(*it)->lines.size(); i++){
      if (doWeCompute[i]){
        // compute the distance from the geometry to the mesh
        dist[i] = MLineGEdgeDistance ( (*it)->lines[i] , *it );
        maxDistCAD = std::max(maxDistCAD,dist[i]);
        distCAD += dist [i] * factor;
      }
    }
    // be clever to compute the derivative : iterate on all
    // Distance = \sum_{lines} Distance (line, GEdge)
    // For a high order vertex, we compute the derivative only by
    // recomputing the distance to one only line
    const double eps = 1.e-6;
    for (unsigned int i=0;i<(*it)->lines.size(); i++){
      if (doWeCompute[i]){
        for (int j=2 ; j<(*it)->lines[i]->getNumVertices()  ; j++){
          MVertex *v = (*it)->lines[i]->getVertex(j);
          int index = mesh.getFreeVertexStartIndex(v);
          //	printf("%d %d (%d %d)\n",v->getNum(),index,v->onWhat()->tag(),v->onWhat()->dim());
          if (index >= 0){
            double t;
            v->getParameter(0,t);
            SPoint3 pp (v->x(),v->y(),v->z());
            GPoint gp = (*it)->point(t+eps);
            v->setParameter(0,t+eps);
            v->setXYZ(gp.x(),gp.y(),gp.z());
            double dist2 = MLineGEdgeDistance ( (*it)->lines[i] , *it );
            double deriv = (dist2 - dist[i])/eps;
            v->setXYZ(pp.x(),pp.y(),pp.z());
            v->setParameter(0,t);
            //	  printf("%g %g %g\n",dist[i],dist2, MLineGEdgeDistance ( (*it)->lines[i] , *it ));
            // get the index of the vertex
            gradObj[index] += deriv * factor;
          }
        }
      }
      //    printf("done\n");
      // For a low order vertex classified on the GEdge, we recompute
    // two distances for the two MLines connected to the vertex
      for (unsigned int i=0;i<(*it)->lines.size()-1; i++){
        MVertex *v =  (*it)->lines[i]->getVertex(1);
        int index = mesh.getFreeVertexStartIndex(v);
        if (index >= 0){
          double t;
          v->getParameter(0,t);
          SPoint3 pp (v->x(),v->y(),v->z());
          GPoint gp = (*it)->point(t+eps);
          v->setParameter(0,t+eps);
          v->setXYZ(gp.x(),gp.y(),gp.z());
          MLine *l1 = (*it)->lines[i];
          MLine *l2 = (*it)->lines[i+1];
          //	printf("%d %d -- %d %d\n",l1->getVertex(0)->getNum(),l1->getVertex(1)->getNum(),l2->getVertex(0)->getNum(),l2->getVertex(1)->getNum());
          double deriv =
            (MLineGEdgeDistance ( l1 , *it ) - dist[i])  /eps +
            (MLineGEdgeDistance ( l2 , *it ) - dist[i+1])/eps;
          v->setXYZ(pp.x(),pp.y(),pp.z());
          v->setParameter(0,t);
          gradObj[index] += deriv * factor;
        }
      }
    }
  }
  //  printf("computing distance : 1D part %12.5E\n",distCAD);

  // now the 3D part !

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

void Mesh::approximationErrorAndGradients(int iEl, double &f, std::vector<double> &gradF, double eps,
					  simpleFunction<double> &fct)
{
  std::vector<SPoint3> _xyz_temp;
  for (int iV = 0; iV < nVert(); iV++){
    _xyz_temp.push_back(SPoint3( _vert[iV]->x(), _vert[iV]->y(), _vert[iV]->z()));
    _vert[iV]->setXYZ(_xyz[iV].x(),_xyz[iV].y(),_xyz[iV].z());
  }

  MElement *element = _el[iEl];

  f = approximationError (fct, element);
  // FIME
  //  if (iEl < 1)printf("approx error elem %d = %g\n",iEl,f);
  int currentId = 0;
  // compute the size of the gradient
  // depends on how many dofs exist per vertex (0,1,2 or 3)
  for (size_t i = 0; i < element->getNumVertices(); ++i) {
    if (_el2FV[iEl][i] >= 0) {// some free coordinates
      currentId += _nPCFV[_el2FV[iEl][i]];
    }
  }
  gradF.clear();
  gradF.resize(currentId, 0.);
  currentId = 0;
  for (size_t i = 0; i < element->getNumVertices(); ++i) {
    if (_el2FV[iEl][i] >= 0) {// some free coordinates
      MVertex *v =  element->getVertex(i);
      // vertex classified on a model edge
      if (_nPCFV[_el2FV[iEl][i]] == 1){
	double t = _uvw[_el2FV[iEl][i]].x();
	GEdge *ge = (GEdge*)v->onWhat();
	SPoint3 p (v->x(),v->y(),v->z());
	GPoint d = ge->point(t+eps);
	v->setXYZ(d.x(),d.y(),d.z());
	double f_d = approximationError (fct, element);
	gradF[currentId++] = (f_d-f)/eps;
	if (iEl < 1)printf("df = %g\n",(f_d-f)/eps);
	v->setXYZ(p.x(),p.y(),p.z());
      }
      else if (_nPCFV[_el2FV[iEl][i]] == 2){
	double uu = _uvw[_el2FV[iEl][i]].x();
	double vv = _uvw[_el2FV[iEl][i]].y();
	GFace *gf = (GFace*)v->onWhat();
	SPoint3 p (v->x(),v->y(),v->z());
	GPoint  d = gf->point(uu+eps,vv);
	v->setXYZ(d.x(),d.y(),d.z());
	double f_u = approximationError (fct, element);
	gradF[currentId++] = (f_u-f)/eps;
	d = gf->point(uu,vv+eps);
	v->setXYZ(d.x(),d.y(),d.z());
	double f_v = approximationError (fct, element);
	gradF[currentId++] = (f_v-f)/eps;
	v->setXYZ(p.x(),p.y(),p.z());
	//	if (iEl < 1)printf("df = %g %g\n",(f_u-f)/eps,(f_v-f)/eps);
      }
    }
  }
  for (int iV = 0; iV < nVert(); iV++)
    _vert[iV]->setXYZ(_xyz_temp[iV].x(),_xyz_temp[iV].y(),_xyz_temp[iV].z());

}