C++ Box::Diam方法代码示例

  void Meshing2 :: BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp)
  {
    double filldist = mp.filldist;

    cout << "blockfill local h" << endl;
    cout << "rel filldist = " << filldist << endl;
    PrintMessage (3, "blockfill local h");

    Array<Point<3> > npoints;

    // adfront -> CreateTrees();

    Box<3> bbox ( Box<3>::EMPTY_BOX );
    double maxh = 0;

    for (int i = 0; i < adfront->GetNFL(); i++)
      {
	const FrontLine & line = adfront->GetLine (i);

	const Point<3> & p1 = adfront->GetPoint(line.L().I1());
	const Point<3> & p2 = adfront->GetPoint(line.L().I2());

	double hi = Dist (p1, p2);
	if (hi > maxh) maxh = hi;

	bbox.Add (p1);
	bbox.Add (p2);
      }


    cout << "bbox = " << bbox << endl;


    // Point<3> mpc = bbox.Center();
    bbox.Increase (bbox.Diam()/2);
    Box<3> meshbox = bbox;

    LocalH loch2 (bbox, 1);

    if (mp.maxh < maxh) maxh = mp.maxh;

    bool changed;
    do
      {
	mesh.LocalHFunction().ClearFlags();

	for (int i = 0; i < adfront->GetNFL(); i++)
	  {
	    const FrontLine & line = adfront->GetLine(i);

	    Box<3> bbox (adfront->GetPoint (line.L().I1()));
	    bbox.Add (adfront->GetPoint (line.L().I2()));


	    double filld = filldist * bbox.Diam();
	    bbox.Increase (filld);

	    mesh.LocalHFunction().CutBoundary (bbox);
	  }


	mesh.LocalHFunction().FindInnerBoxes (adfront, NULL);

	npoints.SetSize(0);
	mesh.LocalHFunction().GetInnerPoints (npoints);

	changed = false;
	for (int i = 0; i < npoints.Size(); i++)
	  {
	    if (mesh.LocalHFunction().GetH(npoints[i]) > 1.5 * maxh)
	      {
		mesh.LocalHFunction().SetH (npoints[i], maxh);
		changed = true;
	      }
	  }
      }
    while (changed);

    if (debugparam.slowchecks)
      (*testout) << "Blockfill with points: " << endl;
    *testout << "loch = " << mesh.LocalHFunction() << endl;

    *testout << "npoints = " << endl << npoints << endl;

    for (int i = 1; i <= npoints.Size(); i++)
      {
	if (meshbox.IsIn (npoints.Get(i)))
	  {
	    PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
	    adfront->AddPoint (npoints.Get(i), gpnum);

	    if (debugparam.slowchecks)
	      {
		(*testout) << npoints.Get(i) << endl;

		Point<2> p2d (npoints.Get(i)(0), npoints.Get(i)(1));
		if (!adfront->Inside(p2d))
		  {
		    cout << "add outside point" << endl;
		    (*testout) << "outside" << endl;
//.........这里部分代码省略.........

  void VisualSceneGeometry :: BuildScene (int zoomall)
  {
    VisualScene::BuildScene(zoomall); // setting light ...
    Box<3> box;
    int hasp = 0;
    for (int i = 0; i < geometry->GetNTopLevelObjects(); i++)
      {
	const TriangleApproximation & ta =
	  *geometry->GetTriApprox(i);
	if (!&ta) continue;

	for (int j = 0; j < ta.GetNP(); j++)      
	  {
	    if (hasp)
	      box.Add (ta.GetPoint(j));
	    else
	      {
		hasp = 1;
		box.Set (ta.GetPoint(j));
	      }
	  }
      }
    if (hasp)
      {
	center = box.Center();
	rad = box.Diam() / 2;
      }
    else
      {
	center = Point3d(0,0,0);
	rad = 1;
      }

    CalcTransformationMatrices();

    for (int i = 0; i < trilists.Size(); i++)
      glDeleteLists (trilists[i], 1);
    trilists.SetSize(0);

    for (int i = 0; i < geometry->GetNTopLevelObjects(); i++)
      {
	trilists.Append (glGenLists (1));
	glNewList (trilists.Last(), GL_COMPILE); 
	glEnable (GL_NORMALIZE);
	const TriangleApproximation & ta =
	  *geometry->GetTriApprox(i);
	if (&ta) 
	  {
	    glEnableClientState(GL_VERTEX_ARRAY);
	    glVertexPointer(3, GL_DOUBLE, 0, &ta.GetPoint(0)(0));

	    glEnableClientState(GL_NORMAL_ARRAY);
	    glNormalPointer(GL_DOUBLE, 0, &ta.GetNormal(0)(0));
	    
	    for (int j = 0; j < ta.GetNT(); j++)
	      glDrawElements(GL_TRIANGLES, 3, GL_UNSIGNED_INT, & (ta.GetTriangle(j)[0]));

	    glDisableClientState(GL_VERTEX_ARRAY);
	    glDisableClientState(GL_NORMAL_ARRAY);
            /*
	    for (int j = 0; j < ta.GetNT(); j++)
	      {
		glBegin (GL_TRIANGLES);
		for (int k = 0; k < 3; k++)
		  {
		    int pi = ta.GetTriangle(j)[k];
		    glNormal3dv (ta.GetNormal(pi));
		    glVertex3dv (ta.GetPoint(pi));
                    cout << "v = " << ta.GetPoint(pi) << endl;
		  }
		glEnd ();
	      }
            */
	  }
	glEndList ();
      }
  }

bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
        const TopoDS_Shape& aShape)
{
    netgen::multithread.terminate = 0;
    //netgen::multithread.task = "Surface meshing";

    SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
    SMESH_MesherHelper helper(aMesh);
    helper.SetElementsOnShape( true );

    NETGENPlugin_NetgenLibWrapper ngLib;
    ngLib._isComputeOk = false;

    netgen::Mesh   ngMeshNoLocSize;
#if NETGEN_VERSION < 6
    netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh,  & ngMeshNoLocSize };
#else
    netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh.get(),  & ngMeshNoLocSize };
#endif
    netgen::OCCGeometry occgeoComm;

    // min / max sizes are set as follows:
    // if ( _hypParameters )
    //    min and max are defined by the user
    // else if ( _hypLengthFromEdges )
    //    min = aMesher.GetDefaultMinSize()
    //    max = average segment len of a FACE
    // else if ( _hypMaxElementArea )
    //    min = aMesher.GetDefaultMinSize()
    //    max = f( _hypMaxElementArea )
    // else
    //    min = aMesher.GetDefaultMinSize()
    //    max = max segment len of a FACE

    NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false);
    aMesher.SetParameters( _hypParameters ); // _hypParameters -> netgen::mparam
    const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true;
    if ( _hypMaxElementArea )
    {
        netgen::mparam.maxh = sqrt( 2. * _hypMaxElementArea->GetMaxArea() / sqrt(3.0) );
    }
    if ( _hypQuadranglePreference )
        netgen::mparam.quad = true;

    // local size is common for all FACEs in aShape?
    const bool isCommonLocalSize = ( !_hypLengthFromEdges && !_hypMaxElementArea && netgen::mparam.uselocalh );
    const bool isDefaultHyp = ( !_hypLengthFromEdges && !_hypMaxElementArea && !_hypParameters );

    if ( isCommonLocalSize ) // compute common local size in ngMeshes[0]
    {
        //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm
        aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM );

        // local size set at MESHCONST_ANALYSE step depends on
        // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user
        if ( !_hypParameters || netgen::mparam.minh < DBL_MIN )
        {
            if ( !_hypParameters )
                netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.;
            netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh );
        }
        // set local size depending on curvature and NOT closeness of EDGEs
        netgen::occparam.resthcloseedgeenable = false;
        //netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
        occgeoComm.face_maxh = netgen::mparam.maxh;
        netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] );
        occgeoComm.emap.Clear();
        occgeoComm.vmap.Clear();

        // set local size according to size of existing segments
        const double factor = netgen::occparam.resthcloseedgefac;
        TopTools_IndexedMapOfShape edgeMap;
        TopExp::MapShapes( aMesh.GetShapeToMesh(), TopAbs_EDGE, edgeMap );
        for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
        {
            const TopoDS_Shape& edge = edgeMap( iE );
            if ( SMESH_Algo::isDegenerated( TopoDS::Edge( edge ))/* ||
           helper.IsSubShape( edge, aShape )*/)
                continue;
            SMESHDS_SubMesh* smDS = meshDS->MeshElements( edge );
            if ( !smDS ) continue;
            SMDS_ElemIteratorPtr segIt = smDS->GetElements();
            while ( segIt->more() )
            {
                const SMDS_MeshElement* seg = segIt->next();
                SMESH_TNodeXYZ n1 = seg->GetNode(0);
                SMESH_TNodeXYZ n2 = seg->GetNode(1);
                gp_XYZ p = 0.5 * ( n1 + n2 );
                netgen::Point3d pi(p.X(), p.Y(), p.Z());
                ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() );
            }
        }
    }
    netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization

    // ==================
    // Loop on all FACEs
    // ==================

    vector< const SMDS_MeshNode* > nodeVec;
//.........这里部分代码省略.........

  void MeshFromSpline2D (SplineGeometry2d & geometry,
			 Mesh *& mesh, 
			 MeshingParameters & mp)
  {
    PrintMessage (1, "Generate Mesh from spline geometry");

    double h = mp.maxh;

    Box<2> bbox = geometry.GetBoundingBox ();

    if (bbox.Diam() < h) 
      {
	h = bbox.Diam();
	mp.maxh = h;
      }

    mesh = new Mesh;
    mesh->SetDimension (2);

    geometry.PartitionBoundary (h, *mesh);

    // marks mesh points for hp-refinement
    for (int i = 0; i < geometry.GetNP(); i++)
      if (geometry.GetPoint(i).hpref)
	{
	  double mindist = 1e99;
	  PointIndex mpi(0);
	  Point<2> gp = geometry.GetPoint(i);
	  Point<3> gp3(gp(0), gp(1), 0);
	  for (PointIndex pi = PointIndex::BASE; 
	       pi < mesh->GetNP()+PointIndex::BASE; pi++)
	    if (Dist2(gp3, (*mesh)[pi]) < mindist)
	      {
		mpi = pi;
		mindist = Dist2(gp3, (*mesh)[pi]);
	      }
	  (*mesh)[mpi].Singularity(1.);
	}


    int maxdomnr = 0;
    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
      {
	if ( (*mesh)[si].domin > maxdomnr) maxdomnr = (*mesh)[si].domin;
	if ( (*mesh)[si].domout > maxdomnr) maxdomnr = (*mesh)[si].domout;
      }

    mesh->ClearFaceDescriptors();
    for (int i = 1; i <= maxdomnr; i++)
      mesh->AddFaceDescriptor (FaceDescriptor (i, 0, 0, i));

    // set Array<string*> bcnames... 
    // number of bcnames
    int maxsegmentindex = 0;
    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
      {
	if ( (*mesh)[si].si > maxsegmentindex) maxsegmentindex = (*mesh)[si].si;
      }

    mesh->SetNBCNames(maxsegmentindex);

    for ( int sindex = 0; sindex < maxsegmentindex; sindex++ )
      mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) );

    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
      (*mesh)[si].SetBCName ( (*mesh).GetBCNamePtr( (*mesh)[si].si-1 ) );

    Point3d pmin(bbox.PMin()(0), bbox.PMin()(1), -bbox.Diam());
    Point3d pmax(bbox.PMax()(0), bbox.PMax()(1), bbox.Diam());

    mesh->SetLocalH (pmin, pmax, mparam.grading);
    mesh->SetGlobalH (h);
  
    mesh->CalcLocalH();

    int bnp = mesh->GetNP(); // boundary points

    int hquad = mparam.quad;


    for (int domnr = 1; domnr <= maxdomnr; domnr++)
      if (geometry.GetDomainTensorMeshing (domnr))
        { // tensor product mesh
          
          Array<PointIndex, PointIndex::BASE> nextpi(bnp);
          Array<int, PointIndex::BASE> si1(bnp), si2(bnp);
          PointIndex firstpi;
          
          nextpi = -1;
          si1 = -1;
          si2 = -1;
          for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
            {
              int p1 = -1, p2 = -2;

              if ( (*mesh)[si].domin == domnr)
                { p1 = (*mesh)[si][0]; p2 = (*mesh)[si][1]; }
              if ( (*mesh)[si].domout == domnr)
                { p1 = (*mesh)[si][1]; p2 = (*mesh)[si][0]; }
              
//.........这里部分代码省略.........

/*
double stlgh;
double GetH(const Point3d& p, double x)
{
  return stlgh;//+0.5)*(x+0.5);
}
*/
STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
                         Array<Point3d>& mp, double ghi,
                         class Mesh& mesh) const
{
    static int timer1a = NgProfiler::CreateTimer ("mesh stl-line 1a");
    static int timer1b = NgProfiler::CreateTimer ("mesh stl-line 1b");
    static int timer2 = NgProfiler::CreateTimer ("mesh stl-line 2");
    static int timer3 = NgProfiler::CreateTimer ("mesh stl-line 3");

    NgProfiler::StartTimer (timer1a);

    STLLine* line = new STLLine(geometry);

    //stlgh = ghi; //uebergangsloesung!!!!

    double len = GetLength(ap);
    double inthl = 0; //integral of 1/h
    double dist = 0;
    double h;
    int ind;
    Point3d p;

    Box<3> bbox;
    GetBoundingBox (ap, bbox);
    double diam = bbox.Diam();

    double minh = mesh.LocalHFunction().GetMinH (bbox.PMin(), bbox.PMax());

    double maxseglen = 0;
    for (int i = 1; i <= GetNS(); i++)
        maxseglen = max2 (maxseglen, GetSegLen (ap, i));

    int nph = 10+int(maxseglen / minh); //anzahl der integralauswertungen pro segment

    Array<double> inthi(GetNS()*nph);
    Array<double> curvelen(GetNS()*nph);

    NgProfiler::StopTimer (timer1a);
    NgProfiler::StartTimer (timer1b);


    for (int i = 1; i <= GetNS(); i++)
    {
        //double seglen = GetSegLen(ap,i);
        for (int j = 1; j <= nph; j++)
        {
            p = GetPointInDist(ap,dist,ind);
            //h = GetH(p,dist/len);
            h = mesh.GetH(p);


            dist += GetSegLen(ap,i)/(double)nph;

            inthl += GetSegLen(ap,i)/nph/(h);
            inthi.Elem((i-1)*nph+j) = GetSegLen(ap,i)/nph/h;
            curvelen.Elem((i-1)*nph+j) = GetSegLen(ap,i)/nph;
        }
    }


    int inthlint = int(inthl+1);

    if ( (inthlint < 3) && (StartP() == EndP()))
    {
        inthlint = 3;
    }
    if ( (inthlint == 1) && ShouldSplit())
    {
        inthlint = 2;
    }

    double fact = inthl/(double)inthlint;
    dist = 0;
    int j = 1;


    p = ap.Get(StartP());
    int pn = AddPointIfNotExists(mp, p, 1e-10*diam);

    int segn = 1;
    line->AddPoint(pn);
    line->AddLeftTrig(GetLeftTrig(segn));
    line->AddRightTrig(GetRightTrig(segn));
    line->AddDist(dist);

    NgProfiler::StopTimer (timer1b);
    NgProfiler::StartTimer (timer2);

    inthl = 0; //restart each meshseg
    for (int i = 1; i <= inthlint; i++)
    {
        while (inthl < 1.000000001 && j <= inthi.Size())
        {
//.........这里部分代码省略.........