C++ Mesh::AddPoint方法代码示例

Mesh *StockMeshes::NewSphereMesh(float width, float height, float depth, const Color &color,
                                 int precision, const Texture *texture) {
  Mesh *result = new Mesh(precision*precision+2, 2*precision*precision, true, true, true);
  float w = width/2, h = height/2, d = depth/2;

  // Create points
  result->AddPoint(Point3(0, h, 0), Point3(0, h, 0), color, 0.5f, 0);
  for (int i = 0; i < precision; i++)
    for (int j = 0; j < precision; j++) {
      float v = float(i+1) / (precision+1);
      float u = float(j) / precision;
      float r = sin(v*kPi);
      Point3 p(-cos(2*u*kPi)*w*r, cos(v*kPi)*h, sin(2*u*kPi)*w*d*r);
      result->AddPoint(p, p, color, u, v);
    }
  result->AddPoint(Point3(0, -h, 0), Point3(0, -h, 0), color, 0.5f, 0);

  // Create triangles
  for (int j = 0; j < precision; j++) {
    int j2 = (j+1) % precision;
    result->AddTriangle(0, 1 + j, 1 + j2, texture);
  }
  for (int i = 0; i < precision-1; i++)
    for (int j = 0; j < precision; j++) {
      int j2 = (j+1) % precision;
      result->AddTriangle(1 + i*precision + j, 1 + (i+1)*precision + j, 1 + (i+1)*precision + j2,
                          texture);
      result->AddTriangle(1 + i*precision + j, 1 + (i+1)*precision + j2, 1 + i*precision + j2,
                          texture);
    }
  for (int j = 0; j < precision; j++) {
    int j2 = (j+1) % precision;
    result->AddTriangle(precision*(precision-1) + j + 1, 1 + precision*precision,
                        precision*(precision-1) + j2 + 1, texture);
  }
  return result;
}

  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;
//.........这里部分代码省略.........

   /* ------------------ 2D Meshing Functions ------------------------- */
   DLL_HEADER void Ng_AddPoint_2D (Ng_Mesh * mesh, double * x)
   {
      Mesh * m = (Mesh*)mesh;

      m->AddPoint (Point3d (x[0], x[1], 0));
   }

Mesh *StockMeshes::NewBoxMesh(float width, float height, float depth, const Color &color,
                              const Texture *texture) {
  float x = width/2, y = height/2, z = depth/2;
  Mesh *result = new Mesh(24, 12, true, true, true);

  // Top face
  result->AddPoint(Point3(-x,  y, -z), Point3(0, 1, 0), color, 0, 1);
  result->AddPoint(Point3(-x,  y,  z), Point3(0, 1, 0), color, 0, 0);
  result->AddPoint(Point3( x,  y,  z), Point3(0, 1, 0), color, 1, 0);
  result->AddPoint(Point3( x,  y, -z), Point3(0, 1, 0), color, 1, 1);
  result->AddTriangle(0, 1, 2, texture);
  result->AddTriangle(0, 2, 3, texture);

  // Front face
  result->AddPoint(Point3(-x, -y,  -z), Point3(0, 0, -1), color, 0, 1);
  result->AddPoint(Point3(-x,  y,  -z), Point3(0, 0, -1), color, 0, 0);
  result->AddPoint(Point3( x,  y,  -z), Point3(0, 0, -1), color, 1, 0);
  result->AddPoint(Point3( x, -y,  -z), Point3(0, 0, -1), color, 1, 1);
  result->AddTriangle(4, 5, 6, texture);
  result->AddTriangle(4, 6, 7, texture);

  // Left face
  result->AddPoint(Point3(-x, -y,  z), Point3(-1, 0, 0), color, 0, 1);
  result->AddPoint(Point3(-x,  y,  z), Point3(-1, 0, 0), color, 0, 0);
  result->AddPoint(Point3(-x,  y, -z), Point3(-1, 0, 0), color, 1, 0);
  result->AddPoint(Point3(-x, -y, -z), Point3(-1, 0, 0), color, 1, 1);
  result->AddTriangle(8, 9, 10, texture);
  result->AddTriangle(8, 10, 11, texture);

  // Back face
  result->AddPoint(Point3( x, -y,  z), Point3(0, 0, 1), color, 0, 1);
  result->AddPoint(Point3( x,  y,  z), Point3(0, 0, 1), color, 0, 0);
  result->AddPoint(Point3(-x,  y,  z), Point3(0, 0, 1), color, 1, 0);
  result->AddPoint(Point3(-x, -y,  z), Point3(0, 0, 1), color, 1, 1);
  result->AddTriangle(12, 13, 14, texture);
  result->AddTriangle(12, 14, 15, texture);

  // Right face
  result->AddPoint(Point3( x, -y, -z), Point3(1, 0, 0), color, 0, 1);
  result->AddPoint(Point3( x,  y, -z), Point3(1, 0, 0), color, 0, 0);
  result->AddPoint(Point3( x,  y,  z), Point3(1, 0, 0), color, 1, 0);
  result->AddPoint(Point3( x, -y,  z), Point3(1, 0, 0), color, 1, 1);
  result->AddTriangle(16, 17, 18, texture);
  result->AddTriangle(16, 18, 19, texture);

  // Bottom face
  result->AddPoint(Point3(-x, -y,  z), Point3(0, -1, 0), color, 0, 1);
  result->AddPoint(Point3(-x, -y, -z), Point3(0, -1, 0), color, 0, 0);
  result->AddPoint(Point3( x, -y, -z), Point3(0, -1, 0), color, 1, 0);
  result->AddPoint(Point3( x, -y,  z), Point3(0, -1, 0), color, 1, 1);
  result->AddTriangle(20, 21, 22, texture);
  result->AddTriangle(20, 22, 23, texture);

  return result;
}

  void RefinePrisms (Mesh & mesh, const CSGeometry * geom, 
		     ZRefinementOptions & opt)
  {
    int i, j;
    bool found, change;
    int cnt = 0;


    // markers for z-refinement:  p1, p2, levels  
    // p1-p2 is an edge to be refined
    ARRAY<INDEX_3> ref_uniform;
    ARRAY<INDEX_3> ref_singular;
    ARRAY<INDEX_4 > ref_slices;

    BitArray first_id(geom->identifications.Size());
    first_id.Set();

  
    INDEX_2_HASHTABLE<int> & identpts = 
      mesh.GetIdentifications().GetIdentifiedPoints ();

    if (&identpts)
      {
	for (i = 1; i <= identpts.GetNBags(); i++)
	  for (j = 1; j <= identpts.GetBagSize(i); j++)
	    {
	      INDEX_2 pair;
	      int idnr;
	      identpts.GetData(i, j, pair, idnr);
	      const CloseSurfaceIdentification * csid = 
		dynamic_cast<const CloseSurfaceIdentification*> 
		(geom->identifications.Get(idnr));
	      if (csid)
		{
		  if (!csid->GetSlices().Size())
		    {
		      if (first_id.Test (idnr))
			{
			  first_id.Clear(idnr);
			  ref_uniform.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels()));
			  ref_singular.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels1()));
			  ref_singular.Append (INDEX_3 (pair.I2(), pair.I1(), csid->RefLevels2()));
			}
		    }
		  else
		    {   
		      //const ARRAY<double> & slices = csid->GetSlices();
		      INDEX_4 i4;
		      i4[0] = pair.I1();
		      i4[1] = pair.I2();
		      i4[2] = idnr;
		      i4[3] = csid->GetSlices().Size();
		      ref_slices.Append (i4);
		    }
		}
	    }
      }

  
  
    ARRAY<EdgePointGeomInfo> epgi;

    while (1)
      {
	cnt++;
	PrintMessage (3, "Z-Refinement, level = ", cnt);
	INDEX_2_HASHTABLE<int> refedges(mesh.GetNSE()+1);


	found = 0;
	// mark prisms due to close surface flags:
	int oldsize = ref_uniform.Size();
	for (i = 1; i <= oldsize; i++)
	  {
	    int pi1 = ref_uniform.Get(i).I1();
	    int pi2 = ref_uniform.Get(i).I2();
	    int levels = ref_uniform.Get(i).I3();

	    if (levels > 0)
	      {
		const Point3d & p1 = mesh.Point(pi1);
		const Point3d & p2 = mesh.Point(pi2);
		int npi(0);
	      
		INDEX_2 edge(pi1, pi2);
		edge.Sort();
		if (!refedges.Used(edge))
		  {
		    Point3d np = Center (p1, p2);
		    npi = mesh.AddPoint (np);
		    refedges.Set (edge, npi);
		    found = 1;
		  }

		ref_uniform.Elem(i) = INDEX_3(pi1, npi, levels-1);
		ref_uniform.Append (INDEX_3(pi2, npi, levels-1));
	      }
	  }
	for (i = 1; i <= ref_singular.Size(); i++)
	  {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
               int selNP = sel.GetNP();
               for(j = 1; j <= selNP; j++)
               {
                  // Set the bitarray to indicate that the 
                  // point is part of the required set
                  bndnodes.Set(sel.PNum(j));
		  
                  // Vec3d& surfacenormal = Vec3d();   ????
                  Vec3d surfacenormal;

                  // Calculate the surface normal at the current point 
                  // with respect to the current surface element
                  GetSurfaceNormal(mesh,sel,j,surfacenormal);
                  
                  // Add the surface normal to the already existent one 
                  // (This gives the effective normal direction at corners 
                  //  and curved areas)
                  growthvectors.Elem(sel.PNum(j)) = growthvectors.Elem(sel.PNum(j)) 
                                                    + surfacenormal;
               }
            }
         }

         // Add additional points into the mesh structure in order to 
         // clone the surface elements.
         // Also invert the growth vectors so that they point inwards, 
         // and normalize them
         cout << "Cloning points and calculating growth vectors...." << endl;

         for (i = 1; i <= np; i++)
         {
            if (bndnodes.Test(i))
            {
               mapto.Elem(i) = mesh.AddPoint (mesh.Point (i));

               growthvectors.Elem(i).Normalize();
               growthvectors.Elem(i) *= -1.0;
            }
            else
            {
               mapto.Elem(i) = 0;
               growthvectors.Elem(i) = Vec3d(0,0,0);
            }
         }


         // Add quad surface elements at edges for surfaces which 
         // dont have boundary layers

         // Bit array to keep track of segments already processed
         BitArray segsel(nseg);

         // Set them all to "1" to initially activate all segments
         segsel.Set();

         cout << "Adding 2D Quad elements on required surfaces...." << endl;

         for (i = 1; i <= nseg; i++)
         {
            int seg_p1 = mesh.LineSegment(i)[0];
            int seg_p2 = mesh.LineSegment(i)[1];

            // Only go in if the segment is still active, and if both its 
            // surface index is part of the "hit-list"
            if(segsel.Test(i) && surfid.Contains(mesh.LineSegment(i).si))
            {

void SplineSeg<D> :: Partition (double h, double elto0,
				Mesh & mesh, Point3dTree & searchtree, int segnr) const
{
  int i, j;
  double l, r1, r2, ra;
  double lold, dt, frac;
  int n = 100;
  Point<D> p, pold, mark, oldmark;
  ARRAY<double> curvepoints;
  double edgelength, edgelengthold;
  l = Length();

  r1 = StartPI().refatpoint;
  r2 = EndPI().refatpoint;
  ra = reffak;

  //  cout << "Partition, l = " << l << ", h = " << h << endl;
  CalcPartition (l, h, r1, r2, ra, elto0, curvepoints);
  //  cout << "curvepoints = " << curvepoints << endl;

  dt = 1.0 / n;

  l = 0;
  j = 1;

  pold = GetPoint (0);
  lold = 0;
  oldmark = pold;
  edgelengthold = 0;
  ARRAY<int> locsearch;

  for (i = 1; i <= n; i++)
    {
      p = GetPoint (i*dt);
      l = lold + Dist (p, pold);
      while (j < curvepoints.Size() && (l >= curvepoints[j] || i == n))
	{
	  frac = (curvepoints[j]-lold) / (l-lold);
	  mark = pold + frac * (p-pold);
	  edgelength = i*dt + (frac-1)*dt;
	  {
	    PointIndex pi1 = -1, pi2 = -1;
	  
	    Point3d mark3(mark(0), mark(1), 0);
	    Point3d oldmark3(oldmark(0), oldmark(1), 0);

	    Vec<3> v (1e-4*h, 1e-4*h, 1e-4*h);
	    searchtree.GetIntersecting (oldmark3 - v, oldmark3 + v, locsearch);
	    if (locsearch.Size()) pi1 = locsearch[0];
	      
	    searchtree.GetIntersecting (mark3 - v, mark3 + v, locsearch);
	    if (locsearch.Size()) pi2 = locsearch[0];
	    /*	    
	      for (PointIndex pk = PointIndex::BASE; 
	      pk < mesh.GetNP()+PointIndex::BASE; pk++)
	      {
	      if (Dist (mesh[pk], oldmark3) < 1e-4 * h) pi1 = pk;
	      if (Dist (mesh[pk], mark3) < 1e-4 * h) pi2 = pk;
	      }
	    */
	    

	    //	    cout << "pi1 = " << pi1 << endl;
	    //	    cout << "pi2 = " << pi2 << endl;
	    
	    if (pi1 == -1)
	      {
		pi1 = mesh.AddPoint(oldmark3);
		searchtree.Insert (oldmark3, pi1);
	      }
	    if (pi2 == -1)
	      {
		pi2 = mesh.AddPoint(mark3);
		searchtree.Insert (mark3, pi2);
	      }

	    // cout << "pi1 = " << pi1 << endl;
	    // cout << "pi2 = " << pi2 << endl;
	  
	    Segment seg;
	    seg.edgenr = segnr;
	    seg.si = bc; // segnr;
	    seg.p1 = pi1;
	    seg.p2 = pi2;
	    seg.domin = leftdom;
	    seg.domout = rightdom;
	    seg.epgeominfo[0].edgenr = segnr;
	    seg.epgeominfo[0].dist = edgelengthold;
	    seg.epgeominfo[1].edgenr = segnr;
	    seg.epgeominfo[1].dist = edgelength;
	    seg.singedge_left = hpref_left;
	    seg.singedge_right = hpref_right;
	    mesh.AddSegment (seg);
	  }
	
	  oldmark = mark;
	  edgelengthold = edgelength;
	  j++;
	}
    
//.........这里部分代码省略.........

  void Refinement :: MakeSecondOrder (Mesh & mesh)
  {
    int nseg, nse, ne;

    mesh.ComputeNVertices();
    mesh.SetNP(mesh.GetNV());
  
    INDEX_2_HASHTABLE<int> between(mesh.GetNP() + 5);


    bool thinlayers = 0;
    for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++)
      if (mesh[ei].GetType() == PRISM ||
	  mesh[ei].GetType() == PRISM12)
	thinlayers = 1;
    

    nseg = mesh.GetNSeg();
    for (SegmentIndex si = 0; si < nseg; si++)
      {
	Segment & el = mesh.LineSegment(si);

	INDEX_2 i2 = INDEX_2::Sort (el[0], el[1]);

	if (between.Used(i2))
	  el[2] = between.Get(i2);
	else
	  {
	    Point<3> pb;
	    EdgePointGeomInfo ngi;
            PointBetween (mesh.Point (el[0]),
                          mesh.Point (el[1]), 0.5,
			  el.surfnr1, el.surfnr2,
			  el.epgeominfo[0], el.epgeominfo[1],
			  pb, ngi);
	  
	    el[2] = mesh.AddPoint (pb, mesh.Point(el[0]).GetLayer(), 
				   EDGEPOINT);
	    between.Set (i2, el[2]);
	  }
      }

    // refine surface elements
    nse = mesh.GetNSE();
    for (SurfaceElementIndex sei = 0; sei < nse; sei++)
      {
	int j;
	const Element2d & el = mesh.SurfaceElement(sei);

	int onp(0);
      
	Element2d newel;
	newel.SetIndex (el.GetIndex());

	static int betw_trig[3][3] =
	  { { 1, 2, 3 },
	    { 0, 2, 4 },
	    { 0, 1, 5 } };
	static int betw_quad6[2][3] =
	  { { 0, 1, 4 },
	    { 3, 2, 5 } };
	static int betw_quad8[4][3] =
	  { { 0, 1, 4 },
	    { 3, 2, 5 },
	    { 0, 3, 6 },
	    { 1, 2, 7 } };
	int (*betw)[3] = NULL;
      
	switch (el.GetType())
	  {
	  case TRIG:
	  case TRIG6:
	    {
	      betw = betw_trig;
	      newel.SetType (TRIG6);
	      onp = 3;
	      break;
	    }
	  case QUAD:
	  case QUAD6: 
	  case QUAD8:
	    {
	      if (thinlayers)
		{
		  betw = betw_quad6;
		  newel.SetType (QUAD6);
		}
	      else
		{
		  betw = betw_quad8;
		  newel.SetType (QUAD8);
		}
	      onp = 4;
	      break;
	    }
	  default:
	    PrintSysError ("Unhandled element in secondorder:", int(el.GetType()));
	  }

	for (j = 0; j < onp; j++)
//.........这里部分代码省略.........