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

  MESHING3_RESULT OptimizeVolume (MeshingParameters & mp, 
				  Mesh & mesh3d)
    //				  const CSGeometry * geometry)
  {
    int i;

    PrintMessage (1, "Volume Optimization");

    /*
      if (!mesh3d.PureTetMesh())
      return MESHING3_OK;
    */

    // (*mycout) << "optstring = " << mp.optimize3d << endl;
    /*
      const char * optstr = globflags.GetStringFlag ("optimize3d", "cmh");
      int optsteps = int (globflags.GetNumFlag ("optsteps3d", 2));
    */

    mesh3d.CalcSurfacesOfNode();
    for (i = 1; i <= mp.optsteps3d; i++)
      {
	if (multithread.terminate)
	  break;

	MeshOptimize3d optmesh(mp);

	// teterrpow = mp.opterrpow;
	// for (size_t j = 1; j <= strlen(mp.optimize3d); j++)
        for (size_t j = 1; j <= mp.optimize3d.length(); j++)
	  {
	    if (multithread.terminate)
	      break;

	    switch (mp.optimize3d[j-1])
	      {
	      case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
	      case 'd': optmesh.SplitImprove(mesh3d); break;
	      case 's': optmesh.SwapImprove(mesh3d); break;
                // case 'u': optmesh.SwapImproveSurface(mesh3d); break;
	      case 't': optmesh.SwapImprove2(mesh3d); break;
#ifdef SOLIDGEOM
	      case 'm': mesh3d.ImproveMesh(*geometry); break;
	      case 'M': mesh3d.ImproveMesh(*geometry); break;
#else
	      case 'm': mesh3d.ImproveMesh(mp); break;
	      case 'M': mesh3d.ImproveMesh(mp); break;
#endif
	      case 'j': mesh3d.ImproveMeshJacobian(mp); break;
	      }
	  }
	mesh3d.mglevels = 1;
	MeshQuality3d (mesh3d);
      }
  
    return MESHING3_OK;
  }

  DLL_HEADER void Optimize2d (Mesh & mesh, MeshingParameters & mp)
  {
    static int timer = NgProfiler::CreateTimer ("optimize2d");
    NgProfiler::RegionTimer reg(timer);

    mesh.CalcSurfacesOfNode();

    const char * optstr = mp.optimize2d.c_str();
    int optsteps = mp.optsteps2d;

    for (int i = 1; i <= optsteps; i++)
      for (size_t j = 1; j <= strlen(optstr); j++)
	{
	  if (multithread.terminate) break;
	  switch (optstr[j-1])
	    {
	    case 's': 
	      {  // topological swap
		MeshOptimize2d meshopt;
                meshopt.SetMetricWeight (mp.elsizeweight);
		meshopt.EdgeSwapping (mesh, 0);
		break;
	      }
	    case 'S': 
	      {  // metric swap
		MeshOptimize2d meshopt;
                meshopt.SetMetricWeight (mp.elsizeweight);
		meshopt.EdgeSwapping (mesh, 1);
		break;
	      }
	    case 'm': 
	      {
		MeshOptimize2d meshopt;
                meshopt.SetMetricWeight (mp.elsizeweight);
		meshopt.ImproveMesh(mesh, mp);
		break;
	      }
	    case 'c': 
	      {
		MeshOptimize2d meshopt;
                meshopt.SetMetricWeight (mp.elsizeweight);
		meshopt.CombineImprove(mesh);
		break;
	      }
	    default:
	      cerr << "Optimization code " << optstr[j-1] << " not defined" << endl;
	    }  
	}
  }

   // Mesh the edges and add Face descriptors to prepare for surface meshing
   DLL_HEADER Ng_Result Ng_OCC_GenerateSurfaceMesh (Ng_OCC_Geometry * geom,
                                                    Ng_Mesh * mesh,
                                                    Ng_Meshing_Parameters * mp)
   {
      int numpoints = 0;

      OCCGeometry * occgeom = (OCCGeometry*)geom;
      Mesh * me = (Mesh*)mesh;

      // Set the internal meshing parameters structure from the nglib meshing 
      // parameters structure
      mp->Transfer_Parameters();


      // Only go into surface meshing if the face descriptors have already been added
      if(!me->GetNFD())
         return NG_ERROR;

      numpoints = me->GetNP();

      // Initially set up only for surface meshing without any optimisation
      int perfstepsend = MESHCONST_MESHSURFACE;

      // Check and if required, enable surface mesh optimisation step
      if(mp->optsurfmeshenable)
      {
         perfstepsend = MESHCONST_OPTSURFACE;
      }

      OCCMeshSurface(*occgeom, *me, perfstepsend);

      me->CalcSurfacesOfNode();
      
      if(me->GetNP() <= numpoints)
         return NG_ERROR;

      if(me->GetNSE() <= 0)
         return NG_ERROR;

      return NG_OK;
   }

  void RemoveIllegalElements (Mesh & mesh3d)
  {
    int it = 10;
    int nillegal, oldn;

    PrintMessage (1, "Remove Illegal Elements");
    // return, if non-pure tet-mesh
    /*
      if (!mesh3d.PureTetMesh())
      return;
    */
    mesh3d.CalcSurfacesOfNode();

    nillegal = mesh3d.MarkIllegalElements();

    MeshingParameters dummymp;
    MeshOptimize3d optmesh(dummymp);
    while (nillegal && (it--) > 0)
      {
	if (multithread.terminate)
	  break;

	PrintMessage (5, nillegal, " illegal tets");
	optmesh.SplitImprove (mesh3d, OPT_LEGAL);

	mesh3d.MarkIllegalElements();  // test
	optmesh.SwapImprove (mesh3d, OPT_LEGAL);
	mesh3d.MarkIllegalElements();  // test
	optmesh.SwapImprove2 (mesh3d, OPT_LEGAL);

	oldn = nillegal;
	nillegal = mesh3d.MarkIllegalElements();

	if (oldn != nillegal)
	  it = 10;
      }
    PrintMessage (5, nillegal, " illegal tets");
  }

//.........这里部分代码省略.........
			    lamedge * (*should.Get(i))(j) +
			    (1.-lamedge) * (*can.Get(i))(j);
		      }
		    else
		      mesh.Point(i) = *can.Get(i);
		  }
		if(facokedge < 1.-1e-8)
		  {
		    ostrstr.str("");
		    ostrstr << "worsening: " <<
		      Validate(mesh,bad_elements,pure_badness,max_worsening,uselocalworsening);

		    PrintMessage(5,ostrstr.str());
		  }
		else
		  Validate(mesh,bad_elements,pure_badness,-1,uselocalworsening);


		ostrstr.str("");
		ostrstr << bad_elements.Size() << " bad elements";
		PrintMessage(5,ostrstr.str());
	      }
	    while (bad_elements.Size() > 0 && 
		   cnttrials < maxtrials &&
		   multithread.terminate != 1);
	  }

	if(cnttrials < maxtrials &&
	   multithread.terminate != 1)
	  {
	    facokedge = factryedge;
	    
	    // smooth faces
	    mesh.CalcSurfacesOfNode();
	    
	    MeshingParameters dummymp;
	    mesh.ImproveMeshJacobianOnSurface(dummymp,isworkingboundary,nv,OPT_QUALITY, &idmaps);
	    
	    for (int i = 1; i <= np; i++)
	      *can.Elem(i) = mesh.Point(i);
	    
	    if(optimizer2d)
	      optimizer2d->ProjectBoundaryPoints(surfaceindex,can,should);
	  }


	oldlamface = lamface;
	lamface *= 6;
	if (lamface > 2)
	  lamface = 2;


	if(cnttrials < maxtrials &&
	   multithread.terminate != 1)
	  {

	    do  // move faces
	      {
		lamface *= 0.5;
		cnttrials++;
		if(cnttrials % 10 == 0)
		  max_worsening *= 1.1;
		factryface = lamface + (1.-lamface) * facokface;

		ostrstr.str("");
		ostrstr << "lamface = " << lamface << ", trying: " << factryface;

  // extern double teterrpow; 
  MESHING3_RESULT MeshVolume (MeshingParameters & mp, Mesh& mesh3d)
  {
     int oldne;
     int meshed;

     Array<INDEX_2> connectednodes;

     if (&mesh3d.LocalHFunction() == NULL) mesh3d.CalcLocalH(mp.grading);

     mesh3d.Compress();

     //  mesh3d.PrintMemInfo (cout);

     if (mp.checkoverlappingboundary)
        if (mesh3d.CheckOverlappingBoundary())
           throw NgException ("Stop meshing since boundary mesh is overlapping");

     int nonconsist = 0;
     for (int k = 1; k <= mesh3d.GetNDomains(); k++)
     {
        PrintMessage (3, "Check subdomain ", k, " / ", mesh3d.GetNDomains());

        mesh3d.FindOpenElements(k);

        /*
        bool res = mesh3d.CheckOverlappingBoundary();
        if (res)
        {
        PrintError ("Surface is overlapping !!");
        nonconsist = 1;
        }
        */

        bool res = (mesh3d.CheckConsistentBoundary() != 0);
        if (res)
        {
           PrintError ("Surface mesh not consistent");
           nonconsist = 1;
        }
     }

     if (nonconsist)
     {
        PrintError ("Stop meshing since surface mesh not consistent");
        throw NgException ("Stop meshing since surface mesh not consistent");
     }

     double globmaxh = mp.maxh;

     for (int k = 1; k <= mesh3d.GetNDomains(); k++)
       {
	 if (multithread.terminate)
           break;
	 
	 PrintMessage (2, "");
	 PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
	 (*testout) << "Meshing subdomain " << k << endl;
	 
	 mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
	 
	 mesh3d.CalcSurfacesOfNode();
	 mesh3d.FindOpenElements(k);
	 
	 if (!mesh3d.GetNOpenElements())
           continue;
	 
	 

	 Box<3> domain_bbox( Box<3>::EMPTY_BOX ); 
	 
	 for (SurfaceElementIndex sei = 0; sei < mesh3d.GetNSE(); sei++)
	   {
	     const Element2d & el = mesh3d[sei];
	     if (el.IsDeleted() ) continue;
	     
	     if (mesh3d.GetFaceDescriptor(el.GetIndex()).DomainIn() == k ||
		 mesh3d.GetFaceDescriptor(el.GetIndex()).DomainOut() == k)
	       
	       for (int j = 0; j < el.GetNP(); j++)
		 domain_bbox.Add (mesh3d[el[j]]);
	   }
	 domain_bbox.Increase (0.01 * domain_bbox.Diam());
	 
	
        for (int qstep = 1; qstep <= 3; qstep++)
	  {
	    // cout << "openquads = " << mesh3d.HasOpenQuads() << endl;
	    if (mesh3d.HasOpenQuads())
	      {
		string rulefile = ngdir;
		
		const char ** rulep = NULL;
		switch (qstep)
		  {
		  case 1:
		    rulefile += "/rules/prisms2.rls";
		    rulep = prismrules2;
		    break;
		  case 2: // connect pyramid to triangle
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	    for (int j = 1; j <= sel.GetNP(); j++)
	      boundp.Set(sel.PNum(j));
	  }


	(*testout) << "bpoints:" << endl;
	for (int i = 1; i <= np; i++)
	  if (boundp.Test(i))
	    (*testout) << i << endl;

	double lam = 0.5;

	while (facok < 1-1e-8 && cnttrials > 0)
	  {
	    lam *= 4;
	    if (lam > 2) lam = 2;

	    do
	      {
		//	      cout << "trials: " << cnttrials << endl;
		lam *= 0.5;
		cnttrials--;

		cout << "lam = " << lam << endl;

		factry = lam + (1-lam) * facok;
		cout << "trying: " << factry << endl;

		for (int i = 1; i <= np; i++)
		  if (boundp.Test(i))
		    {
		      for (int j = 0; j < 3; j++)
			mesh.Point(i)(j) = 
			  lam * should.Get(i)(j) +
			  (1-lam) * can.Get(i)(j);
		    }
		  else
		    mesh.Point(i) = Point<3> (can.Get(i));
	      
		//	      (*testout) << "bad els: " << endl;
		wrongels = 0;
		for (int i = 1; i <= ne; i++)
		  {
		    if (!illegalels.Test(i) && 
			mesh.VolumeElement(i).
			CalcJacobianBadness(mesh.Points()) > 1e10)
		      {
			wrongels++;
			Element & el = mesh.VolumeElement(i);
			el.flags.badel = 1;
		     
		      
			if (lam < 1e-4)
			  illegalels.Set(i);
 

			/*
			  (*testout) << i << ": ";
			  for (j = 1; j <= el.GetNP(); j++)
			  (*testout) << el.PNum(j) << " ";
			  (*testout) << endl;
			*/
		      }
		    else
		      mesh.VolumeElement(i).flags.badel = 0;
		  }
		cout << "wrongels = " << wrongels << endl;
	      }
	    while (wrongels && cnttrials > 0);

	    mesh.CalcSurfacesOfNode();
	    MeshingParameters dummymp;
	    mesh.ImproveMeshJacobian (dummymp, OPT_WORSTCASE);	      
	  
	    facok = factry;
	    for (int i = 1; i <= np; i++)
	      can.Elem(i) = mesh.Point(i);
	  }
      }


      
    for (int i = 1; i <= ne; i++)
      {
	if (illegalels.Test(i))
	  {
	    cout << "illegal element: " << i << endl;
	    mesh.VolumeElement(i).flags.badel = 1;	
	  }
	else
	  mesh.VolumeElement(i).flags.badel = 0;	
      }
  
    /*
      if (cnttrials <= 0)
      {
      cerr << "ERROR: Sorry, illegal elements:" << endl;
      }
    */
  }